1 /* Ergo, version 3.8, a program for linear scaling electronic structure 2 * calculations. 3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 4 * and Anastasia Kruchinina. 5 * 6 * This program is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. 18 * 19 * Primary academic reference: 20 * Ergo: An open-source program for linear-scaling electronic structure 21 * calculations, 22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 23 * Kruchinina, 24 * SoftwareX 7, 107 (2018), 25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 26 * 27 * For further information about Ergo, see <http://www.ergoscf.org>. 28 */ 29 30 /** @file MatrixGeneral.h General matrix class 31 * 32 * Copyright(c) Emanuel Rubensson 2006 33 * 34 * @author Emanuel Rubensson @a responsible @a author 35 * @date January 2006 36 * 37 */ 38 #ifndef MAT_MATRIXGENERAL 39 #define MAT_MATRIXGENERAL 40 #include "MatrixBase.h" 41 namespace mat { 42 /** Normal matrix 43 * 44 * This class belongs to the matrix API 45 * 46 * 47 * Treal: Type for real numbers 48 * 49 * Tmatrix: The matrix class 50 * 51 * Tperm: Permutation used in the matrix class 52 * 53 * @see MatrixBase 54 * @see MatrixSymmetric 55 * @see MatrixTriangular 56 * 57 */ 58 template<typename Treal, typename Tmatrix> 59 class MatrixGeneral : public MatrixBase<Treal, Tmatrix> { 60 public: 61 typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType; 62 MatrixGeneral()63 MatrixGeneral() 64 :MatrixBase<Treal, Tmatrix>() {} /**< Default constructor */ MatrixGeneral(const MatrixGeneral<Treal,Tmatrix> & matr)65 explicit MatrixGeneral(const MatrixGeneral<Treal, Tmatrix>& matr) 66 :MatrixBase<Treal, Tmatrix>(matr) {} /**< Copy constructor */ MatrixGeneral(const MatrixSymmetric<Treal,Tmatrix> & symm)67 explicit MatrixGeneral(const MatrixSymmetric<Treal, Tmatrix>& symm) 68 :MatrixBase<Treal, Tmatrix>(symm) { 69 this->matrixPtr->symToNosym(); 70 }/**< Copy from symmetric matrix constructor */ MatrixGeneral(const MatrixTriangular<Treal,Tmatrix> & triang)71 explicit MatrixGeneral(const MatrixTriangular<Treal, Tmatrix>& triang) 72 :MatrixBase<Treal, Tmatrix>(triang) {} 73 /**< Copy from triangular matrix constructor */ 74 75 #if 0 76 template<typename Tfull> 77 inline void assign_from_full 78 (Tfull const* const fullmatrix, 79 const int totnrows, const int totncols) { 80 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols); 81 } 82 inline void assign_from_full 83 (Treal const* const fullmatrix, 84 const int totnrows, const int totncols) { 85 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols); 86 } 87 #endif 88 assignFromFull(std::vector<Treal> const & fullMat)89 inline void assignFromFull 90 (std::vector<Treal> const & fullMat) { 91 assert((int)fullMat.size() == this->get_nrows() * this->get_ncols()); 92 this->matrixPtr->assignFromFull(fullMat); 93 } 94 fullMatrix(std::vector<Treal> & fullMat)95 inline void fullMatrix(std::vector<Treal> & fullMat) const { 96 this->matrixPtr->fullMatrix(fullMat); 97 } 98 fullMatrix(std::vector<Treal> & fullMat,std::vector<int> const & rowInversePermutation,std::vector<int> const & colInversePermutation)99 inline void fullMatrix 100 (std::vector<Treal> & fullMat, 101 std::vector<int> const & rowInversePermutation, 102 std::vector<int> const & colInversePermutation) const { 103 std::vector<int> rowind; 104 std::vector<int> colind; 105 std::vector<Treal> values; 106 get_all_values(rowind, colind, values, 107 rowInversePermutation, 108 colInversePermutation); 109 fullMat.assign(this->get_nrows()*this->get_ncols(),0); 110 for (unsigned int ind = 0; ind < values.size(); ++ind) 111 fullMat[rowind[ind] + this->get_nrows() * colind[ind]] = 112 values[ind]; 113 } 114 /**< Save matrix as full matrix. 115 * Whole matrix is written in columnwise order. 116 * Both lower and upper triangle. 117 * Permutation is used. 118 */ 119 120 assign_from_sparse(std::vector<int> const & rowind,std::vector<int> const & colind,std::vector<Treal> const & values,SizesAndBlocks const & newRows,SizesAndBlocks const & newCols)121 inline void assign_from_sparse 122 (std::vector<int> const & rowind, 123 std::vector<int> const & colind, 124 std::vector<Treal> const & values, 125 SizesAndBlocks const & newRows, 126 SizesAndBlocks const & newCols) { 127 this->resetSizesAndBlocks(newRows, newCols); 128 this->matrixPtr->assignFromSparse(rowind, colind, values); 129 } 130 /**< Assign from sparse matrix given by three arrays. 131 * The arrays contain row indices, column indices and values. 132 * The indices start at zero. 133 * nval is the length of the three arrays. 134 * @warning All indexing start at zero. 135 */ 136 assign_from_sparse(std::vector<int> const & rowind,std::vector<int> const & colind,std::vector<Treal> const & values,std::vector<int> const & rowPermutation,std::vector<int> const & colPermutation)137 inline void assign_from_sparse 138 (std::vector<int> const & rowind, 139 std::vector<int> const & colind, 140 std::vector<Treal> const & values, 141 std::vector<int> const & rowPermutation, 142 std::vector<int> const & colPermutation) { 143 std::vector<int> newRowind; 144 std::vector<int> newColind; 145 MatrixBase<Treal, Tmatrix>:: 146 getPermutedIndexes(rowind, rowPermutation, newRowind); 147 MatrixBase<Treal, Tmatrix>:: 148 getPermutedIndexes(colind, colPermutation, newColind); 149 this->matrixPtr->assignFromSparse(newRowind, newColind, values); 150 } 151 /**< Same as above, except taking two additional arguments 152 * specifying the permutation of rows and columns. 153 * Also assuming that sizes and blocks are already known 154 */ 155 assign_from_sparse(std::vector<int> const & rowind,std::vector<int> const & colind,std::vector<Treal> const & values,SizesAndBlocks const & newRows,SizesAndBlocks const & newCols,std::vector<int> const & rowPermutation,std::vector<int> const & colPermutation)156 inline void assign_from_sparse 157 (std::vector<int> const & rowind, 158 std::vector<int> const & colind, 159 std::vector<Treal> const & values, 160 SizesAndBlocks const & newRows, 161 SizesAndBlocks const & newCols, 162 std::vector<int> const & rowPermutation, 163 std::vector<int> const & colPermutation) { 164 this->resetSizesAndBlocks(newRows, newCols); 165 assign_from_sparse(rowind, colind, values, 166 rowPermutation, colPermutation); 167 } 168 /**< Same as above, except not assuming that sizes and blocks are set. 169 */ 170 171 get_values(std::vector<int> const & rowind,std::vector<int> const & colind,std::vector<Treal> & values)172 inline void get_values 173 (std::vector<int> const & rowind, 174 std::vector<int> const & colind, 175 std::vector<Treal> & values) const { 176 this->matrixPtr->getValues(rowind, colind, values); 177 } 178 /**< Get values given by row and column index lists. 179 * Input arrays contain row and column indices. 180 * The output array contains values for the given indices. 181 * nval is the length of the three arrays. 182 * @warning All indexing start at zero. 183 */ 184 get_values(std::vector<int> const & rowind,std::vector<int> const & colind,std::vector<Treal> & values,std::vector<int> const & rowPermutation,std::vector<int> const & colPermutation)185 inline void get_values 186 (std::vector<int> const & rowind, 187 std::vector<int> const & colind, 188 std::vector<Treal> & values, 189 std::vector<int> const & rowPermutation, 190 std::vector<int> const & colPermutation) const { 191 std::vector<int> newRowind; 192 std::vector<int> newColind; 193 MatrixBase<Treal, Tmatrix>:: 194 getPermutedIndexes(rowind, rowPermutation, newRowind); 195 MatrixBase<Treal, Tmatrix>:: 196 getPermutedIndexes(colind, colPermutation, newColind); 197 this->matrixPtr->getValues(newRowind, newColind, values); 198 } 199 /**< Same as above, except taking two additional arguments 200 * specifying the permutation of rows and columns. 201 */ 202 get_all_values(std::vector<int> & rowind,std::vector<int> & colind,std::vector<Treal> & values)203 inline void get_all_values 204 (std::vector<int> & rowind, 205 std::vector<int> & colind, 206 std::vector<Treal> & values) const { 207 rowind.resize(0); 208 colind.resize(0); 209 values.resize(0); 210 this->matrixPtr->getAllValues(rowind, colind, values); 211 } 212 /**< Get all values and corresponding row and column index lists, 213 * in matrix. 214 * nval is the length of the three arrays and is preferably 215 * computed with nvalues() before hand. 216 * Returns the number of values. 217 * @see nvalues() 218 * @warning All indexing start at zero. 219 */ 220 get_all_values(std::vector<int> & rowind,std::vector<int> & colind,std::vector<Treal> & values,std::vector<int> const & rowInversePermutation,std::vector<int> const & colInversePermutation)221 inline void get_all_values 222 (std::vector<int> & rowind, 223 std::vector<int> & colind, 224 std::vector<Treal> & values, 225 std::vector<int> const & rowInversePermutation, 226 std::vector<int> const & colInversePermutation) const { 227 std::vector<int> tmpRowind; 228 std::vector<int> tmpColind; 229 tmpRowind.reserve(rowind.capacity()); 230 tmpColind.reserve(colind.capacity()); 231 values.resize(0); 232 this->matrixPtr->getAllValues(tmpRowind, tmpColind, values); 233 MatrixBase<Treal, Tmatrix>:: 234 getPermutedIndexes(tmpRowind, rowInversePermutation, rowind); 235 MatrixBase<Treal, Tmatrix>:: 236 getPermutedIndexes(tmpColind, colInversePermutation, colind); 237 238 } 239 /**< Same as above, except taking two additional arguments 240 * specifying the permutation of rows and columns. 241 * Note, however, that this permutation is the inverse 242 * permutation compared to the permutations provided in the 243 * functions "assign_from_sparse", "add_values", and "get_values" 244 * @warning permutation is inverse compared to other functions 245 */ 246 247 248 #if 0 249 inline void fullmatrix(Treal* const full, 250 const int totnrows, const int totncols) const { 251 this->matrixPtr->fullmatrix(full, totnrows, totncols); 252 } 253 #endif 254 MatrixGeneral<Treal, Tmatrix>& 255 operator=(const MatrixGeneral<Treal, Tmatrix>& mat) { 256 MatrixBase<Treal, Tmatrix>::operator=(mat); 257 return *this; 258 } 259 inline MatrixGeneral<Treal, Tmatrix>& 260 operator=(const Xtrans<MatrixGeneral<Treal, Tmatrix> >& mt) { 261 if (mt.tA) 262 MatrixBase<Treal, Tmatrix>::operator=(transpose(mt.A)); 263 else 264 MatrixBase<Treal, Tmatrix>::operator=(mt.A); 265 return *this; 266 } 267 268 MatrixGeneral<Treal, Tmatrix>& 269 operator=(const MatrixSymmetric<Treal, Tmatrix>& symm) { 270 MatrixBase<Treal, Tmatrix>::operator=(symm); 271 this->matrixPtr->symToNosym(); 272 return *this; 273 } 274 MatrixGeneral<Treal, Tmatrix>& 275 operator=(const MatrixTriangular<Treal, Tmatrix>& triang) { 276 MatrixBase<Treal, Tmatrix>::operator=(triang); 277 return *this; 278 } 279 280 inline MatrixGeneral<Treal, Tmatrix>& operator=(int const k) { 281 *this->matrixPtr = k; 282 return *this; 283 } frob()284 inline Treal frob() const { 285 return this->matrixPtr->frob(); 286 } frob_diff(const MatrixGeneral<Treal,Tmatrix> & A,const MatrixGeneral<Treal,Tmatrix> & B)287 static inline Treal frob_diff 288 (const MatrixGeneral<Treal, Tmatrix>& A, 289 const MatrixGeneral<Treal, Tmatrix>& B) { 290 return Tmatrix::frobDiff(*A.matrixPtr, *B.matrixPtr); 291 } 292 Treal eucl(Treal const requestedAccuracy, 293 int maxIter = -1) const; 294 295 296 void thresh(Treal const threshold, 297 normType const norm); 298 frob_thresh(Treal threshold)299 inline void frob_thresh(Treal threshold) { 300 this->matrixPtr->frob_thresh(threshold); 301 } 302 Treal eucl_thresh(Treal const threshold); 303 gershgorin(Treal & lmin,Treal & lmax)304 inline void gershgorin(Treal& lmin, Treal& lmax) { 305 this->matrixPtr->gershgorin(lmin, lmax); 306 } trace_ab(const MatrixGeneral<Treal,Tmatrix> & A,const MatrixGeneral<Treal,Tmatrix> & B)307 static inline Treal trace_ab 308 (const MatrixGeneral<Treal, Tmatrix>& A, 309 const MatrixGeneral<Treal, Tmatrix>& B) { 310 return Tmatrix::trace_ab(*A.matrixPtr, *B.matrixPtr); 311 } trace_aTb(const MatrixGeneral<Treal,Tmatrix> & A,const MatrixGeneral<Treal,Tmatrix> & B)312 static inline Treal trace_aTb 313 (const MatrixGeneral<Treal, Tmatrix>& A, 314 const MatrixGeneral<Treal, Tmatrix>& B) { 315 return Tmatrix::trace_aTb(*A.matrixPtr, *B.matrixPtr); 316 } nnz()317 inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */ 318 return this->matrixPtr->nnz(); 319 } nvalues()320 inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */ 321 return this->matrixPtr->nvalues(); 322 } 323 write_to_buffer(void * buffer,const int n_bytes)324 inline void write_to_buffer(void* buffer, const int n_bytes) const { 325 this->write_to_buffer_base(buffer, n_bytes, matrix_matr); 326 } read_from_buffer(void * buffer,const int n_bytes)327 inline void read_from_buffer(void* buffer, const int n_bytes) { 328 this->read_from_buffer_base(buffer, n_bytes, matrix_matr); 329 } 330 331 /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */ 332 /** C = alpha * op(A) * op(B) */ 333 MatrixGeneral<Treal, Tmatrix>& operator= 334 (const XYZ<Treal, 335 MatrixGeneral<Treal, Tmatrix>, 336 MatrixGeneral<Treal, Tmatrix> >& smm); 337 338 /** C = op(A) * op(B) */ 339 MatrixGeneral<Treal, Tmatrix>& operator= 340 (const XY<MatrixGeneral<Treal, Tmatrix>, 341 MatrixGeneral<Treal, Tmatrix> >& mm); 342 343 /** C += alpha * op(A) * op(B) */ 344 MatrixGeneral<Treal, Tmatrix>& operator+= 345 (const XYZ<Treal, 346 MatrixGeneral<Treal, Tmatrix>, 347 MatrixGeneral<Treal, Tmatrix> >& smm); 348 349 /** C = alpha * op(A) * op(B) + beta * C */ 350 MatrixGeneral<Treal, Tmatrix>& operator= 351 (const XYZpUV<Treal, 352 MatrixGeneral<Treal, Tmatrix>, 353 MatrixGeneral<Treal, Tmatrix>, 354 Treal, 355 MatrixGeneral<Treal, Tmatrix> >& smmpsm); 356 357 /** C = A + B */ 358 MatrixGeneral<Treal, Tmatrix>& operator= 359 (XpY<MatrixGeneral<Treal, Tmatrix>, 360 MatrixGeneral<Treal, Tmatrix> > const & mpm); 361 /** B += A */ 362 MatrixGeneral<Treal, Tmatrix>& operator+= 363 (MatrixGeneral<Treal, Tmatrix> const & A); 364 MatrixGeneral<Treal, Tmatrix>& operator-= 365 (MatrixGeneral<Treal, Tmatrix> const & A); 366 /** B += alpha * A */ 367 MatrixGeneral<Treal, Tmatrix>& operator+= 368 (XY<Treal, MatrixGeneral<Treal, Tmatrix> > const & sm); 369 370 371 /* OPERATIONS INVOLVING SYMMETRIC MATRICES */ 372 /** C = alpha * A * B : A is symmetric */ 373 MatrixGeneral<Treal, Tmatrix>& operator= 374 (const XYZ<Treal, 375 MatrixSymmetric<Treal, Tmatrix>, 376 MatrixGeneral<Treal, Tmatrix> >& smm); 377 /** C = A * B : A is symmetric */ 378 MatrixGeneral<Treal, Tmatrix>& operator= 379 (const XY<MatrixSymmetric<Treal, Tmatrix>, 380 MatrixGeneral<Treal, Tmatrix> >& mm); 381 /** C += alpha * A * B : A is symmetric */ 382 MatrixGeneral<Treal, Tmatrix>& operator+= 383 (const XYZ<Treal, 384 MatrixSymmetric<Treal, Tmatrix>, 385 MatrixGeneral<Treal, Tmatrix> >& smm); 386 /** C = alpha * A * B + beta * C : A is symmetric */ 387 MatrixGeneral<Treal, Tmatrix>& operator= 388 (const XYZpUV<Treal, 389 MatrixSymmetric<Treal, Tmatrix>, 390 MatrixGeneral<Treal, Tmatrix>, 391 Treal, 392 MatrixGeneral<Treal, Tmatrix> >& smmpsm); 393 /** C = alpha * B * A : A is symmetric */ 394 MatrixGeneral<Treal, Tmatrix>& operator= 395 (const XYZ<Treal, 396 MatrixGeneral<Treal, Tmatrix>, 397 MatrixSymmetric<Treal, Tmatrix> >& smm); 398 /** C = B * A : A is symmetric */ 399 MatrixGeneral<Treal, Tmatrix>& operator= 400 (const XY<MatrixGeneral<Treal, Tmatrix>, 401 MatrixSymmetric<Treal, Tmatrix> >& mm); 402 /** C += alpha * B * A : A is symmetric */ 403 MatrixGeneral<Treal, Tmatrix>& operator+= 404 (const XYZ<Treal, 405 MatrixGeneral<Treal, Tmatrix>, 406 MatrixSymmetric<Treal, Tmatrix> >& smm); 407 /** C = alpha * B * A + beta * C : A is symmetric */ 408 MatrixGeneral<Treal, Tmatrix>& operator= 409 (const XYZpUV<Treal, 410 MatrixGeneral<Treal, Tmatrix>, 411 MatrixSymmetric<Treal, Tmatrix>, 412 Treal, 413 MatrixGeneral<Treal, Tmatrix> >& smmpsm); 414 /** C = alpha * A * B : A and B are symmetric */ 415 MatrixGeneral<Treal, Tmatrix>& operator= 416 (const XYZ<Treal, 417 MatrixSymmetric<Treal, Tmatrix>, 418 MatrixSymmetric<Treal, Tmatrix> >& smm); 419 /** C = A * B : A and B are symmetric */ 420 MatrixGeneral<Treal, Tmatrix>& operator= 421 (const XY<MatrixSymmetric<Treal, Tmatrix>, 422 MatrixSymmetric<Treal, Tmatrix> >& mm); 423 /** C += alpha * A * B : A and B are symmetric */ 424 MatrixGeneral<Treal, Tmatrix>& operator+= 425 (const XYZ<Treal, 426 MatrixSymmetric<Treal, Tmatrix>, 427 MatrixSymmetric<Treal, Tmatrix> >& smm); 428 /** C = alpha * A * B + beta * C : A and B are symmetric */ 429 MatrixGeneral<Treal, Tmatrix>& operator= 430 (const XYZpUV<Treal, 431 MatrixSymmetric<Treal, Tmatrix>, 432 MatrixSymmetric<Treal, Tmatrix>, 433 Treal, 434 MatrixGeneral<Treal, Tmatrix> >& smmpsm); 435 436 /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */ 437 /** B = alpha * op(A) * B : A is upper triangular */ 438 MatrixGeneral<Treal, Tmatrix>& operator= 439 (const XYZ<Treal, 440 MatrixTriangular<Treal, Tmatrix>, 441 MatrixGeneral<Treal, Tmatrix> >& smm); 442 /** B = alpha * B * op(A) : A is upper triangular */ 443 MatrixGeneral<Treal, Tmatrix>& operator= 444 (const XYZ<Treal, 445 MatrixGeneral<Treal, Tmatrix>, 446 MatrixTriangular<Treal, Tmatrix> >& smm); 447 random()448 void random() { 449 this->matrixPtr->random(); 450 } 451 randomZeroStructure(Treal probabilityBeingZero)452 void randomZeroStructure(Treal probabilityBeingZero) { 453 this->matrixPtr->randomZeroStructure(probabilityBeingZero); 454 } 455 456 template<typename TRule> setElementsByRule(TRule & rule)457 void setElementsByRule(TRule & rule) { 458 this->matrixPtr->setElementsByRule(rule); 459 return; 460 } 461 obj_type_id()462 std::string obj_type_id() const {return "MatrixGeneral";} 463 protected: writeToFileProt(std::ofstream & file)464 inline void writeToFileProt(std::ofstream & file) const { 465 this->writeToFileBase(file, matrix_matr); 466 } readFromFileProt(std::ifstream & file)467 inline void readFromFileProt(std::ifstream & file) { 468 this->readFromFileBase(file, matrix_matr); 469 } 470 private: 471 472 }; 473 474 template<typename Treal, typename Tmatrix> 475 Treal MatrixGeneral<Treal, Tmatrix>:: eucl(Treal const requestedAccuracy,int maxIter)476 eucl(Treal const requestedAccuracy, 477 int maxIter) const { 478 VectorType guess; 479 SizesAndBlocks cols; 480 this->getCols(cols); 481 guess.resetSizesAndBlocks(cols); 482 guess.rand(); 483 mat::ATAMatrix<MatrixGeneral<Treal, Tmatrix>, Treal> ata(*this); 484 if (maxIter < 0) 485 maxIter = this->get_nrows() * 100; 486 arn::LanczosLargestMagnitudeEig 487 <Treal, ATAMatrix<MatrixGeneral<Treal, Tmatrix>, Treal>, VectorType> 488 lan(ata, guess, maxIter); 489 lan.setRelTol( 100 * mat::getMachineEpsilon<Treal>() ); 490 lan.run(); 491 Treal eVal = 0; 492 Treal acc = 0; 493 lan.getLargestMagnitudeEig(eVal, acc); 494 Interval<Treal> euclInt( sqrt(eVal-acc), 495 sqrt(eVal+acc) ); 496 if ( euclInt.low() < 0 ) 497 euclInt = Interval<Treal>( 0, sqrt(eVal+acc) ); 498 if ( euclInt.length() / 2.0 > requestedAccuracy ) { 499 std::cout << "req: " << requestedAccuracy 500 << " obt: " << euclInt.length() / 2.0 << std::endl; 501 throw std::runtime_error("Desired accuracy not obtained in Lanczos."); 502 } 503 return euclInt.midPoint(); 504 } 505 506 507 template<typename Treal, typename Tmatrix> 508 void MatrixGeneral<Treal, Tmatrix>:: thresh(Treal const threshold,normType const norm)509 thresh(Treal const threshold, 510 normType const norm) { 511 switch (norm) { 512 case frobNorm: 513 this->frob_thresh(threshold); 514 break; 515 default: 516 throw Failure("MatrixGeneral<Treal, Tmatrix>::" 517 "thresh(Treal const, " 518 "normType const): " 519 "Thresholding not imlpemented for selected norm"); 520 } 521 } 522 523 template<typename Treal, typename Tmatrix> 524 Treal MatrixGeneral<Treal, Tmatrix>:: eucl_thresh(Treal const threshold)525 eucl_thresh(Treal const threshold) { 526 EuclTruncationGeneral<MatrixGeneral<Treal, Tmatrix>, Treal> TruncObj( *this ); 527 return TruncObj.run( threshold ); 528 } 529 530 531 /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */ 532 /* C = alpha * op(A) * op(B) */ 533 template<typename Treal, typename Tmatrix> 534 MatrixGeneral<Treal, Tmatrix>& 535 MatrixGeneral<Treal, Tmatrix>::operator= 536 (const XYZ<Treal, 537 MatrixGeneral<Treal, Tmatrix>, 538 MatrixGeneral<Treal, Tmatrix> >& smm) { 539 assert(this != &smm.B && this != &smm.C); 540 this->matrixPtr.haveDataStructureSet(true); 541 Tmatrix::gemm(smm.tB, smm.tC, smm.A, 542 *smm.B.matrixPtr, *smm.C.matrixPtr, 543 0, *this->matrixPtr); 544 return *this; 545 } 546 547 /* C = op(A) * op(B) */ 548 template<typename Treal, typename Tmatrix> 549 MatrixGeneral<Treal, Tmatrix>& 550 MatrixGeneral<Treal, Tmatrix>::operator= 551 (const XY<MatrixGeneral<Treal, Tmatrix>, 552 MatrixGeneral<Treal, Tmatrix> >& mm) { 553 assert(this != &mm.A && this != &mm.B); 554 Tmatrix::gemm(mm.tA, mm.tB, 1.0, 555 *mm.A.matrixPtr, *mm.B.matrixPtr, 556 0, *this->matrixPtr); 557 return *this; 558 } 559 560 /* C += alpha * op(A) * op(B) */ 561 template<typename Treal, typename Tmatrix> 562 MatrixGeneral<Treal, Tmatrix>& 563 MatrixGeneral<Treal, Tmatrix>::operator+= 564 (const XYZ<Treal, 565 MatrixGeneral<Treal, Tmatrix>, 566 MatrixGeneral<Treal, Tmatrix> >& smm) { 567 assert(this != &smm.B && this != &smm.C); 568 Tmatrix::gemm(smm.tB, smm.tC, smm.A, 569 *smm.B.matrixPtr, *smm.C.matrixPtr, 570 1, *this->matrixPtr); 571 return *this; 572 } 573 574 /* C = alpha * op(A) * op(B) + beta * C */ 575 template<typename Treal, typename Tmatrix> 576 MatrixGeneral<Treal, Tmatrix>& 577 MatrixGeneral<Treal, Tmatrix>::operator= 578 (const XYZpUV<Treal, 579 MatrixGeneral<Treal, Tmatrix>, 580 MatrixGeneral<Treal, Tmatrix>, 581 Treal, 582 MatrixGeneral<Treal, Tmatrix> >& smmpsm) { 583 assert(this != &smmpsm.B && this != &smmpsm.C); 584 assert(!smmpsm.tE); 585 if (this == &smmpsm.E) 586 Tmatrix::gemm(smmpsm.tB, smmpsm.tC, smmpsm.A, 587 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr, 588 smmpsm.D, *this->matrixPtr); 589 else 590 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator=" 591 "(const XYZpUV<Treal, MatrixGeneral" 592 "<Treal, Tmatrix> >&) : D = alpha " 593 "* op(A) * op(B) + beta * C not supported for C != D"); 594 return *this; 595 } 596 597 598 /* C = A + B */ 599 template<typename Treal, typename Tmatrix> 600 inline MatrixGeneral<Treal, Tmatrix>& 601 MatrixGeneral<Treal, Tmatrix>::operator= 602 (const XpY<MatrixGeneral<Treal, Tmatrix>, 603 MatrixGeneral<Treal, Tmatrix> >& mpm) { 604 assert(this != &mpm.A); 605 (*this) = mpm.B; 606 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr); 607 return *this; 608 } 609 /* B += A */ 610 template<typename Treal, typename Tmatrix> 611 inline MatrixGeneral<Treal, Tmatrix>& 612 MatrixGeneral<Treal, Tmatrix>::operator+= 613 (MatrixGeneral<Treal, Tmatrix> const & A) { 614 Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr); 615 return *this; 616 } 617 /* B -= A */ 618 template<typename Treal, typename Tmatrix> 619 inline MatrixGeneral<Treal, Tmatrix>& 620 MatrixGeneral<Treal, Tmatrix>::operator-= 621 (MatrixGeneral<Treal, Tmatrix> const & A) { 622 Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr); 623 return *this; 624 } 625 /* B += alpha * A */ 626 template<typename Treal, typename Tmatrix> 627 inline MatrixGeneral<Treal, Tmatrix>& 628 MatrixGeneral<Treal, Tmatrix>::operator+= 629 (XY<Treal, MatrixGeneral<Treal, Tmatrix> > const & sm) { 630 assert(!sm.tB); 631 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr); 632 return *this; 633 } 634 635 636 /* OPERATIONS INVOLVING SYMMETRIC MATRICES */ 637 /* C = alpha * A * B : A is symmetric */ 638 template<typename Treal, typename Tmatrix> 639 MatrixGeneral<Treal, Tmatrix>& 640 MatrixGeneral<Treal, Tmatrix>::operator= 641 (const XYZ<Treal, 642 MatrixSymmetric<Treal, Tmatrix>, 643 MatrixGeneral<Treal, Tmatrix> >& smm) { 644 assert(this != &smm.C); 645 assert(!smm.tB && !smm.tC); 646 this->matrixPtr.haveDataStructureSet(true); 647 Tmatrix::symm('L', 'U', smm.A, 648 *smm.B.matrixPtr, *smm.C.matrixPtr, 649 0, *this->matrixPtr); 650 return *this; 651 } 652 653 /* C = A * B : A is symmetric */ 654 template<typename Treal, typename Tmatrix> 655 MatrixGeneral<Treal, Tmatrix>& 656 MatrixGeneral<Treal, Tmatrix>::operator= 657 (const XY<MatrixSymmetric<Treal, Tmatrix>, 658 MatrixGeneral<Treal, Tmatrix> >& mm) { 659 assert(this != &mm.B); 660 assert(!mm.tB); 661 Tmatrix::symm('L', 'U', 1.0, 662 *mm.A.matrixPtr, *mm.B.matrixPtr, 663 0, *this->matrixPtr); 664 return *this; 665 } 666 667 /* C += alpha * A * B : A is symmetric */ 668 template<typename Treal, typename Tmatrix> 669 MatrixGeneral<Treal, Tmatrix>& 670 MatrixGeneral<Treal, Tmatrix>::operator+= 671 (const XYZ<Treal, 672 MatrixSymmetric<Treal, Tmatrix>, 673 MatrixGeneral<Treal, Tmatrix> >& smm) { 674 assert(this != &smm.C); 675 assert(!smm.tB && !smm.tC); 676 Tmatrix::symm('L', 'U', smm.A, 677 *smm.B.matrixPtr, *smm.C.matrixPtr, 678 1, *this->matrixPtr); 679 return *this; 680 } 681 682 /* C = alpha * A * B + beta * C : A is symmetric */ 683 template<typename Treal, typename Tmatrix> 684 MatrixGeneral<Treal, Tmatrix>& 685 MatrixGeneral<Treal, Tmatrix>::operator= 686 (const XYZpUV<Treal, 687 MatrixSymmetric<Treal, Tmatrix>, 688 MatrixGeneral<Treal, Tmatrix>, 689 Treal, 690 MatrixGeneral<Treal, Tmatrix> >& smmpsm) { 691 assert(this != &smmpsm.C); 692 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE); 693 if (this == &smmpsm.E) 694 Tmatrix::symm('L', 'U', smmpsm.A, 695 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr, 696 smmpsm.D, *this->matrixPtr); 697 else 698 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator=" 699 "(const XYZpUV<Treal, MatrixGeneral" 700 "<Treal, Tmatrix>, MatrixSymmetric<Treal, " 701 "Tmatrix>, Treal, MatrixGeneral" 702 "<Treal, Tmatrix> >&) " 703 ": D = alpha * A * B + beta * C (with A symmetric)" 704 " not supported for C != D"); 705 return *this; 706 } 707 708 /* C = alpha * B * A : A is symmetric */ 709 template<typename Treal, typename Tmatrix> 710 MatrixGeneral<Treal, Tmatrix>& 711 MatrixGeneral<Treal, Tmatrix>::operator= 712 (const XYZ<Treal, 713 MatrixGeneral<Treal, Tmatrix>, 714 MatrixSymmetric<Treal, Tmatrix> >& smm) { 715 assert(this != &smm.B); 716 assert(!smm.tB && !smm.tC); 717 this->matrixPtr.haveDataStructureSet(true); 718 Tmatrix::symm('R', 'U', smm.A, 719 *smm.C.matrixPtr, *smm.B.matrixPtr, 720 0, *this->matrixPtr); 721 return *this; 722 } 723 724 /* C = B * A : A is symmetric */ 725 template<typename Treal, typename Tmatrix> 726 MatrixGeneral<Treal, Tmatrix>& 727 MatrixGeneral<Treal, Tmatrix>::operator= 728 (const XY<MatrixGeneral<Treal, Tmatrix>, 729 MatrixSymmetric<Treal, Tmatrix> >& mm) { 730 assert(this != &mm.A); 731 assert(!mm.tA && !mm.tB); 732 Tmatrix::symm('R', 'U', 1.0, 733 *mm.B.matrixPtr, *mm.A.matrixPtr, 734 0, *this->matrixPtr); 735 return *this; 736 } 737 738 /* C += alpha * B * A : A is symmetric */ 739 template<typename Treal, typename Tmatrix> 740 MatrixGeneral<Treal, Tmatrix>& 741 MatrixGeneral<Treal, Tmatrix>::operator+= 742 (const XYZ<Treal, 743 MatrixGeneral<Treal, Tmatrix>, 744 MatrixSymmetric<Treal, Tmatrix> >& smm) { 745 assert(this != &smm.B); 746 assert(!smm.tB && !smm.tC); 747 Tmatrix::symm('R', 'U', smm.A, 748 *smm.C.matrixPtr, *smm.B.matrixPtr, 749 1, *this->matrixPtr); 750 return *this; 751 } 752 753 /* C = alpha * B * A + beta * C : A is symmetric */ 754 template<typename Treal, typename Tmatrix> 755 MatrixGeneral<Treal, Tmatrix>& 756 MatrixGeneral<Treal, Tmatrix>::operator= 757 (const XYZpUV<Treal, 758 MatrixGeneral<Treal, Tmatrix>, 759 MatrixSymmetric<Treal, Tmatrix>, 760 Treal, 761 MatrixGeneral<Treal, Tmatrix> >& smmpsm) { 762 assert(this != &smmpsm.B); 763 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE); 764 if (this == &smmpsm.E) 765 Tmatrix::symm('R', 'U', smmpsm.A, 766 *smmpsm.C.matrixPtr, *smmpsm.B.matrixPtr, 767 smmpsm.D, *this->matrixPtr); 768 else 769 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator=" 770 "(const XYZpUV<Treal, MatrixSymmetric" 771 "<Treal, Tmatrix>, MatrixGeneral<Treal, " 772 "Tmatrix>, Treal, MatrixGeneral" 773 "<Treal, Tmatrix> >&) " 774 ": D = alpha * B * A + beta * C (with A symmetric)" 775 " not supported for C != D"); 776 return *this; 777 } 778 779 780 /** C = alpha * A * B : A and B are symmetric */ 781 template<typename Treal, typename Tmatrix> 782 MatrixGeneral<Treal, Tmatrix>& 783 MatrixGeneral<Treal, Tmatrix>::operator= 784 (const XYZ<Treal, 785 MatrixSymmetric<Treal, Tmatrix>, 786 MatrixSymmetric<Treal, Tmatrix> >& smm) { 787 assert(!smm.tB && !smm.tC); 788 this->matrixPtr.haveDataStructureSet(true); 789 Tmatrix::ssmm(smm.A, 790 *smm.B.matrixPtr, 791 *smm.C.matrixPtr, 792 0, 793 *this->matrixPtr); 794 return *this; 795 } 796 797 /** C = A * B : A and B are symmetric */ 798 template<typename Treal, typename Tmatrix> 799 MatrixGeneral<Treal, Tmatrix>& 800 MatrixGeneral<Treal, Tmatrix>::operator= 801 (const XY<MatrixSymmetric<Treal, Tmatrix>, 802 MatrixSymmetric<Treal, Tmatrix> >& mm) { 803 assert(!mm.tB); 804 Tmatrix::ssmm(1.0, 805 *mm.A.matrixPtr, 806 *mm.B.matrixPtr, 807 0, 808 *this->matrixPtr); 809 return *this; 810 } 811 812 /** C += alpha * A * B : A and B are symmetric */ 813 template<typename Treal, typename Tmatrix> 814 MatrixGeneral<Treal, Tmatrix>& 815 MatrixGeneral<Treal, Tmatrix>::operator+= 816 (const XYZ<Treal, 817 MatrixSymmetric<Treal, Tmatrix>, 818 MatrixSymmetric<Treal, Tmatrix> >& smm) { 819 assert(!smm.tB && !smm.tC); 820 Tmatrix::ssmm(smm.A, 821 *smm.B.matrixPtr, 822 *smm.C.matrixPtr, 823 1, 824 *this->matrixPtr); 825 return *this; 826 } 827 828 829 /** C = alpha * A * B + beta * C : A and B are symmetric */ 830 template<typename Treal, typename Tmatrix> 831 MatrixGeneral<Treal, Tmatrix>& 832 MatrixGeneral<Treal, Tmatrix>::operator= 833 (const XYZpUV<Treal, 834 MatrixSymmetric<Treal, Tmatrix>, 835 MatrixSymmetric<Treal, Tmatrix>, 836 Treal, 837 MatrixGeneral<Treal, Tmatrix> >& smmpsm) { 838 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE); 839 if (this == &smmpsm.E) 840 Tmatrix::ssmm(smmpsm.A, 841 *smmpsm.B.matrixPtr, 842 *smmpsm.C.matrixPtr, 843 smmpsm.D, 844 *this->matrixPtr); 845 else 846 throw Failure("MatrixGeneral<Treal, Tmatrix>::" 847 "operator=(const XYZpUV<" 848 "Treal, MatrixSymmetric<Treal, Tmatrix>, " 849 "MatrixSymmetric<Treal, Tmatrix>, Treal, " 850 "MatrixGeneral<Treal, Tmatrix> >&) " 851 ": D = alpha * A * B + beta * C (with A and B symmetric)" 852 " not supported for C != D"); 853 return *this; 854 } 855 856 857 858 /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */ 859 860 /* B = alpha * op(A) * B : A is upper triangular */ 861 template<typename Treal, typename Tmatrix> 862 MatrixGeneral<Treal, Tmatrix>& 863 MatrixGeneral<Treal, Tmatrix>::operator= 864 (const XYZ<Treal, 865 MatrixTriangular<Treal, Tmatrix>, 866 MatrixGeneral<Treal, Tmatrix> >& smm) { 867 assert(!smm.tC); 868 if (this == &smm.C) 869 Tmatrix::trmm('L', 'U', smm.tB, smm.A, 870 *smm.B.matrixPtr, *this->matrixPtr); 871 else 872 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator=" 873 "(const XYZ<Treal, MatrixTriangular" 874 "<Treal, Tmatrix>, MatrixGeneral<Treal," 875 " Tmatrix> >& : D = alpha * op(A) * B (with" 876 " A upper triangular) not supported for B != D"); 877 return *this; 878 } 879 880 881 /* A = alpha * A * op(B) : B is upper triangular */ 882 template<typename Treal, typename Tmatrix> 883 MatrixGeneral<Treal, Tmatrix>& 884 MatrixGeneral<Treal, Tmatrix>::operator= 885 (const XYZ<Treal, 886 MatrixGeneral<Treal, Tmatrix>, 887 MatrixTriangular<Treal, Tmatrix> >& smm) { 888 assert(!smm.tB); 889 if (this == &smm.B) 890 Tmatrix::trmm('R', 'U', smm.tC, smm.A, 891 *smm.C.matrixPtr, *this->matrixPtr); 892 else 893 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator=" 894 "(const XYZ<Treal, MatrixGeneral" 895 "<Treal, Tmatrix>, MatrixTriangular<Treal," 896 " Tmatrix> >& : D = alpha * A * op(B) (with" 897 " B upper triangular) not supported for A != D"); 898 return *this; 899 } 900 901 902 /******* FUNCTIONS DECLARED OUTSIDE CLASS */ 903 template<typename Treal, typename Tmatrix> trace(const XYZ<Treal,MatrixGeneral<Treal,Tmatrix>,MatrixGeneral<Treal,Tmatrix>> & smm)904 Treal trace(const XYZ<Treal, 905 MatrixGeneral<Treal, Tmatrix>, 906 MatrixGeneral<Treal, Tmatrix> >& smm) { 907 if ((!smm.tB && !smm.tC) || (smm.tB && smm.tC)) 908 return smm.A * MatrixGeneral<Treal, Tmatrix>:: 909 trace_ab(smm.B, smm.C); 910 else if (smm.tB) 911 return smm.A * MatrixGeneral<Treal, Tmatrix>:: 912 trace_aTb(smm.B, smm.C); 913 else 914 return smm.A * MatrixGeneral<Treal, Tmatrix>:: 915 trace_aTb(smm.C, smm.B); 916 } 917 918 919 } /* end namespace mat */ 920 #endif 921 922 923