Changeset 75

Show
Ignore:
Timestamp:
04/18/08 14:00:29 (16 years ago)
Author:
smidl
Message:

oprava operaci na ld

Location:
bdm/math
Files:
2 modified

Legend:

Unmodified
Added
Removed
  • bdm/math/libDC.cpp

    r37 r75  
    162162        double x = 0.0, sum; 
    163163        int i,j; 
     164        vec s(v.length()); 
     165        vec S=L*v; 
    164166 
    165167        for ( i=0; i<D.length(); i++ ) { //rows of L 
    166168                sum = 0.0; 
    167                 for ( j=i; j<=i; j++ ){sum+=L( i,j )*v( i );} 
     169                for ( j=0; j<=i; j++ ){sum+=L( i,j )*v( j );} 
    168170                x +=D( i )*sum*sum; 
     171                s(i)=sum; 
     172        }; 
     173        return x; 
     174} 
     175 
     176double ldmat::invqform( const vec &v ) const { 
     177        double x = 0.0; 
     178        int i; 
     179        vec pom(v.length()); 
     180         
     181        backward_substitution(L.T(),v,pom); 
     182         
     183        for ( i=0; i<D.length(); i++ ) { //rows of L 
     184                x +=pom(i)*pom(i)/D(i); 
    169185        }; 
    170186        return x; 
     
    298314                        j = i + k; //change in .m 1+1=2, here 0+0+1=1 
    299315                        s = L( j, i ); 
    300                         for ( m = i + 1; m < ( j - 1 ); m++ ) { 
     316                        for ( m = i + 1; m < ( j ); m++ ) { 
    301317                                s += L( m, i ) * Il( j, m ); 
    302318                        } 
  • bdm/math/libDC.h

    r51 r75  
    3333        public: 
    3434                /*! 
    35                  * Perfroms a rank-1 update by outer product of vectors: $V = V + w v v'$. 
     35                 * Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$. 
    3636                 * @param  v Vector forming the outer product to be added 
    3737                 * @param  w weight of updating; can be negative 
     
    4646                virtual mat to_mat() =0; 
    4747 
    48                 /*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$ 
     48                /*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$ 
    4949                @param C multiplying matrix, 
    5050                */ 
    5151                virtual void mult_sym ( const mat &C ) =0; 
    5252                 
    53                 /*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix $C$, i.e. $V = C'*V*C$ 
     53                /*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'*V*C\f$ 
    5454                @param C multiplying matrix, 
    5555                */ 
     
    6464 
    6565                /*! 
    66                 \brief Multiplies square root of $V$ by vector $x$. 
     66                \brief Multiplies square root of \f$V\f$ by vector \f$x\f$. 
    6767 
    6868                Used e.g. in generating normal samples. 
     
    7171 
    7272                /*! 
    73                 \brief Evaluates quadratic form $x= v'*V*v$; 
     73                \brief Evaluates quadratic form \f$x= v'*V*v\f$; 
    7474 
    7575                */ 
    7676                virtual double qform (const vec &v ) const =0; 
     77 
     78                /*! 
     79                \brief Evaluates quadratic form \f$x= v'*inv(V)*v\f$; 
     80 
     81                */ 
     82                virtual double invqform (const vec &v ) const =0; 
    7783 
    7884//      //! easy version of the 
     
    112118                void mult_sym ( const mat &C); 
    113119                void mult_sym_t ( const mat &C); 
    114                 //! store result of \c mult_sym in external matrix $U$ 
     120                //! store result of \c mult_sym in external matrix \f$U\f$ 
    115121                void mult_sym ( const mat &C, fsqmat &U) const; 
    116                 //! store result of \c mult_sym_t in external matrix $U$ 
     122                //! store result of \c mult_sym_t in external matrix \f$U\f$ 
    117123                void mult_sym_t ( const mat &C, fsqmat &U) const; 
    118124                void clear(); 
     
    140146                double logdet() const {return log ( det ( M ) );}; 
    141147                double qform (const  vec &v ) const {return ( v* ( M*v ) );}; 
     148                double invqform (const  vec &v ) const {return ( v* ( itpp::inv(M)*v ) );}; 
    142149                vec sqrt_mult (const vec &v ) const {mat Ch=chol(M); return Ch*v;}; 
    143150 
     
    163170/*! \brief Matrix stored in LD form, (typically known as UD)  
    164171 
    165 Matrix is decomposed as follows: \f[M = L'DL\f] where only $L$ and $D$ matrices are stored. 
     172Matrix is decomposed as follows: \f[M = L'DL\f] where only \f$L\f$ and \f$D\f$ matrices are stored. 
    166173All inplace operations modifies only these and the need to compose and decompose the matrix is avoided. 
    167174*/ 
     
    194201                double logdet() const; 
    195202                double qform (const vec &v ) const; 
     203                double invqform (const vec &v ) const; 
    196204//      sqmat& operator -= ( const sqmat & ld2 ); 
    197205                void clear(); 
     
    205213                virtual void inv ( ldmat &Inv ) const; 
    206214 
    207                 /*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class. 
     215                /*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class. 
    208216                @param C matrix to multiply with 
    209217                @param U a space where the inverse is stored. 
     
    211219                void mult_sym ( const mat &C, ldmat &U) const; 
    212220 
    213                 /*! \brief Symmetric multiplication of $U$ by a transpose of a general matrix $C$, result of which is stored in the current class. 
     221                /*! \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. 
    214222                @param C matrix to multiply with 
    215223                @param U a space where the inverse is stored. 
     
    218226 
    219227 
    220                 /*! \brief Transforms general $A'D0 A$ into pure $L'DL$ 
    221  
    222                 The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$ 
     228                /*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$ 
     229 
     230                The new decomposition fullfills: \f$A'*diag(D)*A = self.L'*diag(self.D)*self.L\f$ 
    223231                @param A general matrix 
    224232                @param D0 general vector 
     
    245253 
    246254        protected: 
    247                 //! Positive vector $D$ 
     255                //! Positive vector \f$D\f$ 
    248256                vec D; 
    249                 //! Lower-triangular matrix $L$ 
     257                //! Lower-triangular matrix \f$L\f$ 
    250258                mat L; 
    251259