Changeset 21 for bdm/math

Timestamp:

02/18/08 09:04:46 (16 years ago)

Author:

mido

Message:

patch ldformu

Files:

• bdm/math/libDC.cpp (modified) (3 diffs)

bdm/math/libDC.cpp

@@ -42 +42 @@
         D0 = ones(dim);
         
-
-        using std::cout;
-        cout<<V;
         chol(V,F);
         // L and D will be allocated by ldform()
@@ -205 +202 @@
 }
 
-void ldmat::ldform( mat &A,vec &D0 ) {
+void ldmat::ldform( mat &A,vec &D0 )
+{
         int m = A.rows();
         int n = A.cols();
@@ -215 +213 @@
         L=concat_vertical( zeros( n,n ), diag( sqrt( D0 ) )*A );
         D=zeros( n+m );
-        
-        //unnecessary big L and D will be made smaller at the end of file
-        
-        vec w=zeros( n );
-        vec v=zeros(n);
+
+        //unnecessary big L and D will be made smaller at the end of file       
+        vec w=zeros( n+m );
+   
         double sum, beta, pom;
 
         int cc=0;
-        int i=n; // +1 in .m
+        int i=n; // indexovani o 1 niz, nez v matlabu
         int ii,j,jj;
-        while (( i> ( n-mn+1-cc ) )&&( i>1 ) ) {
+        while ( (i>n-mn-cc) && (i>0) ) 
+        {
                 i--;
                 sum = 0.0;
-                v.set_size( m+i-( n-cc ) ); //prepare v
-                for ( ii=n-cc;ii<m+i;i++ ) {
+
+                int last_v = m+i-n+cc+1;
+        
+                vec v = zeros( last_v + 1 ); //prepare v
+                for ( ii=n-cc-1;ii<m+i+1;ii++ ) 
+                {
                         sum+= L( ii,i )*L( ii,i );
-                        v( ii )=L( ii,i ); //assign v
-                }
-                if ( L( m+i,i )==0 ) {
-                        beta = sqrt( sum );
-                } else {
-                        beta = L( m+i,i )+sign( L( m+i,i ) )*sqrt( sum ) ;
-                }
-                if ( std::fabs( beta )<eps ) {
+                        v( ii-n+cc+1 )=L( ii,i ); //assign v
+                }
+
+                if ( L( m+i,i )==0 ) 
+                        beta = sqrt( sum );             
+                else
+                        beta = L( m+i,i )+sign( L( m+i,i ) )*sqrt( sum );               
+     
+                if ( std::fabs( beta )<eps )
+                {
                         cc++;
-                        L.set_row( n+1-cc, L.get_row( m+i ) );
+                        L.set_row( n-cc, L.get_row( m+i ) );
                         L.set_row( m+i,0 );
                         D( m+i )=0; L( m+i,i )=1;
-                        L.set_submatrix( n+1-cc,m+i,i,i,0 );
+                        L.set_submatrix( n-cc,m+i-1,i,i,0 );
                         continue;
                 }
 
-                sum-=v( v.length()-1 )*v( v.length()-1 ); //
+                sum-=v(last_v)*v(last_v);
                 sum/=beta*beta;
                 sum++;
 
                 v/=beta;
-                v( v.length()-1 )=1;
-
-                pom=-2/sum;
-                for ( j=i;i>=0;i-- ) {
-                        w( j )=0.0;
-                        for ( ii=n-cc;ii<m+i;ii++ ) {
-                                w( j )+= v( ii )*L( ii,j );
-                        }
-                        w( j )*=pom;
-                }
-
-                for ( ii=n-cc;ii<m+i;ii++ ) {
-                        for ( jj=0;jj<i-1;jj++ ) {
-                                L( ii,jj )+= v( ii )*w( jj );
-                        }
-                }
-                for ( ii=n-cc;ii<m+i;ii++ ) {
+                v(last_v)=1;
+
+                pom=-2.0/sum;
+                // echo to venca   
+
+                for ( j=i;j>=0;j-- ) 
+                {
+                        double w_elem = 0;                      
+                        for ( ii=n-     cc;ii<=m+i+1;ii++ ) 
+                                w_elem+= v( ii-n+cc )*L( ii-1,j );                      
+                        w(j)=w_elem*pom;
+                }
+
+                for ( ii=n-cc-1;ii<=m+i;ii++ ) 
+                        for ( jj=0;jj<i;jj++ ) 
+                                L( ii,jj )+= v( ii-n+cc+1)*w( jj );
+
+                for ( ii=n-cc-1;ii<m+i;ii++ )
                         L( ii,i )= 0;
-                }
+
                 L( m+i,i )+=w( i );
-
                 D( m+i )=L( m+i,i )*L( m+i,i );
-                for ( ii=0;ii<i;ii++ ) {
-                        L( m+i,ii )/=L( m+i,i );
-                }
-        }
-        if ( i>0 ) {
-                for ( ii=0;ii<i-1;ii++ ) {
+
+                for ( ii=0;ii<=i;ii++ ) 
+                        L( m+i,ii )/=L( m+i,i );                
+        }
+
+        if ( i>=0 )
+                for ( ii=0;ii<i;ii++ ) 
+                {
                         jj = D.length()-1-n+ii;
+                        D(jj) = 0;
                         L.set_row(jj,0);
                         L(jj,jj)=1;
                 }
-        }
-
-        //cut-out L and D;
+
         L.del_rows(0,m-1);
         D.del(0,m-1);