Changeset 22 for bdm/math

Show
Ignore:
Timestamp:
02/18/08 17:50:37 (16 years ago)
Author:
smidl
Message:

upravy Kalmana

Location:
bdm/math
Files:
2 modified

Legend:

Unmodified
Added
Removed

  • bdm/math/libDC.cpp

    r21 r22  
    1313mat ltuinv( const mat &L ); 
    1414 
     15void fsqmat::opupdt ( const vec &v, double w ) {M+=outer_product ( v,v*w );}; 
     16mat fsqmat::to_mat()  {return M;}; 
     17void fsqmat::mult_sym ( const mat &C, bool trans ) {M=C *M*C.T();}; 
     18void fsqmat::mult_sym ( const mat &C, fsqmat &U, bool trans ) 
     19{ 
     20        mat UM= ( C *(M*C.T()) ); 
     21        U=UM; 
     22}; 
     23void fsqmat::inv ( fsqmat &Inv ) {mat IM = itpp::inv ( M ); Inv=IM;}; 
     24void fsqmat::clear() {M.clear();}; 
     25fsqmat::fsqmat ( const mat &M0 ) 
     26{ 
     27        it_assert_debug ( ( M0.cols() ==M0.rows() ),"M0 must be square" ); 
     28        M=M0;dim=M0.cols(); 
     29}; 
     30fsqmat::fsqmat() {}; 
    1531 
    1632 

  • bdm/math/libDC.h

    r19 r22  
    2222All operations defined on this class should be optimized for the chosed decomposition. 
    2323*/ 
    24 class sqmat { 
    25 public: 
    26         /*! 
    27          * Perfroms a rank-1 update by outer product of vectors: $V = V + w v v'$. 
    28          * @param  v Vector forming the outer product to be added 
    29          * @param  w weight of updating; can be negative 
     24class sqmat 
     25{ 
     26        public: 
     27                /*! 
     28                 * Perfroms a rank-1 update by outer product of vectors: $V = V + w v v'$. 
     29                 * @param  v Vector forming the outer product to be added 
     30                 * @param  w weight of updating; can be negative 
    3031 
    31         BLAS-2b operation. 
    32         */ 
    33         virtual void opupdt( const vec &v, double w ) =0; 
     32                BLAS-2b operation. 
     33                */ 
     34                virtual void opupdt ( const vec &v, double w ) =0; 
    3435 
    35         /*! \brief Conversion to full matrix. 
    36         */ 
     36                /*! \brief Conversion to full matrix. 
     37                */ 
    3738 
    38         virtual mat to_mat() =0; 
     39                virtual mat to_mat() =0; 
    3940 
    40         /*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$ 
    41         @param C multiplying matrix, 
    42         @param trans if true, product $V = C'*V*C$ will be computed instead; 
    43         */ 
    44         virtual void mult_sym( const mat &C, bool trans=true ) =0; 
     41                /*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$ 
     42                @param C multiplying matrix, 
     43                @param trans if true, product $V = C'*V*C$ will be computed instead; 
     44                */ 
     45                virtual void mult_sym ( const mat &C, bool trans=true ) =0; 
    4546 
    46          
    47         /*! 
    48         \brief Logarithm of a determinant. 
    49          
    50         */ 
    51         virtual double logdet() =0; 
    5247 
    53         /*! 
    54         \brief Multiplies square root of $V$ by vector $x$. 
    55          
    56         Used e.g. in generating normal samples. 
    57         */ 
    58         virtual vec sqrt_mult(vec &v) =0; 
    59          
    60         /*! 
    61         \brief Evaluates quadratic form $x= v'*V*v$; 
    62          
    63         */ 
    64         virtual double qform(vec &v) =0; 
     48                /*! 
     49                \brief Logarithm of a determinant. 
    6550 
    66 //      //! easy version of the  
     51                */ 
     52                virtual double logdet() =0; 
     53 
     54                /*! 
     55                \brief Multiplies square root of $V$ by vector $x$. 
     56 
     57                Used e.g. in generating normal samples. 
     58                */ 
     59                virtual vec sqrt_mult ( vec &v ) =0; 
     60 
     61                /*! 
     62                \brief Evaluates quadratic form $x= v'*V*v$; 
     63 
     64                */ 
     65                virtual double qform ( vec &v ) =0; 
     66 
     67//      //! easy version of the 
    6768//      sqmat inv(); 
    6869 
    69         //! Clearing matrix so that it corresponds to zeros. 
    70         virtual void clear() =0; 
    71          
    72         //! Reimplementing common functions of mat: cols(). 
    73         virtual int cols() =0; 
     70                //! Clearing matrix so that it corresponds to zeros. 
     71                virtual void clear() =0; 
    7472 
    75         //! Reimplementing common functions of mat: cols(). 
    76         virtual int rows() =0; 
     73                //! Reimplementing common functions of mat: cols(). 
     74                int cols() const {return dim;}; 
    7775 
    78 protected: 
    79         int dim; 
     76                //! Reimplementing common functions of mat: cols(). 
     77                int rows() const {return dim;}; 
     78 
     79        protected: 
     80                int dim; 
    8081}; 
    8182 
    … …  
    8586This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices; 
    8687*/ 
    87 class fsqmat: sqmat { 
    88         void opupdt( const vec &v, double w ); 
    89         mat to_mat(); 
    90         void mult_sym( const mat &C, bool trans=false ); 
    91    void mult_sym( const mat &C, fsqmat &U, bool trans=false ); 
    92         void inv(fsqmat &Inv); 
    93         void clear(); 
    94          
    95          
    96         //! Constructor 
    97         fsqmat(const mat &M); 
    98          
    99         /*! \brief Matrix inversion preserving the chosen form. 
     88class fsqmat: public sqmat 
     89{ 
     90        protected: 
     91                mat M; 
     92        public: 
     93                void opupdt ( const vec &v, double w ); 
     94                mat to_mat() ; 
     95                void mult_sym ( const mat &C, bool trans=false ); 
     96                void mult_sym ( const mat &C, fsqmat &U, bool trans=false ); 
     97                void clear(); 
    10098 
    101         @param Inv a space where the inverse is stored. 
    102          
    103         */ 
    104         virtual void inv(fsqmat* Inv); 
    105          
    10699 
     100                fsqmat(); // mat will be initialized OK 
     101                //! Constructor 
     102                fsqmat ( const mat &M ); 
     103 
     104                /*! \brief Matrix inversion preserving the chosen form. 
     105 
     106                @param Inv a space where the inverse is stored. 
     107 
     108                */ 
     109                virtual void inv ( fsqmat &Inv ); 
     110 
     111                double logdet() {return log ( det ( M ) );}; 
     112                double qform ( vec &v ) {return ( v* ( M*v ) );}; 
     113                vec sqrt_mult ( vec &v ) {it_error ( "not implemented" );return v;}; 
     114 
     115                fsqmat& operator += ( const fsqmat &A ) {M+=A.M;return *this;}; 
     116                fsqmat& operator -= ( const fsqmat &A ) {M-=A.M;return *this;}; 
     117                fsqmat& operator *= ( double x ) {M*=x;return *this;}; 
     118//              fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;}; 
    107119}; 
    108120 
    109 class ldmat: sqmat { 
    110 public: 
     121class ldmat: sqmat 
     122{ 
     123        public: 
    111124 
    112         //! Construct by copy of L and D. 
    113         ldmat( const mat &L, const vec &D ); 
    114         //! Construct by decomposition of full matrix V. 
    115         ldmat( mat V ); 
    116         //! Construct diagonal matrix with diagonal D0 
    117         ldmat( vec D0 ); 
    118         ldmat (); 
     125                //! Construct by copy of L and D. 
     126                ldmat ( const mat &L, const vec &D ); 
     127                //! Construct by decomposition of full matrix V. 
     128                ldmat ( mat V ); 
     129                //! Construct diagonal matrix with diagonal D0 
     130                ldmat ( vec D0 ); 
     131                ldmat (); 
    119132 
    120         // Reimplementation of compulsory operatios 
     133                // Reimplementation of compulsory operatios 
    121134 
    122         void opupdt( const vec &v, double w ); 
    123         mat to_mat(); 
    124    void mult_sym( const mat &C, bool trans=false ); 
    125         void add ( const ldmat &ld2, double w=1.0 ); 
    126         double logdet(); 
    127         double qform(vec &v); 
     135                void opupdt ( const vec &v, double w ); 
     136                mat to_mat(); 
     137                void mult_sym ( const mat &C, bool trans=false ); 
     138                void add ( const ldmat &ld2, double w=1.0 ); 
     139                double logdet(); 
     140                double qform ( vec &v ); 
    128141//      sqmat& operator -= ( const sqmat & ld2 ); 
    129         void clear(); 
    130         int cols(); 
    131         int rows(); 
    132         vec sqrt_mult(vec &v); 
     142                void clear(); 
     143                int cols(); 
     144                int rows(); 
     145                vec sqrt_mult ( vec &v ); 
    133146 
    134         /*! \brief Matrix inversion preserving the chosen form. 
     147                /*! \brief Matrix inversion preserving the chosen form. 
    135148 
    136         @param Inv a space where the inverse is stored. 
    137          
    138         */ 
    139         virtual void inv(ldmat &Inv); 
    140          
    141         /*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class. 
     149                @param Inv a space where the inverse is stored. 
    142150 
    143         @param Inv a space where the inverse is stored. 
    144          
    145         */ 
    146    void mult_sym( const mat &C, ldmat &U, bool trans=false ); 
     151                */ 
     152                virtual void inv ( ldmat &Inv ); 
    147153 
    148         /*! \brief Transforms general $A'D0A$ into pure $L'DL$ 
    149          
    150         The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$ 
     154                /*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class. 
    151155 
    152         @param A general matrix 
    153         @param D0 general vector 
    154          
    155         */ 
    156    void ldform( mat &A, vec &D0 ); 
     156                @param Inv a space where the inverse is stored. 
    157157 
    158         ldmat& operator += (const ldmat &ldA); 
    159         ldmat& operator -= (const ldmat &ldA); 
    160         ldmat& operator *= (double x); 
    161          
    162         friend std::ostream &operator<< ( std::ostream &os, ldmat &sq ); 
     158                */ 
     159                void mult_sym ( const mat &C, ldmat &U, bool trans=false ); 
    163160 
    164 protected: 
    165         vec D; 
    166         mat L; 
     161                /*! \brief Transforms general $A'D0A$ into pure $L'DL$ 
     162 
     163                The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$ 
     164 
     165                @param A general matrix 
     166                @param D0 general vector 
     167 
     168                */ 
     169                void ldform ( mat &A, vec &D0 ); 
     170 
     171                ldmat& operator += ( const ldmat &ldA ); 
     172                ldmat& operator -= ( const ldmat &ldA ); 
     173                ldmat& operator *= ( double x ); 
     174 
     175                friend std::ostream &operator<< ( std::ostream &os, ldmat &sq ); 
     176 
     177        protected: 
     178                vec D; 
     179                mat L; 
    167180 
    168181}; 
    … …  
    170183//////// Operations: 
    171184 
    172 inline ldmat& ldmat::operator += (const ldmat &ldA)  {this->add(ldA);return *this;} 
    173 inline ldmat& ldmat::operator -= (const ldmat &ldA)  {this->add(ldA,-1.0);return *this;} 
    174 inline int ldmat::cols(){return dim;} 
    175 inline int ldmat::rows(){return dim;} 
     185inline ldmat& ldmat::operator += ( const ldmat &ldA )  {this->add ( ldA );return *this;} 
     186inline ldmat& ldmat::operator -= ( const ldmat &ldA )  {this->add ( ldA,-1.0 );return *this;} 
     187inline int ldmat::cols() {return dim;} 
     188inline int ldmat::rows() {return dim;} 
    176189 
    177190#endif // DC_H