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