MatrixGeneral.h

Go to the documentation of this file.
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure
00002  * calculations.
00003  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
00004  * 
00005  * This program is free software: you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation, either version 3 of the License, or
00008  * (at your option) any later version.
00009  * 
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  * 
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00017  * 
00018  * Primary academic reference:
00019  * Kohn−Sham Density Functional Theory Electronic Structure Calculations 
00020  * with Linearly Scaling Computational Time and Memory Usage,
00021  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
00022  * J. Chem. Theory Comput. 7, 340 (2011),
00023  * <http://dx.doi.org/10.1021/ct100611z>
00024  * 
00025  * For further information about Ergo, see <http://www.ergoscf.org>.
00026  */
00027 
00036 #ifndef MAT_MATRIXGENERAL
00037 #define MAT_MATRIXGENERAL
00038 #include "MatrixBase.h"
00039 namespace mat {
00056   template<typename Treal, typename Tmatrix>
00057     class MatrixGeneral : public MatrixBase<Treal, Tmatrix> {
00058   public:
00059     typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType;
00060 
00061     MatrixGeneral()
00062       :MatrixBase<Treal, Tmatrix>() {} 
00063     explicit MatrixGeneral(const MatrixGeneral<Treal, Tmatrix>& matr)
00064       :MatrixBase<Treal, Tmatrix>(matr) {} 
00065     explicit MatrixGeneral(const MatrixSymmetric<Treal, Tmatrix>& symm)
00066       :MatrixBase<Treal, Tmatrix>(symm) { 
00067       this->matrixPtr->symToNosym();
00068     }
00069     explicit MatrixGeneral(const MatrixTriangular<Treal, Tmatrix>& triang)
00070       :MatrixBase<Treal, Tmatrix>(triang) {}
00073 #if 0
00074     template<typename Tfull>
00075       inline void assign_from_full
00076       (Tfull const* const fullmatrix, 
00077        const int totnrows, const int totncols) {
00078       this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);      
00079     }    
00080     inline void assign_from_full
00081       (Treal const* const fullmatrix, 
00082        const int totnrows, const int totncols) {
00083       this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);      
00084     }
00085 #endif
00086 
00087     inline void assignFromFull
00088       (std::vector<Treal> const & fullMat) {
00089       assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
00090       this->matrixPtr->assignFromFull(fullMat);      
00091     }
00092 
00093     inline void fullMatrix(std::vector<Treal> & fullMat) const {
00094       this->matrixPtr->fullMatrix(fullMat);
00095     }    
00096 
00097     inline void fullMatrix
00098       (std::vector<Treal> & fullMat,
00099        std::vector<int> const & rowInversePermutation, 
00100        std::vector<int> const & colInversePermutation) const {
00101       std::vector<int> rowind;
00102       std::vector<int> colind; 
00103       std::vector<Treal> values;
00104       get_all_values(rowind, colind, values, 
00105                      rowInversePermutation, 
00106                      colInversePermutation);
00107       fullMat.assign(this->get_nrows()*this->get_ncols(),0);
00108       for (unsigned int ind = 0; ind < values.size(); ++ind) 
00109         fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
00110           values[ind];
00111     }
00119     inline void assign_from_sparse
00120       (std::vector<int> const & rowind, 
00121        std::vector<int> const & colind, 
00122        std::vector<Treal> const & values, 
00123        SizesAndBlocks const & newRows,
00124        SizesAndBlocks const & newCols) {
00125       this->resetSizesAndBlocks(newRows, newCols);
00126       this->matrixPtr->assignFromSparse(rowind, colind, values);
00127     }
00135     inline void assign_from_sparse
00136       (std::vector<int> const & rowind, 
00137        std::vector<int> const & colind, 
00138        std::vector<Treal> const & values, 
00139        std::vector<int> const & rowPermutation, 
00140        std::vector<int> const & colPermutation) {
00141       std::vector<int> newRowind;
00142       std::vector<int> newColind; 
00143       MatrixBase<Treal, Tmatrix>::
00144 	getPermutedIndexes(rowind, rowPermutation, newRowind);
00145       MatrixBase<Treal, Tmatrix>::
00146 	getPermutedIndexes(colind, colPermutation, newColind);
00147       this->matrixPtr->assignFromSparse(newRowind, newColind, values);
00148     }
00154     inline void assign_from_sparse
00155       (std::vector<int> const & rowind, 
00156        std::vector<int> const & colind, 
00157        std::vector<Treal> const & values, 
00158        SizesAndBlocks const & newRows,
00159        SizesAndBlocks const & newCols,
00160        std::vector<int> const & rowPermutation, 
00161        std::vector<int> const & colPermutation) {
00162       this->resetSizesAndBlocks(newRows, newCols);
00163       assign_from_sparse(rowind, colind, values, 
00164                          rowPermutation, colPermutation);
00165     }
00170     inline void get_values
00171       (std::vector<int> const & rowind, 
00172        std::vector<int> const & colind, 
00173        std::vector<Treal> & values) const {
00174       this->matrixPtr->getValues(rowind, colind, values);
00175     }
00183     inline void get_values
00184       (std::vector<int> const & rowind, 
00185        std::vector<int> const & colind, 
00186        std::vector<Treal> & values,
00187        std::vector<int> const & rowPermutation, 
00188        std::vector<int> const & colPermutation) const {
00189       std::vector<int> newRowind;
00190       std::vector<int> newColind; 
00191       MatrixBase<Treal, Tmatrix>::
00192 	getPermutedIndexes(rowind, rowPermutation, newRowind);
00193       MatrixBase<Treal, Tmatrix>::
00194 	getPermutedIndexes(colind, colPermutation, newColind);
00195       this->matrixPtr->getValues(newRowind, newColind, values);
00196     }
00201     inline void get_all_values
00202       (std::vector<int> & rowind, 
00203        std::vector<int> & colind, 
00204        std::vector<Treal> & values) const {
00205       rowind.resize(0);
00206       colind.resize(0);
00207       values.resize(0);
00208       this->matrixPtr->getAllValues(rowind, colind, values);
00209     }
00219     inline void get_all_values
00220       (std::vector<int> & rowind, 
00221        std::vector<int> & colind, 
00222        std::vector<Treal> & values,
00223        std::vector<int> const & rowInversePermutation, 
00224        std::vector<int> const & colInversePermutation) const {
00225       std::vector<int> tmpRowind;
00226       std::vector<int> tmpColind; 
00227       tmpRowind.reserve(rowind.capacity());
00228       tmpColind.reserve(colind.capacity());
00229       values.resize(0);
00230       this->matrixPtr->getAllValues(tmpRowind, tmpColind, values);
00231       MatrixBase<Treal, Tmatrix>::
00232 	getPermutedIndexes(tmpRowind, rowInversePermutation, rowind);
00233       MatrixBase<Treal, Tmatrix>::
00234 	getPermutedIndexes(tmpColind, colInversePermutation, colind);
00235       
00236     }
00246 #if 0
00247     inline void fullmatrix(Treal* const full, 
00248                            const int totnrows, const int totncols) const {
00249       this->matrixPtr->fullmatrix(full, totnrows, totncols);
00250     }
00251 #endif
00252     MatrixGeneral<Treal, Tmatrix>& 
00253       operator=(const MatrixGeneral<Treal, Tmatrix>& mat) {
00254       MatrixBase<Treal, Tmatrix>::operator=(mat);
00255       return *this;
00256     } 
00257     inline MatrixGeneral<Treal, Tmatrix>& 
00258       operator=(const Xtrans<MatrixGeneral<Treal, Tmatrix> >& mt) {
00259       if (mt.tA)
00260         MatrixBase<Treal, Tmatrix>::operator=(transpose(mt.A));
00261       else
00262         MatrixBase<Treal, Tmatrix>::operator=(mt.A);
00263       return *this;
00264     }
00265 
00266     MatrixGeneral<Treal, Tmatrix>& 
00267       operator=(const MatrixSymmetric<Treal, Tmatrix>& symm) {
00268       MatrixBase<Treal, Tmatrix>::operator=(symm);
00269       this->matrixPtr->symToNosym();
00270       return *this;
00271     } 
00272     MatrixGeneral<Treal, Tmatrix>& 
00273       operator=(const MatrixTriangular<Treal, Tmatrix>& triang) {
00274       MatrixBase<Treal, Tmatrix>::operator=(triang);
00275       return *this;
00276     } 
00277 
00278     inline MatrixGeneral<Treal, Tmatrix>& operator=(int const k) {
00279       *this->matrixPtr = k;
00280       return *this;
00281     }
00282     inline Treal frob() const {
00283       return this->matrixPtr->frob();
00284     }
00285     static inline Treal frob_diff
00286       (const MatrixGeneral<Treal, Tmatrix>& A,
00287        const MatrixGeneral<Treal, Tmatrix>& B) {
00288       return Tmatrix::frobDiff(*A.matrixPtr, *B.matrixPtr);
00289     }
00290     Treal eucl(Treal const requestedAccuracy,
00291                int maxIter = -1) const;
00292 
00293 
00294     void thresh(Treal const threshold,
00295                 normType const norm);
00296 
00297     inline void frob_thresh(Treal threshold) {
00298       this->matrixPtr->frob_thresh(threshold);
00299     }
00300     Treal eucl_thresh(Treal const threshold);
00301     
00302     inline void gersgorin(Treal& lmin, Treal& lmax) {
00303       this->matrixPtr->gersgorin(lmin, lmax);
00304     }    
00305     static inline Treal trace_ab
00306       (const MatrixGeneral<Treal, Tmatrix>& A,
00307        const MatrixGeneral<Treal, Tmatrix>& B) {
00308       return Tmatrix::trace_ab(*A.matrixPtr, *B.matrixPtr);
00309     }
00310     static inline Treal trace_aTb
00311       (const MatrixGeneral<Treal, Tmatrix>& A,
00312        const MatrixGeneral<Treal, Tmatrix>& B) {
00313       return Tmatrix::trace_aTb(*A.matrixPtr, *B.matrixPtr);
00314     }
00315     inline size_t nnz() const {  /* Note: size_t instead of int here to avoid integer overflow. */
00316       return this->matrixPtr->nnz();
00317     }
00318     inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
00319       return this->matrixPtr->nvalues();
00320     }
00321 
00322     inline void write_to_buffer(void* buffer, const int n_bytes) const {
00323       this->write_to_buffer_base(buffer, n_bytes, matrix_matr);
00324     }
00325     inline void read_from_buffer(void* buffer, const int n_bytes) {
00326       this->read_from_buffer_base(buffer, n_bytes, matrix_matr);
00327     }
00328     
00329     /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */
00331     MatrixGeneral<Treal, Tmatrix>& operator= 
00332       (const XYZ<Treal, 
00333        MatrixGeneral<Treal, Tmatrix>, 
00334        MatrixGeneral<Treal, Tmatrix> >& smm);
00335 
00337     MatrixGeneral<Treal, Tmatrix>& operator= 
00338       (const XY<MatrixGeneral<Treal, Tmatrix>, 
00339        MatrixGeneral<Treal, Tmatrix> >& mm);
00340 
00342     MatrixGeneral<Treal, Tmatrix>& operator+=   
00343       (const XYZ<Treal, 
00344        MatrixGeneral<Treal, Tmatrix>, 
00345        MatrixGeneral<Treal, Tmatrix> >& smm);
00346 
00348     MatrixGeneral<Treal, Tmatrix>& operator= 
00349       (const XYZpUV<Treal, 
00350        MatrixGeneral<Treal, Tmatrix>,
00351        MatrixGeneral<Treal, Tmatrix>,
00352        Treal, 
00353        MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00354 
00356     MatrixGeneral<Treal, Tmatrix>& operator= 
00357       (XpY<MatrixGeneral<Treal, Tmatrix>,
00358        MatrixGeneral<Treal, Tmatrix> > const & mpm);
00360     MatrixGeneral<Treal, Tmatrix>& operator+= 
00361       (MatrixGeneral<Treal, Tmatrix> const & A);
00362     MatrixGeneral<Treal, Tmatrix>& operator-= 
00363       (MatrixGeneral<Treal, Tmatrix> const & A);
00365     MatrixGeneral<Treal, Tmatrix>& operator+= 
00366       (XY<Treal, MatrixGeneral<Treal, Tmatrix> > const & sm);
00367         
00368 
00369     /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
00371     MatrixGeneral<Treal, Tmatrix>& operator= 
00372       (const XYZ<Treal, 
00373        MatrixSymmetric<Treal, Tmatrix>, 
00374        MatrixGeneral<Treal, Tmatrix> >& smm);
00376     MatrixGeneral<Treal, Tmatrix>& operator= 
00377       (const XY<MatrixSymmetric<Treal, Tmatrix>, 
00378        MatrixGeneral<Treal, Tmatrix> >& mm);
00380     MatrixGeneral<Treal, Tmatrix>& operator+=   
00381       (const XYZ<Treal, 
00382        MatrixSymmetric<Treal, Tmatrix>, 
00383        MatrixGeneral<Treal, Tmatrix> >& smm);
00385     MatrixGeneral<Treal, Tmatrix>& operator= 
00386       (const XYZpUV<Treal, 
00387        MatrixSymmetric<Treal, Tmatrix>, 
00388        MatrixGeneral<Treal, Tmatrix>, 
00389        Treal,
00390        MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00392     MatrixGeneral<Treal, Tmatrix>& operator= 
00393       (const XYZ<Treal, 
00394        MatrixGeneral<Treal, Tmatrix>,
00395        MatrixSymmetric<Treal, Tmatrix> >& smm);
00397     MatrixGeneral<Treal, Tmatrix>& operator= 
00398       (const XY<MatrixGeneral<Treal, Tmatrix>,
00399        MatrixSymmetric<Treal, Tmatrix> >& mm);
00401     MatrixGeneral<Treal, Tmatrix>& operator+=   
00402       (const XYZ<Treal, 
00403        MatrixGeneral<Treal, Tmatrix>, 
00404        MatrixSymmetric<Treal, Tmatrix> >& smm);
00406     MatrixGeneral<Treal, Tmatrix>& operator= 
00407       (const XYZpUV<Treal, 
00408        MatrixGeneral<Treal, Tmatrix>, 
00409        MatrixSymmetric<Treal, Tmatrix>, 
00410        Treal,
00411        MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00413     MatrixGeneral<Treal, Tmatrix>& operator= 
00414       (const XYZ<Treal, 
00415        MatrixSymmetric<Treal, Tmatrix>, 
00416        MatrixSymmetric<Treal, Tmatrix> >& smm);
00418     MatrixGeneral<Treal, Tmatrix>& operator= 
00419       (const XY<MatrixSymmetric<Treal, Tmatrix>, 
00420        MatrixSymmetric<Treal, Tmatrix> >& mm);
00422     MatrixGeneral<Treal, Tmatrix>& operator+= 
00423       (const XYZ<Treal, 
00424        MatrixSymmetric<Treal, Tmatrix>, 
00425        MatrixSymmetric<Treal, Tmatrix> >& smm);
00427     MatrixGeneral<Treal, Tmatrix>& operator= 
00428       (const XYZpUV<Treal, 
00429        MatrixSymmetric<Treal, Tmatrix>, 
00430        MatrixSymmetric<Treal, Tmatrix>, 
00431        Treal,
00432        MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00433     
00434     /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */    
00436     MatrixGeneral<Treal, Tmatrix>& operator= 
00437       (const XYZ<Treal, 
00438        MatrixTriangular<Treal, Tmatrix>, 
00439        MatrixGeneral<Treal, Tmatrix> >& smm);
00441     MatrixGeneral<Treal, Tmatrix>& operator= 
00442       (const XYZ<Treal, 
00443        MatrixGeneral<Treal, Tmatrix>, 
00444        MatrixTriangular<Treal, Tmatrix> >& smm);
00445     
00446     void random() {
00447       this->matrixPtr->random();
00448     }
00449 
00450     void randomZeroStructure(Treal probabilityBeingZero) {
00451       this->matrixPtr->randomZeroStructure(probabilityBeingZero);
00452     }
00453 
00454     template<typename TRule>
00455       void setElementsByRule(TRule & rule) {
00456       this->matrixPtr->setElementsByRule(rule);
00457       return;
00458     }
00459 
00460     std::string obj_type_id() const {return "MatrixGeneral";}
00461     protected:
00462     inline void writeToFileProt(std::ofstream & file) const {
00463       this->writeToFileBase(file, matrix_matr);
00464     }
00465     inline void readFromFileProt(std::ifstream & file) {
00466       this->readFromFileBase(file, matrix_matr);
00467     }
00468     private:
00469 
00470   };
00471 
00472   template<typename Treal, typename Tmatrix>
00473     Treal MatrixGeneral<Treal, Tmatrix>::
00474     eucl(Treal const requestedAccuracy,
00475          int maxIter) const {
00476     VectorType guess;
00477     SizesAndBlocks cols;
00478     this->getCols(cols);
00479     guess.resetSizesAndBlocks(cols);
00480     guess.rand();
00481     mat::ATAMatrix<MatrixGeneral<Treal, Tmatrix>, Treal> ata(*this);    
00482     if (maxIter < 0)
00483       maxIter = this->get_nrows() * 100;
00484     arn::LanczosLargestMagnitudeEig
00485       <Treal, ATAMatrix<MatrixGeneral<Treal, Tmatrix>, Treal>, VectorType>
00486       lan(ata, guess, maxIter);
00487     lan.setRelTol( 100 * std::numeric_limits<Treal>::epsilon() );
00488     lan.run();
00489     Treal eVal = 0;
00490     Treal acc = 0;
00491     lan.getLargestMagnitudeEig(eVal, acc);
00492     Interval<Treal> euclInt( sqrt(eVal-acc),
00493                              sqrt(eVal+acc) );    
00494     if ( euclInt.low() < 0 )
00495       euclInt = Interval<Treal>( 0, sqrt(eVal+acc) );
00496     if ( euclInt.length() / 2.0 > requestedAccuracy ) {
00497       std::cout << "req: " << requestedAccuracy
00498                 << "  obt: " << euclInt.length() / 2.0 << std::endl;
00499       throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
00500     }
00501     return euclInt.midPoint();
00502   }
00503 
00504 
00505   template<typename Treal, typename Tmatrix>
00506     void MatrixGeneral<Treal, Tmatrix>::
00507     thresh(Treal const threshold,
00508            normType const norm) {
00509     switch (norm) {
00510     case frobNorm:
00511       this->frob_thresh(threshold);
00512       break;
00513     default:
00514       throw Failure("MatrixGeneral<Treal, Tmatrix>::"
00515                     "thresh(Treal const, "
00516                     "normType const): "
00517                     "Thresholding not imlpemented for selected norm");
00518     }
00519   }
00520 
00521   template<typename Treal, typename Tmatrix>
00522     Treal MatrixGeneral<Treal, Tmatrix>::
00523     eucl_thresh(Treal const threshold) {
00524     EuclTruncationGeneral<MatrixGeneral<Treal, Tmatrix>, Treal> TruncObj( *this );
00525     return TruncObj.run( threshold );    
00526   }
00527 
00528 
00529   /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */
00530     /* C = alpha * op(A) * op(B)  */
00531   template<typename Treal, typename Tmatrix>
00532     MatrixGeneral<Treal, Tmatrix>& 
00533     MatrixGeneral<Treal, Tmatrix>::operator= 
00534     (const XYZ<Treal, 
00535      MatrixGeneral<Treal, Tmatrix>, 
00536      MatrixGeneral<Treal, Tmatrix> >& smm) {
00537     assert(this != &smm.B && this != &smm.C);
00538     this->matrixPtr.haveDataStructureSet(true);
00539     Tmatrix::gemm(smm.tB, smm.tC, smm.A, 
00540                   *smm.B.matrixPtr, *smm.C.matrixPtr, 
00541                   0, *this->matrixPtr);
00542     return *this;
00543   }
00544 
00545     /* C = op(A) * op(B)  */
00546   template<typename Treal, typename Tmatrix>
00547     MatrixGeneral<Treal, Tmatrix>& 
00548     MatrixGeneral<Treal, Tmatrix>::operator= 
00549     (const XY<MatrixGeneral<Treal, Tmatrix>, 
00550      MatrixGeneral<Treal, Tmatrix> >& mm) {
00551     assert(this != &mm.A && this != &mm.B);
00552     Tmatrix::gemm(mm.tA, mm.tB, 1.0, 
00553                   *mm.A.matrixPtr, *mm.B.matrixPtr, 
00554                   0, *this->matrixPtr);
00555     return *this;
00556   }
00557    
00558   /* C += alpha * op(A) * op(B) */
00559   template<typename Treal, typename Tmatrix>
00560     MatrixGeneral<Treal, Tmatrix>& 
00561     MatrixGeneral<Treal, Tmatrix>::operator+=   
00562     (const XYZ<Treal, 
00563      MatrixGeneral<Treal, Tmatrix>, 
00564      MatrixGeneral<Treal, Tmatrix> >& smm) {
00565     assert(this != &smm.B && this != &smm.C);
00566     Tmatrix::gemm(smm.tB, smm.tC, smm.A, 
00567                   *smm.B.matrixPtr, *smm.C.matrixPtr, 
00568                   1, *this->matrixPtr);
00569     return *this;
00570   }
00571   
00572   /* C = alpha * op(A) * op(B) + beta * C  */
00573   template<typename Treal, typename Tmatrix>
00574     MatrixGeneral<Treal, Tmatrix>& 
00575     MatrixGeneral<Treal, Tmatrix>::operator= 
00576     (const XYZpUV<Treal, 
00577      MatrixGeneral<Treal, Tmatrix>, 
00578      MatrixGeneral<Treal, Tmatrix>,
00579      Treal,
00580      MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00581     assert(this != &smmpsm.B && this != &smmpsm.C);
00582     assert(!smmpsm.tE);
00583     if (this == &smmpsm.E)
00584       Tmatrix::gemm(smmpsm.tB, smmpsm.tC, smmpsm.A, 
00585                     *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr, 
00586                     smmpsm.D, *this->matrixPtr);
00587     else
00588       throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00589                     "(const XYZpUV<Treal, MatrixGeneral"
00590                     "<Treal, Tmatrix> >&) :  D = alpha "
00591                     "* op(A) * op(B) + beta * C not supported for C != D");
00592     return *this;
00593   }
00594 
00595 
00596   /* C =  A + B   */
00597   template<typename Treal, typename Tmatrix>
00598     inline MatrixGeneral<Treal, Tmatrix>& 
00599     MatrixGeneral<Treal, Tmatrix>::operator= 
00600     (const XpY<MatrixGeneral<Treal, Tmatrix>,
00601      MatrixGeneral<Treal, Tmatrix> >& mpm) {
00602     assert(this != &mpm.A);
00603     (*this) = mpm.B;
00604     Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
00605     return *this;
00606   }
00607   /* B += A */
00608   template<typename Treal, typename Tmatrix>
00609     inline MatrixGeneral<Treal, Tmatrix>& 
00610     MatrixGeneral<Treal, Tmatrix>::operator+= 
00611     (MatrixGeneral<Treal, Tmatrix> const & A) {
00612     Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
00613     return *this;
00614   }
00615   /* B -= A */
00616   template<typename Treal, typename Tmatrix>
00617     inline MatrixGeneral<Treal, Tmatrix>& 
00618     MatrixGeneral<Treal, Tmatrix>::operator-= 
00619     (MatrixGeneral<Treal, Tmatrix> const & A) {
00620     Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
00621     return *this;
00622   }
00623   /* B += alpha * A */
00624   template<typename Treal, typename Tmatrix>
00625     inline MatrixGeneral<Treal, Tmatrix>& 
00626     MatrixGeneral<Treal, Tmatrix>::operator+= 
00627     (XY<Treal, MatrixGeneral<Treal, Tmatrix> > const & sm) {
00628     assert(!sm.tB);
00629     Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
00630     return *this;
00631   }
00632 
00633     
00634   /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
00635   /* C = alpha * A * B                      : A is symmetric */
00636   template<typename Treal, typename Tmatrix>
00637     MatrixGeneral<Treal, Tmatrix>& 
00638     MatrixGeneral<Treal, Tmatrix>::operator= 
00639     (const XYZ<Treal, 
00640      MatrixSymmetric<Treal, Tmatrix>, 
00641      MatrixGeneral<Treal, Tmatrix> >& smm) {
00642     assert(this != &smm.C);
00643     assert(!smm.tB && !smm.tC);
00644     this->matrixPtr.haveDataStructureSet(true);
00645     Tmatrix::symm('L', 'U', smm.A, 
00646                   *smm.B.matrixPtr, *smm.C.matrixPtr,
00647                   0, *this->matrixPtr);
00648     return *this;
00649   }
00650 
00651   /* C = A * B                      : A is symmetric */
00652   template<typename Treal, typename Tmatrix>
00653     MatrixGeneral<Treal, Tmatrix>& 
00654     MatrixGeneral<Treal, Tmatrix>::operator= 
00655     (const XY<MatrixSymmetric<Treal, Tmatrix>, 
00656      MatrixGeneral<Treal, Tmatrix> >& mm) {
00657     assert(this != &mm.B);
00658     assert(!mm.tB);
00659     Tmatrix::symm('L', 'U', 1.0, 
00660                   *mm.A.matrixPtr, *mm.B.matrixPtr,
00661                   0, *this->matrixPtr);
00662     return *this;
00663   }
00664   
00665   /* C += alpha * A * B                     : A is symmetric */
00666   template<typename Treal, typename Tmatrix>
00667     MatrixGeneral<Treal, Tmatrix>& 
00668     MatrixGeneral<Treal, Tmatrix>::operator+=   
00669     (const XYZ<Treal, 
00670      MatrixSymmetric<Treal, Tmatrix>, 
00671      MatrixGeneral<Treal, Tmatrix> >& smm) {
00672     assert(this != &smm.C);
00673     assert(!smm.tB && !smm.tC);
00674     Tmatrix::symm('L', 'U', smm.A, 
00675                   *smm.B.matrixPtr, *smm.C.matrixPtr,
00676                   1, *this->matrixPtr);
00677     return *this;
00678   }
00679     
00680   /* C = alpha * A * B + beta * C           : A is symmetric */
00681   template<typename Treal, typename Tmatrix>
00682     MatrixGeneral<Treal, Tmatrix>& 
00683     MatrixGeneral<Treal, Tmatrix>::operator= 
00684     (const XYZpUV<Treal, 
00685      MatrixSymmetric<Treal, Tmatrix>, 
00686      MatrixGeneral<Treal, Tmatrix>, 
00687      Treal,
00688      MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00689     assert(this != &smmpsm.C);
00690     assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
00691     if (this == &smmpsm.E)
00692       Tmatrix::symm('L', 'U', smmpsm.A, 
00693                     *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
00694                     smmpsm.D, *this->matrixPtr);
00695     else
00696       throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00697                     "(const XYZpUV<Treal, MatrixGeneral"
00698                     "<Treal, Tmatrix>, MatrixSymmetric<Treal, "
00699                     "Tmatrix>, Treal, MatrixGeneral"
00700                     "<Treal, Tmatrix> >&) "
00701                     ":  D = alpha * A * B + beta * C (with A symmetric)"
00702                     " not supported for C != D");
00703     return *this;
00704   }
00705 
00706   /* C = alpha * B * A                      : A is symmetric */
00707   template<typename Treal, typename Tmatrix>
00708   MatrixGeneral<Treal, Tmatrix>& 
00709     MatrixGeneral<Treal, Tmatrix>::operator= 
00710   (const XYZ<Treal, 
00711    MatrixGeneral<Treal, Tmatrix>, 
00712    MatrixSymmetric<Treal, Tmatrix> >& smm) {
00713     assert(this != &smm.B);
00714     assert(!smm.tB && !smm.tC);
00715     this->matrixPtr.haveDataStructureSet(true);
00716     Tmatrix::symm('R', 'U', smm.A, 
00717                   *smm.C.matrixPtr, *smm.B.matrixPtr,
00718                   0, *this->matrixPtr);
00719     return *this;
00720   }
00721 
00722   /* C = B * A                      : A is symmetric */
00723   template<typename Treal, typename Tmatrix>
00724   MatrixGeneral<Treal, Tmatrix>& 
00725     MatrixGeneral<Treal, Tmatrix>::operator= 
00726   (const XY<MatrixGeneral<Treal, Tmatrix>, 
00727    MatrixSymmetric<Treal, Tmatrix> >& mm) {
00728     assert(this != &mm.A);
00729     assert(!mm.tA && !mm.tB);
00730     Tmatrix::symm('R', 'U', 1.0, 
00731                   *mm.B.matrixPtr, *mm.A.matrixPtr,
00732                   0, *this->matrixPtr);
00733     return *this;
00734   }
00735   
00736   /* C += alpha * B * A                      : A is symmetric */
00737   template<typename Treal, typename Tmatrix>
00738     MatrixGeneral<Treal, Tmatrix>& 
00739     MatrixGeneral<Treal, Tmatrix>::operator+=   
00740     (const XYZ<Treal, 
00741      MatrixGeneral<Treal, Tmatrix>, 
00742      MatrixSymmetric<Treal, Tmatrix> >& smm) {
00743     assert(this != &smm.B);
00744     assert(!smm.tB && !smm.tC);
00745     Tmatrix::symm('R', 'U', smm.A, 
00746                   *smm.C.matrixPtr, *smm.B.matrixPtr,
00747                   1, *this->matrixPtr);
00748     return *this;
00749   }
00750     
00751   /* C = alpha * B * A + beta * C           : A is symmetric */
00752   template<typename Treal, typename Tmatrix>
00753     MatrixGeneral<Treal, Tmatrix>& 
00754     MatrixGeneral<Treal, Tmatrix>::operator= 
00755     (const XYZpUV<Treal, 
00756      MatrixGeneral<Treal, Tmatrix>, 
00757      MatrixSymmetric<Treal, Tmatrix>, 
00758      Treal,
00759      MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00760     assert(this != &smmpsm.B);
00761     assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
00762     if (this == &smmpsm.E)
00763       Tmatrix::symm('R', 'U', smmpsm.A, 
00764                     *smmpsm.C.matrixPtr, *smmpsm.B.matrixPtr,
00765                     smmpsm.D, *this->matrixPtr);
00766     else
00767       throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00768                     "(const XYZpUV<Treal, MatrixSymmetric"
00769                     "<Treal, Tmatrix>, MatrixGeneral<Treal, "
00770                     "Tmatrix>, Treal, MatrixGeneral"
00771                     "<Treal, Tmatrix> >&) "
00772                     ":  D = alpha * B * A + beta * C (with A symmetric)"
00773                     " not supported for C != D");
00774     return *this;
00775   }
00776 
00777 
00779   template<typename Treal, typename Tmatrix>
00780     MatrixGeneral<Treal, Tmatrix>& 
00781     MatrixGeneral<Treal, Tmatrix>::operator= 
00782     (const XYZ<Treal, 
00783      MatrixSymmetric<Treal, Tmatrix>, 
00784      MatrixSymmetric<Treal, Tmatrix> >& smm) {
00785     assert(!smm.tB && !smm.tC);
00786     this->matrixPtr.haveDataStructureSet(true);
00787     Tmatrix::ssmm(smm.A, 
00788                   *smm.B.matrixPtr, 
00789                   *smm.C.matrixPtr,
00790                   0, 
00791                   *this->matrixPtr);
00792     return *this;    
00793   }
00794 
00796   template<typename Treal, typename Tmatrix>
00797     MatrixGeneral<Treal, Tmatrix>& 
00798     MatrixGeneral<Treal, Tmatrix>::operator= 
00799     (const XY<MatrixSymmetric<Treal, Tmatrix>, 
00800      MatrixSymmetric<Treal, Tmatrix> >& mm) {
00801     assert(!mm.tB);
00802     Tmatrix::ssmm(1.0, 
00803                   *mm.A.matrixPtr, 
00804                   *mm.B.matrixPtr,
00805                   0, 
00806                   *this->matrixPtr);
00807     return *this;    
00808   }
00809 
00811   template<typename Treal, typename Tmatrix>
00812     MatrixGeneral<Treal, Tmatrix>& 
00813     MatrixGeneral<Treal, Tmatrix>::operator+= 
00814     (const XYZ<Treal, 
00815      MatrixSymmetric<Treal, Tmatrix>, 
00816      MatrixSymmetric<Treal, Tmatrix> >& smm) {
00817     assert(!smm.tB && !smm.tC);
00818     Tmatrix::ssmm(smm.A, 
00819                   *smm.B.matrixPtr, 
00820                   *smm.C.matrixPtr,
00821                   1, 
00822                   *this->matrixPtr);
00823     return *this;    
00824   }
00825 
00826 
00828   template<typename Treal, typename Tmatrix>
00829     MatrixGeneral<Treal, Tmatrix>& 
00830     MatrixGeneral<Treal, Tmatrix>::operator= 
00831     (const XYZpUV<Treal, 
00832      MatrixSymmetric<Treal, Tmatrix>, 
00833      MatrixSymmetric<Treal, Tmatrix>, 
00834      Treal,
00835      MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00836     assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
00837     if (this == &smmpsm.E)
00838       Tmatrix::ssmm(smmpsm.A, 
00839                     *smmpsm.B.matrixPtr, 
00840                     *smmpsm.C.matrixPtr,
00841                     smmpsm.D, 
00842                     *this->matrixPtr);
00843     else
00844       throw Failure("MatrixGeneral<Treal, Tmatrix>::"
00845                     "operator=(const XYZpUV<"
00846                     "Treal, MatrixSymmetric<Treal, Tmatrix>, "
00847                     "MatrixSymmetric<Treal, Tmatrix>, Treal, "
00848                     "MatrixGeneral<Treal, Tmatrix> >&) "
00849                     ":  D = alpha * A * B + beta * C (with A and B symmetric)"
00850                     " not supported for C != D");
00851     return *this;    
00852   }
00853    
00854 
00855     
00856     /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */
00857     
00858     /* B = alpha * op(A) * B                  : A is upper triangular */
00859     template<typename Treal, typename Tmatrix>
00860       MatrixGeneral<Treal, Tmatrix>& 
00861       MatrixGeneral<Treal, Tmatrix>::operator= 
00862       (const XYZ<Treal, 
00863        MatrixTriangular<Treal, Tmatrix>, 
00864        MatrixGeneral<Treal, Tmatrix> >& smm) {
00865       assert(!smm.tC);
00866       if (this == &smm.C)
00867         Tmatrix::trmm('L', 'U', smm.tB, smm.A, 
00868                       *smm.B.matrixPtr, *this->matrixPtr);
00869       else
00870         throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00871                       "(const XYZ<Treal, MatrixTriangular"
00872                       "<Treal, Tmatrix>, MatrixGeneral<Treal,"
00873                       " Tmatrix> >& : D = alpha * op(A) * B (with"
00874                       " A upper triangular) not supported for B != D");
00875       return *this;
00876     }
00877 
00878 
00879     /* A = alpha * A * op(B)                  : B is upper triangular */
00880     template<typename Treal, typename Tmatrix>
00881       MatrixGeneral<Treal, Tmatrix>& 
00882       MatrixGeneral<Treal, Tmatrix>::operator= 
00883       (const XYZ<Treal, 
00884        MatrixGeneral<Treal, Tmatrix>, 
00885        MatrixTriangular<Treal, Tmatrix> >& smm) {
00886       assert(!smm.tB);
00887       if (this == &smm.B)
00888         Tmatrix::trmm('R', 'U', smm.tC, smm.A, 
00889                       *smm.C.matrixPtr, *this->matrixPtr);
00890       else
00891         throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00892                       "(const XYZ<Treal, MatrixGeneral"
00893                       "<Treal, Tmatrix>, MatrixTriangular<Treal,"
00894                       " Tmatrix> >& : D = alpha * A * op(B) (with"
00895                       " B upper triangular) not supported for A != D");
00896       return *this;
00897     }
00898     
00899 
00900     /******* FUNCTIONS DECLARED OUTSIDE CLASS */
00901     template<typename Treal, typename Tmatrix>
00902       Treal trace(const XYZ<Treal, 
00903                   MatrixGeneral<Treal, Tmatrix>, 
00904                   MatrixGeneral<Treal, Tmatrix> >& smm) {
00905       if ((!smm.tB && !smm.tC) || (smm.tB && smm.tC))
00906         return smm.A * MatrixGeneral<Treal, Tmatrix>::
00907 	  trace_ab(smm.B, smm.C);
00908       else if (smm.tB) 
00909         return smm.A * MatrixGeneral<Treal, Tmatrix>::
00910 	  trace_aTb(smm.B, smm.C);
00911       else
00912         return smm.A * MatrixGeneral<Treal, Tmatrix>::
00913 	  trace_aTb(smm.C, smm.B);
00914     }
00915 
00916 
00917 } /* end namespace mat */
00918 #endif
00919 
00920 

Generated on Mon Sep 17 14:30:40 2012 for ergo by  doxygen 1.4.7