Changeset 477 for library/bdm/math

Timestamp:

08/05/09 14:40:03 (15 years ago)

Author:

mido

Message:

panove, vite, jak jsem peclivej na upravu kodu.. snad se vam bude libit:) konfigurace je v souboru /system/astylerc

Location:

library/bdm/math

Files:

• chmat.cpp (modified) (1 diff)
• chmat.h (modified) (3 diffs)
• functions.h (modified) (2 diffs)
• psi.cpp (modified) (1 diff)
• square_mat.cpp (modified) (11 diffs)
• square_mat.h (modified) (4 diffs)

library/bdm/math/chmat.cpp

@@ -9 +9 @@
         mat Z;
         mat R;
-        mat V(1,v.length());
-        V.set_row(0,v*sqrt(w));
-        Z = concat_vertical ( Ch,V );
-        qr ( Z,R );
-        Ch = R ( 0, Ch.rows()-1, 0, Ch.cols()-1 );
+        mat V ( 1, v.length() );
+        V.set_row ( 0, v*sqrt ( w ) );
+        Z = concat_vertical ( Ch, V );
+        qr ( Z, R );
+        Ch = R ( 0, Ch.rows() - 1, 0, Ch.cols() - 1 );
 };
-mat chmat::to_mat() const {mat F=Ch.T() *Ch;return F;};
+mat chmat::to_mat() const {
+        mat F = Ch.T() * Ch;
+        return F;
+};
 void chmat::mult_sym ( const mat &C ) {
-        it_assert_debug(C.rows()==dim, "Wrong dimension of U");
-        if(!qr(Ch*C.T(), Ch)) {it_warning("QR unstable in chmat mult_sym");}
+        it_assert_debug ( C.rows() == dim, "Wrong dimension of U" );
+        if ( !qr ( Ch*C.T(), Ch ) ) {
+                it_warning ( "QR unstable in chmat mult_sym" );
+        }
 };
 void chmat::mult_sym ( const mat &C , chmat &U ) const {
-        it_assert_debug(C.rows()==U.dim, "Wrong dimension of U");
-        if(!qr(Ch*C.T(), U.Ch)) {it_warning("QR unstable in chmat mult_sym");}
+        it_assert_debug ( C.rows() == U.dim, "Wrong dimension of U" );
+        if ( !qr ( Ch*C.T(), U.Ch ) ) {
+                it_warning ( "QR unstable in chmat mult_sym" );
+        }
 };
 void chmat::mult_sym_t ( const mat &C ) {
-        it_assert_debug(C.cols()==dim, "Wrong dimension of U");
-        if(!qr(Ch*C, Ch)) {it_warning("QR unstable in chmat mult_sym");}
+        it_assert_debug ( C.cols() == dim, "Wrong dimension of U" );
+        if ( !qr ( Ch*C, Ch ) ) {
+                it_warning ( "QR unstable in chmat mult_sym" );
+        }
 };
 void chmat::mult_sym_t ( const mat &C, chmat &U ) const {
-        it_assert_debug(C.cols()==U.dim, "Wrong dimension of U");
-        if(!qr(Ch*C, U.Ch)) {it_warning("QR unstable in chmat mult_sym");}
+        it_assert_debug ( C.cols() == U.dim, "Wrong dimension of U" );
+        if ( !qr ( Ch*C, U.Ch ) ) {
+                it_warning ( "QR unstable in chmat mult_sym" );
+        }
 };
 double chmat::logdet() const {
-        double ldet=0.0; int i;
+        double ldet = 0.0;
+        int i;
         //sum of logs of (possibly negative!) diagonal entries
-        for ( i=0;i<Ch.rows();i++ ) {ldet+=log ( std::fabs ( Ch ( i,i ) ) );}
+        for ( i = 0; i < Ch.rows(); i++ ) {
+                ldet += log ( std::fabs ( Ch ( i, i ) ) );
+        }
         return 2*ldet; //compensate for Ch being sqrt()
 };
 //TODO can be done more efficiently using BLAS, see triangular matrices
