MatrixSymmetric.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_MatrixSymmetric
00037 #define MAT_MatrixSymmetric
00038 
00039 #include <algorithm>
00040 
00041 #include "MatrixBase.h"
00042 #include "Interval.h"
00043 #include "LanczosLargestMagnitudeEig.h"
00044 #include "mat_utils.h"
00045 #include "truncation.h"
00046 
00047 namespace mat {
00065   template<typename Treal, typename Tmatrix>
00066     class MatrixSymmetric : public MatrixBase<Treal, Tmatrix> {
00067   public:
00068     typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType;
00069     typedef Treal real;
00070 
00071     MatrixSymmetric()
00072       :MatrixBase<Treal, Tmatrix>() {} 
00073     explicit MatrixSymmetric(const MatrixSymmetric<Treal, Tmatrix>& symm)
00074       :MatrixBase<Treal, Tmatrix>(symm) {} 
00075     explicit MatrixSymmetric(const XY<Treal, MatrixSymmetric<Treal, Tmatrix> >& sm)
00076       :MatrixBase<Treal, Tmatrix>() { *this = sm.A * sm.B; }
00077     explicit MatrixSymmetric(const MatrixGeneral<Treal, Tmatrix>& matr)
00078       :MatrixBase<Treal, Tmatrix>(matr) {
00079       this->matrixPtr->nosymToSym();
00080     } 
00083 #if 0
00084     template<typename Tfull>
00085       inline void assign_from_full
00086       (Tfull const* const fullmatrix, 
00087        int const totnrows, int const totncols) {
00088        assert(totnrows == totncols);
00089       this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);      
00090       this->matrixPtr->nosym_to_sym();
00091     }    
00092     inline void assign_from_full
00093       (Treal const* const fullmatrix, 
00094        int const totnrows, int const totncols) {
00095       assert(totnrows == totncols);
00096       this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);      
00097       this->matrixPtr->nosym_to_sym();
00098     }
00099 #endif
00100 
00101     inline void assignFromFull
00102       (std::vector<Treal> const & fullMat) {
00103       assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
00104       this->matrixPtr->assignFromFull(fullMat);      
00105       this->matrixPtr->nosymToSym();
00106     }
00107     
00108     inline void assignFromFull
00109       (std::vector<Treal> const & fullMat,
00110        std::vector<int> const & rowPermutation, 
00111        std::vector<int> const & colPermutation) {
00112       assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
00113       std::vector<int> rowind(fullMat.size());
00114       std::vector<int> colind(fullMat.size());
00115       int ind = 0;
00116       for (int col = 0; col < this->get_ncols(); ++col)
00117         for (int row = 0; row < this->get_nrows(); ++row) {
00118           rowind[ind] = row;
00119           colind[ind] = col;
00120           ++ind;
00121         }
00122       this->assign_from_sparse(rowind, 
00123                                colind, 
00124                                fullMat, 
00125                                rowPermutation, 
00126                                colPermutation);
00127       this->matrixPtr->nosymToSym();
00128     }
00129 
00130     inline void fullMatrix(std::vector<Treal> & fullMat) const {
00131       this->matrixPtr->syFullMatrix(fullMat);
00132     }
00139     inline void fullMatrix
00140       (std::vector<Treal> & fullMat,
00141        std::vector<int> const & rowInversePermutation, 
00142        std::vector<int> const & colInversePermutation) const {
00143       std::vector<int> rowind;
00144       std::vector<int> colind; 
00145       std::vector<Treal> values;
00146       get_all_values(rowind, colind, values, 
00147                      rowInversePermutation, 
00148                      colInversePermutation);
00149       fullMat.assign(this->get_nrows()*this->get_ncols(),0);
00150       assert(rowind.size() == values.size());
00151       assert(rowind.size() == colind.size());
00152       for (unsigned int ind = 0; ind < values.size(); ++ind) {
00153         assert(rowind[ind] + this->get_nrows() * colind[ind] < 
00154                this->get_nrows()*this->get_ncols());
00155         fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
00156           values[ind];
00157         fullMat[colind[ind] + this->get_nrows() * rowind[ind]] =
00158           values[ind];
00159       }
00160     }
00167     inline void assign_from_sparse
00168       (std::vector<int> const & rowind, 
00169        std::vector<int> const & colind, 
00170        std::vector<Treal> const & values) {
00171       this->matrixPtr->syAssignFromSparse(rowind, colind, values);
00172     }
00184     inline void assign_from_sparse
00185       (std::vector<int> const & rowind, 
00186        std::vector<int> const & colind, 
00187        std::vector<Treal> const & values, 
00188        SizesAndBlocks const & newRows,
00189        SizesAndBlocks const & newCols) {
00190       this->resetSizesAndBlocks(newRows, newCols);
00191       this->matrixPtr->syAssignFromSparse(rowind, colind, values);
00192     }
00202     inline void assign_from_sparse
00203       (std::vector<int> const & rowind, 
00204        std::vector<int> const & colind, 
00205        std::vector<Treal> const & values, 
00206        std::vector<int> const & rowPermutation, 
00207        std::vector<int> const & colPermutation) {
00208       std::vector<int> newRowind;
00209       std::vector<int> newColind; 
00210       this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
00211                                       colind, colPermutation, newColind);
00212       
00213       this->matrixPtr->syAssignFromSparse(newRowind, newColind, values);
00214     }
00220     inline void assign_from_sparse
00221       (std::vector<int> const & rowind, 
00222        std::vector<int> const & colind, 
00223        std::vector<Treal> const & values, 
00224        SizesAndBlocks const & newRows,
00225        SizesAndBlocks const & newCols,
00226        std::vector<int> const & rowPermutation, 
00227        std::vector<int> const & colPermutation) {
00228       this->resetSizesAndBlocks(newRows, newCols);
00229       assign_from_sparse(rowind, colind, values, 
00230                          rowPermutation, colPermutation);
00231     }
00239     inline void add_values
00240       (std::vector<int> const & rowind, 
00241        std::vector<int> const & colind, 
00242        std::vector<Treal> const & values) {
00243       this->matrixPtr->syAddValues(rowind, colind, values);
00244     }
00245 
00249     inline void add_values
00250       (std::vector<int> const & rowind, 
00251        std::vector<int> const & colind, 
00252        std::vector<Treal> const & values,
00253        std::vector<int> const & rowPermutation, 
00254        std::vector<int> const & colPermutation) {
00255       std::vector<int> newRowind;
00256       std::vector<int> newColind; 
00257       this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
00258                                       colind, colPermutation, newColind);
00259       this->matrixPtr->syAddValues(newRowind, newColind, values);
00260     }
00261 
00262 
00263 
00264     inline void get_values
00265       (std::vector<int> const & rowind, 
00266        std::vector<int> const & colind, 
00267        std::vector<Treal> & values) const {
00268       this->matrixPtr->syGetValues(rowind, colind, values);
00269     }
00277     inline void get_values
00278       (std::vector<int> const & rowind, 
00279        std::vector<int> const & colind, 
00280        std::vector<Treal> & values,
00281        std::vector<int> const & rowPermutation, 
00282        std::vector<int> const & colPermutation) const {
00283       std::vector<int> newRowind;
00284       std::vector<int> newColind; 
00285       this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
00286                                       colind, colPermutation, newColind);
00287       this->matrixPtr->syGetValues(newRowind, newColind, values);
00288     }
00294     inline void get_all_values
00295       (std::vector<int> & rowind, 
00296        std::vector<int> & colind, 
00297        std::vector<Treal> & values) const {
00298       rowind.resize(0);
00299       colind.resize(0);
00300       values.resize(0);
00301       this->matrixPtr->syGetAllValues(rowind, colind, values);
00302     }
00308     inline void get_all_values
00309       (std::vector<int> & rowind, 
00310        std::vector<int> & colind, 
00311        std::vector<Treal> & values,
00312        std::vector<int> const & rowInversePermutation, 
00313        std::vector<int> const & colInversePermutation) const {
00314       std::vector<int> tmpRowind;
00315       std::vector<int> tmpColind; 
00316       tmpRowind.reserve(rowind.capacity());
00317       tmpColind.reserve(colind.capacity());
00318       values.resize(0);
00319       this->matrixPtr->syGetAllValues(tmpRowind, tmpColind, values);
00320       this->getPermutedAndSymmetrized(tmpRowind, rowInversePermutation, rowind,
00321                                       tmpColind, colInversePermutation, colind);
00322     }
00333     MatrixSymmetric<Treal, Tmatrix>& 
00334       operator=(const MatrixSymmetric<Treal, Tmatrix>& symm) {
00335       MatrixBase<Treal, Tmatrix>::operator=(symm);
00336       return *this;
00337     } 
00338     MatrixSymmetric<Treal, Tmatrix>& 
00339       operator=(const MatrixGeneral<Treal, Tmatrix>& matr) {
00340       MatrixBase<Treal, Tmatrix>::operator=(matr);
00341       this->matrixPtr->nosymToSym();
00342       return *this;
00343     } 
00344     inline MatrixSymmetric<Treal, Tmatrix>& operator=(int const k) {
00345       *this->matrixPtr = k;
00346       return *this;
00347     }
00348 
00349     inline Treal frob() const {
00350       return this->matrixPtr->syFrob();
00351     }
00352     Treal mixed_norm(Treal const requestedAccuracy,
00353                      int maxIter = -1) const;
00354     Treal eucl(Treal const requestedAccuracy,
00355                int maxIter = -1) const;
00356 
00357     void quickEuclBounds(Treal & euclLowerBound, 
00358                          Treal & euclUpperBound) const {
00359       Treal frobTmp = frob();
00360       euclLowerBound = frobTmp  / template_blas_sqrt( (Treal)this->get_nrows() );
00361       euclUpperBound = frobTmp;
00362     }
00363 
00364 
00370       static Interval<Treal> 
00371         diff(const MatrixSymmetric<Treal, Tmatrix>& A,
00372              const MatrixSymmetric<Treal, Tmatrix>& B,
00373              normType const norm, 
00374              Treal const requestedAccuracy);
00382       static Interval<Treal> 
00383       diffIfSmall(const MatrixSymmetric<Treal, Tmatrix>& A,
00384                   const MatrixSymmetric<Treal, Tmatrix>& B,
00385                   normType const norm, 
00386                   Treal const requestedAccuracy,
00387                   Treal const maxAbsVal);
00391     static inline Treal frob_diff
00392       (const MatrixSymmetric<Treal, Tmatrix>& A,
00393        const MatrixSymmetric<Treal, Tmatrix>& B) {
00394       return Tmatrix::syFrobDiff(*A.matrixPtr, *B.matrixPtr);
00395     }
00396 
00400     static Treal eucl_diff
00401       (const MatrixSymmetric<Treal, Tmatrix>& A,
00402        const MatrixSymmetric<Treal, Tmatrix>& B,
00403        Treal const requestedAccuracy);
00404     
00408     static Treal mixed_diff
00409       (const MatrixSymmetric<Treal, Tmatrix>& A,
00410        const MatrixSymmetric<Treal, Tmatrix>& B,
00411        Treal const requestedAccuracy);
00412     
00419     static Interval<Treal> euclDiffIfSmall
00420       (const MatrixSymmetric<Treal, Tmatrix>& A,
00421        const MatrixSymmetric<Treal, Tmatrix>& B,
00422        Treal const requestedAccuracy,
00423        Treal const maxAbsVal,
00424        VectorType * const eVecPtr = 0);
00425       
00426 
00437     Treal thresh(Treal const threshold,
00438                  normType const norm);
00439     
00446     inline Treal frob_thresh(Treal const threshold) {
00447       return 2.0 * this->matrixPtr->frob_thresh(threshold / 2);
00448     } 
00449     
00450     Treal eucl_thresh(Treal const threshold,
00451                       MatrixTriangular<Treal, Tmatrix> const * const Zptr = NULL); 
00452     
00453     Treal eucl_element_level_thresh(Treal const threshold); 
00454     
00455     void getSizesAndBlocksForFrobNormMat
00456       ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const;
00457 
00458     Treal mixed_norm_thresh(Treal const threshold);
00459     
00460     void simple_blockwise_frob_thresh(Treal const threshold) {
00461       this->matrixPtr->frobThreshLowestLevel(threshold*threshold, 0);
00462     }
00463     
00464     inline void gersgorin(Treal& lmin, Treal& lmax) const {
00465       this->matrixPtr->sy_gersgorin(lmin, lmax);
00466     }
00467     static inline Treal trace_ab
00468       (const MatrixSymmetric<Treal, Tmatrix>& A,
00469        const MatrixSymmetric<Treal, Tmatrix>& B) {
00470       return Tmatrix::sy_trace_ab(*A.matrixPtr, *B.matrixPtr);
00471     }
00472     inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
00473       return this->matrixPtr->sy_nnz();
00474     }
00475     inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
00476       return this->matrixPtr->sy_nvalues();
00477     }
00478     inline void write_to_buffer(void* buffer, const int n_bytes) const {
00479       this->write_to_buffer_base(buffer, n_bytes, matrix_symm);
00480     }
00481     inline void read_from_buffer(void* buffer, const int n_bytes) {
00482       this->read_from_buffer_base(buffer, n_bytes, matrix_symm);
00483     }
00484 
00485 
00487     MatrixSymmetric<Treal, Tmatrix>& operator= 
00488       (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm);
00490     MatrixSymmetric<Treal, Tmatrix>& operator=
00491       (const XYZpUV<Treal, 
00492        MatrixSymmetric<Treal, Tmatrix>, 
00493        MatrixSymmetric<Treal, Tmatrix>,
00494        Treal, 
00495        MatrixSymmetric<Treal, Tmatrix> >& sm2psm);
00497     MatrixSymmetric<Treal, Tmatrix>& operator=
00498       (const XYZ<Treal, 
00499        MatrixSymmetric<Treal, Tmatrix>, 
00500        MatrixSymmetric<Treal, Tmatrix> >& sm2);
00502     MatrixSymmetric<Treal, Tmatrix>& operator+=
00503       (const XYZ<Treal, 
00504        MatrixSymmetric<Treal, Tmatrix>, 
00505        MatrixSymmetric<Treal, Tmatrix> >& sm2);
00507     MatrixSymmetric<Treal, Tmatrix>& operator= 
00508       (const XYZpUV<Treal, 
00509        MatrixGeneral<Treal, Tmatrix>,
00510        MatrixGeneral<Treal, Tmatrix>,
00511        Treal,
00512        MatrixSymmetric<Treal, Tmatrix> >& smmpsm);
00514     MatrixSymmetric<Treal, Tmatrix>& operator= 
00515       (const XYZ<Treal, 
00516        MatrixGeneral<Treal, Tmatrix>, 
00517        MatrixGeneral<Treal, Tmatrix> >& smm);
00519     MatrixSymmetric<Treal, Tmatrix>& operator+= 
00520       (const XYZ<Treal, 
00521        MatrixGeneral<Treal, Tmatrix>, 
00522        MatrixGeneral<Treal, Tmatrix> >& smm);
00523     
00527     MatrixSymmetric<Treal, Tmatrix>& operator= 
00528       (const XYZ<MatrixTriangular<Treal, Tmatrix>, 
00529        MatrixSymmetric<Treal, Tmatrix>, 
00530        MatrixTriangular<Treal, Tmatrix> >& zaz);
00531     
00536     static void ssmmUpperTriangleOnly(const Treal alpha, 
00537                                       const MatrixSymmetric<Treal, Tmatrix>& A, 
00538                                       const MatrixSymmetric<Treal, Tmatrix>& B, 
00539                                       const Treal beta, 
00540                                       MatrixSymmetric<Treal, Tmatrix>& C);
00541     
00542     
00543     /* Addition */    
00545     MatrixSymmetric<Treal, Tmatrix>& operator= 
00546       (XpY<MatrixSymmetric<Treal, Tmatrix>,
00547        MatrixSymmetric<Treal, Tmatrix> > const & mpm);
00549     MatrixSymmetric<Treal, Tmatrix>& operator= 
00550       (XmY<MatrixSymmetric<Treal, Tmatrix>,
00551        MatrixSymmetric<Treal, Tmatrix> > const & mm);
00553     MatrixSymmetric<Treal, Tmatrix>& operator+= 
00554       (MatrixSymmetric<Treal, Tmatrix> const & A);
00556     MatrixSymmetric<Treal, Tmatrix>& operator-=
00557       (MatrixSymmetric<Treal, Tmatrix> const & A);
00558 
00560     MatrixSymmetric<Treal, Tmatrix>& operator+= 
00561       (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm);
00562 
00564     MatrixSymmetric<Treal, Tmatrix>& operator-= 
00565       (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm);
00566     
00567     template<typename Top>
00568       Treal accumulateWith(Top & op) {
00569       return this->matrixPtr->syAccumulateWith(op);
00570     }
00571 
00572     void random() {
00573       this->matrixPtr->syRandom();
00574     }
00575     
00576     void randomZeroStructure(Treal probabilityBeingZero) {
00577       this->matrixPtr->syRandomZeroStructure(probabilityBeingZero);
00578     }
00579     
00584     template<typename TRule>
00585       void setElementsByRule(TRule & rule) {
00586       this->matrixPtr->sySetElementsByRule(rule);
00587       return;
00588     }
00589 
00592     void transfer( MatrixSymmetric<Treal, Tmatrix> & dest ) {
00593       ValidPtr<Tmatrix>::swap( this->matrixPtr, dest.matrixPtr );
00594       // *this now contains previous content of dest
00595       this->clear(); 
00596     }
00597 
00598     template<typename Tvector>
00599       void matVecProd(Tvector & y, Tvector const & x) const {
00600       Treal const ONE = 1.0;
00601       y = (ONE * (*this) * x);
00602     }
00603     
00604     
00605     std::string obj_type_id() const {return "MatrixSymmetric";}
00606   protected:
00607     inline void writeToFileProt(std::ofstream & file) const {
00608       this->writeToFileBase(file, matrix_symm);
00609     }
00610     inline void readFromFileProt(std::ifstream & file) {
00611       this->readFromFileBase(file, matrix_symm);
00612     }
00613     
00619     static void getPermutedAndSymmetrized
00620       (std::vector<int> const & rowind, 
00621        std::vector<int> const & rowPermutation,
00622        std::vector<int> & newRowind,
00623        std::vector<int> const & colind, 
00624        std::vector<int> const & colPermutation,
00625        std::vector<int> & newColind) {
00626       MatrixBase<Treal, Tmatrix>::
00627 	getPermutedIndexes(rowind, rowPermutation, newRowind);
00628       MatrixBase<Treal, Tmatrix>::
00629 	getPermutedIndexes(colind, colPermutation, newColind);
00630       int tmp;
00631       for (unsigned int i = 0; i < newRowind.size(); ++i) {
00632         if (newRowind[i] > newColind[i]) {
00633           tmp = newRowind[i];
00634           newRowind[i] = newColind[i];
00635           newColind[i] = tmp;
00636         }
00637       }
00638     }
00639   private:
00640   };
00641   
00642   template<typename Treal, typename Tmatrix>
00643     Treal MatrixSymmetric<Treal, Tmatrix>::
00644     mixed_norm(Treal const requestedAccuracy,
00645                int maxIter) const {
00646     // Construct SizesAndBlocks for frobNormMat
00647     SizesAndBlocks rows_new;
00648     SizesAndBlocks cols_new;
00649     this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
00650     // Now we can construct an empty matrix where the Frobenius norms
00651     // of lowest level nonzero submatrices will be stored
00652     MatrixSymmetric<Treal, typename Tmatrix::ElementType> frobNormMat;
00653     frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
00654     frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
00655     return frobNormMat.eucl(requestedAccuracy, maxIter);
00656   }
00657 
00658 
00659   template<typename Treal, typename Tmatrix>
00660     Treal MatrixSymmetric<Treal, Tmatrix>::
00661     eucl(Treal const requestedAccuracy,
00662          int maxIter) const {
00663     assert(requestedAccuracy >= 0);
00664     /* Check if norm is really small, in that case quick return */
00665     Treal frobTmp = this->frob(); 
00666     if (frobTmp < requestedAccuracy)
00667       return (Treal)0.0;
00668     if (maxIter < 0)
00669       maxIter = this->get_nrows() * 100;
00670     VectorType guess;
00671     SizesAndBlocks cols;
00672     this->getCols(cols);
00673     guess.resetSizesAndBlocks(cols);
00674     guess.rand();
00675     // Elias note 2010-03-26: changed this back from "new code" to "old code" to reduce memory usage. 
00676 #if 0 // "new code"
00677     MatrixSymmetric<Treal, Tmatrix> Copy(*this);
00678     Copy.frob_thresh(requestedAccuracy / 2.0);
00679     arn::LanczosLargestMagnitudeEig
00680       <Treal, MatrixSymmetric<Treal, Tmatrix>, VectorType>
00681       lan(Copy, guess, maxIter);
00682     lan.setAbsTol( requestedAccuracy / 2.0 );    
00683 #else // "old code"
00684     arn::LanczosLargestMagnitudeEig
00685       <Treal, MatrixSymmetric<Treal, Tmatrix>, VectorType>
00686       lan(*this, guess, maxIter);
00687     lan.setAbsTol( requestedAccuracy );    
00688 #endif
00689     lan.run();
00690     Treal eVal = 0;
00691     Treal acc = 0;
00692     lan.getLargestMagnitudeEig(eVal, acc);
00693     return template_blas_fabs(eVal);
00694   }
00695 
00696   template<typename Treal, typename Tmatrix>
00697     Interval<Treal> MatrixSymmetric<Treal, Tmatrix>::
00698     diff(const MatrixSymmetric<Treal, Tmatrix>& A,
00699          const MatrixSymmetric<Treal, Tmatrix>& B,
00700          normType const norm, Treal const requestedAccuracy) {
00701     Treal diff;
00702     Treal eNMin;
00703     switch (norm) {
00704     case frobNorm:
00705       diff = frob_diff(A, B);
00706       return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
00707       break;
00708     case euclNorm:
00709       diff = eucl_diff(A, B, requestedAccuracy);
00710       eNMin = diff - requestedAccuracy;
00711       eNMin = eNMin >= 0 ? eNMin : 0; 
00712       return Interval<Treal>(eNMin, diff + requestedAccuracy);
00713       break;
00714     default:
00715       throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
00716                     "diff(const MatrixSymmetric<Treal, Tmatrix>&, "
00717                     "const MatrixSymmetric<Treal, Tmatrix>&, "
00718                     "normType const, Treal): "
00719                     "Diff not implemented for selected norm");
00720     }    
00721   }
00722   
00723 #if 1
00724   template<typename Treal, typename Tmatrix>
00725     Interval<Treal> MatrixSymmetric<Treal, Tmatrix>::
00726     diffIfSmall(const MatrixSymmetric<Treal, Tmatrix>& A,
00727                 const MatrixSymmetric<Treal, Tmatrix>& B,
00728                 normType const norm, 
00729                 Treal const requestedAccuracy,
00730                 Treal const maxAbsVal) {
00731     Treal diff;
00732     switch (norm) {
00733     case frobNorm:
00734       {
00735         diff = frob_diff(A, B);
00736         return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
00737       }
00738       break;
00739     case euclNorm:
00740       return euclDiffIfSmall(A, B, requestedAccuracy, maxAbsVal);
00741       break;
00742     default:
00743       throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
00744                     "diffIfSmall"
00745                     "(const MatrixSymmetric<Treal, Tmatrix>&, "
00746                     "const MatrixSymmetric<Treal, Tmatrix>&, "
00747                     "normType const, Treal const, Treal const): "
00748                     "Diff not implemented for selected norm");
00749     }    
00750   }
00751 
00752 #endif
00753 
00754 
00755   template<typename Treal, typename Tmatrix>
00756     Treal MatrixSymmetric<Treal, Tmatrix>::eucl_diff
00757     (const MatrixSymmetric<Treal, Tmatrix>& A,
00758      const MatrixSymmetric<Treal, Tmatrix>& B,
00759      Treal const requestedAccuracy) {
00760     // DiffMatrix is a lightweight proxy object:
00761     mat::DiffMatrix< MatrixSymmetric<Treal, Tmatrix>, Treal> Diff(A, B);
00762     Treal maxAbsVal = 2 * frob_diff(A,B); 
00763     // Now, maxAbsVal should be larger than the Eucl norm
00764     // Note that mat::euclIfSmall lies outside this class 
00765     Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
00766     Interval<Treal> euclInt = 
00767       mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal); 
00768     return euclInt.midPoint();
00769   }
00770 
00771   template<typename Treal, typename Tmatrix>
00772     Treal MatrixSymmetric<Treal, Tmatrix>::mixed_diff
00773     (const MatrixSymmetric<Treal, Tmatrix>& A,
00774      const MatrixSymmetric<Treal, Tmatrix>& B,
00775      Treal const requestedAccuracy) {
00776     MatrixSymmetric<Treal, typename Tmatrix::ElementType> frobNormMat;
00777     {
00778       SizesAndBlocks rows_new;
00779       SizesAndBlocks cols_new;
00780       A.getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
00781       frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
00782       frobNormMat.getMatrix().syAssignDiffFrobNormsLowestLevel(A.getMatrix(),B.getMatrix());
00783     }
00784     return frobNormMat.eucl(requestedAccuracy);
00785   }
00786 
00787 #if 1
00788   template<typename Treal, typename Tmatrix>
00789     Interval<Treal> MatrixSymmetric<Treal, Tmatrix>::euclDiffIfSmall
00790     (const MatrixSymmetric<Treal, Tmatrix>& A,
00791      const MatrixSymmetric<Treal, Tmatrix>& B,
00792      Treal const requestedAccuracy,
00793      Treal const maxAbsVal,
00794      VectorType * const eVecPtr) {
00795     // DiffMatrix is a lightweight proxy object:
00796     mat::DiffMatrix< MatrixSymmetric<Treal, Tmatrix>, Treal> Diff(A, B);
00797     // Note that this function lies outside this class 
00798     Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
00799     Interval<Treal> tmpInterval = mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtr);
00800     // Emanuel note: Ugly fix to make certain tests pass, we expand
00801     // the interval up to the requested accuracy. Note that larger
00802     // intervals may occur if the norm is not 'small'. It happens that
00803     // Lanczos misconverges to for example the second largest
00804     // eigenvalue. This happens in particular when the first and second
00805     // eigenvalues are very close (of the order of the requested
00806     // accuracy). Expanding the interval makes the largest eigenvalue
00807     // (at least for certain cases) end up inside the interval even
00808     // though Lanczos has misconverged.
00809     if ( tmpInterval.length() < 2*requestedAccuracy )
00810       return Interval<Treal>( tmpInterval.midPoint()-requestedAccuracy, tmpInterval.midPoint()+requestedAccuracy );
00811     return tmpInterval;
00812   }
00813   
00814 #endif
00815 
00816 
00817   template<typename Treal, typename Tmatrix>
00818     Treal MatrixSymmetric<Treal, Tmatrix>::
00819     thresh(Treal const threshold,
00820            normType const norm) {
00821     switch (norm) {
00822     case frobNorm:
00823       return this->frob_thresh(threshold);
00824       break;
00825     case euclNorm:
00826       return this->eucl_thresh(threshold);
00827       break;
00828     case mixedNorm:
00829       return this->mixed_norm_thresh(threshold);
00830       break;
00831     default:
00832       throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
00833                     "thresh(Treal const, "
00834                     "normType const): "
00835                     "Thresholding not implemented for selected norm");
00836     }
00837   }
00838   
00839 #if 1
00840 
00841   template<typename Treal, typename Tmatrix>
00842     Treal MatrixSymmetric<Treal, Tmatrix>::
00843     eucl_thresh(Treal const threshold,
00844                 MatrixTriangular<Treal, Tmatrix> const * const Zptr) {    
00845     if ( Zptr == NULL ) {
00846       EuclTruncationSymm<MatrixSymmetric<Treal, Tmatrix>, Treal> TruncObj(*this);
00847       return TruncObj.run( threshold );
00848     }
00849     EuclTruncationSymmWithZ<MatrixSymmetric<Treal, Tmatrix>, MatrixTriangular<Treal, Tmatrix>, Treal> TruncObj(*this, *Zptr);
00850     return TruncObj.run( threshold );
00851   }
00852   
00853 #endif  
00854 
00855 
00856   template<typename Treal, typename Tmatrix>
00857     void MatrixSymmetric<Treal, Tmatrix>::getSizesAndBlocksForFrobNormMat
00858     ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const {
00859     std::vector<int> rows_block_sizes;
00860     std::vector<int> cols_block_sizes;
00861     int n_rows;
00862     int n_cols;
00863     {
00864       SizesAndBlocks rows;
00865       SizesAndBlocks cols;
00866       this->getRows(rows);
00867       this->getCols(cols);
00868       rows.getBlockSizeVector( rows_block_sizes );
00869       cols.getBlockSizeVector( cols_block_sizes );
00870       rows_block_sizes.pop_back(); // Remove the '1' at the end
00871       cols_block_sizes.pop_back(); // Remove the '1' at the end
00872       n_rows = rows.getNTotalScalars();
00873       n_cols = cols.getNTotalScalars();
00874       int factor_rows = rows_block_sizes[rows_block_sizes.size()-1]; 
00875       int factor_cols = cols_block_sizes[cols_block_sizes.size()-1]; 
00876       for (unsigned int ind = 0; ind < rows_block_sizes.size(); ++ind) 
00877         rows_block_sizes[ind] = rows_block_sizes[ind] / factor_rows;
00878       for (unsigned int ind = 0; ind < cols_block_sizes.size(); ++ind) 
00879         cols_block_sizes[ind] = cols_block_sizes[ind] / factor_cols;
00880       // Now set the number of (scalar) rows and cols, should be equal
00881       // to the number of blocks at the lowest level of the original
00882       // matrix
00883       if (n_rows % factor_rows) 
00884         n_rows = n_rows / factor_rows + 1; 
00885       else
00886         n_rows = n_rows / factor_rows;
00887       if (n_cols % factor_cols) 
00888         n_cols = n_cols / factor_cols + 1; 
00889       else
00890         n_cols = n_cols / factor_cols;
00891     }
00892     rows_new = SizesAndBlocks( rows_block_sizes, n_rows );
00893     cols_new = SizesAndBlocks( cols_block_sizes, n_cols );
00894   }
00895 
00896 
00897   template<typename Treal, typename Tmatrix>
00898     Treal MatrixSymmetric<Treal, Tmatrix>::
00899     mixed_norm_thresh(Treal const threshold) {
00900     assert(threshold >= (Treal)0.0);
00901     if (threshold == (Treal)0.0)
00902       return (Treal)0;
00903     // Construct SizesAndBlocks for frobNormMat
00904     SizesAndBlocks rows_new;
00905     SizesAndBlocks cols_new;
00906     this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
00907     // Now we can construct an empty matrix where the Frobenius norms
00908     // of lowest level nonzero submatrices will be stored
00909     MatrixSymmetric<Treal, typename Tmatrix::ElementType> frobNormMat;
00910     frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
00911     // We want the following step to dominate the mixed_norm truncation (this
00912     // is where Frobenius norms of submatrices are computed, i.e. it
00913     // is here we loop over all matrix elements.)
00914     frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
00915     EuclTruncationSymmElementLevel<MatrixSymmetric<Treal, typename Tmatrix::ElementType>, Treal> TruncObj( frobNormMat );
00916     Treal mixed_norm_result = TruncObj.run( threshold );
00917     frobNormMat.getMatrix().truncateAccordingToSparsityPattern(this->getMatrix());
00918     return mixed_norm_result;
00919   }
00920 
00921 
00922   /* B = alpha * A */
00923   template<typename Treal, typename Tmatrix>
00924     inline MatrixSymmetric<Treal, Tmatrix>& 
00925     MatrixSymmetric<Treal, Tmatrix>::operator= 
00926     (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm) {
00927     assert(!sm.tB);
00928     /* Data structure set by assign - therefore set haveDataStructure to true */
00929     this->matrixPtr.haveDataStructureSet(true); 
00930     this->matrixPtr->assign(sm.A, *sm.B.matrixPtr);
00931     return *this;
00932   }
00933   /* C = alpha * A * A + beta * C          : A and C are symmetric    */
00934   template<typename Treal, typename Tmatrix>
00935     MatrixSymmetric<Treal, Tmatrix>& 
00936     MatrixSymmetric<Treal, Tmatrix>::operator=
00937     (const XYZpUV<Treal, 
00938      MatrixSymmetric<Treal, Tmatrix>, 
00939      MatrixSymmetric<Treal, Tmatrix>,
00940      Treal, 
00941      MatrixSymmetric<Treal, Tmatrix> >& sm2psm) {
00942     assert(this != &sm2psm.B);
00943     if (this == &sm2psm.E && &sm2psm.B == &sm2psm.C) { 
00944       /* Operation is C = alpha * A * A + beta * C */
00945       Tmatrix::sysq('U', 
00946                     sm2psm.A, *sm2psm.B.matrixPtr, 
00947                     sm2psm.D, *this->matrixPtr);
00948       return *this;
00949     }
00950     else /* this != &sm2psm.C || &sm2psm.B != &sm2psm.C */
00951       throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
00952                     "(const XYZpUV<Treal, MatrixSymmetric"
00953                     "<Treal, Tmatrix> >& sm2psm) :  "
00954                     "D = alpha * A * B + beta * C not supported for C != D"
00955                     " and for symmetric matrices not for A != B since this "
00956                     "generally will result in a nonsymmetric matrix");
00957   }
00958 
00959   /* C = alpha * A * A                     : A and C are symmetric    */
00960   template<typename Treal, typename Tmatrix>
00961     MatrixSymmetric<Treal, Tmatrix>& 
00962     MatrixSymmetric<Treal, Tmatrix>::operator=
00963     (const XYZ<Treal, 
00964      MatrixSymmetric<Treal, Tmatrix>, 
00965      MatrixSymmetric<Treal, Tmatrix> >& sm2) {
00966     assert(this != &sm2.B);
00967     if (&sm2.B == &sm2.C) { 
00968       this->matrixPtr.haveDataStructureSet(true); 
00969       Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 0, *this->matrixPtr);
00970       return *this;
00971     }
00972     else
00973       throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
00974                     "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
00975                     " MatrixSymmetric<Treal, Tmatrix> >& sm2) :  "
00976                     "Operation C = alpha * A * B with only symmetric "
00977                     "matrices not supported for A != B");
00978   }
00979   
00980   /* C += alpha * A * A                    : A and C are symmetric    */
00981   template<typename Treal, typename Tmatrix>
00982     MatrixSymmetric<Treal, Tmatrix>& 
00983     MatrixSymmetric<Treal, Tmatrix>::operator+=
00984     (const XYZ<Treal, 
00985      MatrixSymmetric<Treal, Tmatrix>, 
00986      MatrixSymmetric<Treal, Tmatrix> >& sm2) {
00987     assert(this != &sm2.B);
00988     if (&sm2.B == &sm2.C) { 
00989       Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 1, *this->matrixPtr);
00990       return *this;
00991     }
00992     else
00993       throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
00994                     "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
00995                     " MatrixSymmetric<Treal, Tmatrix> >& sm2) :  "
00996                     "Operation C += alpha * A * B with only symmetric "
00997                     "matrices not supported for A != B");
00998   }
00999 
01000   /* C = alpha * A * transpose(A) + beta * C        : C is symmetric  */
01001   template<typename Treal, typename Tmatrix>
01002     MatrixSymmetric<Treal, Tmatrix>& 
01003     MatrixSymmetric<Treal, Tmatrix>::operator= 
01004     (const XYZpUV<Treal, 
01005      MatrixGeneral<Treal, Tmatrix>, 
01006      MatrixGeneral<Treal, Tmatrix>, 
01007      Treal,
01008      MatrixSymmetric<Treal, Tmatrix> >& smmpsm) {
01009     if (this == &smmpsm.E)
01010       if (&smmpsm.B == &smmpsm.C)
01011         if (!smmpsm.tB && smmpsm.tC) {
01012           Tmatrix::syrk('U', false, 
01013                         smmpsm.A, *smmpsm.B.matrixPtr, 
01014                         smmpsm.D, *this->matrixPtr);
01015         }
01016         else 
01017           throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01018                         "(const XYZpUV<Treal, MatrixGeneral"
01019                         "<Treal, Tmatrix>, "
01020                         "MatrixGeneral<Treal, Tmatrix>, Treal, "
01021                         "MatrixSymmetric<Treal, Tmatrix> >&) : "
01022                         "C = alpha * A' * A + beta * C, not implemented"
01023                         " only C = alpha * A * A' + beta * C");
01024       else
01025         throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01026                       "(const XYZpUV<"
01027                       "Treal, MatrixGeneral<Treal, Tmatrix>, "
01028                       "MatrixGeneral<Treal, Tmatrix>, Treal, "
01029                       "MatrixSymmetric<Treal, Tmatrix> >&) : "
01030                       "You are trying to call C = alpha * A * A' + beta * C"
01031                       " with  A and A' being different objects");
01032     else
01033       throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01034                     "(const XYZpUV<"
01035                     "Treal, MatrixGeneral<Treal, Tmatrix>, "
01036                     "MatrixGeneral<Treal, Tmatrix>, Treal, "
01037                     "MatrixSymmetric<Treal, Tmatrix> >&) :  "
01038                     "D = alpha * A * A' + beta * C not supported for C != D");
01039     return *this;
01040   }
01041 
01042   /* C = alpha * A * transpose(A)                   : C is symmetric  */
01043   template<typename Treal, typename Tmatrix>
01044     MatrixSymmetric<Treal, Tmatrix>& 
01045     MatrixSymmetric<Treal, Tmatrix>::operator= 
01046     (const XYZ<Treal, 
01047      MatrixGeneral<Treal, Tmatrix>, 
01048      MatrixGeneral<Treal, Tmatrix> >& smm) {
01049     if (&smm.B == &smm.C)
01050       if (!smm.tB && smm.tC) {
01051         Tmatrix::syrk('U', false, 
01052                       smm.A, *smm.B.matrixPtr, 
01053                       0, *this->matrixPtr);
01054       }
01055       else 
01056         throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01057                       "(const XYZ<"
01058                       "Treal, MatrixGeneral<Treal, Tmatrix>, "
01059                       "MatrixGeneral<Treal, Tmatrix> >&) : "
01060                       "C = alpha * A' * A, not implemented "
01061                       "only C = alpha * A * A'");
01062     else
01063       throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01064                     "(const XYZ<"
01065                     "Treal, MatrixGeneral<Treal, Tmatrix>, "
01066                     "MatrixGeneral<Treal, Tmatrix> >&) : "
01067                     "You are trying to call C = alpha * A * A' "
01068                     "with A and A' being different objects");
01069     return *this;
01070   }
01071   
01072   /* C += alpha * A * transpose(A)                   : C is symmetric  */
01073   template<typename Treal, typename Tmatrix>
01074     MatrixSymmetric<Treal, Tmatrix>& 
01075     MatrixSymmetric<Treal, Tmatrix>::operator+= 
01076     (const XYZ<Treal, 
01077      MatrixGeneral<Treal, Tmatrix>, 
01078      MatrixGeneral<Treal, Tmatrix> >& smm) {
01079     if (&smm.B == &smm.C)
01080       if (!smm.tB && smm.tC) {
01081         Tmatrix::syrk('U', false, 
01082                       smm.A, *smm.B.matrixPtr, 
01083                       1, *this->matrixPtr);
01084       }
01085       else 
01086         throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
01087                       "(const XYZ<"
01088                       "Treal, MatrixGeneral<Treal, Tmatrix>, "
01089                       "MatrixGeneral<Treal, Tmatrix> >&) : "
01090                       "C += alpha * A' * A, not implemented "
01091                       "only C += alpha * A * A'");
01092     else
01093       throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
01094                     "(const XYZ<"
01095                     "Treal, MatrixGeneral<Treal, Tmatrix>, "
01096                     "MatrixGeneral<Treal, Tmatrix> >&) : "
01097                     "You are trying to call C += alpha * A * A' "
01098                     "with  A and A' being different objects");
01099     return *this;
01100   }
01101   
01102 #if 1
01103   /* A = op1(Z) * A * op2(Z)   : Z is upper triangular and A is symmetric */
01104   /* Either op1() or op2() is the transpose operation. */
01105   template<typename Treal, typename Tmatrix>
01106     MatrixSymmetric<Treal, Tmatrix>& 
01107     MatrixSymmetric<Treal, Tmatrix>::operator= 
01108     (const XYZ<MatrixTriangular<Treal, Tmatrix>, 
01109      MatrixSymmetric<Treal, Tmatrix>, 
01110      MatrixTriangular<Treal, Tmatrix> >& zaz) {
01111     if (this == &zaz.B) {
01112       if (&zaz.A == &zaz.C) {
01113         if (zaz.tA && !zaz.tC) {
01114           Tmatrix::trsytriplemm('R', *zaz.A.matrixPtr, *this->matrixPtr); 
01115         }
01116         else if (!zaz.tA && zaz.tC) {
01117           Tmatrix::trsytriplemm('L', *zaz.A.matrixPtr, *this->matrixPtr); 
01118         }
01119         else
01120           throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator=" 
01121                         "(const XYZ<MatrixTriangular<Treal, Tmatrix>," 
01122                         "MatrixSymmetric<Treal, Tmatrix>," 
01123                         "MatrixTriangular<Treal, Tmatrix> >&) : "
01124                         "A = op1(Z) * A * op2(Z) : Either op1 xor op2 must be "
01125                         "the transpose operation.");
01126       }
01127       else
01128         throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator=" 
01129                       "(const XYZ<MatrixTriangular<Treal, Tmatrix>," 
01130                       "MatrixSymmetric<Treal, Tmatrix>," 
01131                       "MatrixTriangular<Treal, Tmatrix> >&) : "
01132                       "A = op1(Z1) * A * op2(Z2) : Z1 and Z2 must be the same "
01133                       "object");
01134     }
01135     else
01136       throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator=" 
01137                     "(const XYZ<MatrixTriangular<Treal, Tmatrix>," 
01138                     "MatrixSymmetric<Treal, Tmatrix>," 
01139                     "MatrixTriangular<Treal, Tmatrix> >&) : "
01140                     "C = op1(Z) * A * op2(Z) : A and C must be the same "
01141                     "object");
01142     return *this;
01143   }
01144 
01145 #endif
01146 
01147 
01152   template<typename Treal, typename Tmatrix>
01153     void MatrixSymmetric<Treal, Tmatrix>::
01154     ssmmUpperTriangleOnly(const Treal alpha, 
01155                           const MatrixSymmetric<Treal, Tmatrix>& A, 
01156                           const MatrixSymmetric<Treal, Tmatrix>& B, 
01157                           const Treal beta, 
01158                           MatrixSymmetric<Treal, Tmatrix>& C) {
01159     Tmatrix::ssmm_upper_tr_only(alpha, *A.matrixPtr, *B.matrixPtr,
01160                                 beta, *C.matrixPtr);
01161   }
01162 
01163 
01164 
01165   /* Addition */  
01166   /* C =  A + B   */
01167   template<typename Treal, typename Tmatrix>
01168     inline MatrixSymmetric<Treal, Tmatrix>& 
01169     MatrixSymmetric<Treal, Tmatrix>::operator= 
01170     (const XpY<MatrixSymmetric<Treal, Tmatrix>,
01171      MatrixSymmetric<Treal, Tmatrix> >& mpm) {
01172     assert(this != &mpm.A);
01173     (*this) = mpm.B;
01174     Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
01175     return *this;
01176   }
01177   /* C =  A - B   */
01178   template<typename Treal, typename Tmatrix>
01179     inline MatrixSymmetric<Treal, Tmatrix>& 
01180     MatrixSymmetric<Treal, Tmatrix>::operator= 
01181     (const XmY<MatrixSymmetric<Treal, Tmatrix>,
01182      MatrixSymmetric<Treal, Tmatrix> >& mmm) {
01183     assert(this != &mmm.B);
01184     (*this) = mmm.A;
01185     Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
01186     return *this;
01187   }
01188   /* B += A */
01189   template<typename Treal, typename Tmatrix>
01190     inline MatrixSymmetric<Treal, Tmatrix>& 
01191     MatrixSymmetric<Treal, Tmatrix>::operator+= 
01192     (MatrixSymmetric<Treal, Tmatrix> const & A) {
01193     Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
01194     return *this;
01195   }
01196   /* B -= A */
01197   template<typename Treal, typename Tmatrix>
01198     inline MatrixSymmetric<Treal, Tmatrix>& 
01199     MatrixSymmetric<Treal, Tmatrix>::operator-=
01200     (MatrixSymmetric<Treal, Tmatrix> const & A) {
01201     Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
01202     return *this;
01203   }
01204   /* B += alpha * A */
01205   template<typename Treal, typename Tmatrix>
01206     inline MatrixSymmetric<Treal, Tmatrix>& 
01207     MatrixSymmetric<Treal, Tmatrix>::operator+= 
01208     (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm) {
01209     assert(!sm.tB);
01210     Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
01211     return *this;
01212   }
01213 
01214   /* B -= alpha * A */
01215   template<typename Treal, typename Tmatrix>
01216     inline MatrixSymmetric<Treal, Tmatrix>& 
01217     MatrixSymmetric<Treal, Tmatrix>::operator-= 
01218     (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm) {
01219     assert(!sm.tB);
01220     Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
01221     return *this;
01222   }
01223 
01228   template<typename Treal, typename Tmatrix, typename Top>
01229     Treal accumulate(MatrixSymmetric<Treal, Tmatrix> & A, 
01230                      Top & op) {
01231     return A.accumulateWith(op);
01232   }
01233 
01234 
01235   
01236 } /* end namespace mat */
01237 #endif
01238 

Generated on Wed Nov 21 09:32:11 2012 for ergo by  doxygen 1.4.7