Changeset 477 for library/bdm/math/chmat.h
- Timestamp:
- 08/05/09 14:40:03 (15 years ago)
- Files:
-
- 1 modified
Legend:
- Unmodified
- Added
- Removed
-
library/bdm/math/chmat.h
r437 r477 40 40 void clear(); 41 41 //! add another chmat \c A2 with weight \c w. 42 void add ( const chmat &A2, double w=1.0 ) { 43 it_assert_debug(dim==A2.dim,"Matrices of unequal dimension"); 44 mat pre=concat_vertical(Ch,sqrt(w)*A2.Ch); 45 mat post=zeros(pre.rows(),pre.cols()); 46 if(!qr(pre,post)) {it_warning("Unstable QR in chmat add");} 47 Ch=post(0,dim-1,0,dim-1); 42 void add ( const chmat &A2, double w = 1.0 ) { 43 it_assert_debug ( dim == A2.dim, "Matrices of unequal dimension" ); 44 mat pre = concat_vertical ( Ch, sqrt ( w ) * A2.Ch ); 45 mat post = zeros ( pre.rows(), pre.cols() ); 46 if ( !qr ( pre, post ) ) { 47 it_warning ( "Unstable QR in chmat add" ); 48 } 49 Ch = post ( 0, dim - 1, 0, dim - 1 ); 48 50 }; 49 51 //!Inversion in the same form, i.e. cholesky 50 void inv ( chmat &Inv ) const { ( Inv.Ch ) = itpp::inv ( Ch ).T();}; //Fixme: can be more efficient 52 void inv ( chmat &Inv ) const { 53 ( Inv.Ch ) = itpp::inv ( Ch ).T(); 54 }; //Fixme: can be more efficient 51 55 ; 52 56 // void inv ( mat &Inv ); … … 55 59 virtual ~chmat() {}; 56 60 //! 57 chmat ( ) : sqmat (), Ch ( ) {};61 chmat ( ) : sqmat (), Ch ( ) {}; 58 62 //! Default constructor 59 chmat ( const int dim0 ) : sqmat ( dim0 ), Ch ( dim0,dim0 ) {};63 chmat ( const int dim0 ) : sqmat ( dim0 ), Ch ( dim0, dim0 ) {}; 60 64 //! Default constructor 61 chmat ( const vec &v ) : sqmat ( v.length() ),Ch ( diag(sqrt(v)) ) {};65 chmat ( const vec &v ) : sqmat ( v.length() ), Ch ( diag ( sqrt ( v ) ) ) {}; 62 66 //! Copy constructor 63 chmat ( const chmat &Ch0) : sqmat ( Ch0.dim),Ch ( Ch0.dim,Ch0.dim ) {Ch=Ch0.Ch;}; 67 chmat ( const chmat &Ch0 ) : sqmat ( Ch0.dim ), Ch ( Ch0.dim, Ch0.dim ) { 68 Ch = Ch0.Ch; 69 }; 64 70 //! Default constructor (m3k:cholform) 65 chmat ( const mat &M ) : sqmat ( M.rows() ), Ch ( M.rows(),M.cols() ) {71 chmat ( const mat &M ) : sqmat ( M.rows() ), Ch ( M.rows(), M.cols() ) { 66 72 mat Q; 67 it_assert_debug ( M.rows() == M.cols(),"chmat:: input matrix must be square!" );68 Ch =chol ( M );73 it_assert_debug ( M.rows() == M.cols(), "chmat:: input matrix must be square!" ); 74 Ch = chol ( M ); 69 75 }; 70 76 //! Constructor 71 chmat ( const chmat &M, const ivec &perm ):sqmat(M.rows()){it_error("not implemneted");}; 77 chmat ( const chmat &M, const ivec &perm ) : sqmat ( M.rows() ) { 78 it_error ( "not implemneted" ); 79 }; 72 80 //! Access function 73 mat & _Ch() {return Ch;} 81 mat & _Ch() { 82 return Ch; 83 } 74 84 //! Access function 75 const mat & _Ch() const {return Ch;} 85 const mat & _Ch() const { 86 return Ch; 87 } 76 88 //! Access functions 77 void setD ( const vec &nD ) {Ch=diag ( sqrt(nD) );} 89 void setD ( const vec &nD ) { 90 Ch = diag ( sqrt ( nD ) ); 91 } 78 92 //! Access functions 79 93 void setCh ( const vec &chQ ) { 80 it_assert_debug (chQ.length()==dim*dim,"");81 copy_vector (dim*dim, chQ._data(), Ch._data());94 it_assert_debug ( chQ.length() == dim*dim, "" ); 95 copy_vector ( dim*dim, chQ._data(), Ch._data() ); 82 96 } 83 97 //! Access functions 84 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 98 void setD ( const vec &nD, int i ) { 99 for ( int j = i; j < nD.length(); j++ ) { 100 Ch ( j, j ) = sqrt ( nD ( j - i ) ); //Fixme can be more general 101 } 102 } 85 103 86 104 //! Operators 87 105 chmat& operator += ( const chmat &A2 ); 88 106 chmat& operator -= ( const chmat &A2 ); 89 chmat& operator * ( const double &d ){Ch*sqrt(d); return *this;}; 90 chmat& operator = ( const chmat &A2 ){Ch=A2.Ch;dim=A2.dim;return *this;} 91 chmat& operator *= ( double x ){Ch*=sqrt(x); return *this;}; 107 chmat& operator * ( const double &d ) { 108 Ch*sqrt ( d ); 109 return *this; 110 }; 111 chmat& operator = ( const chmat &A2 ) { 112 Ch = A2.Ch; 113 dim = A2.dim; 114 return *this; 115 } 116 chmat& operator *= ( double x ) { 117 Ch *= sqrt ( x ); 118 return *this; 119 }; 92 120 }; 93 121 … … 95 123 //////// Operations: 96 124 //!mapping of add operation to operators 97 inline chmat& chmat::operator += ( const chmat &A2 ) {this->add ( A2 );return *this;} 125 inline chmat& chmat::operator += ( const chmat & A2 ) { 126 this->add ( A2 ); 127 return *this; 128 } 98 129 //!mapping of negative add operation to operators 99 inline chmat& chmat::operator -= ( const chmat &A2 ) {this->add ( A2,-1.0 );return *this;} 130 inline chmat& chmat::operator -= ( const chmat & A2 ) { 131 this->add ( A2, -1.0 ); 132 return *this; 133 } 100 134 101 135 #endif // CHMAT_H