-vec chmat::sqrt_mult ( const vec &v ) const {vec pom; pom = Ch*v; return pom;};
-double chmat::qform ( const vec &v ) const {vec pom; pom = Ch*v; return pom*pom;};
-double chmat::invqform ( const vec &v ) const {
-        vec pom(v.length());
-        forward_substitution(Ch.T(),v,pom);
+vec chmat::sqrt_mult ( const vec &v ) const {
+        vec pom;
+        pom = Ch * v;
+        return pom;
+};
+double chmat::qform ( const vec &v ) const {
+        vec pom;
+        pom = Ch * v;
         return pom*pom;
 };
-void chmat::clear() {Ch.clear();};
+double chmat::invqform ( const vec &v ) const {
+        vec pom ( v.length() );
+        forward_substitution ( Ch.T(), v, pom );
+        return pom*pom;
+};
+void chmat::clear() {
+        Ch.clear();
+};

library/bdm/math/chmat.h

@@ -40 +40 @@
         void clear();
         //! add another chmat \c A2 with weight \c w.
-        void add ( const chmat &A2, double w=1.0 ) {
-                it_assert_debug(dim==A2.dim,"Matrices of unequal dimension");
-                mat pre=concat_vertical(Ch,sqrt(w)*A2.Ch); 
-                mat post=zeros(pre.rows(),pre.cols()); 
-                if(!qr(pre,post)) {it_warning("Unstable QR in chmat add");}
-                Ch=post(0,dim-1,0,dim-1);
+        void add ( const chmat &A2, double w = 1.0 ) {
+                it_assert_debug ( dim == A2.dim, "Matrices of unequal dimension" );
+                mat pre = concat_vertical ( Ch, sqrt ( w ) * A2.Ch );
+                mat post = zeros ( pre.rows(), pre.cols() );
+                if ( !qr ( pre, post ) ) {
+                        it_warning ( "Unstable QR in chmat add" );
+                }
+                Ch = post ( 0, dim - 1, 0, dim - 1 );
         };
         //!Inversion in the same form, i.e. cholesky
-        void inv ( chmat &Inv ) const   { ( Inv.Ch ) = itpp::inv ( Ch ).T();}; //Fixme: can be more efficient
+        void inv ( chmat &Inv ) const   {
+                ( Inv.Ch ) = itpp::inv ( Ch ).T();
+        }; //Fixme: can be more efficient
         ;
 //      void inv ( mat &Inv );
@@ -55 +59 @@
         virtual ~chmat() {};
         //!
-        chmat ( ) : sqmat (),Ch ( ) {};
+        chmat ( ) : sqmat (), Ch ( ) {};
         //! Default constructor
-        chmat ( const int dim0 ) : sqmat ( dim0 ),Ch ( dim0,dim0 ) {};
+        chmat ( const int dim0 ) : sqmat ( dim0 ), Ch ( dim0, dim0 ) {};
         //! Default constructor
-        chmat ( const vec &v) : sqmat ( v.length() ),Ch ( diag(sqrt(v)) ) {};
+        chmat ( const vec &v ) : sqmat ( v.length() ), Ch ( diag ( sqrt ( v ) ) ) {};
         //! Copy constructor
-        chmat ( const chmat &Ch0) : sqmat ( Ch0.dim),Ch ( Ch0.dim,Ch0.dim ) {Ch=Ch0.Ch;};
+        chmat ( const chmat &Ch0 ) : sqmat ( Ch0.dim ), Ch ( Ch0.dim, Ch0.dim ) {
+                Ch = Ch0.Ch;
+        };
         //! Default constructor (m3k:cholform)
-        chmat ( const mat &M ) : sqmat ( M.rows() ),Ch ( M.rows(),M.cols() ) {
+        chmat ( const mat &M ) : sqmat ( M.rows() ), Ch ( M.rows(), M.cols() ) {
                 mat Q;
-                it_assert_debug ( M.rows() ==M.cols(),"chmat:: input matrix must be square!" );
-                Ch=chol ( M );
+                it_assert_debug ( M.rows() == M.cols(), "chmat:: input matrix must be square!" );
+                Ch = chol ( M );
         };
         //! Constructor
-        chmat ( const chmat &M, const ivec &perm ):sqmat(M.rows()){it_error("not implemneted");};
+        chmat ( const chmat &M, const ivec &perm ) : sqmat ( M.rows() ) {
+                it_error ( "not implemneted" );
+        };
         //! Access function
-        mat & _Ch() {return Ch;}
+        mat & _Ch() {
+                return Ch;
+        }
         //! Access function
-        const mat & _Ch() const {return Ch;}
+        const mat & _Ch() const {
+                return Ch;
+        }
         //! Access functions
-        void setD ( const vec &nD ) {Ch=diag ( sqrt(nD) );}
+        void setD ( const vec &nD ) {
+                Ch = diag ( sqrt ( nD ) );
+        }
         //! Access functions
         void setCh ( const vec &chQ ) {
-                it_assert_debug(chQ.length()==dim*dim,""); 
-                copy_vector(dim*dim, chQ._data(), Ch._data()); 
+                it_assert_debug ( chQ.length() == dim*dim, "" );
+                copy_vector ( dim*dim, chQ._data(), Ch._data() );
         }
         //! Access functions
-        void setD ( const vec &nD, int i ) {for ( int j=i;j<nD.length();j++ ) {Ch( j,j ) =sqrt(nD ( j-i ));}} //Fixme can be more general
+        void setD ( const vec &nD, int i ) {
+                for ( int j = i; j < nD.length(); j++ ) {
+                        Ch ( j, j ) = sqrt ( nD ( j - i ) );    //Fixme can be more general
+                }
+        }
 
         //! Operators
         chmat& operator += ( const chmat &A2 );
         chmat& operator -= ( const chmat &A2 );
-        chmat& operator * ( const double &d ){Ch*sqrt(d); return *this;};
-        chmat& operator = ( const chmat &A2 ){Ch=A2.Ch;dim=A2.dim;return *this;}
-        chmat& operator *= ( double x ){Ch*=sqrt(x); return *this;};
+        chmat& operator * ( const double &d ) {
+                Ch*sqrt ( d );
+                return *this;
+        };
+        chmat& operator = ( const chmat &A2 ) {
+                Ch = A2.Ch;
+                dim = A2.dim;
+                return *this;
+        }
+        chmat& operator *= ( double x ) {
+                Ch *= sqrt ( x );
+                return *this;
+        };
 };
 
@@ -95 +123 @@
 //////// Operations:
 //!mapping of add operation to operators
-inline chmat& chmat::operator += ( const chmat &A2 )  {this->add ( A2 );return *this;}
+inline chmat& chmat::operator += ( const chmat & A2 )  {
+        this->add ( A2 );
+        return *this;
+}
 //!mapping of negative add operation to operators
-inline chmat& chmat::operator -= ( const chmat &A2 )  {this->add ( A2,-1.0 );return *this;}
+inline chmat& chmat::operator -= ( const chmat & A2 )  {
+        this->add ( A2, -1.0 );
+        return *this;
+}
 
 #endif // CHMAT_H

library/bdm/math/functions.h

@@ -15 +15 @@
 #include "../base/bdmbase.h"
 
-namespace bdm{
+namespace bdm {
 
 //! class representing function \f$f(x) = a\f$, here \c rv is empty
-class constfn : public fnc
-{
-                //! value of the function
-                vec val;
+class constfn : public fnc {
+        //! value of the function
+        vec val;
 
-        public:
-                //vec eval() {return val;};
-                //! inherited
-                vec eval ( const vec &cond ) {return val;};
-                //!Default constructor
-                constfn ( const vec &val0 ) :fnc(), val ( val0 ) {dimy=val.length();};
+public:
+        //vec eval() {return val;};
+        //! inherited
+        vec eval ( const vec &cond ) {
+                return val;
+        };
+        //!Default constructor
+        constfn ( const vec &val0 ) : fnc(), val ( val0 ) {
+                dimy = val.length();
+        };
 };
 
 //! Class representing function \f$f(x) = Ax+B\f$
-class linfn: public fnc
-{
-                //! Identification of \f$x\f$
-                RV rv;
-                //! Matrix A
-                mat A;
-                //! vector B
-                vec B;
-        public :
-                vec eval (const vec &cond ) {it_assert_debug ( cond.length() ==A.cols(), "linfn::eval Wrong cond." );return A*cond+B;};
+class linfn: public fnc {
+        //! Identification of \f$x\f$
+        RV rv;
+        //! Matrix A
+        mat A;
+        //! vector B
+        vec B;
+public :
+        vec eval ( const vec &cond ) {
+                it_assert_debug ( cond.length() == A.cols(), "linfn::eval Wrong cond." );
+                return A*cond + B;
+        };
 
 //              linfn evalsome ( ivec &rvind );
-                //!default constructor
-                linfn ( ) : fnc(), A ( ),B () { };
-                //! Set values of \c A and \c B
-                void set_parameters ( const mat &A0 , const vec &B0 ) {A=A0; B=B0; dimy=A.rows();};
+        //!default constructor
+        linfn ( ) : fnc(), A ( ), B () { };
+        //! Set values of \c A and \c B
+        void set_parameters ( const mat &A0 , const vec &B0 ) {
+                A = A0;
+                B = B0;
+                dimy = A.rows();
+        };
 };
 
@@ -60 +69 @@
 2) It could be generalized into multivariate form, (which was original meaning of \c fnc ).
 */
-class diffbifn: public fnc
-{
-        protected:
-                //! Indentifier of the first rv.
-                RV rvx;
-                //! Indentifier of the second rv.
-                RV rvu;
-                //! cache for rvx.count()
-                int dimx;
-                //! cache for rvu.count()
-                int dimu;
-        public:
-                //! Evaluates \f$f(x0,u0)\f$ (VS: Do we really need common eval? )
-                vec eval ( const vec &cond )
-                {
-                        it_assert_debug ( cond.length() == ( dimx+dimu ), "linfn::eval Wrong cond." );
-                        return eval ( cond ( 0,dimx-1 ),cond ( dimx,dimx+dimu-1 ) );//-1 = end (in matlab)
-                };
+class diffbifn: public fnc {
+protected:
+        //! Indentifier of the first rv.
+        RV rvx;
+        //! Indentifier of the second rv.
+        RV rvu;
+        //! cache for rvx.count()
+        int dimx;
+        //! cache for rvu.count()
+        int dimu;
+public:
+        //! Evaluates \f$f(x0,u0)\f$ (VS: Do we really need common eval? )
+        vec eval ( const vec &cond ) {
+                it_assert_debug ( cond.length() == ( dimx + dimu ), "linfn::eval Wrong cond." );
+                return eval ( cond ( 0, dimx - 1 ), cond ( dimx, dimx + dimu - 1 ) );//-1 = end (in matlab)
+        };
 
-                //! Evaluates \f$f(x0,u0)\f$
-                virtual vec eval ( const vec &x0, const vec &u0 ) {return zeros ( dimy );};
-                //! Evaluates \f$A=\frac{d}{dx}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed. @param x0 numeric value of \f$x\f$, @param u0 numeric value of \f$u\f$ @param A a place where the result will be stored.
-                virtual void dfdx_cond ( const vec &x0, const vec &u0, mat &A , bool full=true ) {};
-                //! Evaluates \f$A=\frac{d}{du}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed.        @param x0 numeric value of \f$x\f$, @param u0 numeric value of \f$u\f$ @param A a place where the result will be stored.
-                virtual void dfdu_cond ( const vec &x0, const vec &u0, mat &A, bool full=true ) {};
-                //!Default constructor (dimy is not set!)
-                diffbifn () : fnc() {};
-                //! access function
-                int _dimx() const{return dimx;}
-                //! access function
-                int _dimu() const{return dimu;}
+        //! Evaluates \f$f(x0,u0)\f$
+        virtual vec eval ( const vec &x0, const vec &u0 ) {
+                return zeros ( dimy );
+        };
+        //! Evaluates \f$A=\frac{d}{dx}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed. @param x0 numeric value of \f$x\f$, @param u0 numeric value of \f$u\f$ @param A a place where the result will be stored.
+        virtual void dfdx_cond ( const vec &x0, const vec &u0, mat &A , bool full = true ) {};
+        //! Evaluates \f$A=\frac{d}{du}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed.        @param x0 numeric value of \f$x\f$, @param u0 numeric value of \f$u\f$ @param A a place where the result will be stored.
+        virtual void dfdu_cond ( const vec &x0, const vec &u0, mat &A, bool full = true ) {};
+        //!Default constructor (dimy is not set!)
+        diffbifn () : fnc() {};
+        //! access function
+        int _dimx() const {
+                return dimx;
+        }
+        //! access function
+        int _dimu() const {
+                return dimu;
+        }
 };
 
 //! Class representing function \f$f(x,u) = Ax+Bu\f$
 //TODO can be generalized into multilinear form!
-class bilinfn: public diffbifn
-{
-                mat A;
-                mat B;
-        public :
-                //!\name Constructors
-                //!@{
-                
-                bilinfn () : diffbifn () ,A() ,B()      {};
-                bilinfn (const mat A0, const mat B0) {set_parameters(A0,B0);};
-                //! Alternative constructor
-                void set_parameters(const mat A0, const mat B0){
-                        it_assert_debug(A0.rows()==B0.rows(),"");
-                        A=A0;B=B0;
-                        dimy=A.rows();
-                        dimx=A.cols();
-                        dimu=B.cols();
-                }
-                //!@}
-                
-                //!\name Mathematical operations
-                //!@{
-                inline vec eval ( const  vec &x0, const vec &u0 )
-                {
-                        it_assert_debug ( x0.length() ==dimx, "linfn::eval Wrong xcond." );
-                        it_assert_debug ( u0.length() ==dimu, "linfn::eval Wrong ucond." );
-                        return A*x0+B*u0;
-                }
+class bilinfn: public diffbifn {
+        mat A;
+        mat B;
+public :
+        //!\name Constructors
+        //!@{
 
-                void dfdx_cond ( const vec &x0, const vec &u0, mat &F, bool full )
-                {
-                        it_assert_debug ( ( F.cols() ==A.cols() ) & ( F.rows() ==A.rows() ),"Allocated F is not compatible." );
-                        if ( full ) F=A;        //else : nothing has changed no need to regenerate
-                }
-                //!
-                void dfdu_cond ( const vec &x0, const vec &u0, mat &F,  bool full=true )
-                {
-                        it_assert_debug ( ( F.cols() ==B.cols() ) & ( F.rows() ==B.rows() ),"Allocated F is not compatible." );
-                        if ( full ) F=B;        //else : nothing has changed no need to regenerate
-                }
-                //!@}
+        bilinfn () : diffbifn () , A() , B()    {};
+        bilinfn ( const mat A0, const mat B0 ) {
+                set_parameters ( A0, B0 );
+        };
+        //! Alternative constructor
+        void set_parameters ( const mat A0, const mat B0 ) {
+                it_assert_debug ( A0.rows() == B0.rows(), "" );
+                A = A0;
+                B = B0;
+                dimy = A.rows();
+                dimx = A.cols();
+                dimu = B.cols();
+        }
+        //!@}
+
+        //!\name Mathematical operations
+        //!@{
+        inline vec eval ( const  vec &x0, const vec &u0 ) {
+                it_assert_debug ( x0.length() == dimx, "linfn::eval Wrong xcond." );
+                it_assert_debug ( u0.length() == dimu, "linfn::eval Wrong ucond." );
+                return A*x0 + B*u0;
+        }
+
+        void dfdx_cond ( const vec &x0, const vec &u0, mat &F, bool full ) {
+                it_assert_debug ( ( F.cols() == A.cols() ) & ( F.rows() == A.rows() ), "Allocated F is not compatible." );
+                if ( full ) F = A;      //else : nothing has changed no need to regenerate
+        }
+        //!
+        void dfdu_cond ( const vec &x0, const vec &u0, mat &F,  bool full = true ) {
+                it_assert_debug ( ( F.cols() == B.cols() ) & ( F.rows() == B.rows() ), "Allocated F is not compatible." );
+                if ( full ) F = B;      //else : nothing has changed no need to regenerate
+        }
+        //!@}
 };
 

library/bdm/math/psi.cpp

@@ -12 +12 @@
 #define el 0.5772156649015329
 
-double psi(double x)
-{
-    double s,ps,xa,x2;
-    int n,k;
-    static double a[] = {
-        -0.8333333333333e-01,
-         0.83333333333333333e-02,
-        -0.39682539682539683e-02,
-         0.41666666666666667e-02,
-        -0.75757575757575758e-02,
-         0.21092796092796093e-01,
-        -0.83333333333333333e-01,
-         0.4432598039215686};
+double psi ( double x ) {
+        double s, ps, xa, x2;
+        int n, k;
+        static double a[] = {
+                -0.8333333333333e-01,
+                0.83333333333333333e-02,
+                -0.39682539682539683e-02,
+                0.41666666666666667e-02,
+                -0.75757575757575758e-02,
+                0.21092796092796093e-01,
+                -0.83333333333333333e-01,
+                0.4432598039215686
+        };
 
-    xa = fabs(x);
-    s = 0.0;
-    if ((x == (int)x) && (x <= 0.0)) {
-        ps = 1e308;
-        return ps;
-    }
-    if (xa == (int)xa) {
-        n = xa;
-        for (k=1;k<n;k++) {
-            s += 1.0/k;
-        }
-        ps =  s-el;
-    }
-    else if ((xa+0.5) == ((int)(xa+0.5))) {
-        n = xa-0.5;
-        for (k=1;k<=n;k++) {
-            s += 1.0/(2.0*k-1.0);
-        }
-        ps = 2.0*s-el-1.386294361119891;
-    }
-    else {
-        if (xa < 10.0) {
-            n = 10-(int)xa;
-            for (k=0;k<n;k++) {
-                s += 1.0/(xa+k);
-            }
-            xa += n;
-        }
-        x2 = 1.0/(xa*xa);
-        ps = log(xa)-0.5/xa+x2*(((((((a[7]*x2+a[6])*x2+a[5])*x2+
-            a[4])*x2+a[3])*x2+a[2])*x2+a[1])*x2+a[0]);
-        ps -= s;
-    }
-    if (x < 0.0)
-        ps = ps - M_PI*cos(M_PI*x)/sin(M_PI*x)-1.0/x;
-        return ps;
+        xa = fabs ( x );
+        s = 0.0;
+        if ( ( x == ( int ) x ) && ( x <= 0.0 ) ) {
+                ps = 1e308;
+                return ps;
+        }
+        if ( xa == ( int ) xa ) {
+                n = xa;
+                for ( k = 1; k < n; k++ ) {
+                        s += 1.0 / k;
+                }
+                ps =  s - el;
+        } else if ( ( xa + 0.5 ) == ( ( int ) ( xa + 0.5 ) ) ) {
+                n = xa - 0.5;
+                for ( k = 1; k <= n; k++ ) {
+                        s += 1.0 / ( 2.0 * k - 1.0 );
+                }
+                ps = 2.0 * s - el - 1.386294361119891;
+        } else {
+                if ( xa < 10.0 ) {
+                        n = 10 - ( int ) xa;
+                        for ( k = 0; k < n; k++ ) {
+                                s += 1.0 / ( xa + k );
+                        }
+                        xa += n;
+                }
+                x2 = 1.0 / ( xa * xa );
+                ps = log ( xa ) - 0.5 / xa + x2 * ( ( ( ( ( ( ( a[7] * x2 + a[6] ) * x2 + a[5] ) * x2 +
+                                                        a[4] ) * x2 + a[3] ) * x2 + a[2] ) * x2 + a[1] ) * x2 + a[0] );
+                ps -= s;
+        }
+        if ( x < 0.0 )
+                ps = ps - M_PI * cos ( M_PI * x ) / sin ( M_PI * x ) - 1.0 / x;
+        return ps;
 }

library/bdm/math/square_mat.cpp

@@ -6 +6 @@
 using std::endl;
 
-void fsqmat::opupdt ( const vec &v, double w ) {M+=outer_product ( v,v*w );};
-mat fsqmat::to_mat() const {return M;};
-void fsqmat::mult_sym ( const mat &C) {M=C *M*C.T();};
-void fsqmat::mult_sym_t ( const mat &C) {M=C.T() *M*C;};
-void fsqmat::mult_sym ( const mat &C, fsqmat &U) const { U.M = ( C *(M*C.T()) );};
-void fsqmat::mult_sym_t ( const mat &C, fsqmat &U) const { U.M = ( C.T() *(M*C) );};
-void fsqmat::inv ( fsqmat &Inv ) const {mat IM = itpp::inv ( M ); Inv=IM;};
-void fsqmat::clear() {M.clear();};
-fsqmat::fsqmat ( const mat &M0 ) : sqmat(M0.cols())
-{
-        it_assert_debug ( ( M0.cols() ==M0.rows() ),"M0 must be square" );
-        M=M0;
+void fsqmat::opupdt ( const vec &v, double w ) {
+        M += outer_product ( v, v * w );
+};
+mat fsqmat::to_mat() const {
+        return M;
+};
+void fsqmat::mult_sym ( const mat &C ) {
+        M = C * M * C.T();
+};
+void fsqmat::mult_sym_t ( const mat &C ) {
+        M = C.T() * M * C;
+};
+void fsqmat::mult_sym ( const mat &C, fsqmat &U ) const {
+        U.M = ( C * ( M * C.T() ) );
+};
+void fsqmat::mult_sym_t ( const mat &C, fsqmat &U ) const {
+        U.M = ( C.T() * ( M * C ) );
+};
+void fsqmat::inv ( fsqmat &Inv ) const {
+        mat IM = itpp::inv ( M );
+        Inv = IM;
+};
+void fsqmat::clear() {
+        M.clear();
+};
+fsqmat::fsqmat ( const mat &M0 ) : sqmat ( M0.cols() ) {
+        it_assert_debug ( ( M0.cols() == M0.rows() ), "M0 must be square" );
+        M = M0;
 };
 
 //fsqmat::fsqmat() {};
 
-fsqmat::fsqmat(const int dim0): sqmat(dim0), M(dim0,dim0) {};
+fsqmat::fsqmat ( const int dim0 ) : sqmat ( dim0 ), M ( dim0, dim0 ) {};
 
 std::ostream &operator<< ( std::ostream &os, const fsqmat &ld ) {
@@ -30 +46 @@
 
 
-ldmat::ldmat( const mat &exL, const vec &exD ) : sqmat(exD.length()) {
+ldmat::ldmat ( const mat &exL, const vec &exD ) : sqmat ( exD.length() ) {
         D = exD;
         L = exL;
 }
 
-ldmat::ldmat() :sqmat(0) {}
-
-ldmat::ldmat(const int dim0): sqmat(dim0), D(dim0),L(dim0,dim0) {}
-
-ldmat::ldmat(const vec D0):sqmat(D0.length()) {
+ldmat::ldmat() : sqmat ( 0 ) {}
+
+ldmat::ldmat ( const int dim0 ) : sqmat ( dim0 ), D ( dim0 ), L ( dim0, dim0 ) {}
+
+ldmat::ldmat ( const vec D0 ) : sqmat ( D0.length() ) {
         D = D0;
-        L = eye(dim);
-}
-
-ldmat::ldmat( const mat &V ):sqmat(V.cols()) {
+        L = eye ( dim );
+}
+
+ldmat::ldmat ( const mat &V ) : sqmat ( V.cols() ) {
 //TODO check if correct!! Based on heuristic observation of lu()
 
-        it_assert_debug( dim == V.rows(),"ldmat::ldmat matrix V is not square!" );
-                
+        it_assert_debug ( dim == V.rows(), "ldmat::ldmat matrix V is not square!" );
+
         // L and D will be allocated by ldform()
-        
+
         //Chol is unstable
-        this->ldform(chol(V),ones(dim));
+        this->ldform ( chol ( V ), ones ( dim ) );
 //      this->ldform(ul(V),ones(dim));
 }
 
-void ldmat::opupdt( const vec &v,  double w ) {
+void ldmat::opupdt ( const vec &v,  double w ) {
         int dim = D.length();
         double kr;
@@ -65 +81 @@
         double *rraw = r._data();
 
-        it_assert_debug( v.length() == dim, "LD::ldupdt vector v is not compatible with this ld." );
+        it_assert_debug ( v.length() == dim, "LD::ldupdt vector v is not compatible with this ld." );
 
         for ( int i = dim - 1; i >= 0; i-- ) {
-                dydr( rraw, Lraw + i, &w, Draw + i, rraw + i, 0, i, &kr, 1, dim );
+                dydr ( rraw, Lraw + i, &w, Draw + i, rraw + i, 0, i, &kr, 1, dim );
         }
 }
@@ -80 +96 @@
 mat ldmat::to_mat() const {
         int dim = D.length();
-        mat V( dim, dim );
+        mat V ( dim, dim );
         double sum;
         int r, c, cc;
 
-        for ( r = 0;r < dim;r++ ) { //row cycle
-                for ( c = r;c < dim;c++ ) {
+        for ( r = 0; r < dim; r++ ) { //row cycle
+                for ( c = r; c < dim; c++ ) {
                         //column cycle, using symmetricity => c=r!
                         sum = 0.0;
-                        for ( cc = c;cc < dim;cc++ ) { //cycle over the remaining part of the vector
-                                sum += L( cc, r ) * D( cc ) * L( cc, c );
+                        for ( cc = c; cc < dim; cc++ ) { //cycle over the remaining part of the vector
+                                sum += L ( cc, r ) * D ( cc ) * L ( cc, c );
                                 //here L(cc,r) = L(r,cc)';
                         }
-                        V( r, c ) = sum;
+                        V ( r, c ) = sum;
                         // symmetricity
-                        if ( r != c ) {V( c, r ) = sum;};
-                }
-        }
-        mat V2 = L.transpose()*diag( D )*L;
+                        if ( r != c ) {
+                                V ( c, r ) = sum;
+                        };
+                }
+        }
+        mat V2 = L.transpose() * diag ( D ) * L;
         return V2;
 }
 
 
-void ldmat::add( const ldmat &ld2, double w ) {
+void ldmat::add ( const ldmat &ld2, double w ) {
         int dim = D.length();
 
-        it_assert_debug( ld2.D.length() == dim, "LD.add() incompatible sizes of LDs;" );
+        it_assert_debug ( ld2.D.length() == dim, "LD.add() incompatible sizes of LDs;" );
 
         //Fixme can be done more efficiently either via dydr or ldform
         for ( int r = 0; r < dim; r++ ) {
                 // Add columns of ld2.L' (i.e. rows of ld2.L) as dyads weighted by ld2.D
-                this->opupdt( ld2.L.get_row( r ), w*ld2.D( r ) );
-        }
-}
-
-void ldmat::clear(){L.clear(); for ( int i=0;i<L.cols();i++ ){L( i,i )=1;}; D.clear();}
-
-void ldmat::inv( ldmat &Inv ) const {
+                this->opupdt ( ld2.L.get_row ( r ), w*ld2.D ( r ) );
+        }
+}
+
+void ldmat::clear() {
+        L.clear();
+        for ( int i = 0; i < L.cols(); i++ ) {
+                L ( i, i ) = 1;
+        };
+        D.clear();
+}
+
+void ldmat::inv ( ldmat &Inv ) const {
         Inv.clear();   //Inv = zero in LD
-        mat U = ltuinv( L );
-
-        Inv.ldform( U.transpose(), 1.0 / D );
-}
-
-void ldmat::mult_sym( const mat &C) {
-        mat A = L*C.T();
-        this->ldform(A,D);
-}
-
-void ldmat::mult_sym_t( const mat &C) {
-        mat A = L*C;
-        this->ldform(A,D);
-}
-
-void ldmat::mult_sym( const mat &C, ldmat &U) const {
-        mat A=L*C.T(); //could be done more efficiently using BLAS
-        U.ldform(A,D);
-}
-
-void ldmat::mult_sym_t( const mat &C, ldmat &U) const {
-        mat A=L*C;
-/*      vec nD=zeros(U.rows());
-        nD.replace_mid(0, D); //I case that D < nD*/
-        U.ldform(A,D);
+        mat U = ltuinv ( L );
+
+        Inv.ldform ( U.transpose(), 1.0 / D );
+}
+
+void ldmat::mult_sym ( const mat &C ) {
+        mat A = L * C.T();
+        this->ldform ( A, D );
+}
+
+void ldmat::mult_sym_t ( const mat &C ) {
+        mat A = L * C;
+        this->ldform ( A, D );
+}
+
+void ldmat::mult_sym ( const mat &C, ldmat &U ) const {
+        mat A = L * C.T(); //could be done more efficiently using BLAS
+        U.ldform ( A, D );
+}
+
+void ldmat::mult_sym_t ( const mat &C, ldmat &U ) const {
+        mat A = L * C;
+        /*      vec nD=zeros(U.rows());
+                nD.replace_mid(0, D); //I case that D < nD*/
+        U.ldform ( A, D );
 }
 
@@ -150 +174 @@
         int i;
 // sum logarithms of diagobal elements
-        for ( i=0; i<D.length(); i++ ){ldet+=log( D( i ) );};
+        for ( i = 0; i < D.length(); i++ ) {
+                ldet += log ( D ( i ) );
+        };
         return ldet;
 }
 
-double ldmat::qform( const vec &v ) const {
+double ldmat::qform ( const vec &v ) const {
         double x = 0.0, sum;
-        int i,j;
-
-        for ( i=0; i<D.length(); i++ ) { //rows of L
+        int i, j;
+
+        for ( i = 0; i < D.length(); i++ ) { //rows of L
                 sum = 0.0;
-                for ( j=0; j<=i; j++ ){sum+=L( i,j )*v( j );}
-                x +=D( i )*sum*sum;
+                for ( j = 0; j <= i; j++ ) {
+                        sum += L ( i, j ) * v ( j );
+                }
+                x += D ( i ) * sum * sum;
         };
         return x;
 }
 
-double ldmat::invqform( const vec &v ) const {
+double ldmat::invqform ( const vec &v ) const {
         double x = 0.0;
         int i;
-        vec pom(v.length());
-        
-        backward_substitution(L.T(),v,pom);
-        
-        for ( i=0; i<D.length(); i++ ) { //rows of L
-                x +=pom(i)*pom(i)/D(i);
+        vec pom ( v.length() );
+
+        backward_substitution ( L.T(), v, pom );
+
+        for ( i = 0; i < D.length(); i++ ) { //rows of L
+                x += pom ( i ) * pom ( i ) / D ( i );
         };
         return x;
@@ -180 +208 @@
 
 ldmat& ldmat::operator *= ( double x ) {
-        D*=x;
+        D *= x;
         return *this;
 }
 
-vec ldmat::sqrt_mult( const vec &x ) const {
-        int i,j;
-        vec res( dim );
+vec ldmat::sqrt_mult ( const vec &x ) const {
+        int i, j;
+        vec res ( dim );
         //double sum;
-        for ( i=0;i<dim;i++ ) {//for each element of result
-                res( i ) = 0.0;
-                for ( j=i;j<dim;j++ ) {//sum D(j)*L(:,i).*x
-                        res( i ) += sqrt( D( j ) )*L( j,i )*x( j );
+        for ( i = 0; i < dim; i++ ) {//for each element of result
+                res ( i ) = 0.0;
+                for ( j = i; j < dim; j++ ) {//sum D(j)*L(:,i).*x
+                        res ( i ) += sqrt ( D ( j ) ) * L ( j, i ) * x ( j );
                 }
         }
@@ -198 +226 @@
 }
 
-void ldmat::ldform(const mat &A,const vec &D0 )
-{
+void ldmat::ldform ( const mat &A, const vec &D0 ) {
         int m = A.rows();
         int n = A.cols();
-        int mn = (m<n) ? m :n ;
+        int mn = ( m < n ) ? m : n ;
 
 //      it_assert_debug( A.cols()==dim,"ldmat::ldform A is not compatible" );
-        it_assert_debug( D0.length()==A.rows(),"ldmat::ldform Vector D must have the length as row count of A" );
-
-        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+m );
-   
+        it_assert_debug ( D0.length() == A.rows(), "ldmat::ldform Vector D must have the length as row count of A" );
+
+        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 + m );
+
         double sum, beta, pom;
 
-        int cc=0;
-        int i=n; // indexovani o 1 niz, nez v matlabu
-        int ii,j,jj;
-        while ( (i>n-mn-cc) && (i>0) ) 
-        {
+        int cc = 0;
+        int i = n; // indexovani o 1 niz, nez v matlabu
+        int ii, j, jj;
+        while ( ( i > n - mn - cc ) && ( i > 0 ) ) {
                 i--;
                 sum = 0.0;
 
-                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-n+cc+1 )=L( ii,i ); //assign v
-                }
-
-                if ( L( m+i,i )==0 ) 
-                        beta = sqrt( sum );             
+                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 - 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 )
-                {
+                        beta = L ( m + i, i ) + sign ( L ( m + i, i ) ) * sqrt ( sum );
+
+                if ( std::fabs ( beta ) < eps ) {
                         cc++;
-                        L.set_row( n-cc, L.get_row( m+i ) );
-                        L.set_row( m+i,zeros(L.cols()) );
-                        D( m+i )=0; L( m+i,i )=1;
-                        L.set_submatrix( n-cc,m+i-1,i,i,0 );
+                        L.set_row ( n - cc, L.get_row ( m + i ) );
+                        L.set_row ( m + i, zeros ( L.cols() ) );
+                        D ( m + i ) = 0;
+                        L ( m + i, i ) = 1;
+                        L.set_submatrix ( n - cc, m + i - 1, i, i, 0 );
                         continue;
                 }
 
-                sum-=v(last_v)*v(last_v);
-                sum/=beta*beta;
+                sum -= v ( last_v ) * v ( last_v );
+                sum /= beta * beta;
                 sum++;
 
-                v/=beta;
-                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;ii++ ) 
-                {
-                        jj = D.length()-1-n+ii;
-                        D(jj) = 0;
-                        L.set_row(jj,zeros(L.cols())); //TODO: set_row accepts Num_T
-                        L(jj,jj)=1;
-                }
-
-        L.del_rows(0,m-1);
-        D.del(0,m-1);
-        
+                v /= beta;
+                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; ii++ ) {
+                        jj = D.length() - 1 - n + ii;
+                        D ( jj ) = 0;
+                        L.set_row ( jj, zeros ( L.cols() ) ); //TODO: set_row accepts Num_T
+                        L ( jj, jj ) = 1;
+                }
+
+        L.del_rows ( 0, m - 1 );
+        D.del ( 0, m - 1 );
+
         dim = L.rows();
 }
@@ -296 +319 @@
 //////// Auxiliary Functions
 
-mat ltuinv( const mat &L ) {
+mat ltuinv ( const mat &L ) {
         int dim = L.cols();
-        mat Il = eye( dim );
+        mat Il = eye ( dim );
         int i, j, k, m;
         double s;
 
 //Fixme blind transcription of ltuinv.m
-        for ( k = 1; k < ( dim );k++ ) {
-                for ( i = 0; i < ( dim - k );i++ ) {
+        for ( k = 1; k < ( dim ); k++ ) {
+                for ( i = 0; i < ( dim - k ); i++ ) {
                         j = i + k; //change in .m 1+1=2, here 0+0+1=1
-                        s = L( j, i );
+                        s = L ( j, i );
                         for ( m = i + 1; m < ( j ); m++ ) {
-                                s += L( m, i ) * Il( j, m );
+                                s += L ( m, i ) * Il ( j, m );
                         }
-                        Il( j, i ) = -s;
+                        Il ( j, i ) = -s;
                 }
         }
@@ -317 +340 @@
 }
 
-void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx )
+void dydr ( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx )
 /********************************************************************
 
@@ -353 +376 @@
         double threshold = 1e-4;
 
-        if ( fabs( *Dr ) < mzero ) *Dr = 0;
+        if ( fabs ( *Dr ) < mzero ) *Dr = 0;
         r0 = *R;
         *R = 0.0;
@@ -366 +389 @@
                 *kr = 0.0;
                 if ( *Df < -threshold ) {
-                        it_warning( "Problem in dydr: subraction of dyad results in negative definitness. Likely mistake in calling function." );
+                        it_warning ( "Problem in dydr: subraction of dyad results in negative definitness. Likely mistake in calling function." );
                 }
                 *Df = 0.0;

library/bdm/math/square_mat.h

@@ -20 +20 @@
 
 //! Auxiliary function dydr; dyadic reduction
-void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx );
+void dydr ( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx );
 
 //! Auxiliary function ltuinv; inversion of a triangular matrix;
 //TODO can be done via: dtrtri.f from lapack
-mat ltuinv( const mat &L );
+mat ltuinv ( const mat &L );
 
 /*! \brief Abstract class for representation of double symmetric matrices in square-root form.
@@ -30 +30 @@
 All operations defined on this class should be optimized for the chosen decomposition.
 */
-class sqmat
-{
-        public:
-                /*!
-                 * Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$.
-                 * @param  v Vector forming the outer product to be added
-                 * @param  w weight of updating; can be negative
-
-                 BLAS-2b operation.
-                 */
-                virtual void opupdt ( const vec &v, double w ) =0;
-
-                /*! \brief Conversion to full matrix.
-                */
-
-                virtual mat to_mat() const =0;
-
-                /*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$
-                @param C multiplying matrix,
-                */
-                virtual void mult_sym ( const mat &C ) =0;
-                
-                /*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'*V*C\f$
-                @param C multiplying matrix,
-                */
-                virtual void mult_sym_t ( const mat &C ) =0;
-
-
-                /*!
-                \brief Logarithm of a determinant.
-
-                */
-                virtual double logdet() const =0;
-
-                /*!
-                \brief Multiplies square root of \f$V\f$ by vector \f$x\f$.
-
-                Used e.g. in generating normal samples.
-                */
-                virtual vec sqrt_mult (const vec &v ) const =0;
-
-                /*!
-                \brief Evaluates quadratic form \f$x= v'*V*v\f$;
-
-                */
-                virtual double qform (const vec &v ) const =0;
-
-                /*!
-                \brief Evaluates quadratic form \f$x= v'*inv(V)*v\f$;
-
-                */
-                virtual double invqform (const vec &v ) const =0;
+class sqmat {
+public:
+        /*!
+         * Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$.
+         * @param  v Vector forming the outer product to be added
+         * @param  w weight of updating; can be negative
+
+         BLAS-2b operation.
+         */
+        virtual void opupdt ( const vec &v, double w ) = 0;
+
+        /*! \brief Conversion to full matrix.
+        */
+
+        virtual mat to_mat() const = 0;
+
+        /*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$
+        @param C multiplying matrix,
+        */
+        virtual void mult_sym ( const mat &C ) = 0;
+
+        /*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'*V*C\f$
+        @param C multiplying matrix,
+        */
+        virtual void mult_sym_t ( const mat &C ) = 0;
+
+
+        /*!
+        \brief Logarithm of a determinant.
+
+        */
+        virtual double logdet() const = 0;
+
+        /*!
+        \brief Multiplies square root of \f$V\f$ by vector \f$x\f$.
+
+        Used e.g. in generating normal samples.
+        */
+        virtual vec sqrt_mult ( const vec &v ) const = 0;
+
+        /*!
+        \brief Evaluates quadratic form \f$x= v'*V*v\f$;
+
+        */
+        virtual double qform ( const vec &v ) const = 0;
+
+        /*!
+        \brief Evaluates quadratic form \f$x= v'*inv(V)*v\f$;
+
+        */
+        virtual double invqform ( const vec &v ) const = 0;
 
 //      //! easy version of the
 //      sqmat inv();
 
-                //! Clearing matrix so that it corresponds to zeros.
-                virtual void clear() =0;
-
-                //! Reimplementing common functions of mat: cols().
-                int cols() const {return dim;};
-
-                //! Reimplementing common functions of mat: rows().
-                int rows() const {return dim;};
-
-                //! Destructor for future use;
-                virtual ~sqmat(){};
-                //! Default constructor
-                sqmat(const int dim0): dim(dim0){};
-                //! Default constructor
-                sqmat(): dim(0){};
-        protected:
-                //! dimension of the square matrix
-                int dim;
+        //! Clearing matrix so that it corresponds to zeros.
+        virtual void clear() = 0;
+
+        //! Reimplementing common functions of mat: cols().
+        int cols() const {
+                return dim;
+        };
+
+        //! Reimplementing common functions of mat: rows().
+        int rows() const {
+                return dim;
+        };
+
+        //! Destructor for future use;
+        virtual ~sqmat() {};
+        //! Default constructor
+        sqmat ( const int dim0 ) : dim ( dim0 ) {};
+        //! Default constructor
+        sqmat() : dim ( 0 ) {};
+protected:
+        //! dimension of the square matrix
+        int dim;
 };
 
@@ -111 +114 @@
 This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices;
 */
-class fsqmat: public sqmat
-{
-        protected:
-                //! Full matrix on which the operations are performed
-                mat M;
-        public:
-                void opupdt ( const vec &v, double w );
-                mat to_mat() const;
-                void mult_sym ( const mat &C);
-                void mult_sym_t ( const mat &C);
-                //! store result of \c mult_sym in external matrix \f$U\f$
-                void mult_sym ( const mat &C, fsqmat &U) const;
-                //! store result of \c mult_sym_t in external matrix \f$U\f$
-                void mult_sym_t ( const mat &C, fsqmat &U) const;
-                void clear();
-
-                //! Default initialization
-                fsqmat(){}; // mat will be initialized OK
-                //! Default initialization with proper size
-                fsqmat(const int dim0); // mat will be initialized OK
-                //! Constructor
-                fsqmat ( const mat &M );
-                //! Constructor
-                fsqmat ( const fsqmat &M, const ivec &perm ):sqmat(M.rows()){it_error("not implemneted");};
-                //! Constructor
-                fsqmat ( const vec &d ):sqmat(d.length()){M=diag(d);};
-
-                //! Destructor for future use;
-                virtual ~fsqmat(){};
-
-
-                /*! \brief Matrix inversion preserving the chosen form.
-
-                @param Inv a space where the inverse is stored.
-
-                */
-                void inv ( fsqmat &Inv ) const;
-
-                double logdet() const {return log ( det ( M ) );};
-                double qform (const  vec &v ) const {return ( v* ( M*v ) );};
-                double invqform (const  vec &v ) const {return ( v* ( itpp::inv(M)*v ) );};
-                vec sqrt_mult (const vec &v ) const {mat Ch=chol(M); return Ch*v;};
-
-                //! Add another matrix in fsq form with weight w
-                void add ( const fsqmat &fsq2, double w=1.0 ){M+=fsq2.M;};
-
-                //! Access functions
-                void setD (const vec &nD){M=diag(nD);}
-                //! Access functions
-                vec getD (){return diag(M);}
-                //! Access functions
-                void setD (const vec &nD, int i){for(int j=i;j<nD.length();j++){M(j,j)=nD(j-i);}} //Fixme can be more general
-
-
-                //! add another fsqmat matrix
-                fsqmat& operator += ( const fsqmat &A ) {M+=A.M;return *this;};
-                //! subtrack another fsqmat matrix
-                fsqmat& operator -= ( const fsqmat &A ) {M-=A.M;return *this;};
-                //! multiply by a scalar
-                fsqmat& operator *= ( double x ) {M*=x;return *this;};
+class fsqmat: public sqmat {
+protected:
+        //! Full matrix on which the operations are performed
+        mat M;
+public:
+        void opupdt ( const vec &v, double w );
+        mat to_mat() const;
+        void mult_sym ( const mat &C );
+        void mult_sym_t ( const mat &C );
+        //! store result of \c mult_sym in external matrix \f$U\f$
+        void mult_sym ( const mat &C, fsqmat &U ) const;
+        //! store result of \c mult_sym_t in external matrix \f$U\f$
+        void mult_sym_t ( const mat &C, fsqmat &U ) const;
+        void clear();
+
+        //! Default initialization
+        fsqmat() {}; // mat will be initialized OK
+        //! Default initialization with proper size
+        fsqmat ( const int dim0 ); // mat will be initialized OK
+        //! Constructor
+        fsqmat ( const mat &M );
+        //! Constructor
+        fsqmat ( const fsqmat &M, const ivec &perm ) : sqmat ( M.rows() ) {
+                it_error ( "not implemneted" );
+        };
+        //! Constructor
+        fsqmat ( const vec &d ) : sqmat ( d.length() ) {
+                M = diag ( d );
+        };
+
+        //! Destructor for future use;
+        virtual ~fsqmat() {};
+
+
+        /*! \brief Matrix inversion preserving the chosen form.
+
+        @param Inv a space where the inverse is stored.
+
+        */
+        void inv ( fsqmat &Inv ) const;
+
+        double logdet() const {
+                return log ( det ( M ) );
+        };
+        double qform ( const  vec &v ) const {
+                return ( v* ( M*v ) );
+        };
+        double invqform ( const  vec &v ) const {
+                return ( v* ( itpp::inv ( M ) *v ) );
+        };
+        vec sqrt_mult ( const vec &v ) const {
+                mat Ch = chol ( M );
+                return Ch*v;
+        };
+
+        //! Add another matrix in fsq form with weight w
+        void add ( const fsqmat &fsq2, double w = 1.0 ) {
+                M += fsq2.M;
+        };
+
+        //! Access functions
+        void setD ( const vec &nD ) {
+                M = diag ( nD );
+        }
+        //! Access functions
+        vec getD () {
+                return diag ( M );
+        }
+        //! Access functions
+        void setD ( const vec &nD, int i ) {
+                for ( int j = i; j < nD.length(); j++ ) {
+                        M ( j, j ) = nD ( j - i );    //Fixme can be more general
+                }
+        }
+
+
+        //! add another fsqmat matrix
+        fsqmat& operator += ( const fsqmat &A ) {
+                M += A.M;
+                return *this;
+        };
+        //! subtrack another fsqmat matrix
+        fsqmat& operator -= ( const fsqmat &A ) {
+                M -= A.M;
+                return *this;
+        };
+        //! multiply by a scalar
+        fsqmat& operator *= ( double x ) {
+                M *= x;
+                return *this;
+        };
 //              fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;};
-                //! print full matrix
-                friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq );
+        //! print full matrix
+        friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq );
 
 };
 
-/*! \brief Matrix stored in LD form, (commonly known as UD) 
+/*! \brief Matrix stored in LD form, (commonly known as UD)
 
 Matrix is decomposed as follows: \f[M = L'DL\f] where only \f$L\f$ and \f$D\f$ matrices are stored.
 All inplace operations modifies only these and the need to compose and decompose the matrix is avoided.
 */
-class ldmat: public sqmat
-{
-        public:
-                //! Construct by copy of L and D.
-                ldmat ( const mat &L, const vec &D );
-                //! Construct by decomposition of full matrix V.
-                ldmat (const mat &V );
-                //! Construct by restructuring of V0 accordint to permutation vector perm.
-                ldmat (const ldmat &V0, const ivec &perm):sqmat(V0.rows()){     ldform(V0.L.get_cols(perm), V0.D);};
-                //! Construct diagonal matrix with diagonal D0
-                ldmat ( vec D0 );
-                //!Default constructor
-                ldmat ();
-                //! Default initialization with proper size
-                ldmat(const int dim0);
-                
-                //! Destructor for future use;
-                virtual ~ldmat(){};
-
-                // Reimplementation of compulsory operatios
-
-                void opupdt ( const vec &v, double w );
-                mat to_mat() const;
-                void mult_sym ( const mat &C);
-                void mult_sym_t ( const mat &C);
-                //! Add another matrix in LD form with weight w
-                void add ( const ldmat &ld2, double w=1.0 );
-                double logdet() const;
-                double qform (const vec &v ) const;
-                double invqform (const vec &v ) const;
-                void clear();
-                vec sqrt_mult ( const vec &v ) const;
-
-                
-                /*! \brief Matrix inversion preserving the chosen form.
-                @param Inv a space where the inverse is stored.
-                */
-                void inv ( ldmat &Inv ) const;
-
-                /*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class.
-                @param C matrix to multiply with
-                @param U a space where the inverse is stored.
-                */
-                void mult_sym ( const mat &C, ldmat &U) const;
-
-                /*! \brief Symmetric multiplication of \f$U\f$ by a transpose of a general matrix \f$C\f$, result of which is stored in the current class.
-                @param C matrix to multiply with
-                @param U a space where the inverse is stored.
-                */
-                void mult_sym_t ( const mat &C, ldmat &U) const;
-
-
-                /*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$
-
-                The new decomposition fullfills: \f$A'*diag(D)*A = self.L'*diag(self.D)*self.L\f$
-                @param A general matrix
-                @param D0 general vector
-                */
-                void ldform (const mat &A,const vec &D0 );
-
-                //! Access functions
-                void setD (const vec &nD){D=nD;}
-                //! Access functions
-                void setD (const vec &nD, int i){D.replace_mid(i,nD);} //Fixme can be more general
-                //! Access functions
-                void setL (const vec &nL){L=nL;}
-
-                //! Access functions
-                const vec& _D() const {return D;}
-                //! Access functions
-                const mat& _L() const {return L;}
-
-                //! add another ldmat matrix
-                ldmat& operator += ( const ldmat &ldA );
-                //! subtract another ldmat matrix
-                ldmat& operator -= ( const ldmat &ldA );
-                //! multiply by a scalar
-                ldmat& operator *= ( double x );
-
-                //! print both \c L and \c D
-                friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq );
-
-        protected:
-                //! Positive vector \f$D\f$
-                vec D;
-                //! Lower-triangular matrix \f$L\f$
-                mat L;
+class ldmat: public sqmat {
+public:
+        //! Construct by copy of L and D.
+        ldmat ( const mat &L, const vec &D );
+        //! Construct by decomposition of full matrix V.
+        ldmat ( const mat &V );
+        //! Construct by restructuring of V0 accordint to permutation vector perm.
+        ldmat ( const ldmat &V0, const ivec &perm ) : sqmat ( V0.rows() ) {
+                ldform ( V0.L.get_cols ( perm ), V0.D );
+        };
+        //! Construct diagonal matrix with diagonal D0
+        ldmat ( vec D0 );
+        //!Default constructor
+        ldmat ();
+        //! Default initialization with proper size
+        ldmat ( const int dim0 );
+
+        //! Destructor for future use;
+        virtual ~ldmat() {};
+
+        // Reimplementation of compulsory operatios
+
+        void opupdt ( const vec &v, double w );
+        mat to_mat() const;
+        void mult_sym ( const mat &C );
+        void mult_sym_t ( const mat &C );
+        //! Add another matrix in LD form with weight w
+        void add ( const ldmat &ld2, double w = 1.0 );
+        double logdet() const;
+        double qform ( const vec &v ) const;
+        double invqform ( const vec &v ) const;
+        void clear();
+        vec sqrt_mult ( const vec &v ) const;
+
+
+        /*! \brief Matrix inversion preserving the chosen form.
+        @param Inv a space where the inverse is stored.
+        */
+        void inv ( ldmat &Inv ) const;
+
+        /*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class.
+        @param C matrix to multiply with
+        @param U a space where the inverse is stored.
+        */
+        void mult_sym ( const mat &C, ldmat &U ) const;
+
+        /*! \brief Symmetric multiplication of \f$U\f$ by a transpose of a general matrix \f$C\f$, result of which is stored in the current class.
+        @param C matrix to multiply with
+        @param U a space where the inverse is stored.
+        */
+        void mult_sym_t ( const mat &C, ldmat &U ) const;
+
+
+        /*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$
+
+        The new decomposition fullfills: \f$A'*diag(D)*A = self.L'*diag(self.D)*self.L\f$
+        @param A general matrix
+        @param D0 general vector
+        */
+        void ldform ( const mat &A, const vec &D0 );
+
+        //! Access functions
+        void setD ( const vec &nD ) {
+                D = nD;
+        }
+        //! Access functions
+        void setD ( const vec &nD, int i ) {
+                D.replace_mid ( i, nD );    //Fixme can be more general
+        }
+        //! Access functions
+        void setL ( const vec &nL ) {
+                L = nL;
+        }
+
+        //! Access functions
+        const vec& _D() const {
+                return D;
+        }
+        //! Access functions
+        const mat& _L() const {
+                return L;
+        }
+
+        //! add another ldmat matrix
+        ldmat& operator += ( const ldmat &ldA );
+        //! subtract another ldmat matrix
+        ldmat& operator -= ( const ldmat &ldA );
+        //! multiply by a scalar
+        ldmat& operator *= ( double x );
+
+        //! print both \c L and \c D
+        friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq );
+
+protected:
+        //! Positive vector \f$D\f$
+        vec D;
+        //! Lower-triangular matrix \f$L\f$
+        mat L;
 
 };
@@ -274 +319 @@
 //////// Operations:
 //!mapping of add operation to operators
-inline ldmat& ldmat::operator += ( const ldmat &ldA )  {this->add ( ldA );return *this;}
+inline ldmat& ldmat::operator += ( const ldmat & ldA )  {
+        this->add ( ldA );
+        return *this;
+}
 //!mapping of negative add operation to operators
-inline ldmat& ldmat::operator -= ( const ldmat &ldA )  {this->add ( ldA,-1.0 );return *this;}
+inline ldmat& ldmat::operator -= ( const ldmat & ldA )  {
+        this->add ( ldA, -1.0 );
+        return *this;
+}
 
 #endif // DC_H