ergo
|
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