Show
Ignore:
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

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • library/bdm/math/square_mat.h

    r471 r477  
    2020 
    2121//! Auxiliary function dydr; dyadic reduction 
    22 void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx ); 
     22void dydr ( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx ); 
    2323 
    2424//! Auxiliary function ltuinv; inversion of a triangular matrix; 
    2525//TODO can be done via: dtrtri.f from lapack 
    26 mat ltuinv( const mat &L ); 
     26mat ltuinv ( const mat &L ); 
    2727 
    2828/*! \brief Abstract class for representation of double symmetric matrices in square-root form. 
    … …  
    3030All operations defined on this class should be optimized for the chosen decomposition. 
    3131*/ 
    32 class sqmat 
    33 { 
    34         public: 
    35                 /*! 
    36                  * Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$. 
    37                  * @param  v Vector forming the outer product to be added 
    38                  * @param  w weight of updating; can be negative 
    39  
    40                  BLAS-2b operation. 
    41                  */ 
    42                 virtual void opupdt ( const vec &v, double w ) =0; 
    43  
    44                 /*! \brief Conversion to full matrix. 
    45                 */ 
    46  
    47                 virtual mat to_mat() const =0; 
    48  
    49                 /*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$ 
    50                 @param C multiplying matrix, 
    51                 */ 
    52                 virtual void mult_sym ( const mat &C ) =0; 
    53                  
    54                 /*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'*V*C\f$ 
    55                 @param C multiplying matrix, 
    56                 */ 
    57                 virtual void mult_sym_t ( const mat &C ) =0; 
    58  
    59  
    60                 /*! 
    61                 \brief Logarithm of a determinant. 
    62  
    63                 */ 
    64                 virtual double logdet() const =0; 
    65  
    66                 /*! 
    67                 \brief Multiplies square root of \f$V\f$ by vector \f$x\f$. 
    68  
    69                 Used e.g. in generating normal samples. 
    70                 */ 
    71                 virtual vec sqrt_mult (const vec &v ) const =0; 
    72  
    73                 /*! 
    74                 \brief Evaluates quadratic form \f$x= v'*V*v\f$; 
    75  
    76                 */ 
    77                 virtual double qform (const vec &v ) const =0; 
    78  
    79                 /*! 
    80                 \brief Evaluates quadratic form \f$x= v'*inv(V)*v\f$; 
    81  
    82                 */ 
    83                 virtual double invqform (const vec &v ) const =0; 
     32class sqmat { 
     33public: 
     34        /*! 
     35         * Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$. 
     36         * @param  v Vector forming the outer product to be added 
     37         * @param  w weight of updating; can be negative 
     38 
     39         BLAS-2b operation. 
     40         */ 
     41        virtual void opupdt ( const vec &v, double w ) = 0; 
     42 
     43        /*! \brief Conversion to full matrix. 
     44        */ 
     45 
     46        virtual mat to_mat() const = 0; 
     47 
     48        /*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$ 
     49        @param C multiplying matrix, 
     50        */ 
     51        virtual void mult_sym ( const mat &C ) = 0; 
     52 
     53        /*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'*V*C\f$ 
     54        @param C multiplying matrix, 
     55        */ 
     56        virtual void mult_sym_t ( const mat &C ) = 0; 
     57 
     58 
     59        /*! 
     60        \brief Logarithm of a determinant. 
     61 
     62        */ 
     63        virtual double logdet() const = 0; 
     64 
     65        /*! 
     66        \brief Multiplies square root of \f$V\f$ by vector \f$x\f$. 
     67 
     68        Used e.g. in generating normal samples. 
     69        */ 
     70        virtual vec sqrt_mult ( const vec &v ) const = 0; 
     71 
     72        /*! 
     73        \brief Evaluates quadratic form \f$x= v'*V*v\f$; 
     74 
     75        */ 
     76        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; 
    8483 
    8584//      //! easy version of the 
    8685//      sqmat inv(); 
    8786 
    88                 //! Clearing matrix so that it corresponds to zeros. 
    89                 virtual void clear() =0; 
    90  
    91                 //! Reimplementing common functions of mat: cols(). 
    92                 int cols() const {return dim;}; 
    93  
    94                 //! Reimplementing common functions of mat: rows(). 
    95                 int rows() const {return dim;}; 
    96  
    97                 //! Destructor for future use; 
    98                 virtual ~sqmat(){}; 
    99                 //! Default constructor 
    100                 sqmat(const int dim0): dim(dim0){}; 
    101                 //! Default constructor 
    102                 sqmat(): dim(0){}; 
    103         protected: 
    104                 //! dimension of the square matrix 
    105                 int dim; 
     87        //! Clearing matrix so that it corresponds to zeros. 
     88        virtual void clear() = 0; 
     89 
     90        //! Reimplementing common functions of mat: cols(). 
     91        int cols() const { 
     92                return dim; 
     93        }; 
     94 
     95        //! Reimplementing common functions of mat: rows(). 
     96        int rows() const { 
     97                return dim; 
     98        }; 
     99 
     100        //! Destructor for future use; 
     101        virtual ~sqmat() {}; 
     102        //! Default constructor 
     103        sqmat ( const int dim0 ) : dim ( dim0 ) {}; 
     104        //! Default constructor 
     105        sqmat() : dim ( 0 ) {}; 
     106protected: 
     107        //! dimension of the square matrix 
     108        int dim; 
    106109}; 
    107110 
    … …  
    111114This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices; 
    112115*/ 
    113 class fsqmat: public sqmat 
    114 { 
    115         protected: 
    116                 //! Full matrix on which the operations are performed 
    117                 mat M; 
    118         public: 
    119                 void opupdt ( const vec &v, double w ); 
    120                 mat to_mat() const; 
    121                 void mult_sym ( const mat &C); 
    122                 void mult_sym_t ( const mat &C); 
    123                 //! store result of \c mult_sym in external matrix \f$U\f$ 
    124                 void mult_sym ( const mat &C, fsqmat &U) const; 
    125                 //! store result of \c mult_sym_t in external matrix \f$U\f$ 
    126                 void mult_sym_t ( const mat &C, fsqmat &U) const; 
    127                 void clear(); 
    128  
    129                 //! Default initialization 
    130                 fsqmat(){}; // mat will be initialized OK 
    131                 //! Default initialization with proper size 
    132                 fsqmat(const int dim0); // mat will be initialized OK 
    133                 //! Constructor 
    134                 fsqmat ( const mat &M ); 
    135                 //! Constructor 
    136                 fsqmat ( const fsqmat &M, const ivec &perm ):sqmat(M.rows()){it_error("not implemneted");}; 
    137                 //! Constructor 
    138                 fsqmat ( const vec &d ):sqmat(d.length()){M=diag(d);}; 
    139  
    140                 //! Destructor for future use; 
    141                 virtual ~fsqmat(){}; 
    142  
    143  
    144                 /*! \brief Matrix inversion preserving the chosen form. 
    145  
    146                 @param Inv a space where the inverse is stored. 
    147  
    148                 */ 
    149                 void inv ( fsqmat &Inv ) const; 
    150  
    151                 double logdet() const {return log ( det ( M ) );}; 
    152                 double qform (const  vec &v ) const {return ( v* ( M*v ) );}; 
    153                 double invqform (const  vec &v ) const {return ( v* ( itpp::inv(M)*v ) );}; 
    154                 vec sqrt_mult (const vec &v ) const {mat Ch=chol(M); return Ch*v;}; 
    155  
    156                 //! Add another matrix in fsq form with weight w 
    157                 void add ( const fsqmat &fsq2, double w=1.0 ){M+=fsq2.M;}; 
    158  
    159                 //! Access functions 
    160                 void setD (const vec &nD){M=diag(nD);} 
    161                 //! Access functions 
    162                 vec getD (){return diag(M);} 
    163                 //! Access functions 
    164                 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 
    165  
    166  
    167                 //! add another fsqmat matrix 
    168                 fsqmat& operator += ( const fsqmat &A ) {M+=A.M;return *this;}; 
    169                 //! subtrack another fsqmat matrix 
    170                 fsqmat& operator -= ( const fsqmat &A ) {M-=A.M;return *this;}; 
    171                 //! multiply by a scalar 
    172                 fsqmat& operator *= ( double x ) {M*=x;return *this;}; 
     116class fsqmat: public sqmat { 
     117protected: 
     118        //! Full matrix on which the operations are performed 
     119        mat M; 
     120public: 
     121        void opupdt ( const vec &v, double w ); 
     122        mat to_mat() const; 
     123        void mult_sym ( const mat &C ); 
     124        void mult_sym_t ( const mat &C ); 
     125        //! store result of \c mult_sym in external matrix \f$U\f$ 
     126        void mult_sym ( const mat &C, fsqmat &U ) const; 
     127        //! store result of \c mult_sym_t in external matrix \f$U\f$ 
     128        void mult_sym_t ( const mat &C, fsqmat &U ) const; 
     129        void clear(); 
     130 
     131        //! Default initialization 
     132        fsqmat() {}; // mat will be initialized OK 
     133        //! Default initialization with proper size 
     134        fsqmat ( const int dim0 ); // mat will be initialized OK 
     135        //! Constructor 
     136        fsqmat ( const mat &M ); 
     137        //! Constructor 
     138        fsqmat ( const fsqmat &M, const ivec &perm ) : sqmat ( M.rows() ) { 
     139                it_error ( "not implemneted" ); 
     140        }; 
     141        //! Constructor 
     142        fsqmat ( const vec &d ) : sqmat ( d.length() ) { 
     143                M = diag ( d ); 
     144        }; 
     145 
     146        //! Destructor for future use; 
     147        virtual ~fsqmat() {}; 
     148 
     149 
     150        /*! \brief Matrix inversion preserving the chosen form. 
     151 
     152        @param Inv a space where the inverse is stored. 
     153 
     154        */ 
     155        void inv ( fsqmat &Inv ) const; 
     156 
     157        double logdet() const { 
     158                return log ( det ( M ) ); 
     159        }; 
     160        double qform ( const  vec &v ) const { 
     161                return ( v* ( M*v ) ); 
     162        }; 
     163        double invqform ( const  vec &v ) const { 
     164                return ( v* ( itpp::inv ( M ) *v ) ); 
     165        }; 
     166        vec sqrt_mult ( const vec &v ) const { 
     167                mat Ch = chol ( M ); 
     168                return Ch*v; 
     169        }; 
     170 
     171        //! Add another matrix in fsq form with weight w 
     172        void add ( const fsqmat &fsq2, double w = 1.0 ) { 
     173                M += fsq2.M; 
     174        }; 
     175 
     176        //! Access functions 
     177        void setD ( const vec &nD ) { 
     178                M = diag ( nD ); 
     179        } 
     180        //! Access functions 
     181        vec getD () { 
     182                return diag ( M ); 
     183        } 
     184        //! Access functions 
     185        void setD ( const vec &nD, int i ) { 
     186                for ( int j = i; j < nD.length(); j++ ) { 
     187                        M ( j, j ) = nD ( j - i );    //Fixme can be more general 
     188                } 
     189        } 
     190 
     191 
     192        //! add another fsqmat matrix 
     193        fsqmat& operator += ( const fsqmat &A ) { 
     194                M += A.M; 
     195                return *this; 
     196        }; 
     197        //! subtrack another fsqmat matrix 
     198        fsqmat& operator -= ( const fsqmat &A ) { 
     199                M -= A.M; 
     200                return *this; 
     201        }; 
     202        //! multiply by a scalar 
     203        fsqmat& operator *= ( double x ) { 
     204                M *= x; 
     205                return *this; 
     206        }; 
    173207//              fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;}; 
    174                 //! print full matrix 
    175                 friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq ); 
     208        //! print full matrix 
     209        friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq ); 
    176210 
    177211}; 
    178212 
    179 /*! \brief Matrix stored in LD form, (commonly known as UD)  
     213/*! \brief Matrix stored in LD form, (commonly known as UD) 
    180214 
    181215Matrix is decomposed as follows: \f[M = L'DL\f] where only \f$L\f$ and \f$D\f$ matrices are stored. 
    182216All inplace operations modifies only these and the need to compose and decompose the matrix is avoided. 
    183217*/ 
    184 class ldmat: public sqmat 
    185 { 
    186         public: 
    187                 //! Construct by copy of L and D. 
    188                 ldmat ( const mat &L, const vec &D ); 
    189                 //! Construct by decomposition of full matrix V. 
    190                 ldmat (const mat &V ); 
    191                 //! Construct by restructuring of V0 accordint to permutation vector perm. 
    192                 ldmat (const ldmat &V0, const ivec &perm):sqmat(V0.rows()){     ldform(V0.L.get_cols(perm), V0.D);}; 
    193                 //! Construct diagonal matrix with diagonal D0 
    194                 ldmat ( vec D0 ); 
    195                 //!Default constructor 
    196                 ldmat (); 
    197                 //! Default initialization with proper size 
    198                 ldmat(const int dim0); 
    199                  
    200                 //! Destructor for future use; 
    201                 virtual ~ldmat(){}; 
    202  
    203                 // Reimplementation of compulsory operatios 
    204  
    205                 void opupdt ( const vec &v, double w ); 
    206                 mat to_mat() const; 
    207                 void mult_sym ( const mat &C); 
    208                 void mult_sym_t ( const mat &C); 
    209                 //! Add another matrix in LD form with weight w 
    210                 void add ( const ldmat &ld2, double w=1.0 ); 
    211                 double logdet() const; 
    212                 double qform (const vec &v ) const; 
    213                 double invqform (const vec &v ) const; 
    214                 void clear(); 
    215                 vec sqrt_mult ( const vec &v ) const; 
    216  
    217                  
    218                 /*! \brief Matrix inversion preserving the chosen form. 
    219                 @param Inv a space where the inverse is stored. 
    220                 */ 
    221                 void inv ( ldmat &Inv ) const; 
    222  
    223                 /*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class. 
    224                 @param C matrix to multiply with 
    225                 @param U a space where the inverse is stored. 
    226                 */ 
    227                 void mult_sym ( const mat &C, ldmat &U) const; 
    228  
    229                 /*! \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. 
    230                 @param C matrix to multiply with 
    231                 @param U a space where the inverse is stored. 
    232                 */ 
    233                 void mult_sym_t ( const mat &C, ldmat &U) const; 
    234  
    235  
    236                 /*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$ 
    237  
    238                 The new decomposition fullfills: \f$A'*diag(D)*A = self.L'*diag(self.D)*self.L\f$ 
    239                 @param A general matrix 
    240                 @param D0 general vector 
    241                 */ 
    242                 void ldform (const mat &A,const vec &D0 ); 
    243  
    244                 //! Access functions 
    245                 void setD (const vec &nD){D=nD;} 
    246                 //! Access functions 
    247                 void setD (const vec &nD, int i){D.replace_mid(i,nD);} //Fixme can be more general 
    248                 //! Access functions 
    249                 void setL (const vec &nL){L=nL;} 
    250  
    251                 //! Access functions 
    252                 const vec& _D() const {return D;} 
    253                 //! Access functions 
    254                 const mat& _L() const {return L;} 
    255  
    256                 //! add another ldmat matrix 
    257                 ldmat& operator += ( const ldmat &ldA ); 
    258                 //! subtract another ldmat matrix 
    259                 ldmat& operator -= ( const ldmat &ldA ); 
    260                 //! multiply by a scalar 
    261                 ldmat& operator *= ( double x ); 
    262  
    263                 //! print both \c L and \c D 
    264                 friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq ); 
    265  
    266         protected: 
    267                 //! Positive vector \f$D\f$ 
    268                 vec D; 
    269                 //! Lower-triangular matrix \f$L\f$ 
    270                 mat L; 
     218class ldmat: public sqmat { 
     219public: 
     220        //! Construct by copy of L and D. 
     221        ldmat ( const mat &L, const vec &D ); 
     222        //! Construct by decomposition of full matrix V. 
     223        ldmat ( const mat &V ); 
     224        //! Construct by restructuring of V0 accordint to permutation vector perm. 
     225        ldmat ( const ldmat &V0, const ivec &perm ) : sqmat ( V0.rows() ) { 
     226                ldform ( V0.L.get_cols ( perm ), V0.D ); 
     227        }; 
     228        //! Construct diagonal matrix with diagonal D0 
     229        ldmat ( vec D0 ); 
     230        //!Default constructor 
     231        ldmat (); 
     232        //! Default initialization with proper size 
     233        ldmat ( const int dim0 ); 
     234 
     235        //! Destructor for future use; 
     236        virtual ~ldmat() {}; 
     237 
     238        // Reimplementation of compulsory operatios 
     239 
     240        void opupdt ( const vec &v, double w ); 
     241        mat to_mat() const; 
     242        void mult_sym ( const mat &C ); 
     243        void mult_sym_t ( const mat &C ); 
     244        //! Add another matrix in LD form with weight w 
     245        void add ( const ldmat &ld2, double w = 1.0 ); 
     246        double logdet() const; 
     247        double qform ( const vec &v ) const; 
     248        double invqform ( const vec &v ) const; 
     249        void clear(); 
     250        vec sqrt_mult ( const vec &v ) const; 
     251 
     252 
     253        /*! \brief Matrix inversion preserving the chosen form. 
     254        @param Inv a space where the inverse is stored. 
     255        */ 
     256        void inv ( ldmat &Inv ) const; 
     257 
     258        /*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class. 
     259        @param C matrix to multiply with 
     260        @param U a space where the inverse is stored. 
     261        */ 
     262        void mult_sym ( const mat &C, ldmat &U ) const; 
     263 
     264        /*! \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. 
     265        @param C matrix to multiply with 
     266        @param U a space where the inverse is stored. 
     267        */ 
     268        void mult_sym_t ( const mat &C, ldmat &U ) const; 
     269 
     270 
     271        /*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$ 
     272 
     273        The new decomposition fullfills: \f$A'*diag(D)*A = self.L'*diag(self.D)*self.L\f$ 
     274        @param A general matrix 
     275        @param D0 general vector 
     276        */ 
     277        void ldform ( const mat &A, const vec &D0 ); 
     278 
     279        //! Access functions 
     280        void setD ( const vec &nD ) { 
     281                D = nD; 
     282        } 
     283        //! Access functions 
     284        void setD ( const vec &nD, int i ) { 
     285                D.replace_mid ( i, nD );    //Fixme can be more general 
     286        } 
     287        //! Access functions 
     288        void setL ( const vec &nL ) { 
     289                L = nL; 
     290        } 
     291 
     292        //! Access functions 
     293        const vec& _D() const { 
     294                return D; 
     295        } 
     296        //! Access functions 
     297        const mat& _L() const { 
     298                return L; 
     299        } 
     300 
     301        //! add another ldmat matrix 
     302        ldmat& operator += ( const ldmat &ldA ); 
     303        //! subtract another ldmat matrix 
     304        ldmat& operator -= ( const ldmat &ldA ); 
     305        //! multiply by a scalar 
     306        ldmat& operator *= ( double x ); 
     307 
     308        //! print both \c L and \c D 
     309        friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq ); 
     310 
     311protected: 
     312        //! Positive vector \f$D\f$ 
     313        vec D; 
     314        //! Lower-triangular matrix \f$L\f$ 
     315        mat L; 
    271316 
    272317}; 
    … …  
    274319//////// Operations: 
    275320//!mapping of add operation to operators 
    276 inline ldmat& ldmat::operator += ( const ldmat &ldA )  {this->add ( ldA );return *this;} 
     321inline ldmat& ldmat::operator += ( const ldmat & ldA )  { 
     322        this->add ( ldA ); 
     323        return *this; 
     324} 
    277325//!mapping of negative add operation to operators 
    278 inline ldmat& ldmat::operator -= ( const ldmat &ldA )  {this->add ( ldA,-1.0 );return *this;} 
     326inline ldmat& ldmat::operator -= ( const ldmat & ldA )  { 
     327        this->add ( ldA, -1.0 ); 
     328        return *this; 
     329} 
    279330 
    280331#endif // DC_H