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 00058 #ifndef MAT_MATRIX 00059 #define MAT_MATRIX 00060 00061 #include <math.h> 00062 #include <cstdlib> 00063 #include <algorithm> 00064 00065 #include "MatrixHierarchicBase.h" 00066 #include "matrix_proxy.h" 00067 #include "gblas.h" 00068 #include "sort.h" 00069 //#define MAT_NAIVE 00070 00071 namespace mat{ 00072 enum side {left, right}; 00073 enum inchversion {unstable, stable}; 00074 template<class Treal, class Telement> 00075 class Vector; 00076 template<typename Tperm> 00077 struct AccessMap; 00078 00079 class SingletonForTimings { 00080 private: 00081 double accumulatedTime; 00082 public: 00083 void reset() { accumulatedTime = 0; } 00084 double getAccumulatedTime() { return accumulatedTime; } 00085 void addTime(double timeToAdd) { 00086 #ifdef _OPENMP 00087 #pragma omp critical 00088 #endif 00089 { 00090 accumulatedTime += timeToAdd; 00091 } 00092 } 00093 static SingletonForTimings & instance() { 00094 static SingletonForTimings theInstance; 00095 return theInstance; 00096 } 00097 private: 00098 SingletonForTimings(const SingletonForTimings & other); 00099 SingletonForTimings() : accumulatedTime(0) { } 00100 }; 00101 00102 00111 template<class Treal, class Telement = Treal> 00112 class Matrix: public MatrixHierarchicBase<Treal, Telement> { 00113 public: 00114 typedef Telement ElementType; 00115 typedef Vector<Treal, typename ElementType::VectorType> VectorType; 00116 00117 friend class Vector<Treal, Telement>; 00118 Matrix():MatrixHierarchicBase<Treal, Telement>(){} 00119 00120 00121 void allocate() { 00122 assert(!this->is_empty()); 00123 assert(this->is_zero()); 00124 //#define MAT_USE_ALLOC_TIMER 00125 #ifdef MAT_USE_ALLOC_TIMER 00126 mat::Time theTimer; theTimer.tic(); 00127 #endif 00128 this->elements = new Telement[this->nElements()]; 00129 #ifdef MAT_USE_ALLOC_TIMER 00130 SingletonForTimings::instance().addTime(theTimer.toc()); 00131 #endif 00132 SizesAndBlocks colSAB; 00133 SizesAndBlocks rowSAB; 00134 for (int col = 0; col < this->cols.getNBlocks(); col++) { 00135 colSAB = this->cols.getSizesAndBlocksForLowerLevel(col); 00136 for (int row = 0; row < this->rows.getNBlocks(); row++) { 00137 /* This could be improved - precompute rowSAB as well as colSAB */ 00138 rowSAB = this->rows.getSizesAndBlocksForLowerLevel(row); 00139 (*this)(row,col).resetRows(rowSAB); 00140 (*this)(row,col).resetCols(colSAB); 00141 } 00142 } 00143 } 00144 00145 /* Full matrix assigns etc */ 00146 void assignFromFull(std::vector<Treal> const & fullMat); 00147 void fullMatrix(std::vector<Treal> & fullMat) const; 00148 void syFullMatrix(std::vector<Treal> & fullMat) const; 00149 void syUpTriFullMatrix(std::vector<Treal> & fullMat) const; 00150 00151 /* Sparse matrix assigns etc */ 00152 void assignFromSparse(std::vector<int> const & rowind, 00153 std::vector<int> const & colind, 00154 std::vector<Treal> const & values); 00155 void assignFromSparse(std::vector<int> const & rowind, 00156 std::vector<int> const & colind, 00157 std::vector<Treal> const & values, 00158 std::vector<int> const & indexes); 00159 /* Adds values (+=) to elements */ 00160 void addValues(std::vector<int> const & rowind, 00161 std::vector<int> const & colind, 00162 std::vector<Treal> const & values); 00163 void addValues(std::vector<int> const & rowind, 00164 std::vector<int> const & colind, 00165 std::vector<Treal> const & values, 00166 std::vector<int> const & indexes); 00167 00168 void syAssignFromSparse(std::vector<int> const & rowind, 00169 std::vector<int> const & colind, 00170 std::vector<Treal> const & values); 00171 00172 void syAddValues(std::vector<int> const & rowind, 00173 std::vector<int> const & colind, 00174 std::vector<Treal> const & values); 00175 00176 void getValues(std::vector<int> const & rowind, 00177 std::vector<int> const & colind, 00178 std::vector<Treal> & values) const; 00179 void getValues(std::vector<int> const & rowind, 00180 std::vector<int> const & colind, 00181 std::vector<Treal> &, 00182 std::vector<int> const & indexes) const; 00183 void syGetValues(std::vector<int> const & rowind, 00184 std::vector<int> const & colind, 00185 std::vector<Treal> & values) const; 00186 void getAllValues(std::vector<int> & rowind, 00187 std::vector<int> & colind, 00188 std::vector<Treal> &) const; 00189 void syGetAllValues(std::vector<int> & rowind, 00190 std::vector<int> & colind, 00191 std::vector<Treal> &) const; 00192 00193 00194 Matrix<Treal, Telement>& 00195 operator=(const Matrix<Treal, Telement>& mat) { 00196 MatrixHierarchicBase<Treal, Telement>::operator=(mat); 00197 return *this; 00198 } 00199 00200 00201 void clear(); 00202 ~Matrix() { 00203 clear(); 00204 } 00205 void writeToFile(std::ofstream & file) const; 00206 void readFromFile(std::ifstream & file); 00207 00208 void random(); 00209 void syRandom(); 00213 void randomZeroStructure(Treal probabilityBeingZero); 00214 void syRandomZeroStructure(Treal probabilityBeingZero); 00215 00216 template<typename TRule> 00217 void setElementsByRule(TRule & rule); 00218 template<typename TRule> 00219 void sySetElementsByRule(TRule & rule); 00220 template<typename TRule> 00221 void trSetElementsByRule(TRule & rule) { 00222 // Setting elements for triangular is the same as for symmetric 00223 sySetElementsByRule(rule); 00224 } 00225 00226 void addIdentity(Treal alpha); /* C += alpha * I */ 00227 00228 static void transpose(Matrix<Treal, Telement> const & A, 00229 Matrix<Treal, Telement> & AT); 00230 00231 void symToNosym(); 00232 void nosymToSym(); 00233 00234 /* Basic linear algebra routines */ 00235 00236 /* Set matrix to zero (k = 0) or identity (k = 1) */ 00237 Matrix<Treal, Telement>& operator=(int const k); 00238 00239 Matrix<Treal, Telement>& operator*=(const Treal alpha); 00240 00241 static void gemm(const bool tA, const bool tB, const Treal alpha, 00242 const Matrix<Treal, Telement>& A, 00243 const Matrix<Treal, Telement>& B, 00244 const Treal beta, 00245 Matrix<Treal, Telement>& C); 00246 static void symm(const char side, const char uplo, const Treal alpha, 00247 const Matrix<Treal, Telement>& A, 00248 const Matrix<Treal, Telement>& B, 00249 const Treal beta, 00250 Matrix<Treal, Telement>& C); 00251 static void syrk(const char uplo, const bool tA, const Treal alpha, 00252 const Matrix<Treal, Telement>& A, 00253 const Treal beta, 00254 Matrix<Treal, Telement>& C); 00255 /* C = alpha * A * A + beta * C where A and C are symmetric */ 00256 static void sysq(const char uplo, const Treal alpha, 00257 const Matrix<Treal, Telement>& A, 00258 const Treal beta, 00259 Matrix<Treal, Telement>& C); 00260 /* C = alpha * A * B + beta * C where A and B are symmetric */ 00261 static void ssmm(const Treal alpha, 00262 const Matrix<Treal, Telement>& A, 00263 const Matrix<Treal, Telement>& B, 00264 const Treal beta, 00265 Matrix<Treal, Telement>& C); 00266 /* C = alpha * A * B + beta * C where A and B are symmetric 00267 * and only the upper triangle of C is computed. 00268 */ 00269 static void ssmm_upper_tr_only(const Treal alpha, 00270 const Matrix<Treal, Telement>& A, 00271 const Matrix<Treal, Telement>& B, 00272 const Treal beta, 00273 Matrix<Treal, Telement>& C); 00274 00275 static void trmm(const char side, const char uplo, const bool tA, 00276 const Treal alpha, 00277 const Matrix<Treal, Telement>& A, 00278 Matrix<Treal, Telement>& B); 00279 00280 /* Frobenius norms */ 00281 00282 /* Returns the Frobenius norm of the matrix. */ 00283 inline Treal frob() const { 00284 return template_blas_sqrt(this->frobSquared()); 00285 } 00286 /* Returns the squared Frobenius norm */ 00287 Treal frobSquared() const; 00288 inline Treal syFrob() const { 00289 return template_blas_sqrt(this->syFrobSquared()); 00290 } 00291 Treal syFrobSquared() const; 00292 00293 inline static Treal frobDiff 00294 (const Matrix<Treal, Telement>& A, 00295 const Matrix<Treal, Telement>& B) { 00296 return template_blas_sqrt(frobSquaredDiff(A, B)); 00297 } 00298 static Treal frobSquaredDiff 00299 (const Matrix<Treal, Telement>& A, 00300 const Matrix<Treal, Telement>& B); 00301 00302 inline static Treal syFrobDiff 00303 (const Matrix<Treal, Telement>& A, 00304 const Matrix<Treal, Telement>& B) { 00305 return template_blas_sqrt(syFrobSquaredDiff(A, B)); 00306 } 00307 static Treal syFrobSquaredDiff 00308 (const Matrix<Treal, Telement>& A, 00309 const Matrix<Treal, Telement>& B); 00310 00311 /* Traces */ 00312 Treal trace() const; 00313 static Treal trace_ab(const Matrix<Treal, Telement>& A, 00314 const Matrix<Treal, Telement>& B); 00315 static Treal trace_aTb(const Matrix<Treal, Telement>& A, 00316 const Matrix<Treal, Telement>& B); 00317 static Treal sy_trace_ab(const Matrix<Treal, Telement>& A, 00318 const Matrix<Treal, Telement>& B); 00319 00320 static void add(const Treal alpha, /* B += alpha * A */ 00321 const Matrix<Treal, Telement>& A, 00322 Matrix<Treal, Telement>& B); 00323 void assign(Treal const alpha, /* *this = alpha * A */ 00324 Matrix<Treal, Telement> const & A); 00325 00326 00327 /********** Help functions for thresholding */ 00328 // int getnnzBlocksLowestLevel() const; 00329 void getFrobSqLowestLevel(std::vector<Treal> & frobsq) const; 00330 void frobThreshLowestLevel 00331 (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix); 00332 00333 void getFrobSqElementLevel(std::vector<Treal> & frobsq) const; 00334 void frobThreshElementLevel 00335 (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix); 00336 00337 00338 #if 0 00339 inline void frobThreshLowestLevel 00340 (Treal const threshold, 00341 Matrix<Treal, Telement> * ErrorMatrix) { 00342 bool a,b; 00343 frobThreshLowestLevel(threshold, ErrorMatrix, a, b); 00344 } 00345 #endif 00346 00349 void assignFrobNormsLowestLevel 00350 ( Matrix<Treal, Matrix<Treal, Telement> > const & A ); 00352 void syAssignFrobNormsLowestLevel 00353 ( Matrix<Treal, Matrix<Treal, Telement> > const & A ); 00354 00358 void assignDiffFrobNormsLowestLevel 00359 ( Matrix<Treal, Matrix<Treal, Telement> > const & A, 00360 Matrix<Treal, Matrix<Treal, Telement> > const & B ); 00364 void syAssignDiffFrobNormsLowestLevel 00365 ( Matrix<Treal, Matrix<Treal, Telement> > const & A, 00366 Matrix<Treal, Matrix<Treal, Telement> > const & B ); 00367 00370 void truncateAccordingToSparsityPattern( Matrix<Treal, Matrix<Treal, Telement> > & A ) const; 00371 00372 00373 /********** End of help functions for thresholding */ 00374 00375 static void gemm_upper_tr_only(const bool tA, const bool tB, 00376 const Treal alpha, 00377 const Matrix<Treal, Telement>& A, 00378 const Matrix<Treal, Telement>& B, 00379 const Treal beta, 00380 Matrix<Treal, Telement>& C); 00381 static void sytr_upper_tr_only(char const side, const Treal alpha, 00382 Matrix<Treal, Telement>& A, 00383 const Matrix<Treal, Telement>& Z); 00384 static void trmm_upper_tr_only(const char side, const char uplo, 00385 const bool tA, const Treal alpha, 00386 const Matrix<Treal, Telement>& A, 00387 Matrix<Treal, Telement>& B); 00388 static void trsytriplemm(char const side, 00389 const Matrix<Treal, Telement>& Z, 00390 Matrix<Treal, Telement>& A); 00391 00392 inline Treal frob_thresh 00393 (Treal const threshold, 00394 Matrix<Treal, Telement> * ErrorMatrix = 0) { 00395 return template_blas_sqrt(frob_squared_thresh(threshold * threshold, 00396 ErrorMatrix)); 00397 } 00403 Treal frob_squared_thresh 00404 (Treal const threshold, 00405 Matrix<Treal, Telement> * ErrorMatrix = 0); 00411 static void syInch(const Matrix<Treal, Telement>& A, 00412 Matrix<Treal, Telement>& Z, 00413 const Treal threshold = 0, 00414 const side looking = left, 00415 const inchversion version = unstable); 00416 00417 void gersgorin(Treal& lmin, Treal& lmax) const; /* Computes bounds for*/ 00418 /* real part of eigenvalues. The matrix must be quadratic (of course) */ 00419 void sy_gersgorin(Treal& lmin, Treal& lmax) const { 00420 Matrix<Treal, Telement> tmp = (*this); 00421 tmp.symToNosym(); 00422 tmp.gersgorin(lmin, lmax); 00423 return; 00424 } 00425 00426 void add_abs_col_sums(Treal* abscolsums) const; /* Adds the absolute */ 00427 /* column sums to the abscolsums array. */ 00428 /* abscolsums(i) += sum(abs(C(:,i))) for all i (C == *this) */ 00429 /* Used by e.g. gersgorin eigenvalue inclusion */ 00430 void get_diagonal(Treal* diag) const; /*Copy diagonal to the diag array*/ 00431 00432 size_t memory_usage() const; /* Returns memory usage in bytes */ 00433 00434 /* Note: we use size_t instead of int for nnz and nvalues to avoid 00435 integer overflow. */ 00436 size_t nnz() const; 00437 size_t sy_nnz() const; 00441 inline size_t nvalues() const { 00442 return nnz(); 00443 } 00446 size_t sy_nvalues() const; 00453 template<typename Top> 00454 Treal syAccumulateWith(Top & op) { 00455 Treal res = 0; 00456 if (!this->is_zero()) { 00457 for (int col = 0; col < this->ncols(); col++) { 00458 for (int row = 0; row < col; row++) 00459 res += 2 * (*this)(row, col).geAccumulateWith(op); 00460 res += (*this)(col, col).syAccumulateWith(op); 00461 } 00462 } 00463 return res; 00464 } 00466 template<typename Top> 00467 Treal geAccumulateWith(Top & op) { 00468 Treal res = 0; 00469 if (!this->is_zero()) { 00470 for (int col = 0; col < this->ncols(); col++) 00471 for (int row = 0; row < this->nrows(); row++) 00472 res += (*this)(row, col).geAccumulateWith(op); 00473 } 00474 return res; 00475 } 00476 00477 static inline unsigned int level() {return Telement::level() + 1;} 00478 00479 Treal maxAbsValue() const { 00480 if (this->is_zero()) 00481 return 0; 00482 else { 00483 Treal maxAbsGlobal = 0; 00484 Treal maxAbsLocal = 0; 00485 for (int ind = 0; ind < this->nElements(); ++ind) { 00486 maxAbsLocal = this->elements[ind].maxAbsValue(); 00487 maxAbsGlobal = maxAbsGlobal > maxAbsLocal ? 00488 maxAbsGlobal : maxAbsLocal; 00489 } /* end for */ 00490 return maxAbsGlobal; 00491 } 00492 } 00493 00494 protected: 00495 private: 00496 }; /* end class */ 00497 00498 00499 #if 1 00500 /* Full matrix assigns etc */ 00501 template<class Treal, class Telement> 00502 void Matrix<Treal, Telement>:: 00503 assignFromFull(std::vector<Treal> const & fullMat) { 00504 int nTotalRows = this->rows.getNTotalScalars(); 00505 int nTotalCols = this->cols.getNTotalScalars(); 00506 assert((int)fullMat.size() == nTotalRows * nTotalCols); 00507 if (this->is_zero()) 00508 allocate(); 00509 for (int col = 0; col < this->ncols(); col++) 00510 for (int row = 0; row < this->nrows(); row++) 00511 (*this)(row, col).assignFromFull(fullMat); 00512 } 00513 00514 template<class Treal, class Telement> 00515 void Matrix<Treal, Telement>:: 00516 fullMatrix(std::vector<Treal> & fullMat) const { 00517 int nTotalRows = this->rows.getNTotalScalars(); 00518 int nTotalCols = this->cols.getNTotalScalars(); 00519 fullMat.resize(nTotalRows * nTotalCols); 00520 if (this->is_zero()) { 00521 int rowOffset = this->rows.getOffset(); 00522 int colOffset = this->cols.getOffset(); 00523 for (int col = 0; col < this->nScalarsCols(); col++) 00524 for (int row = 0; row < this->nScalarsRows(); row++) 00525 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0; 00526 } 00527 else { 00528 for (int col = 0; col < this->ncols(); col++) 00529 for (int row = 0; row < this->nrows(); row++) 00530 (*this)(row, col).fullMatrix(fullMat); 00531 } 00532 } 00533 00534 template<class Treal, class Telement> 00535 void Matrix<Treal, Telement>:: 00536 syFullMatrix(std::vector<Treal> & fullMat) const { 00537 int nTotalRows = this->rows.getNTotalScalars(); 00538 int nTotalCols = this->cols.getNTotalScalars(); 00539 fullMat.resize(nTotalRows * nTotalCols); 00540 if (this->is_zero()) { 00541 int rowOffset = this->rows.getOffset(); 00542 int colOffset = this->cols.getOffset(); 00543 for (int col = 0; col < this->nScalarsCols(); col++) 00544 for (int row = 0; row <= col; row++) { 00545 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0; 00546 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0; 00547 } 00548 } 00549 else { 00550 for (int col = 0; col < this->ncols(); col++) { 00551 for (int row = 0; row < col; row++) 00552 (*this)(row, col).syUpTriFullMatrix(fullMat); 00553 (*this)(col, col).syFullMatrix(fullMat); 00554 } 00555 } 00556 } 00557 00558 template<class Treal, class Telement> 00559 void Matrix<Treal, Telement>:: 00560 syUpTriFullMatrix(std::vector<Treal> & fullMat) const { 00561 int nTotalRows = this->rows.getNTotalScalars(); 00562 int nTotalCols = this->cols.getNTotalScalars(); 00563 fullMat.resize(nTotalRows * nTotalCols); 00564 if (this->is_zero()) { 00565 int rowOffset = this->rows.getOffset(); 00566 int colOffset = this->cols.getOffset(); 00567 for (int col = 0; col < this->nScalarsCols(); col++) 00568 for (int row = 0; row <= this->nScalarsRows(); row++) { 00569 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0; 00570 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0; 00571 } 00572 } 00573 else { 00574 for (int col = 0; col < this->ncols(); col++) 00575 for (int row = 0; row < this->nrows(); row++) 00576 (*this)(row, col).syUpTriFullMatrix(fullMat); 00577 } 00578 } 00579 00580 #endif 00581 00582 00583 template<class Treal, class Telement> 00584 void Matrix<Treal, Telement>:: 00585 assignFromSparse(std::vector<int> const & rowind, 00586 std::vector<int> const & colind, 00587 std::vector<Treal> const & values) { 00588 assert(rowind.size() == colind.size() && 00589 rowind.size() == values.size()); 00590 std::vector<int> indexes(values.size()); 00591 for (unsigned int ind = 0; ind < values.size(); ++ind) 00592 indexes[ind] = ind; 00593 assignFromSparse(rowind, colind, values, indexes); 00594 } 00595 00596 template<class Treal, class Telement> 00597 void Matrix<Treal, Telement>:: 00598 assignFromSparse(std::vector<int> const & rowind, 00599 std::vector<int> const & colind, 00600 std::vector<Treal> const & values, 00601 std::vector<int> const & indexes) { 00602 if (indexes.empty()) { 00603 this->clear(); 00604 return; 00605 } 00606 if (this->is_zero()) 00607 allocate(); 00608 00609 std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks()); 00610 std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks()); 00611 int currentInd = 0; 00612 00613 00614 std::vector<int>::const_iterator it; 00615 for ( it = indexes.begin(); it < indexes.end(); it++ ) 00616 columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it); 00617 00618 /* Go through all column buckets. */ 00619 for (int col = 0; col < this->ncols(); col++) { 00620 /* Go through current column bucket and distribute to row buckets. */ 00621 while (!columnBuckets[col].empty()) { 00622 currentInd = columnBuckets[col].back(); 00623 columnBuckets[col].pop_back(); 00624 rowBuckets[ this->rows.whichBlock 00625 ( rowind[currentInd] ) ].push_back (currentInd); 00626 } 00627 /* Make calls to lower level for every row bucket. */ 00628 for (int row = 0; row < this->nrows(); row++) { 00629 (*this)(row,col).assignFromSparse(rowind, colind, values, rowBuckets[row]); 00630 rowBuckets[row].clear(); 00631 } /* end row loop */ 00632 } /* end column loop */ 00633 } 00634 00635 template<class Treal, class Telement> 00636 void Matrix<Treal, Telement>:: 00637 addValues(std::vector<int> const & rowind, 00638 std::vector<int> const & colind, 00639 std::vector<Treal> const & values) { 00640 assert(rowind.size() == colind.size() && 00641 rowind.size() == values.size()); 00642 std::vector<int> indexes(values.size()); 00643 for (unsigned int ind = 0; ind < values.size(); ++ind) 00644 indexes[ind] = ind; 00645 addValues(rowind, colind, values, indexes); 00646 } 00647 00648 template<class Treal, class Telement> 00649 void Matrix<Treal, Telement>:: 00650 addValues(std::vector<int> const & rowind, 00651 std::vector<int> const & colind, 00652 std::vector<Treal> const & values, 00653 std::vector<int> const & indexes) { 00654 if (indexes.empty()) 00655 return; 00656 if (this->is_zero()) 00657 allocate(); 00658 00659 std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks()); 00660 std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks()); 00661 int currentInd = 0; 00662 00663 std::vector<int>::const_iterator it; 00664 for ( it = indexes.begin(); it < indexes.end(); it++ ) 00665 columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it); 00666 00667 /* Go through all column buckets. */ 00668 for (int col = 0; col < this->ncols(); col++) { 00669 /* Go through current column bucket and distribute to row buckets. */ 00670 while (!columnBuckets[col].empty()) { 00671 currentInd = columnBuckets[col].back(); 00672 columnBuckets[col].pop_back(); 00673 rowBuckets[ this->rows.whichBlock 00674 ( rowind[currentInd] ) ].push_back (currentInd); 00675 } 00676 /* Make calls to lower level for every row bucket. */ 00677 for (int row = 0; row < this->nrows(); row++) { 00678 (*this)(row,col).addValues(rowind, colind, values, rowBuckets[row]); 00679 rowBuckets[row].clear(); 00680 } /* end row loop */ 00681 } /* end column loop */ 00682 } 00683 00684 template<class Treal, class Telement> 00685 void Matrix<Treal, Telement>:: 00686 syAssignFromSparse(std::vector<int> const & rowind, 00687 std::vector<int> const & colind, 00688 std::vector<Treal> const & values) { 00689 assert(rowind.size() == colind.size() && 00690 rowind.size() == values.size()); 00691 bool upperTriangleOnly = true; 00692 for (unsigned int ind = 0; ind < values.size(); ++ind) { 00693 upperTriangleOnly = 00694 upperTriangleOnly && colind[ind] >= rowind[ind]; 00695 } 00696 if (!upperTriangleOnly) 00697 throw Failure("Matrix<Treal, Telement>::" 00698 "syAddValues(std::vector<int> const &, " 00699 "std::vector<int> const &, " 00700 "std::vector<Treal> const &, int const): " 00701 "Only upper triangle can contain elements when assigning " 00702 "symmetric or triangular matrix "); 00703 assignFromSparse(rowind, colind, values); 00704 } 00705 00706 template<class Treal, class Telement> 00707 void Matrix<Treal, Telement>:: 00708 syAddValues(std::vector<int> const & rowind, 00709 std::vector<int> const & colind, 00710 std::vector<Treal> const & values) { 00711 assert(rowind.size() == colind.size() && 00712 rowind.size() == values.size()); 00713 bool upperTriangleOnly = true; 00714 for (unsigned int ind = 0; ind < values.size(); ++ind) { 00715 upperTriangleOnly = 00716 upperTriangleOnly && colind[ind] >= rowind[ind]; 00717 } 00718 if (!upperTriangleOnly) 00719 throw Failure("Matrix<Treal, Telement>::" 00720 "syAddValues(std::vector<int> const &, " 00721 "std::vector<int> const &, " 00722 "std::vector<Treal> const &, int const): " 00723 "Only upper triangle can contain elements when assigning " 00724 "symmetric or triangular matrix "); 00725 addValues(rowind, colind, values); 00726 } 00727 00728 template<class Treal, class Telement> 00729 void Matrix<Treal, Telement>:: 00730 getValues(std::vector<int> const & rowind, 00731 std::vector<int> const & colind, 00732 std::vector<Treal> & values) const { 00733 assert(rowind.size() == colind.size()); 00734 values.resize(rowind.size(),0); 00735 std::vector<int> indexes(rowind.size()); 00736 for (unsigned int ind = 0; ind < rowind.size(); ++ind) 00737 indexes[ind] = ind; 00738 getValues(rowind, colind, values, indexes); 00739 } 00740 00741 template<class Treal, class Telement> 00742 void Matrix<Treal, Telement>:: 00743 getValues(std::vector<int> const & rowind, 00744 std::vector<int> const & colind, 00745 std::vector<Treal> & values, 00746 std::vector<int> const & indexes) const { 00747 assert(!this->is_empty()); 00748 if (indexes.empty()) 00749 return; 00750 std::vector<int>::const_iterator it; 00751 if (this->is_zero()) { 00752 for ( it = indexes.begin(); it < indexes.end(); it++ ) 00753 values[*it] = 0; 00754 return; 00755 } 00756 00757 std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks()); 00758 std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks()); 00759 int currentInd = 0; 00760 00761 for ( it = indexes.begin(); it < indexes.end(); it++ ) 00762 columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it); 00763 00764 /* Go through all column buckets. */ 00765 for (int col = 0; col < this->ncols(); col++) { 00766 /* Go through current column bucket and distribute to row buckets. */ 00767 while (!columnBuckets[col].empty()) { 00768 currentInd = columnBuckets[col].back(); 00769 columnBuckets[col].pop_back(); 00770 rowBuckets[ this->rows.whichBlock 00771 ( rowind[currentInd] ) ].push_back (currentInd); 00772 } 00773 /* Make calls to lower level for every row bucket. */ 00774 for (int row = 0; row < this->nrows(); row++) { 00775 (*this)(row,col).getValues(rowind, colind, values, rowBuckets[row]); 00776 rowBuckets[row].clear(); 00777 } /* end row loop */ 00778 } /* end column loop */ 00779 } 00780 00781 template<class Treal, class Telement> 00782 void Matrix<Treal, Telement>:: 00783 syGetValues(std::vector<int> const & rowind, 00784 std::vector<int> const & colind, 00785 std::vector<Treal> & values) const { 00786 assert(rowind.size() == colind.size()); 00787 bool upperTriangleOnly = true; 00788 for (int unsigned ind = 0; ind < rowind.size(); ++ind) { 00789 upperTriangleOnly = 00790 upperTriangleOnly && colind[ind] >= rowind[ind]; 00791 } 00792 if (!upperTriangleOnly) 00793 throw Failure("Matrix<Treal, Telement>::" 00794 "syGetValues(std::vector<int> const &, " 00795 "std::vector<int> const &, " 00796 "std::vector<Treal> const &, int const): " 00797 "Only upper triangle when retrieving elements from " 00798 "symmetric or triangular matrix "); 00799 getValues(rowind, colind, values); 00800 } 00801 00802 00803 template<class Treal, class Telement> 00804 void Matrix<Treal, Telement>:: 00805 getAllValues(std::vector<int> & rowind, 00806 std::vector<int> & colind, 00807 std::vector<Treal> & values) const { 00808 assert(rowind.size() == colind.size() && 00809 rowind.size() == values.size()); 00810 if (!this->is_zero()) 00811 for (int ind = 0; ind < this->nElements(); ++ind) 00812 this->elements[ind].getAllValues(rowind, colind, values); 00813 } 00814 00815 template<class Treal, class Telement> 00816 void Matrix<Treal, Telement>:: 00817 syGetAllValues(std::vector<int> & rowind, 00818 std::vector<int> & colind, 00819 std::vector<Treal> & values) const { 00820 assert(rowind.size() == colind.size() && 00821 rowind.size() == values.size()); 00822 if (!this->is_zero()) 00823 for (int col = 0; col < this->ncols(); ++col) { 00824 for (int row = 0; row < col; ++row) 00825 (*this)(row, col).getAllValues(rowind, colind, values); 00826 (*this)(col, col).syGetAllValues(rowind, colind, values); 00827 } 00828 } 00829 00830 00831 template<class Treal, class Telement> 00832 void Matrix<Treal, Telement>::clear() { 00833 if (!this->is_zero()) 00834 for (int i = 0; i < this->nElements(); i++) 00835 this->elements[i].clear(); 00836 delete[] this->elements; 00837 this->elements = 0; 00838 } 00839 00840 template<class Treal, class Telement> 00841 void Matrix<Treal, Telement>:: 00842 writeToFile(std::ofstream & file) const { 00843 int const ZERO = 0; 00844 int const ONE = 1; 00845 if (this->is_zero()) { 00846 char * tmp = (char*)&ZERO; 00847 file.write(tmp,sizeof(int)); 00848 } 00849 else { 00850 char * tmp = (char*)&ONE; 00851 file.write(tmp,sizeof(int)); 00852 for (int i = 0; i < this->nElements(); i++) 00853 this->elements[i].writeToFile(file); 00854 } 00855 } 00856 template<class Treal, class Telement> 00857 void Matrix<Treal, Telement>:: 00858 readFromFile(std::ifstream & file) { 00859 int const ZERO = 0; 00860 int const ONE = 1; 00861 char tmp[sizeof(int)]; 00862 file.read(tmp, (std::ifstream::pos_type)sizeof(int)); 00863 switch ((int)*tmp) { 00864 case ZERO: 00865 this->clear(); 00866 break; 00867 case ONE: 00868 if (this->is_zero()) 00869 allocate(); 00870 for (int i = 0; i < this->nElements(); i++) 00871 this->elements[i].readFromFile(file); 00872 break; 00873 default: 00874 throw Failure("Matrix<Treal, Telement>::" 00875 "readFromFile(std::ifstream & file):" 00876 "File corruption int value not 0 or 1"); 00877 } 00878 } 00879 00880 template<class Treal, class Telement> 00881 void Matrix<Treal, Telement>::random() { 00882 if (this->is_zero()) 00883 allocate(); 00884 for (int ind = 0; ind < this->nElements(); ind++) 00885 this->elements[ind].random(); 00886 } 00887 00888 template<class Treal, class Telement> 00889 void Matrix<Treal, Telement>::syRandom() { 00890 if (this->is_zero()) 00891 allocate(); 00892 /* Above diagonal */ 00893 for (int col = 1; col < this->ncols(); col++) 00894 for (int row = 0; row < col; row++) 00895 (*this)(row, col).random(); 00896 /* Diagonal */ 00897 for (int rc = 0; rc < this->nrows(); rc++) 00898 (*this)(rc,rc).syRandom(); 00899 } 00900 00901 template<class Treal, class Telement> 00902 void Matrix<Treal, Telement>:: 00903 randomZeroStructure(Treal probabilityBeingZero) { 00904 if (!this->highestLevel() && 00905 probabilityBeingZero > rand() / (Treal)RAND_MAX) 00906 this->clear(); 00907 else { 00908 if (this->is_zero()) 00909 allocate(); 00910 for (int ind = 0; ind < this->nElements(); ind++) 00911 this->elements[ind].randomZeroStructure(probabilityBeingZero); 00912 } 00913 } 00914 00915 template<class Treal, class Telement> 00916 void Matrix<Treal, Telement>:: 00917 syRandomZeroStructure(Treal probabilityBeingZero) { 00918 if (!this->highestLevel() && 00919 probabilityBeingZero > rand() / (Treal)RAND_MAX) 00920 this->clear(); 00921 else { 00922 if (this->is_zero()) 00923 allocate(); 00924 /* Above diagonal */ 00925 for (int col = 1; col < this->ncols(); col++) 00926 for (int row = 0; row < col; row++) 00927 (*this)(row, col).randomZeroStructure(probabilityBeingZero); 00928 /* Diagonal */ 00929 for (int rc = 0; rc < this->nrows(); rc++) 00930 (*this)(rc,rc).syRandomZeroStructure(probabilityBeingZero); 00931 } 00932 } 00933 00934 template<class Treal, class Telement> 00935 template<typename TRule> 00936 void Matrix<Treal, Telement>:: 00937 setElementsByRule(TRule & rule) { 00938 if (this->is_zero()) 00939 allocate(); 00940 for (int ind = 0; ind < this->nElements(); ind++) 00941 this->elements[ind].setElementsByRule(rule); 00942 } 00943 template<class Treal, class Telement> 00944 template<typename TRule> 00945 void Matrix<Treal, Telement>:: 00946 sySetElementsByRule(TRule & rule) { 00947 if (this->is_zero()) 00948 allocate(); 00949 /* Above diagonal */ 00950 for (int col = 1; col < this->ncols(); col++) 00951 for (int row = 0; row < col; row++) 00952 (*this)(row, col).setElementsByRule(rule); 00953 /* Diagonal */ 00954 for (int rc = 0; rc < this->nrows(); rc++) 00955 (*this)(rc,rc).sySetElementsByRule(rule); 00956 } 00957 00958 00959 template<class Treal, class Telement> 00960 void Matrix<Treal, Telement>:: 00961 addIdentity(Treal alpha) { 00962 if (this->is_empty()) 00963 throw Failure("Matrix<Treal, Telement>::addIdentity(Treal): " 00964 "Cannot add identity to empty matrix."); 00965 if (this->ncols() != this->nrows()) 00966 throw Failure("Matrix<Treal, Telement>::addIdentity(Treal): " 00967 "Matrix must be square to add identity"); 00968 if (this->is_zero()) 00969 allocate(); 00970 for (int ind = 0; ind < this->ncols(); ind++) 00971 (*this)(ind,ind).addIdentity(alpha); 00972 } 00973 00974 template<class Treal, class Telement> 00975 void Matrix<Treal, Telement>:: 00976 transpose(Matrix<Treal, Telement> const & A, 00977 Matrix<Treal, Telement> & AT) { 00978 if (A.is_zero()) { /* Condition also matches empty matrices. */ 00979 AT.rows = A.cols; 00980 AT.cols = A.rows; 00981 delete[] AT.elements; 00982 AT.elements = 0; 00983 return; 00984 } 00985 if (AT.is_zero() || (AT.nElements() != A.nElements())) { 00986 delete[] AT.elements; 00987 AT.elements = new Telement[A.nElements()]; 00988 } 00989 AT.cols = A.rows; 00990 AT.rows = A.cols; 00991 for (int row = 0; row < AT.nrows(); row++) 00992 for (int col = 0; col < AT.ncols(); col++) 00993 Telement::transpose(A(col,row), AT(row,col)); 00994 } 00995 00996 template<class Treal, class Telement> 00997 void Matrix<Treal, Telement>:: 00998 symToNosym() { 00999 try { 01000 if (this->nrows() == this->ncols()) { 01001 if (!this->is_zero()) { 01002 /* Fix the diagonal: */ 01003 for (int rc = 0; rc < this->ncols(); rc++) 01004 (*this)(rc, rc).symToNosym(); 01005 /* Fix the lower triangle */ 01006 for (int row = 1; row < this->nrows(); row++) 01007 for (int col = 0; col < row; col++) 01008 Telement::transpose((*this)(col, row), (*this)(row, col)); 01009 } 01010 } 01011 else 01012 throw Failure("Matrix<Treal, Telement>::symToNosym(): " 01013 "Only quadratic matrices can be symmetric"); 01014 } 01015 catch(Failure& f) { 01016 std::cout<<"Failure caught:Matrix<Treal, Telement>::symToNosym()" 01017 <<std::endl; 01018 throw f; 01019 } 01020 } 01021 01022 template<class Treal, class Telement> 01023 void Matrix<Treal, Telement>:: 01024 nosymToSym() { 01025 if (this->nrows() == this->ncols()) { 01026 if (!this->is_zero()) { 01027 /* Fix the diagonal: */ 01028 for (int rc = 0; rc < this->ncols(); rc++) 01029 (*this)(rc, rc).nosymToSym(); 01030 /* Fix the lower triangle */ 01031 for (int col = 0; col < this->ncols() - 1; col++) 01032 for (int row = col + 1; row < this->nrows(); row++) 01033 (*this)(row, col).clear(); 01034 } 01035 } 01036 else 01037 throw Failure("Matrix<Treal, Telement>::nosymToSym(): " 01038 "Only quadratic matrices can be symmetric"); 01039 } 01040 01041 /* Basic linear algebra routines */ 01042 01043 template<class Treal, class Telement> 01044 Matrix<Treal, Telement>& 01045 Matrix<Treal, Telement>::operator=(int const k) { 01046 switch (k) { 01047 case 0: 01048 this->clear(); 01049 break; 01050 case 1: 01051 if (this->ncols() != this->nrows()) 01052 throw Failure 01053 ("Matrix::operator=(int k = 1): " 01054 "Matrix must be quadratic to become identity matrix."); 01055 this->clear(); 01056 this->allocate(); 01057 for (int rc = 0; rc < this->ncols(); rc++) /*Set diagonal to identity*/ 01058 (*this)(rc,rc) = 1; 01059 break; 01060 default: 01061 throw Failure("Matrix::operator=(int k) only " 01062 "implemented for k = 0, k = 1"); 01063 } 01064 return *this; 01065 } 01066 01067 01068 template<class Treal, class Telement> 01069 Matrix<Treal, Telement>& Matrix<Treal, Telement>:: 01070 operator*=(const Treal alpha) { 01071 if (!this->is_zero() && alpha != 1) { 01072 for (int ind = 0; ind < this->nElements(); ind++) 01073 this->elements[ind] *= alpha; 01074 } 01075 return *this; 01076 } 01077 01078 /* C = beta * C + alpha * A * B */ 01079 template<class Treal, class Telement> 01080 void Matrix<Treal, Telement>:: 01081 gemm(const bool tA, const bool tB, const Treal alpha, 01082 const Matrix<Treal, Telement>& A, 01083 const Matrix<Treal, Telement>& B, 01084 const Treal beta, 01085 Matrix<Treal, Telement>& C) { 01086 assert(!A.is_empty()); 01087 assert(!B.is_empty()); 01088 if (C.is_empty()) { 01089 assert(beta == 0); 01090 if (tA) 01091 C.resetRows(A.cols); 01092 else 01093 C.resetRows(A.rows); 01094 if (tB) 01095 C.resetCols(B.rows); 01096 else 01097 C.resetCols(B.cols); 01098 } 01099 01100 if ((A.is_zero() || B.is_zero() || alpha == 0) && 01101 (C.is_zero() || beta == 0)) 01102 C = 0; 01103 else { 01104 Treal beta_tmp = beta; 01105 if (C.is_zero()) { 01106 C.allocate(); 01107 beta_tmp = 0; 01108 } 01109 if (!A.is_zero() && !B.is_zero() && alpha != 0) { 01110 MAT_OMP_INIT; 01111 if (!tA && !tB) 01112 if (A.ncols() == B.nrows() && 01113 A.nrows() == C.nrows() && 01114 B.ncols() == C.ncols()) { 01115 #ifdef _OPENMP 01116 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 01117 #endif 01118 for (int col = 0; col < C.ncols(); col++) { 01119 MAT_OMP_START; 01120 for (int row = 0; row < C.nrows(); row++) { 01121 Telement::gemm(tA, tB, alpha, 01122 A(row, 0), B(0, col), 01123 beta_tmp, 01124 C(row, col)); 01125 for(int cola = 1; cola < A.ncols(); cola++) 01126 Telement::gemm(tA, tB, alpha, 01127 A(row, cola), B(cola, col), 01128 1.0, 01129 C(row, col)); 01130 } 01131 MAT_OMP_END; 01132 } /* end omp for */ 01133 } 01134 else 01135 throw Failure("Matrix<class Treal, class Telement>::" 01136 "gemm(bool, bool, Treal, " 01137 "Matrix<Treal, Telement>, " 01138 "Matrix<Treal, Telement>, Treal, " 01139 "Matrix<Treal, Telement>): " 01140 "Incorrect matrixdimensions for multiplication"); 01141 else if (tA && !tB) 01142 if (A.nrows() == B.nrows() && 01143 A.ncols() == C.nrows() && 01144 B.ncols() == C.ncols()) { 01145 #ifdef _OPENMP 01146 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 01147 #endif 01148 for (int col = 0; col < C.ncols(); col++) { 01149 MAT_OMP_START; 01150 for (int row = 0; row < C.nrows(); row++) { 01151 Telement::gemm(tA, tB, alpha, 01152 A(0,row), B(0,col), 01153 beta_tmp, 01154 C(row,col)); 01155 for(int cola = 1; cola < A.nrows(); cola++) 01156 Telement::gemm(tA, tB, alpha, 01157 A(cola, row), B(cola, col), 01158 1.0, 01159 C(row,col)); 01160 } 01161 MAT_OMP_END; 01162 } /* end omp for */ 01163 } 01164 else 01165 throw Failure("Matrix<class Treal, class Telement>::" 01166 "gemm(bool, bool, Treal, " 01167 "Matrix<Treal, Telement>, " 01168 "Matrix<Treal, Telement>, Treal, " 01169 "Matrix<Treal, Telement>): " 01170 "Incorrect matrixdimensions for multiplication"); 01171 else if (!tA && tB) 01172 if (A.ncols() == B.ncols() && 01173 A.nrows() == C.nrows() && 01174 B.nrows() == C.ncols()) { 01175 #ifdef _OPENMP 01176 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 01177 #endif 01178 for (int col = 0; col < C.ncols(); col++) { 01179 MAT_OMP_START; 01180 for (int row = 0; row < C.nrows(); row++) { 01181 Telement::gemm(tA, tB, alpha, 01182 A(row, 0), B(col, 0), 01183 beta_tmp, 01184 C(row,col)); 01185 for(int cola = 1; cola < A.ncols(); cola++) 01186 Telement::gemm(tA, tB, alpha, 01187 A(row, cola), B(col, cola), 01188 1.0, 01189 C(row,col)); 01190 } 01191 MAT_OMP_END; 01192 } /* end omp for */ 01193 } 01194 else 01195 throw Failure("Matrix<class Treal, class Telement>::" 01196 "gemm(bool, bool, Treal, " 01197 "Matrix<Treal, Telement>, " 01198 "Matrix<Treal, Telement>, Treal, " 01199 "Matrix<Treal, Telement>): " 01200 "Incorrect matrixdimensions for multiplication"); 01201 else if (tA && tB) 01202 if (A.nrows() == B.ncols() && 01203 A.ncols() == C.nrows() && 01204 B.nrows() == C.ncols()) { 01205 #ifdef _OPENMP 01206 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 01207 #endif 01208 for (int col = 0; col < C.ncols(); col++) { 01209 MAT_OMP_START; 01210 for (int row = 0; row < C.nrows(); row++) { 01211 Telement::gemm(tA, tB, alpha, 01212 A(0, row), B(col, 0), 01213 beta_tmp, 01214 C(row,col)); 01215 for(int cola = 1; cola < A.nrows(); cola++) 01216 Telement::gemm(tA, tB, alpha, 01217 A(cola, row), B(col, cola), 01218 1.0, 01219 C(row,col)); 01220 } 01221 MAT_OMP_END; 01222 } /* end omp for */ 01223 } 01224 else 01225 throw Failure("Matrix<class Treal, class Telement>::" 01226 "gemm(bool, bool, Treal, " 01227 "Matrix<Treal, Telement>, " 01228 "Matrix<Treal, Telement>, " 01229 "Treal, Matrix<Treal, Telement>): " 01230 "Incorrect matrixdimensions for multiplication"); 01231 else throw Failure("Matrix<class Treal, class Telement>::" 01232 "gemm(bool, bool, Treal, " 01233 "Matrix<Treal, Telement>, " 01234 "Matrix<Treal, Telement>, Treal, " 01235 "Matrix<Treal, Telement>):" 01236 "Very strange error!!"); 01237 MAT_OMP_FINALIZE; 01238 } 01239 else 01240 C *= beta; 01241 } 01242 } 01243 01244 template<class Treal, class Telement> 01245 void Matrix<Treal, Telement>:: 01246 symm(const char side, const char uplo, const Treal alpha, 01247 const Matrix<Treal, Telement>& A, 01248 const Matrix<Treal, Telement>& B, 01249 const Treal beta, 01250 Matrix<Treal, Telement>& C) { 01251 assert(!A.is_empty()); 01252 assert(!B.is_empty()); 01253 assert(A.nrows() == A.ncols()); 01254 int dimA = A.nrows(); 01255 if (C.is_empty()) { 01256 assert(beta == 0); 01257 if (side =='L') { 01258 C.resetRows(A.rows); 01259 C.resetCols(B.cols); 01260 } 01261 else { 01262 assert(side == 'R'); 01263 C.resetRows(B.rows); 01264 C.resetCols(A.cols); 01265 } 01266 } 01267 if ((A.is_zero() || B.is_zero() || alpha == 0) && 01268 (C.is_zero() || beta == 0)) 01269 C = 0; 01270 else { 01271 if (uplo == 'U') { 01272 Treal beta_tmp = beta; 01273 if (C.is_zero()) { 01274 C.allocate(); 01275 beta_tmp = 0; 01276 } 01277 if (!A.is_zero() && !B.is_zero() && alpha != 0) { 01278 MAT_OMP_INIT; 01279 if (side =='L') 01280 if (dimA == B.nrows() && 01281 dimA == C.nrows() && 01282 B.ncols() == C.ncols()) { 01283 #ifdef _OPENMP 01284 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 01285 #endif 01286 for (int col = 0; col < C.ncols(); col++) { 01287 MAT_OMP_START; 01288 for (int row = 0; row < C.nrows(); row++) { 01289 /* Diagonal element in matrix A */ 01290 Telement::symm(side, uplo, alpha, A(row, row), 01291 B(row, col), beta_tmp, C(row, col)); 01292 /* Elements below diagonal in A */ 01293 for(int ind = 0; ind < row; ind++) 01294 Telement::gemm(true, false, alpha, 01295 A(ind, row), B(ind, col), 01296 1.0, C(row,col)); 01297 /* Elements above diagonal in A */ 01298 for(int ind = row + 1; ind < dimA; ind++) 01299 Telement::gemm(false, false, alpha, 01300 A(row, ind), B(ind, col), 01301 1.0, C(row,col)); 01302 } 01303 MAT_OMP_END; 01304 } /* end omp for */ 01305 } 01306 else 01307 throw Failure("Matrix<class Treal, class Telement>" 01308 "::symm(bool, bool, Treal, Matrix<Treal, " 01309 "Telement>, Matrix<Treal, Telement>, " 01310 "Treal, Matrix<Treal, Telement>): " 01311 "Incorrect matrixdimensions for multiplication"); 01312 else { /* side == 'R' */ 01313 assert(side == 'R'); 01314 if (B.ncols() == dimA && 01315 B.nrows() == C.nrows() && 01316 dimA == C.ncols()) { 01317 #ifdef _OPENMP 01318 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 01319 #endif 01320 for (int col = 0; col < C.ncols(); col++) { 01321 MAT_OMP_START; 01322 for (int row = 0; row < C.nrows(); row++) { 01323 /* Diagonal element in matrix A */ 01324 Telement::symm(side, uplo, alpha, A(col, col), 01325 B(row, col), beta_tmp, C(row, col)); 01326 /* Elements below diagonal in A */ 01327 for(int ind = col + 1; ind < dimA; ind++) 01328 Telement::gemm(false, true, alpha, 01329 B(row, ind), A(col, ind), 01330 1.0, C(row,col)); 01331 /* Elements above diagonal in A */ 01332 for(int ind = 0; ind < col; ind++) 01333 Telement::gemm(false, false, alpha, 01334 B(row, ind), A(ind, col), 01335 1.0, C(row,col)); 01336 } 01337 MAT_OMP_END; 01338 } /* end omp for */ 01339 } 01340 else 01341 throw Failure("Matrix<class Treal, class Telement>" 01342 "::symm(bool, bool, Treal, Matrix<Treal, " 01343 "Telement>, Matrix<Treal, Telement>, " 01344 "Treal, Matrix<Treal, Telement>): " 01345 "Incorrect matrixdimensions for multiplication"); 01346 } 01347 MAT_OMP_FINALIZE; 01348 } 01349 else 01350 C *= beta; 01351 } 01352 else 01353 throw Failure("Matrix<class Treal, class Telement>::" 01354 "symm only implemented for symmetric matrices in " 01355 "upper triangular storage"); 01356 } 01357 } 01358 01359 template<class Treal, class Telement> 01360 void Matrix<Treal, Telement>:: 01361 syrk(const char uplo, const bool tA, const Treal alpha, 01362 const Matrix<Treal, Telement>& A, 01363 const Treal beta, 01364 Matrix<Treal, Telement>& C) { 01365 assert(!A.is_empty()); 01366 if (C.is_empty()) { 01367 assert(beta == 0); 01368 if (tA) { 01369 C.resetRows(A.cols); 01370 C.resetCols(A.cols); 01371 } 01372 else { 01373 C.resetRows(A.rows); 01374 C.resetCols(A.rows); 01375 } 01376 } 01377 01378 if (C.nrows() == C.ncols() && 01379 ((!tA && A.nrows() == C.nrows()) || (tA && A.ncols() == C.nrows()))) 01380 if (alpha != 0 && !A.is_zero()) { 01381 Treal beta_tmp = beta; 01382 if (C.is_zero()) { 01383 C.allocate(); 01384 beta_tmp = 0; 01385 } 01386 MAT_OMP_INIT; 01387 if (!tA && uplo == 'U') { /* C = alpha * A * A' + beta * C */ 01388 #ifdef _OPENMP 01389 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) 01390 #endif 01391 { 01392 #ifdef _OPENMP 01393 #pragma omp for schedule(dynamic) nowait 01394 #endif 01395 for (int rc = 0; rc < C.ncols(); rc++) { 01396 MAT_OMP_START; 01397 Telement::syrk(uplo, tA, alpha, A(rc, 0), beta_tmp, C(rc, rc)); 01398 for (int cola = 1; cola < A.ncols(); cola++) 01399 Telement::syrk(uplo, tA, alpha, A(rc, cola), 1.0, C(rc, rc)); 01400 MAT_OMP_END; 01401 } 01402 #ifdef _OPENMP 01403 #pragma omp for schedule(dynamic) nowait 01404 #endif 01405 for (int row = 0; row < C.nrows() - 1; row++) { 01406 MAT_OMP_START; 01407 for (int col = row + 1; col < C.ncols(); col++) { 01408 Telement::gemm(tA, !tA, alpha, 01409 A(row, 0), A(col,0), beta_tmp, C(row,col)); 01410 for (int ind = 1; ind < A.ncols(); ind++) 01411 Telement::gemm(tA, !tA, alpha, 01412 A(row, ind), A(col,ind), 1.0, C(row,col)); 01413 } 01414 MAT_OMP_END; 01415 } 01416 } /* end omp parallel */ 01417 } /* end if (!tA) */ 01418 else if (tA && uplo == 'U') { /* C = alpha * A' * A + beta * C */ 01419 #ifdef _OPENMP 01420 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) 01421 #endif 01422 { 01423 #ifdef _OPENMP 01424 #pragma omp for schedule(dynamic) nowait 01425 #endif 01426 for (int rc = 0; rc < C.ncols(); rc++) { 01427 MAT_OMP_START; 01428 Telement::syrk(uplo, tA, alpha, A(0, rc), beta_tmp, C(rc, rc)); 01429 for (int rowa = 1; rowa < A.nrows(); rowa++) 01430 Telement::syrk(uplo, tA, alpha, A(rowa, rc), 1.0, C(rc, rc)); 01431 MAT_OMP_END; 01432 } 01433 #ifdef _OPENMP 01434 #pragma omp for schedule(dynamic) nowait 01435 #endif 01436 for (int row = 0; row < C.nrows() - 1; row++) { 01437 MAT_OMP_START; 01438 for (int col = row + 1; col < C.ncols(); col++) { 01439 Telement::gemm(tA, !tA, alpha, 01440 A(0, row), A(0, col), beta_tmp, C(row,col)); 01441 for (int ind = 1; ind < A.nrows(); ind++) 01442 Telement::gemm(tA, !tA, alpha, 01443 A(ind, row), A(ind, col), 1.0, C(row,col)); 01444 } 01445 MAT_OMP_END; 01446 } 01447 } /* end omp parallel */ 01448 } /* end if (tA) */ 01449 else 01450 throw Failure("Matrix<class Treal, class Telement>::" 01451 "syrk not implemented for lower triangular storage"); 01452 MAT_OMP_FINALIZE; 01453 } 01454 else { 01455 C *= beta; 01456 } 01457 else 01458 throw Failure("Matrix<class Treal, class Telement>::syrk: " 01459 "Incorrect matrix dimensions for symmetric rank-k update"); 01460 } 01461 01462 template<class Treal, class Telement> 01463 void Matrix<Treal, Telement>:: 01464 sysq(const char uplo, const Treal alpha, 01465 const Matrix<Treal, Telement>& A, 01466 const Treal beta, 01467 Matrix<Treal, Telement>& C) { 01468 assert(!A.is_empty()); 01469 if (C.is_empty()) { 01470 assert(beta == 0); 01471 C.resetRows(A.rows); 01472 C.resetCols(A.cols); 01473 } 01474 if (C.nrows() == C.ncols() && 01475 A.nrows() == C.nrows() && A.nrows() == A.ncols()) 01476 if (alpha != 0 && !A.is_zero()) { 01477 if (uplo == 'U') { 01478 Treal beta_tmp = beta; 01479 if (C.is_zero()) { 01480 C.allocate(); 01481 beta_tmp = 0; 01482 } 01483 MAT_OMP_INIT; 01484 #ifdef _OPENMP 01485 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) 01486 #endif 01487 { 01488 #ifdef _OPENMP 01489 #pragma omp for schedule(dynamic) nowait 01490 #endif 01491 for (int rc = 0; rc < C.ncols(); rc++) { 01492 MAT_OMP_START; 01493 Telement::sysq(uplo, alpha, A(rc, rc), beta_tmp, C(rc, rc)); 01494 for (int cola = 0; cola < rc; cola++) 01495 Telement::syrk(uplo, true, alpha, A(cola, rc), 1.0, C(rc, rc)); 01496 for (int cola = rc + 1; cola < A.ncols(); cola++) 01497 Telement::syrk(uplo, false, alpha, A(rc, cola), 1.0, C(rc, rc)); 01498 MAT_OMP_END; 01499 } 01500 /* Maste anvanda symm? */ 01501 #ifdef _OPENMP 01502 #pragma omp for schedule(dynamic) nowait 01503 #endif 01504 for (int row = 0; row < C.nrows() - 1; row++) { 01505 MAT_OMP_START; 01506 for (int col = row + 1; col < C.ncols(); col++) { 01507 /* First the two operations involving diagonal elements */ 01508 Telement::symm('L', 'U', alpha, A(row, row), A(row, col), 01509 beta_tmp, C(row, col)); 01510 Telement::symm('R', 'U', alpha, A(col, col), A(row, col), 01511 1.0, C(row, col)); 01512 /* First element in product is below the diagonal */ 01513 for (int ind = 0; ind < row; ind++) 01514 Telement::gemm(true, false, alpha, 01515 A(ind, row), A(ind, col), 1.0, C(row, col)); 01516 /* None of the elements are below the diagonal */ 01517 for (int ind = row + 1; ind < col; ind++) 01518 Telement::gemm(false, false, alpha, 01519 A(row, ind), A(ind, col), 1.0, C(row, col)); 01520 /* Second element is below the diagonal */ 01521 for (int ind = col + 1; ind < A.ncols(); ind++) 01522 Telement::gemm(false, true, alpha, 01523 A(row, ind), A(col, ind), 1.0, C(row, col)); 01524 } 01525 MAT_OMP_END; 01526 } /* end omp for */ 01527 } /* end omp parallel */ 01528 MAT_OMP_FINALIZE; 01529 } 01530 else 01531 throw Failure("Matrix<class Treal, class Telement>::" 01532 "sysq only implemented for symmetric matrices in " 01533 "upper triangular storage");; 01534 } 01535 else { 01536 C *= beta; 01537 } 01538 else 01539 throw Failure("Matrix<class Treal, class Telement>::sysq: " 01540 "Incorrect matrix dimensions for symmetric square " 01541 "operation"); 01542 } 01543 01544 /* C = alpha * A * B + beta * C where A and B are symmetric */ 01545 template<class Treal, class Telement> 01546 void Matrix<Treal, Telement>:: 01547 ssmm(const Treal alpha, 01548 const Matrix<Treal, Telement>& A, 01549 const Matrix<Treal, Telement>& B, 01550 const Treal beta, 01551 Matrix<Treal, Telement>& C) { 01552 assert(!A.is_empty()); 01553 assert(!B.is_empty()); 01554 if (C.is_empty()) { 01555 assert(beta == 0); 01556 C.resetRows(A.rows); 01557 C.resetCols(B.cols); 01558 } 01559 if (A.ncols() != B.nrows() || 01560 A.nrows() != C.nrows() || 01561 B.ncols() != C.ncols() || 01562 A.nrows() != A.ncols() || 01563 B.nrows() != B.ncols()) { 01564 throw Failure("Matrix<class Treal, class Telement>::" 01565 "ssmm(Treal, " 01566 "Matrix<Treal, Telement>, " 01567 "Matrix<Treal, Telement>, Treal, " 01568 "Matrix<Treal, Telement>): " 01569 "Incorrect matrixdimensions for ssmm multiplication"); 01570 } 01571 if ((A.is_zero() || B.is_zero() || alpha == 0) && 01572 (C.is_zero() || beta == 0)) { 01573 C = 0; 01574 return; 01575 } 01576 Treal beta_tmp = beta; 01577 if (C.is_zero()) { 01578 C.allocate(); 01579 beta_tmp = 0; 01580 } 01581 if (A.is_zero() || B.is_zero() || alpha == 0) { 01582 C *= beta; 01583 return; 01584 } 01585 MAT_OMP_INIT; 01586 #ifdef _OPENMP 01587 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 01588 #endif 01589 for (int col = 0; col < C.ncols(); col++) { 01590 MAT_OMP_START; 01591 /* Diagonal */ 01592 /* C(col, col) = alpha * A(col, col) * B(col, col) + beta * C(col, col)*/ 01593 Telement::ssmm(alpha, A(col,col), B(col, col), 01594 beta_tmp, C(col,col)); 01595 for (int ind = 0; ind < col; ind++) 01596 /* C(col, col) += alpha * A(ind, col)' * B(ind, col) */ 01597 Telement::gemm(true, false, 01598 alpha, A(ind,col), B(ind, col), 01599 1.0, C(col,col)); 01600 for (int ind = col + 1; ind < A.ncols(); ind++) 01601 /* C(col, col) += alpha * A(col, ind) * B(col, ind)' */ 01602 Telement::gemm(false, true, 01603 alpha, A(col, ind), B(col, ind), 01604 1.0, C(col,col)); 01605 /* Upper half */ 01606 for (int row = 0; row < col; row++) { 01607 /* C(row, col) = alpha * A(row, row) * B(row, col) + 01608 * beta * C(row, col) 01609 */ 01610 Telement::symm('L', 'U', 01611 alpha, A(row, row), B(row, col), 01612 beta_tmp, C(row, col)); 01613 /* C(row, col) += alpha * A(row, col) * B(col, col) */ 01614 Telement::symm('R', 'U', 01615 alpha, B(col, col), A(row, col), 01616 1.0, C(row, col)); 01617 for (int ind = 0; ind < row; ind++) 01618 /* C(row, col) += alpha * A(ind, row)' * B(ind, col) */ 01619 Telement::gemm(true, false, 01620 alpha, A(ind, row), B(ind, col), 01621 1.0, C(row,col)); 01622 for (int ind = row + 1; ind < col; ind++) 01623 /* C(row, col) += alpha * A(row, ind) * B(ind, col) */ 01624 Telement::gemm(false, false, 01625 alpha, A(row, ind), B(ind, col), 01626 1.0, C(row,col)); 01627 for (int ind = col + 1; ind < A.ncols(); ind++) 01628 /* C(row, col) += alpha * A(row, ind) * B(col, ind)' */ 01629 Telement::gemm(false, true, 01630 alpha, A(row, ind), B(col, ind), 01631 1.0, C(row,col)); 01632 } 01633 /* Lower half */ 01634 Telement tmpSubMat; 01635 for (int row = col + 1; row < C.nrows(); row++) { 01636 Telement::transpose(C(row, col), tmpSubMat); 01637 /* tmpSubMat = alpha * B(col, col) * A(col, row) + 01638 * beta * tmpSubMat 01639 */ 01640 Telement::symm('L', 'U', 01641 alpha, B(col, col), A(col, row), 01642 beta_tmp, tmpSubMat); 01643 /* tmpSubMat += alpha * B(col, row) * A(row, row) */ 01644 Telement::symm('R', 'U', 01645 alpha, A(row, row), B(col, row), 01646 1.0, tmpSubMat); 01647 for (int ind = 0; ind < col; ind++) 01648 /* tmpSubMat += alpha * B(ind, col)' * A(ind, row) */ 01649 Telement::gemm(true, false, 01650 alpha, B(ind, col), A(ind, row), 01651 1.0, tmpSubMat); 01652 for (int ind = col + 1; ind < row; ind++) 01653 /* tmpSubMat += alpha * B(col, ind) * A(ind, row) */ 01654 Telement::gemm(false, false, 01655 alpha, B(col, ind), A(ind, row), 01656 1.0, tmpSubMat); 01657 for (int ind = row + 1; ind < B.nrows(); ind++) 01658 /* tmpSubMat += alpha * B(col, ind) * A(row, ind)' */ 01659 Telement::gemm(false, true, 01660 alpha, B(col, ind), A(row, ind), 01661 1.0, tmpSubMat); 01662 Telement::transpose(tmpSubMat, C(row, col)); 01663 } 01664 MAT_OMP_END; 01665 } 01666 MAT_OMP_FINALIZE; 01667 } /* end ssmm */ 01668 01669 01670 /* C = alpha * A * B + beta * C where A and B are symmetric 01671 * and only the upper triangle of C is computed. 01672 */ 01673 template<class Treal, class Telement> 01674 void Matrix<Treal, Telement>:: 01675 ssmm_upper_tr_only(const Treal alpha, 01676 const Matrix<Treal, Telement>& A, 01677 const Matrix<Treal, Telement>& B, 01678 const Treal beta, 01679 Matrix<Treal, Telement>& C) { 01680 assert(!A.is_empty()); 01681 assert(!B.is_empty()); 01682 if (C.is_empty()) { 01683 assert(beta == 0); 01684 C.resetRows(A.rows); 01685 C.resetCols(B.cols); 01686 } 01687 if (A.ncols() != B.nrows() || 01688 A.nrows() != C.nrows() || 01689 B.ncols() != C.ncols() || 01690 A.nrows() != A.ncols() || 01691 B.nrows() != B.ncols()) { 01692 throw Failure("Matrix<class Treal, class Telement>::" 01693 "ssmm_upper_tr_only(Treal, " 01694 "Matrix<Treal, Telement>, " 01695 "Matrix<Treal, Telement>, Treal, " 01696 "Matrix<Treal, Telement>): " 01697 "Incorrect matrixdimensions for ssmm_upper_tr_only " 01698 "multiplication"); 01699 } 01700 if ((A.is_zero() || B.is_zero() || alpha == 0) && 01701 (C.is_zero() || beta == 0)) { 01702 C = 0; 01703 return; 01704 } 01705 Treal beta_tmp = beta; 01706 if (C.is_zero()) { 01707 C.allocate(); 01708 beta_tmp = 0; 01709 } 01710 if (A.is_zero() || B.is_zero() || alpha == 0) { 01711 C *= beta; 01712 return; 01713 } 01714 MAT_OMP_INIT; 01715 #ifdef _OPENMP 01716 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 01717 #endif 01718 for (int col = 0; col < C.ncols(); col++) { 01719 MAT_OMP_START; 01720 /* Diagonal */ 01721 /* C(col, col) = alpha * A(col, col) * B(col, col) + beta * C(col, col)*/ 01722 Telement::ssmm_upper_tr_only(alpha, A(col,col), B(col, col), 01723 beta_tmp, C(col,col)); 01724 for (int ind = 0; ind < col; ind++) 01725 /* C(col, col) += alpha * A(ind, col)' * B(ind, col) */ 01726 Telement::gemm_upper_tr_only(true, false, 01727 alpha, A(ind,col), B(ind, col), 01728 1.0, C(col,col)); 01729 for (int ind = col + 1; ind < A.ncols(); ind++) 01730 /* C(col, col) += alpha * A(col, ind) * B(col, ind)' */ 01731 Telement::gemm_upper_tr_only(false, true, 01732 alpha, A(col, ind), B(col, ind), 01733 1.0, C(col,col)); 01734 /* Upper half */ 01735 for (int row = 0; row < col; row++) { 01736 /* C(row, col) = alpha * A(row, row) * B(row, col) + 01737 * beta * C(row, col) 01738 */ 01739 Telement::symm('L', 'U', 01740 alpha, A(row, row), B(row, col), 01741 beta_tmp, C(row, col)); 01742 /* C(row, col) += alpha * A(row, col) * B(col, col) */ 01743 Telement::symm('R', 'U', 01744 alpha, B(col, col), A(row, col), 01745 1.0, C(row, col)); 01746 for (int ind = 0; ind < row; ind++) 01747 /* C(row, col) += alpha * A(ind, row)' * B(ind, col) */ 01748 Telement::gemm(true, false, 01749 alpha, A(ind, row), B(ind, col), 01750 1.0, C(row,col)); 01751 for (int ind = row + 1; ind < col; ind++) 01752 /* C(row, col) += alpha * A(row, ind) * B(ind, col) */ 01753 Telement::gemm(false, false, 01754 alpha, A(row, ind), B(ind, col), 01755 1.0, C(row,col)); 01756 for (int ind = col + 1; ind < A.ncols(); ind++) 01757 /* C(row, col) += alpha * A(row, ind) * B(col, ind)' */ 01758 Telement::gemm(false, true, 01759 alpha, A(row, ind), B(col, ind), 01760 1.0, C(row,col)); 01761 } 01762 MAT_OMP_END; 01763 } 01764 MAT_OMP_FINALIZE; 01765 } /* end ssmm_upper_tr_only */ 01766 01767 01768 01769 template<class Treal, class Telement> 01770 void Matrix<Treal, Telement>:: 01771 trmm(const char side, const char uplo, const bool tA, const Treal alpha, 01772 const Matrix<Treal, Telement>& A, 01773 Matrix<Treal, Telement>& B) { 01774 assert(!A.is_empty()); 01775 assert(!B.is_empty()); 01776 if (alpha != 0 && !A.is_zero() && !B.is_zero()) 01777 if (((side == 'R' && B.ncols() == A.nrows()) || 01778 (side == 'L' && A.ncols() == B.nrows())) && 01779 A.nrows() == A.ncols()) 01780 if (uplo == 'U') 01781 if (!tA) { 01782 if (side == 'R') { 01783 /* Last column must be calculated first */ 01784 for (int col = B.ncols() - 1; col >= 0; col--) 01785 for (int row = 0; row < B.nrows(); row++) { 01786 /* Use first the diagonal element in A */ 01787 /* Otherwise the faster call to trmm cannot be utilized */ 01788 Telement::trmm(side, uplo, tA, alpha, 01789 A(col, col), B(row,col)); 01790 /* And the rest: */ 01791 for (int ind = 0; ind < col; ind++) 01792 Telement::gemm(false, tA, alpha, 01793 B(row,ind), A(ind, col), 01794 1.0, B(row,col)); 01795 } 01796 } /* end if (side == 'R')*/ 01797 else { 01798 assert(side == 'L'); 01799 /* First row must be calculated first */ 01800 for (int row = 0; row < B.nrows(); row++ ) 01801 for (int col = 0; col < B.ncols(); col++) { 01802 Telement::trmm(side, uplo, tA, alpha, 01803 A(row,row), B(row,col)); 01804 for (int ind = row + 1 ; ind < B.nrows(); ind++) 01805 Telement::gemm(tA, false, alpha, 01806 A(row,ind), B(ind,col), 01807 1.0, B(row,col)); 01808 } 01809 } /* end else (side == 'L')*/ 01810 } /* end if (tA == false) */ 01811 else { 01812 assert(tA == true); 01813 if (side == 'R') { 01814 /* First column must be calculated first */ 01815 for (int col = 0; col < B.ncols(); col++) 01816 for (int row = 0; row < B.nrows(); row++) { 01817 Telement::trmm(side, uplo, tA, alpha, 01818 A(col,col), B(row,col)); 01819 for (int ind = col + 1; ind < A.ncols(); ind++) 01820 Telement::gemm(false, tA, alpha, 01821 B(row,ind), A(col,ind), 01822 1.0, B(row,col)); 01823 } 01824 } /* end if (side == 'R')*/ 01825 else { 01826 assert(side == 'L'); 01827 /* Last row must be calculated first */ 01828 for (int row = B.nrows() - 1; row >= 0; row--) 01829 for (int col = 0; col < B.ncols(); col++) { 01830 Telement::trmm(side, uplo, tA, alpha, 01831 A(row,row), B(row,col)); 01832 for (int ind = 0; ind < row; ind++) 01833 Telement::gemm(tA, false, alpha, 01834 A(ind,row), B(ind,col), 01835 1.0, B(row,col)); 01836 } 01837 } /* end else (side == 'L')*/ 01838 } /* end else (tA == true)*/ 01839 else 01840 throw Failure("Matrix<class Treal, class Telement>::" 01841 "trmm not implemented for lower triangular matrices"); 01842 else 01843 throw Failure("Matrix<class Treal, class Telement>::trmm" 01844 ": Incorrect matrix dimensions for multiplication"); 01845 else 01846 B = 0; 01847 } 01848 01849 01850 template<class Treal, class Telement> 01851 Treal Matrix<Treal, Telement>::frobSquared() const { 01852 assert(!this->is_empty()); 01853 if (this->is_zero()) 01854 return 0; 01855 else { 01856 Treal sum(0); 01857 for (int i = 0; i < this->nElements(); i++) 01858 sum += this->elements[i].frobSquared(); 01859 return sum; 01860 } 01861 } 01862 01863 template<class Treal, class Telement> 01864 Treal Matrix<Treal, Telement>:: 01865 syFrobSquared() const { 01866 assert(!this->is_empty()); 01867 if (this->nrows() != this->ncols()) 01868 throw Failure("Matrix<Treal, Telement>::syFrobSquared(): " 01869 "Matrix must be have equal number of rows " 01870 "and cols to be symmetric"); 01871 Treal sum(0); 01872 if (!this->is_zero()) { 01873 for (int col = 1; col < this->ncols(); col++) 01874 for (int row = 0; row < col; row++) 01875 sum += 2 * (*this)(row, col).frobSquared(); 01876 for (int rc = 0; rc < this->ncols(); rc++) 01877 sum += (*this)(rc, rc).syFrobSquared(); 01878 } 01879 return sum; 01880 } 01881 01882 template<class Treal, class Telement> 01883 Treal Matrix<Treal, Telement>:: 01884 frobSquaredDiff 01885 (const Matrix<Treal, Telement>& A, 01886 const Matrix<Treal, Telement>& B) { 01887 assert(!A.is_empty()); 01888 assert(!B.is_empty()); 01889 if (A.nrows() != B.nrows() || A.ncols() != B.ncols()) 01890 throw Failure("Matrix<Treal, Telement>::" 01891 "frobSquaredDiff: Incorrect matrix dimensions"); 01892 Treal sum(0); 01893 if (!A.is_zero() && !B.is_zero()) 01894 for (int i = 0; i < A.nElements(); i++) 01895 sum += Telement::frobSquaredDiff(A.elements[i],B.elements[i]); 01896 else if (A.is_zero() && !B.is_zero()) 01897 sum = B.frobSquared(); 01898 else if (!A.is_zero() && B.is_zero()) 01899 sum = A.frobSquared(); 01900 /* sum is already zero if A.is_zero() && B.is_zero() */ 01901 return sum; 01902 } 01903 01904 template<class Treal, class Telement> 01905 Treal Matrix<Treal, Telement>:: 01906 syFrobSquaredDiff 01907 (const Matrix<Treal, Telement>& A, 01908 const Matrix<Treal, Telement>& B) { 01909 assert(!A.is_empty()); 01910 assert(!B.is_empty()); 01911 if (A.nrows() != B.nrows() || 01912 A.ncols() != B.ncols() || 01913 A.nrows() != A.ncols()) 01914 throw Failure("Matrix<Treal, Telement>::syFrobSquaredDiff: " 01915 "Incorrect matrix dimensions"); 01916 Treal sum(0); 01917 if (!A.is_zero() && !B.is_zero()) { 01918 for (int col = 1; col < A.ncols(); col++) 01919 for (int row = 0; row < col; row++) 01920 sum += 2 * Telement::frobSquaredDiff(A(row, col), B(row, col)); 01921 for (int rc = 0; rc < A.ncols(); rc++) 01922 sum += Telement::syFrobSquaredDiff(A(rc, rc),B(rc, rc)); 01923 } 01924 else if (A.is_zero() && !B.is_zero()) 01925 sum = B.syFrobSquared(); 01926 else if (!A.is_zero() && B.is_zero()) 01927 sum = A.syFrobSquared(); 01928 /* sum is already zero if A.is_zero() && B.is_zero() */ 01929 return sum; 01930 } 01931 01932 template<class Treal, class Telement> 01933 Treal Matrix<Treal, Telement>:: 01934 trace() const { 01935 assert(!this->is_empty()); 01936 if (this->nrows() != this->ncols()) 01937 throw Failure("Matrix<Treal, Telement>::trace(): " 01938 "Matrix must be quadratic"); 01939 if (this->is_zero()) 01940 return 0; 01941 else { 01942 Treal sum = 0; 01943 for (int rc = 0; rc < this->ncols(); rc++) 01944 sum += (*this)(rc,rc).trace(); 01945 return sum; 01946 } 01947 } 01948 01949 template<class Treal, class Telement> 01950 Treal Matrix<Treal, Telement>:: 01951 trace_ab(const Matrix<Treal, Telement>& A, 01952 const Matrix<Treal, Telement>& B) { 01953 assert(!A.is_empty()); 01954 assert(!B.is_empty()); 01955 if (A.nrows() != B.ncols() || A.ncols() != B.nrows()) 01956 throw Failure("Matrix<Treal, Telement>::trace_ab: " 01957 "Wrong matrix dimensions for trace(A * B)"); 01958 Treal tr = 0; 01959 if (!A.is_zero() && !B.is_zero()) 01960 for (int rc = 0; rc < A.nrows(); rc++) 01961 for (int colA = 0; colA < A.ncols(); colA++) 01962 tr += Telement::trace_ab(A(rc,colA), B(colA, rc)); 01963 return tr; 01964 } 01965 01966 template<class Treal, class Telement> 01967 Treal Matrix<Treal, Telement>:: 01968 trace_aTb(const Matrix<Treal, Telement>& A, 01969 const Matrix<Treal, Telement>& B) { 01970 assert(!A.is_empty()); 01971 assert(!B.is_empty()); 01972 if (A.ncols() != B.ncols() || A.nrows() != B.nrows()) 01973 throw Failure("Matrix<Treal, Telement>::trace_aTb: " 01974 "Wrong matrix dimensions for trace(A' * B)"); 01975 Treal tr = 0; 01976 if (!A.is_zero() && !B.is_zero()) 01977 for (int rc = 0; rc < A.ncols(); rc++) 01978 for (int rowB = 0; rowB < B.nrows(); rowB++) 01979 tr += Telement::trace_aTb(A(rowB,rc), B(rowB, rc)); 01980 return tr; 01981 } 01982 01983 template<class Treal, class Telement> 01984 Treal Matrix<Treal, Telement>:: 01985 sy_trace_ab(const Matrix<Treal, Telement>& A, 01986 const Matrix<Treal, Telement>& B) { 01987 assert(!A.is_empty()); 01988 assert(!B.is_empty()); 01989 if (A.nrows() != B.ncols() || A.ncols() != B.nrows() || 01990 A.nrows() != A.ncols()) 01991 throw Failure("Matrix<Treal, Telement>::sy_trace_ab: " 01992 "Wrong matrix dimensions for symmetric trace(A * B)"); 01993 Treal tr = 0; 01994 if (!A.is_zero() && !B.is_zero()) { 01995 /* Diagonal first */ 01996 for (int rc = 0; rc < A.nrows(); rc++) 01997 tr += Telement::sy_trace_ab(A(rc,rc), B(rc, rc)); 01998 /* Using that trace of transpose is equal to that without transpose: */ 01999 for (int colA = 1; colA < A.ncols(); colA++) 02000 for (int rowA = 0; rowA < colA; rowA++) 02001 tr += 2 * Telement::trace_aTb(A(rowA, colA), B(rowA, colA)); 02002 } 02003 return tr; 02004 } 02005 02006 template<class Treal, class Telement> 02007 void Matrix<Treal, Telement>:: 02008 add(const Treal alpha, /* B += alpha * A */ 02009 const Matrix<Treal, Telement>& A, 02010 Matrix<Treal, Telement>& B) { 02011 assert(!A.is_empty()); 02012 assert(!B.is_empty()); 02013 if (A.nrows() != B.nrows() || A.ncols() != B.ncols()) 02014 throw Failure("Matrix<Treal, Telement>::add: " 02015 "Wrong matrix dimensions for addition"); 02016 if (!A.is_zero() && alpha != 0) { 02017 if (B.is_zero()) 02018 B.allocate(); 02019 for (int ind = 0; ind < A.nElements(); ind++) 02020 Telement::add(alpha, A.elements[ind], B.elements[ind]); 02021 } 02022 } 02023 02024 template<class Treal, class Telement> 02025 void Matrix<Treal, Telement>:: 02026 assign(Treal const alpha, /* *this = alpha * A */ 02027 Matrix<Treal, Telement> const & A) { 02028 assert(!A.is_empty()); 02029 if (this->is_empty()) { 02030 this->resetRows(A.rows); 02031 this->resetCols(A.cols); 02032 } 02033 *this = 0; 02034 Matrix<Treal, Telement>:: 02035 add(alpha, A, *this); 02036 } 02037 02038 02039 /********** Help functions for thresholding */ 02040 02041 02042 template<class Treal, class Telement> 02043 void Matrix<Treal, Telement>:: 02044 getFrobSqLowestLevel(std::vector<Treal> & frobsq) const { 02045 if (!this->is_zero()) 02046 for (int ind = 0; ind < this->nElements(); ind++) 02047 this->elements[ind].getFrobSqLowestLevel(frobsq); 02048 } 02049 02050 template<class Treal, class Telement> 02051 void Matrix<Treal, Telement>:: 02052 frobThreshLowestLevel 02053 (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix) { 02054 assert(!this->is_empty()); 02055 bool thisMatIsZero = true; 02056 if (ErrorMatrix) { 02057 assert(!ErrorMatrix->is_empty()); 02058 bool EMatIsZero = true; 02059 if (!ErrorMatrix->is_zero() || !this->is_zero()) { 02060 if (ErrorMatrix->is_zero()) 02061 ErrorMatrix->allocate(); 02062 if (this->is_zero()) 02063 this->allocate(); 02064 for (int ind = 0; ind < this->nElements(); ind++) { 02065 this->elements[ind]. 02066 frobThreshLowestLevel(threshold, &ErrorMatrix->elements[ind]); 02067 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero(); 02068 EMatIsZero = EMatIsZero && ErrorMatrix->elements[ind].is_zero(); 02069 } 02070 if (thisMatIsZero) 02071 this->clear(); 02072 if (EMatIsZero) 02073 ErrorMatrix->clear(); 02074 } 02075 } 02076 else 02077 if (!this->is_zero()) { 02078 for (int ind = 0; ind < this->nElements(); ind++) { 02079 this->elements[ind].frobThreshLowestLevel(threshold, 0); 02080 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero(); 02081 } 02082 if (thisMatIsZero) 02083 this->clear(); 02084 } 02085 } 02086 02087 template<class Treal, class Telement> 02088 void Matrix<Treal, Telement>:: 02089 getFrobSqElementLevel(std::vector<Treal> & frobsq) const { 02090 if (!this->is_zero()) 02091 for (int ind = 0; ind < this->nElements(); ind++) 02092 this->elements[ind].getFrobSqElementLevel(frobsq); 02093 } 02094 02095 template<class Treal, class Telement> 02096 void Matrix<Treal, Telement>:: 02097 frobThreshElementLevel 02098 (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix) { 02099 assert(!this->is_empty()); 02100 bool thisMatIsZero = true; 02101 if (ErrorMatrix) { 02102 assert(!ErrorMatrix->is_empty()); 02103 bool EMatIsZero = true; 02104 if (!ErrorMatrix->is_zero() || !this->is_zero()) { 02105 if (ErrorMatrix->is_zero()) 02106 ErrorMatrix->allocate(); 02107 if (this->is_zero()) 02108 this->allocate(); 02109 for (int ind = 0; ind < this->nElements(); ind++) { 02110 this->elements[ind]. 02111 frobThreshElementLevel(threshold, &ErrorMatrix->elements[ind]); 02112 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero(); 02113 EMatIsZero = EMatIsZero && ErrorMatrix->elements[ind].is_zero(); 02114 } 02115 if (thisMatIsZero) 02116 this->clear(); 02117 if (EMatIsZero) 02118 ErrorMatrix->clear(); 02119 } 02120 } 02121 else 02122 if (!this->is_zero()) { 02123 for (int ind = 0; ind < this->nElements(); ind++) { 02124 this->elements[ind].frobThreshElementLevel(threshold, 0); 02125 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero(); 02126 } 02127 if (thisMatIsZero) 02128 this->clear(); 02129 } 02130 } 02131 02132 02133 02134 template<class Treal, class Telement> 02135 void Matrix<Treal, Telement>::assignFrobNormsLowestLevel 02136 ( Matrix<Treal, Matrix<Treal, Telement> > const & A ) { 02137 if (!A.is_zero()) { 02138 if ( this->is_zero() ) 02139 this->allocate(); 02140 assert( this->nElements() == A.nElements() ); 02141 for (int ind = 0; ind < this->nElements(); ind++) 02142 this->elements[ind].assignFrobNormsLowestLevel(A[ind]); 02143 } 02144 else 02145 this->clear(); 02146 } 02147 02148 template<class Treal, class Telement> 02149 void Matrix<Treal, Telement>::syAssignFrobNormsLowestLevel 02150 ( Matrix<Treal, Matrix<Treal, Telement> > const & A ) { 02151 assert(!A.is_empty()); 02152 if (A.nrows() != A.ncols()) 02153 throw Failure("Matrix<Treal, Telement>::syAssignFrobNormsLowestLevel(...): " 02154 "Matrix must be have equal number of rows " 02155 "and cols to be symmetric"); 02156 if (!A.is_zero()) { 02157 if ( this->is_zero() ) 02158 this->allocate(); 02159 assert( this->nElements() == A.nElements() ); 02160 for (int col = 1; col < this->ncols(); col++) 02161 for (int row = 0; row < col; row++) 02162 (*this)(row, col).assignFrobNormsLowestLevel( A(row,col) ); 02163 for (int rc = 0; rc < this->ncols(); rc++) 02164 (*this)(rc, rc).syAssignFrobNormsLowestLevel( A(rc,rc) ); 02165 } 02166 else 02167 this->clear(); 02168 } 02169 02170 template<class Treal, class Telement> 02171 void Matrix<Treal, Telement>::assignDiffFrobNormsLowestLevel 02172 ( Matrix<Treal, Matrix<Treal, Telement> > const & A, 02173 Matrix<Treal, Matrix<Treal, Telement> > const & B ) { 02174 if ( A.is_zero() && B.is_zero() ) { 02175 // Both A and B are zero 02176 this->clear(); 02177 return; 02178 } 02179 // At least one of A and B is nonzero 02180 if ( this->is_zero() ) 02181 this->allocate(); 02182 if ( !A.is_zero() && !B.is_zero() ) { 02183 // Both are nonzero 02184 assert( this->nElements() == A.nElements() ); 02185 assert( this->nElements() == B.nElements() ); 02186 for (int ind = 0; ind < this->nElements(); ind++) 02187 this->elements[ind].assignDiffFrobNormsLowestLevel( A[ind], B[ind] ); 02188 return; 02189 } 02190 if ( !A.is_zero() ) { 02191 // A is nonzero 02192 this->assignFrobNormsLowestLevel( A ); 02193 return; 02194 } 02195 if ( !B.is_zero() ) { 02196 // B is nonzero 02197 this->assignFrobNormsLowestLevel( B ); 02198 return; 02199 } 02200 } 02201 template<class Treal, class Telement> 02202 void Matrix<Treal, Telement>::syAssignDiffFrobNormsLowestLevel 02203 ( Matrix<Treal, Matrix<Treal, Telement> > const & A, 02204 Matrix<Treal, Matrix<Treal, Telement> > const & B ) { 02205 if ( A.is_zero() && B.is_zero() ) { 02206 // Both A and B are zero 02207 this->clear(); 02208 return; 02209 } 02210 // At least one of A and B is nonzero 02211 if ( this->is_zero() ) 02212 this->allocate(); 02213 if ( !A.is_zero() && !B.is_zero() ) { 02214 // Both are nonzero 02215 assert( this->nElements() == A.nElements() ); 02216 assert( this->nElements() == B.nElements() ); 02217 for (int col = 1; col < this->ncols(); col++) 02218 for (int row = 0; row < col; row++) 02219 (*this)(row, col).assignDiffFrobNormsLowestLevel( A(row,col), B(row,col) ); 02220 for (int rc = 0; rc < this->ncols(); rc++) 02221 (*this)(rc, rc).syAssignDiffFrobNormsLowestLevel( A(rc,rc), B(rc,rc) ); 02222 return; 02223 } 02224 if ( !A.is_zero() ) { 02225 // A is nonzero 02226 this->syAssignFrobNormsLowestLevel( A ); 02227 return; 02228 } 02229 if ( !B.is_zero() ) { 02230 // B is nonzero 02231 this->syAssignFrobNormsLowestLevel( B ); 02232 return; 02233 } 02234 } 02235 02236 02237 02238 template<class Treal, class Telement> 02239 void Matrix<Treal, Telement>:: 02240 truncateAccordingToSparsityPattern( Matrix<Treal, Matrix<Treal, Telement> > & A ) const { 02241 if ( this->is_zero() ) { 02242 A.clear(); 02243 } 02244 else { 02245 assert( !A.is_zero() ); 02246 assert( this->nElements() == A.nElements() ); 02247 for (int ind = 0; ind < this->nElements(); ind++) 02248 this->elements[ind].truncateAccordingToSparsityPattern( A[ind] ); 02249 } 02250 } 02251 02252 02253 02254 /********** End of help functions for thresholding */ 02255 02256 /* C = beta * C + alpha * A * B where only the upper triangle of C is */ 02257 /* referenced and updated */ 02258 template<class Treal, class Telement> 02259 void Matrix<Treal, Telement>:: 02260 gemm_upper_tr_only(const bool tA, const bool tB, const Treal alpha, 02261 const Matrix<Treal, Telement>& A, 02262 const Matrix<Treal, Telement>& B, 02263 const Treal beta, 02264 Matrix<Treal, Telement>& C) { 02265 assert(!A.is_empty()); 02266 assert(!B.is_empty()); 02267 if (C.is_empty()) { 02268 assert(beta == 0); 02269 if (tA) 02270 C.resetRows(A.cols); 02271 else 02272 C.resetRows(A.rows); 02273 if (tB) 02274 C.resetCols(B.rows); 02275 else 02276 C.resetCols(B.cols); 02277 } 02278 if ((A.is_zero() || B.is_zero() || alpha == 0) && 02279 (C.is_zero() || beta == 0)) 02280 C = 0; 02281 else { 02282 Treal beta_tmp = beta; 02283 if (C.is_zero()) { 02284 C.allocate(); 02285 beta_tmp = 0; 02286 } 02287 if (!A.is_zero() && !B.is_zero() && alpha != 0) { 02288 if (!tA && !tB) 02289 if (A.ncols() == B.nrows() && 02290 A.nrows() == C.nrows() && 02291 B.ncols() == C.ncols()) { 02292 for (int col = 0; col < C.ncols(); col++) { 02293 Telement::gemm_upper_tr_only(tA, tB, alpha, 02294 A(col, 0), B(0, col), 02295 beta_tmp, 02296 C(col, col)); 02297 for(int cola = 1; cola < A.ncols(); cola++) 02298 Telement::gemm_upper_tr_only(tA, tB, alpha, 02299 A(col, cola), B(cola, col), 02300 1.0, 02301 C(col,col)); 02302 for (int row = 0; row < col; row++) { 02303 Telement::gemm(tA, tB, alpha, 02304 A(row, 0), B(0, col), 02305 beta_tmp, 02306 C(row,col)); 02307 for(int cola = 1; cola < A.ncols(); cola++) 02308 Telement::gemm(tA, tB, alpha, 02309 A(row, cola), B(cola, col), 02310 1.0, 02311 C(row,col)); 02312 } 02313 } 02314 } 02315 else 02316 throw Failure("Matrix<class Treal, class Telement>::" 02317 "gemm_upper_tr_only(bool, bool, Treal, " 02318 "Matrix<Treal, Telement>, " 02319 "Matrix<Treal, Telement>, " 02320 "Treal, Matrix<Treal, Telement>): " 02321 "Incorrect matrixdimensions for multiplication"); 02322 else if (tA && !tB) 02323 if (A.nrows() == B.nrows() && 02324 A.ncols() == C.nrows() && 02325 B.ncols() == C.ncols()) { 02326 for (int col = 0; col < C.ncols(); col++) { 02327 Telement::gemm_upper_tr_only(tA, tB, alpha, 02328 A(0,col), B(0,col), 02329 beta_tmp, 02330 C(col,col)); 02331 for(int cola = 1; cola < A.nrows(); cola++) 02332 Telement::gemm_upper_tr_only(tA, tB, alpha, 02333 A(cola, col), B(cola, col), 02334 1.0, 02335 C(col, col)); 02336 for (int row = 0; row < col; row++) { 02337 Telement::gemm(tA, tB, alpha, 02338 A(0,row), B(0,col), 02339 beta_tmp, 02340 C(row,col)); 02341 for(int cola = 1; cola < A.nrows(); cola++) 02342 Telement::gemm(tA, tB, alpha, 02343 A(cola, row), B(cola, col), 02344 1.0, 02345 C(row,col)); 02346 } 02347 } 02348 } 02349 else 02350 throw Failure("Matrix<class Treal, class Telement>::" 02351 "gemm_upper_tr_only(bool, bool, Treal, " 02352 "Matrix<Treal, Telement>, " 02353 "Matrix<Treal, Telement>, " 02354 "Treal, Matrix<Treal, Telement>): " 02355 "Incorrect matrixdimensions for multiplication"); 02356 else if (!tA && tB) 02357 if (A.ncols() == B.ncols() && 02358 A.nrows() == C.nrows() && 02359 B.nrows() == C.ncols()) { 02360 for (int col = 0; col < C.ncols(); col++) { 02361 Telement::gemm_upper_tr_only(tA, tB, alpha, 02362 A(col, 0), B(col, 0), 02363 beta_tmp, 02364 C(col,col)); 02365 for(int cola = 1; cola < A.ncols(); cola++) 02366 Telement::gemm_upper_tr_only(tA, tB, alpha, 02367 A(col, cola), B(col, cola), 02368 1.0, 02369 C(col,col)); 02370 for (int row = 0; row < col; row++) { 02371 Telement::gemm(tA, tB, alpha, 02372 A(row, 0), B(col, 0), 02373 beta_tmp, 02374 C(row,col)); 02375 for(int cola = 1; cola < A.ncols(); cola++) 02376 Telement::gemm(tA, tB, alpha, 02377 A(row, cola), B(col, cola), 02378 1.0, 02379 C(row,col)); 02380 } 02381 } 02382 } 02383 else 02384 throw Failure("Matrix<class Treal, class Telement>::" 02385 "gemm_upper_tr_only(bool, bool, Treal, " 02386 "Matrix<Treal, Telement>, " 02387 "Matrix<Treal, Telement>, " 02388 "Treal, Matrix<Treal, Telement>): " 02389 "Incorrect matrixdimensions for multiplication"); 02390 else if (tA && tB) 02391 if (A.nrows() == B.ncols() && 02392 A.ncols() == C.nrows() && 02393 B.nrows() == C.ncols()) { 02394 for (int col = 0; col < C.ncols(); col++) { 02395 Telement::gemm_upper_tr_only(tA, tB, alpha, 02396 A(0, col), B(col, 0), 02397 beta_tmp, 02398 C(col,col)); 02399 for(int cola = 1; cola < A.nrows(); cola++) 02400 Telement::gemm_upper_tr_only(tA, tB, alpha, 02401 A(cola, col), B(col, cola), 02402 1.0, 02403 C(col,col)); 02404 for (int row = 0; row < col; row++) { 02405 Telement::gemm(tA, tB, alpha, 02406 A(0, row), B(col, 0), 02407 beta_tmp, 02408 C(row,col)); 02409 for(int cola = 1; cola < A.nrows(); cola++) 02410 Telement::gemm(tA, tB, alpha, 02411 A(cola, row), B(col, cola), 02412 1.0, 02413 C(row,col)); 02414 } 02415 } 02416 } 02417 else 02418 throw Failure("Matrix<class Treal, class Telement>::" 02419 "gemm_upper_tr_only(bool, bool, Treal, " 02420 "Matrix<Treal, Telement>, " 02421 "Matrix<Treal, Telement>, Treal, " 02422 "Matrix<Treal, Telement>): " 02423 "Incorrect matrixdimensions for multiplication"); 02424 else throw Failure("Matrix<class Treal, class Telement>::" 02425 "gemm_upper_tr_only(bool, bool, Treal, " 02426 "Matrix<Treal, Telement>, " 02427 "Matrix<Treal, Telement>, Treal, " 02428 "Matrix<Treal, Telement>):" 02429 "Very strange error!!"); 02430 } 02431 else 02432 C *= beta; 02433 } 02434 } 02435 02436 /* A = alpha * A * Z or A = alpha * Z * A where A is symmetric, */ 02437 /* Z is upper triangular and */ 02438 /* only the upper triangle of the result is calculated */ 02439 /* side == 'R' for A = alpha * A * Z */ 02440 /* side == 'L' for A = alpha * Z * A */ 02441 template<class Treal, class Telement> 02442 void Matrix<Treal, Telement>:: 02443 sytr_upper_tr_only(char const side, const Treal alpha, 02444 Matrix<Treal, Telement>& A, 02445 const Matrix<Treal, Telement>& Z) { 02446 assert(!A.is_empty()); 02447 assert(!Z.is_empty()); 02448 if (alpha != 0 && !A.is_zero() && !Z.is_zero()) { 02449 if (A.nrows() == A.ncols() && 02450 Z.nrows() == Z.ncols() && 02451 A.nrows() == Z.nrows()) { 02452 if (side == 'R') { 02453 /* Last column must be calculated first */ 02454 for (int col = A.ncols() - 1; col >= 0; col--) { 02455 // A(col, col) = alpha * A(col, col) * Z(col, col) 02456 Telement::sytr_upper_tr_only(side, alpha, 02457 A(col, col), Z(col, col)); 02458 for (int ind = 0; ind < col; ind++) { 02459 // A(col, col) += alpha * A(ind, col)' * Z(ind, col) 02460 Telement::gemm_upper_tr_only(true, false, alpha, A(ind, col), 02461 Z(ind, col), 1.0, A(col, col)); 02462 } 02463 /* Last row must be calculated first */ 02464 for (int row = col - 1; row >= 0; row--) { 02465 // A(row, col) = alpha * A(row, col) * Z(col, col); 02466 Telement::trmm(side, 'U', false, 02467 alpha, Z(col, col), A(row, col)); 02468 // A(row, col) += alpha * A(row, row) * Z(row, col); 02469 Telement::symm('L', 'U', alpha, A(row, row), Z(row, col), 02470 1.0, A(row, col)); 02471 for (int ind = 0; ind < row; ind++) { 02472 // A(row, col) += alpha * A(ind, row)' * Z(ind, col); 02473 Telement::gemm(true, false, alpha, A(ind, row), Z(ind, col), 02474 1.0, A(row, col)); 02475 } 02476 for (int ind = row + 1; ind < col; ind++) { 02477 // A(row, col) += alpha * A(row, ind) * Z(ind, col); 02478 Telement::gemm(false, false, alpha, A(row, ind), Z(ind, col), 02479 1.0, A(row, col)); 02480 } 02481 } 02482 } 02483 } 02484 else { /* side == 'L' */ 02485 assert(side == 'L'); 02486 /* First column must be calculated first */ 02487 for (int col = 0; col < A.ncols(); col++) { 02488 /* First row must be calculated first */ 02489 for (int row = 0; row < col; row++) { 02490 // A(row, col) = alpha * Z(row, row) * A(row, col) 02491 Telement::trmm(side, 'U', false, alpha, 02492 Z(row, row), A(row, col)); 02493 // A(row, col) += alpha * Z(row, col) * A(col, col) 02494 Telement::symm('R', 'U', alpha, A(col, col), Z(row, col), 02495 1.0, A(row, col)); 02496 for (int ind = row + 1; ind < col; ind++) 02497 // A(row, col) += alpha * Z(row, ind) * A(ind, col) 02498 Telement::gemm(false, false, alpha, Z(row, ind), A(ind, col), 02499 1.0, A(row, col)); 02500 for (int ind = col + 1; ind < A.nrows(); ind++) 02501 // A(row, col) += alpha * Z(row, ind) * A(col, ind)' 02502 Telement::gemm(false, true, alpha, Z(row, ind), A(col, ind), 02503 1.0, A(row, col)); 02504 } 02505 Telement::sytr_upper_tr_only(side, alpha, 02506 A(col, col), Z(col, col)); 02507 for (int ind = col + 1; ind < A.ncols(); ind++) 02508 Telement::gemm_upper_tr_only(false, true, alpha, Z(col, ind), 02509 A(col, ind), 1.0, A(col, col)); 02510 } 02511 } 02512 } 02513 else 02514 throw Failure("Matrix<class Treal, class Telement>::" 02515 "sytr_upper_tr_only: Incorrect matrix dimensions " 02516 "for symmetric-triangular multiplication"); 02517 } 02518 else 02519 A = 0; 02520 } 02521 02522 /* The result B is assumed to be symmetric, i.e. only upper triangle */ 02523 /* calculated and hence only upper triangle of input B is used. */ 02524 template<class Treal, class Telement> 02525 void Matrix<Treal, Telement>:: 02526 trmm_upper_tr_only(const char side, const char uplo, 02527 const bool tA, const Treal alpha, 02528 const Matrix<Treal, Telement>& A, 02529 Matrix<Treal, Telement>& B) { 02530 assert(!A.is_empty()); 02531 assert(!B.is_empty()); 02532 if (alpha != 0 && !A.is_zero() && !B.is_zero()) 02533 if (((side == 'R' && B.ncols() == A.nrows()) || 02534 (side == 'L' && A.ncols() == B.nrows())) && 02535 A.nrows() == A.ncols()) 02536 if (uplo == 'U') 02537 if (!tA) { 02538 throw Failure("Matrix<Treal, class Telement>::" 02539 "trmm_upper_tr_only : " 02540 "not possible for upper triangular not transposed " 02541 "matrices A since lower triangle of B is needed"); 02542 } /* end if (tA == false) */ 02543 else { 02544 assert(tA == true); 02545 if (side == 'R') { 02546 /* First column must be calculated first */ 02547 for (int col = 0; col < B.ncols(); col++) { 02548 Telement::trmm_upper_tr_only(side, uplo, tA, alpha, 02549 A(col,col), B(col,col)); 02550 for (int ind = col + 1; ind < A.ncols(); ind++) 02551 Telement::gemm_upper_tr_only(false, tA, alpha, 02552 B(col,ind), A(col,ind), 02553 1.0, B(col,col)); 02554 for (int row = 0; row < col; row++) { 02555 Telement::trmm(side, uplo, tA, alpha, 02556 A(col,col), B(row,col)); 02557 for (int ind = col + 1; ind < A.ncols(); ind++) 02558 Telement::gemm(false, tA, alpha, 02559 B(row,ind), A(col,ind), 02560 1.0, B(row,col)); 02561 } 02562 } 02563 } /* end if (side == 'R')*/ 02564 else { 02565 assert(side == 'L'); 02566 /* Last row must be calculated first */ 02567 for (int row = B.nrows() - 1; row >= 0; row--) { 02568 Telement::trmm_upper_tr_only(side, uplo, tA, alpha, 02569 A(row,row), B(row,row)); 02570 for (int ind = 0; ind < row; ind++) 02571 Telement::gemm_upper_tr_only(tA, false, alpha, 02572 A(ind,row), B(ind,row), 02573 1.0, B(row,row)); 02574 for (int col = row + 1; col < B.ncols(); col++) { 02575 Telement::trmm(side, uplo, tA, alpha, 02576 A(row,row), B(row,col)); 02577 for (int ind = 0; ind < row; ind++) 02578 Telement::gemm(tA, false, alpha, 02579 A(ind,row), B(ind,col), 02580 1.0, B(row,col)); 02581 } 02582 } 02583 } /* end else (side == 'L')*/ 02584 } /* end else (tA == true)*/ 02585 else 02586 throw Failure("Matrix<class Treal, class Telement>::" 02587 "trmm_upper_tr_only not implemented for lower " 02588 "triangular matrices"); 02589 else 02590 throw Failure("Matrix<class Treal, class Telement>::" 02591 "trmm_upper_tr_only: Incorrect matrix dimensions " 02592 "for multiplication"); 02593 else 02594 B = 0; 02595 } 02596 02597 /* A = Z' * A * Z or A = Z * A * Z' */ 02598 /* where A is symmetric and Z is (nonunit) upper triangular */ 02599 /* side == 'R' for A = Z' * A * Z */ 02600 /* side == 'L' for A = Z * A * Z' */ 02601 template<class Treal, class Telement> 02602 void Matrix<Treal, Telement>:: 02603 trsytriplemm(char const side, 02604 const Matrix<Treal, Telement>& Z, 02605 Matrix<Treal, Telement>& A) { 02606 if (side == 'R') { 02607 Matrix<Treal, Telement>:: 02608 sytr_upper_tr_only('R', 1.0, A, Z); 02609 Matrix<Treal, Telement>:: 02610 trmm_upper_tr_only('L', 'U', true, 1.0, Z, A); 02611 } 02612 else { 02613 assert(side == 'L'); 02614 Matrix<Treal, Telement>:: 02615 sytr_upper_tr_only('L', 1.0, A, Z); 02616 Matrix<Treal, Telement>:: 02617 trmm_upper_tr_only('R', 'U', true, 1.0, Z, A); 02618 } 02619 } 02620 02621 template<class Treal, class Telement> 02622 Treal Matrix<Treal, Telement>:: 02623 frob_squared_thresh(Treal const threshold, 02624 Matrix<Treal, Telement> * ErrorMatrix) { 02625 assert(!this->is_empty()); 02626 if (ErrorMatrix && ErrorMatrix->is_empty()) { 02627 ErrorMatrix->resetRows(this->rows); 02628 ErrorMatrix->resetCols(this->cols); 02629 } 02630 assert(threshold >= (Treal)0.0); 02631 if (threshold == (Treal)0.0) 02632 return 0; 02633 02634 std::vector<Treal> frobsq_norms; 02635 getFrobSqLowestLevel(frobsq_norms); 02636 /* Sort in ascending order */ 02637 std::sort(frobsq_norms.begin(), frobsq_norms.end()); 02638 int lastRemoved = -1; 02639 Treal frobsqSum = 0; 02640 int nnzBlocks = frobsq_norms.size(); 02641 while(lastRemoved < nnzBlocks - 1 && frobsqSum < threshold) { 02642 ++lastRemoved; 02643 frobsqSum += frobsq_norms[lastRemoved]; 02644 } 02645 02646 /* Check if entire matrix will be removed */ 02647 if (lastRemoved == nnzBlocks - 1 && frobsqSum < threshold) { 02648 if (ErrorMatrix) 02649 Matrix<Treal, Telement>::swap(*this, *ErrorMatrix); 02650 else 02651 *this = 0; 02652 } 02653 else { 02654 frobsqSum -= frobsq_norms[lastRemoved]; 02655 --lastRemoved; 02656 while(lastRemoved > -1 && frobsq_norms[lastRemoved] == 02657 frobsq_norms[lastRemoved + 1]) { 02658 frobsqSum -= frobsq_norms[lastRemoved]; 02659 --lastRemoved; 02660 } 02661 if (lastRemoved > -1) { 02662 Treal threshLowestLevel = 02663 (frobsq_norms[lastRemoved + 1] + 02664 frobsq_norms[lastRemoved]) / 2; 02665 this->frobThreshLowestLevel(threshLowestLevel, ErrorMatrix); 02666 } 02667 } 02668 return frobsqSum; 02669 } 02670 02671 template<class Treal, class Telement> 02672 void Matrix<Treal, Telement>:: 02673 syInch(const Matrix<Treal, Telement>& A, 02674 Matrix<Treal, Telement>& Z, 02675 const Treal threshold, const side looking, 02676 const inchversion version) { 02677 assert(!A.is_empty()); 02678 if (A.nrows() != A.ncols()) 02679 throw Failure("Matrix<Treal, Telement>::sy_inch: " 02680 "Matrix must be quadratic!"); 02681 if (A.is_zero()) 02682 throw Failure("Matrix<Treal>::sy_inch: Matrix is zero! " 02683 "Not possible to compute inverse cholesky."); 02684 if (version == stable) /* STABILIZED */ 02685 throw Failure("Matrix<Treal>::sy_inch: Only unstable " 02686 "version of sy_inch implemented."); 02687 Treal myThresh = 0; 02688 if (threshold != 0) 02689 myThresh = threshold / (Z.ncols() * Z.nrows()); 02690 Z.resetRows(A.rows); 02691 Z.resetCols(A.cols); 02692 Z.allocate(); 02693 02694 if (looking == left) { /* LEFT-LOOKING INCH */ 02695 if (threshold != 0) 02696 throw Failure("Matrix<Treal, Telement>::syInch: " 02697 "Thresholding not implemented for left-looking inch."); 02698 /* Left and unstable */ 02699 Telement::syInch(A(0,0), Z(0,0), threshold, looking, version); 02700 Telement Ptmp;//, tmp; 02701 for (int i = 1; i < Z.ncols(); i++) { 02702 for (int j = 0; j < i; j++) { 02703 /* Z(k,i) is nonzero for k = 0, 1, ...,j - 2, j - 1, i */ 02704 /* and Z(i,i) = I (yes it is i ^) */ 02705 Ptmp = A(j,i); /* (Z(i,i) == I) */ 02706 for (int k = 0; k < j; k++) /* Ptmp = A(j,:) * Z(:,i) */ 02707 Telement::gemm(true, false, 1.0, /* SYMMETRY USED */ 02708 A(k,j), Z(k,i), 1.0, Ptmp); 02709 Telement::trmm('L', 'U', true, 1.0, Z(j,j), Ptmp); 02710 02711 for (int k = 0; k < j; k++) /* Z(:,i) -= Z(:,j) * Ptmp */ 02712 Telement::gemm(false, false, -1.0, 02713 Z(k,j), Ptmp, 1.0, Z(k,i)); 02714 /* Z(j,j) is triangular: */ 02715 Telement::trmm('L', 'U', false, -1.0, Z(j,j), Ptmp); 02716 Telement::add(1.0, Ptmp, Z(j,i)); 02717 } 02718 Ptmp = A(i,i); /* Z(i,i) == I */ 02719 for (int k = 0; k < i; k++) /* SYMMETRY USED */ 02720 Telement::gemm_upper_tr_only(true, false, 1.0, 02721 A(k,i), Z(k,i), 1.0, Ptmp); 02722 /* Z(i,i) == I !! */ 02723 /* Z(:,i) *= INCH(Ptmp) */ 02724 Telement::syInch(Ptmp, Z(i,i), threshold, looking, version); 02725 for (int k = 0; k < i; k++) { 02726 Telement::trmm('R', 'U', false, 1.0, Z(i,i), Z(k,i)); 02727 /* INCH(Ptmp) is upper triangular*/ 02728 } 02729 } 02730 } /* end if left-looking inch */ 02731 else { /* RIGHT-LOOKING INCH */ 02732 assert(looking == right); /* right and unstable */ 02733 Telement Ptmp; 02734 Treal newThresh = 0; 02735 if (myThresh != 0 && Z.ncols() > 1) 02736 newThresh = myThresh * 0.0001; 02737 else 02738 newThresh = myThresh; 02739 02740 for (int i = 0; i < Z.ncols(); i++) { 02741 /* Ptmp = A(i,:) * Z(:,i) */ 02742 Ptmp = A(i,i); /* Z(i,i) == I */ 02743 for (int k = 0; k < i; k++) 02744 Telement::gemm_upper_tr_only(true, false, 1.0, /* SYMMETRY USED */ 02745 A(k,i), Z(k,i), 1.0, Ptmp); 02746 02747 /* Z(:,i) *= INCH(Ptmp) */ 02748 Telement::syInch(Ptmp, Z(i,i), newThresh, looking, version); 02749 for (int k = 0; k < i; k++) { 02750 Telement::trmm('R', 'U', false, 1.0, Z(i,i), Z(k,i)); 02751 /* INCH(Ptmp) is upper triangular */ 02752 } 02753 02754 for (int j = i + 1; j < Z.ncols(); j++) { 02755 /* Compute Ptmp = Z(i,i)^T * A(i,:) * Z(:,j) */ 02756 /* Z(k,j) is nonzero for k = 0, 1, ...,i - 2, i - 1, j */ 02757 /* and Z(j,j) = I */ 02758 Ptmp = A(i,j); /* (Z(j,j) == I) */ 02759 for (int k = 0; k < i; k++) /* Ptmp = A(i,:) * Z(:,j) */ 02760 Telement::gemm(true, false, 1.0, /* SYMMETRY USED */ 02761 A(k,i), Z(k,j), 1.0, Ptmp); 02762 Telement::trmm('L', 'U', true, 1.0, Z(i,i), Ptmp); 02763 02764 for (int k = 0; k < i; k++) /* Z(:,j) -= Z(:,i) * Ptmp */ 02765 Telement::gemm(false, false, -1.0, 02766 Z(k,i), Ptmp, 1.0, Z(k,j)); 02767 /* Z(i,i) is triangular: */ 02768 Telement::trmm('L', 'U', false, -1.0, Z(i,i), Ptmp); 02769 Telement::add(1.0, Ptmp, Z(i,j)); 02770 } /* end for j */ 02771 02772 /* Truncation starts here */ 02773 if (threshold != 0) { 02774 for (int k = 0; k < i; k++) 02775 Z(k,i).frob_thresh(myThresh); 02776 } 02777 } /* end for i */ 02778 } /* end else right-looking inch */ 02779 } 02780 02781 template<class Treal, class Telement> 02782 void Matrix<Treal, Telement>:: 02783 gersgorin(Treal& lmin, Treal& lmax) const { 02784 assert(!this->is_empty()); 02785 if (this->nScalarsRows() == this->nScalarsCols()) { 02786 int size = this->nScalarsCols(); 02787 Treal* colsums = new Treal[size]; 02788 Treal* diag = new Treal[size]; 02789 for (int ind = 0; ind < size; ind++) 02790 colsums[ind] = 0; 02791 this->add_abs_col_sums(colsums); 02792 this->get_diagonal(diag); 02793 Treal tmp1 = colsums[0] - template_blas_fabs(diag[0]); 02794 Treal tmp2; 02795 lmin = diag[0] - tmp1; 02796 lmax = diag[0] + tmp1; 02797 for (int col = 1; col < size; col++) { 02798 tmp1 = colsums[col] - template_blas_fabs(diag[col]); 02799 tmp2 = diag[col] - tmp1; 02800 lmin = (tmp2 < lmin ? tmp2 : lmin); 02801 tmp2 = diag[col] + tmp1; 02802 lmax = (tmp2 > lmax ? tmp2 : lmax); 02803 } 02804 delete[] diag; 02805 delete[] colsums; 02806 } 02807 else 02808 throw Failure("Matrix<Treal, Telement, int>::gersgorin(Treal, Treal): " 02809 "Matrix must be quadratic"); 02810 } 02811 02812 02813 template<class Treal, class Telement> 02814 void Matrix<Treal, Telement>:: 02815 add_abs_col_sums(Treal* sums) const { 02816 assert(sums); 02817 if (!this->is_zero()) { 02818 int offset = 0; 02819 for (int col = 0; col < this->ncols(); col++) { 02820 for (int row = 0; row < this->nrows(); row++) { 02821 (*this)(row,col).add_abs_col_sums(&sums[offset]); 02822 } 02823 offset += (*this)(0,col).nScalarsCols(); 02824 } 02825 } 02826 } 02827 02828 template<class Treal, class Telement> 02829 void Matrix<Treal, Telement>:: 02830 get_diagonal(Treal* diag) const { 02831 assert(diag); 02832 assert(this->nrows() == this->ncols()); 02833 if (this->is_zero()) 02834 for (int rc = 0; rc < this->nScalarsCols(); rc++) 02835 diag[rc] = 0; 02836 else { 02837 int offset = 0; 02838 for (int rc = 0; rc < this->ncols(); rc++) { 02839 (*this)(rc,rc).get_diagonal(&diag[offset]); 02840 offset += (*this)(rc,rc).nScalarsCols(); 02841 } 02842 } 02843 } 02844 02845 template<class Treal, class Telement> 02846 size_t Matrix<Treal, Telement>::memory_usage() const { 02847 assert(!this->is_empty()); 02848 size_t sum = 0; 02849 if (this->is_zero()) 02850 return (size_t)0; 02851 for (int ind = 0; ind < this->nElements(); ind++) 02852 sum += this->elements[ind].memory_usage(); 02853 return sum; 02854 } 02855 02856 template<class Treal, class Telement> 02857 size_t Matrix<Treal, Telement>::nnz() const { 02858 size_t sum = 0; 02859 if (!this->is_zero()) { 02860 for (int ind = 0; ind < this->nElements(); ind++) 02861 sum += this->elements[ind].nnz(); 02862 } 02863 return sum; 02864 } 02865 template<class Treal, class Telement> 02866 size_t Matrix<Treal, Telement>::sy_nnz() const { 02867 size_t sum = 0; 02868 if (!this->is_zero()) { 02869 /* Above diagonal */ 02870 for (int col = 1; col < this->ncols(); col++) 02871 for (int row = 0; row < col; row++) 02872 sum += (*this)(row, col).nnz(); 02873 /* Below diagonal */ 02874 sum *= 2; 02875 /* Diagonal */ 02876 for (int rc = 0; rc < this->nrows(); rc++) 02877 sum += (*this)(rc,rc).sy_nnz(); 02878 } 02879 return sum; 02880 } 02881 02882 template<class Treal, class Telement> 02883 size_t Matrix<Treal, Telement>::sy_nvalues() const { 02884 size_t sum = 0; 02885 if (!this->is_zero()) { 02886 /* Above diagonal */ 02887 for (int col = 1; col < this->ncols(); col++) 02888 for (int row = 0; row < col; row++) 02889 sum += (*this)(row, col).nvalues(); 02890 /* Diagonal */ 02891 for (int rc = 0; rc < this->nrows(); rc++) 02892 sum += (*this)(rc,rc).sy_nvalues(); 02893 } 02894 return sum; 02895 } 02896 02897 02898 02899 02900 02901 /***************************************************************************/ 02902 /***************************************************************************/ 02903 /* Specialization for Telement = Treal */ 02904 /***************************************************************************/ 02905 /***************************************************************************/ 02906 02907 template<class Treal> 02908 class Matrix<Treal>: public MatrixHierarchicBase<Treal> { 02909 02910 public: 02911 typedef Vector<Treal, Treal> VectorType; 02912 friend class Vector<Treal, Treal>; 02913 02914 Matrix() 02915 :MatrixHierarchicBase<Treal>() { 02916 } 02917 /* Matrix(const Atomblock<Treal>& row_atoms, 02918 const Atomblock<Treal>& col_atoms, 02919 bool z = true, int nr = 0, int nc = 0, char tr = 'N') 02920 :MatrixHierarchicBase<Treal>(row_atoms, col_atoms, z, nr, nc,tr) {} 02921 */ 02922 02923 void allocate() { 02924 assert(!this->is_empty()); 02925 assert(this->is_zero()); 02926 this->elements = new Treal[this->nElements()]; 02927 for (int ind = 0; ind < this->nElements(); ++ind) 02928 this->elements[ind] = 0; 02929 } 02930 02931 /* Full matrix assigns etc */ 02932 void assignFromFull(std::vector<Treal> const & fullMat); 02933 void fullMatrix(std::vector<Treal> & fullMat) const; 02934 void syFullMatrix(std::vector<Treal> & fullMat) const; 02935 void syUpTriFullMatrix(std::vector<Treal> & fullMat) const; 02936 02937 /* Sparse matrix assigns etc */ 02938 void assignFromSparse(std::vector<int> const & rowind, 02939 std::vector<int> const & colind, 02940 std::vector<Treal> const & values); 02941 void assignFromSparse(std::vector<int> const & rowind, 02942 std::vector<int> const & colind, 02943 std::vector<Treal> const & values, 02944 std::vector<int> const & indexes); 02945 02946 /* Adds values (+=) to elements */ 02947 void addValues(std::vector<int> const & rowind, 02948 std::vector<int> const & colind, 02949 std::vector<Treal> const & values); 02950 void addValues(std::vector<int> const & rowind, 02951 std::vector<int> const & colind, 02952 std::vector<Treal> const & values, 02953 std::vector<int> const & indexes); 02954 02955 void syAssignFromSparse(std::vector<int> const & rowind, 02956 std::vector<int> const & colind, 02957 std::vector<Treal> const & values); 02958 02959 void syAddValues(std::vector<int> const & rowind, 02960 std::vector<int> const & colind, 02961 std::vector<Treal> const & values); 02962 02963 void getValues(std::vector<int> const & rowind, 02964 std::vector<int> const & colind, 02965 std::vector<Treal> & values) const; 02966 void getValues(std::vector<int> const & rowind, 02967 std::vector<int> const & colind, 02968 std::vector<Treal> & values, 02969 std::vector<int> const & indexes) const; 02970 void syGetValues(std::vector<int> const & rowind, 02971 std::vector<int> const & colind, 02972 std::vector<Treal> & values) const; 02973 02974 void getAllValues(std::vector<int> & rowind, 02975 std::vector<int> & colind, 02976 std::vector<Treal> & values) const; 02977 void syGetAllValues(std::vector<int> & rowind, 02978 std::vector<int> & colind, 02979 std::vector<Treal> & values) const; 02980 02981 Matrix<Treal>& 02982 operator=(const Matrix<Treal>& mat) { 02983 MatrixHierarchicBase<Treal>::operator=(mat); 02984 return *this; 02985 } 02986 02987 void clear(); 02988 ~Matrix() { 02989 clear(); 02990 } 02991 02992 void writeToFile(std::ofstream & file) const; 02993 void readFromFile(std::ifstream & file); 02994 02995 void random(); 02996 void syRandom(); 02997 void randomZeroStructure(Treal probabilityBeingZero); 02998 void syRandomZeroStructure(Treal probabilityBeingZero); 02999 template<typename TRule> 03000 void setElementsByRule(TRule & rule); 03001 template<typename TRule> 03002 void sySetElementsByRule(TRule & rule); 03003 03004 03005 void addIdentity(Treal alpha); /* C += alpha * I */ 03006 03007 static void transpose(Matrix<Treal> const & A, 03008 Matrix<Treal> & AT); 03009 03010 void symToNosym(); 03011 void nosymToSym(); 03012 03013 /* Set matrix to zero (k = 0) or identity (k = 1) */ 03014 Matrix<Treal>& operator=(int const k); 03015 03016 Matrix<Treal>& operator*=(const Treal alpha); 03017 03018 static void gemm(const bool tA, const bool tB, const Treal alpha, 03019 const Matrix<Treal>& A, 03020 const Matrix<Treal>& B, 03021 const Treal beta, 03022 Matrix<Treal>& C); 03023 static void symm(const char side, const char uplo, const Treal alpha, 03024 const Matrix<Treal>& A, 03025 const Matrix<Treal>& B, 03026 const Treal beta, 03027 Matrix<Treal>& C); 03028 static void syrk(const char uplo, const bool tA, const Treal alpha, 03029 const Matrix<Treal>& A, 03030 const Treal beta, 03031 Matrix<Treal>& C); 03032 /* C = beta * C + alpha * A * A where A and C are symmetric */ 03033 static void sysq(const char uplo, const Treal alpha, 03034 const Matrix<Treal>& A, 03035 const Treal beta, 03036 Matrix<Treal>& C); 03037 /* C = alpha * A * B + beta * C where A and B are symmetric */ 03038 static void ssmm(const Treal alpha, 03039 const Matrix<Treal>& A, 03040 const Matrix<Treal>& B, 03041 const Treal beta, 03042 Matrix<Treal>& C); 03043 /* C = alpha * A * B + beta * C where A and B are symmetric 03044 * and only the upper triangle of C is computed. 03045 */ 03046 static void ssmm_upper_tr_only(const Treal alpha, 03047 const Matrix<Treal>& A, 03048 const Matrix<Treal>& B, 03049 const Treal beta, 03050 Matrix<Treal>& C); 03051 03052 static void trmm(const char side, const char uplo, const bool tA, 03053 const Treal alpha, 03054 const Matrix<Treal>& A, 03055 Matrix<Treal>& B); 03056 03057 /* Frobenius norms */ 03058 inline Treal frob() const {return template_blas_sqrt(frobSquared());} 03059 Treal frobSquared() const; 03060 inline Treal syFrob() const { 03061 return template_blas_sqrt(this->syFrobSquared()); 03062 } 03063 Treal syFrobSquared() const; 03064 03065 inline static Treal frobDiff 03066 (const Matrix<Treal>& A, 03067 const Matrix<Treal>& B) { 03068 return template_blas_sqrt(frobSquaredDiff(A, B)); 03069 } 03070 static Treal frobSquaredDiff 03071 (const Matrix<Treal>& A, 03072 const Matrix<Treal>& B); 03073 03074 inline static Treal syFrobDiff 03075 (const Matrix<Treal>& A, 03076 const Matrix<Treal>& B) { 03077 return template_blas_sqrt(syFrobSquaredDiff(A, B)); 03078 } 03079 static Treal syFrobSquaredDiff 03080 (const Matrix<Treal>& A, 03081 const Matrix<Treal>& B); 03082 03083 Treal trace() const; 03084 static Treal trace_ab(const Matrix<Treal>& A, 03085 const Matrix<Treal>& B); 03086 static Treal trace_aTb(const Matrix<Treal>& A, 03087 const Matrix<Treal>& B); 03088 static Treal sy_trace_ab(const Matrix<Treal>& A, 03089 const Matrix<Treal>& B); 03090 03091 static void add(const Treal alpha, /* B += alpha * A */ 03092 const Matrix<Treal>& A, 03093 Matrix<Treal>& B); 03094 void assign(Treal const alpha, /* *this = alpha * A */ 03095 Matrix<Treal> const & A); 03096 03097 03098 /********** Help functions for thresholding */ 03099 // int getnnzBlocksLowestLevel() const; 03100 void getFrobSqLowestLevel(std::vector<Treal> & frobsq) const; 03101 void frobThreshLowestLevel 03102 (Treal const threshold, Matrix<Treal> * ErrorMatrix); 03103 03104 void getFrobSqElementLevel(std::vector<Treal> & frobsq) const; 03105 void frobThreshElementLevel 03106 (Treal const threshold, Matrix<Treal> * ErrorMatrix); 03107 03108 03109 #if 0 03110 inline void frobThreshLowestLevel 03111 (Treal const threshold, 03112 Matrix<Treal> * ErrorMatrix) { 03113 bool a,b; 03114 frobThreshLowestLevel(threshold, ErrorMatrix, a, b); 03115 } 03116 #endif 03117 03118 void assignFrobNormsLowestLevel 03119 ( Matrix<Treal, Matrix<Treal> > const & A ) { 03120 if (!A.is_zero()) { 03121 if ( this->is_zero() ) 03122 this->allocate(); 03123 assert( this->nElements() == A.nElements() ); 03124 for (int ind = 0; ind < this->nElements(); ind++) 03125 this->elements[ind] = A[ind].frob(); 03126 } 03127 else 03128 this->clear(); 03129 } 03130 03131 void syAssignFrobNormsLowestLevel( Matrix<Treal, Matrix<Treal> > const & A ) { 03132 if (!A.is_zero()) { 03133 if ( this->is_zero() ) 03134 this->allocate(); 03135 assert( this->nElements() == A.nElements() ); 03136 for (int col = 1; col < this->ncols(); col++) 03137 for (int row = 0; row < col; row++) 03138 (*this)(row,col) = A(row,col).frob(); 03139 for (int rc = 0; rc < this->nrows(); rc++) 03140 (*this)(rc,rc) = A(rc,rc).syFrob(); 03141 } 03142 else 03143 this->clear(); 03144 } 03145 03146 void assignDiffFrobNormsLowestLevel( Matrix<Treal, Matrix<Treal> > const & A, 03147 Matrix<Treal, Matrix<Treal> > const & B ) { 03148 if ( A.is_zero() && B.is_zero() ) { 03149 // Both A and B are zero 03150 this->clear(); 03151 return; 03152 } 03153 // At least one of A and B is nonzero 03154 if ( this->is_zero() ) 03155 this->allocate(); 03156 if ( !A.is_zero() && !B.is_zero() ) { 03157 // Both are nonzero 03158 assert( this->nElements() == A.nElements() ); 03159 assert( this->nElements() == B.nElements() ); 03160 for (int ind = 0; ind < this->nElements(); ind++) { 03161 Matrix<Treal> Diff(A[ind]); 03162 Matrix<Treal>::add( -1.0, B[ind], Diff ); 03163 this->elements[ind] = Diff.frob(); 03164 } 03165 return; 03166 } 03167 if ( !A.is_zero() ) { 03168 // A is nonzero 03169 this->assignFrobNormsLowestLevel( A ); 03170 return; 03171 } 03172 if ( !B.is_zero() ) { 03173 // B is nonzero 03174 this->assignFrobNormsLowestLevel( B ); 03175 return; 03176 } 03177 } 03178 void syAssignDiffFrobNormsLowestLevel( Matrix<Treal, Matrix<Treal> > const & A, 03179 Matrix<Treal, Matrix<Treal> > const & B ) { 03180 if ( A.is_zero() && B.is_zero() ) { 03181 // Both A and B are zero 03182 this->clear(); 03183 return; 03184 } 03185 // At least one of A and B is nonzero 03186 if ( this->is_zero() ) 03187 this->allocate(); 03188 if ( !A.is_zero() && !B.is_zero() ) { 03189 // Both are nonzero 03190 assert( this->nElements() == A.nElements() ); 03191 assert( this->nElements() == B.nElements() ); 03192 for (int col = 1; col < this->ncols(); col++) 03193 for (int row = 0; row < col; row++) { 03194 Matrix<Treal> Diff(A(row,col)); 03195 Matrix<Treal>::add( -1.0, B(row,col), Diff ); 03196 (*this)(row, col) = Diff.frob(); 03197 } 03198 for (int rc = 0; rc < this->ncols(); rc++) { 03199 Matrix<Treal> Diff( A(rc,rc) ); 03200 Matrix<Treal>::add( -1.0, B(rc,rc), Diff ); 03201 (*this)(rc, rc) = Diff.syFrob(); 03202 } 03203 return; 03204 } 03205 if ( !A.is_zero() ) { 03206 // A is nonzero 03207 this->syAssignFrobNormsLowestLevel( A ); 03208 return; 03209 } 03210 if ( !B.is_zero() ) { 03211 // B is nonzero 03212 this->syAssignFrobNormsLowestLevel( B ); 03213 return; 03214 } 03215 } 03216 03217 03218 void truncateAccordingToSparsityPattern( Matrix<Treal, Matrix<Treal> > & A ) const { 03219 if ( this->is_zero() ) 03220 A.clear(); 03221 else { 03222 assert( !A.is_zero() ); 03223 assert( this->nElements() == A.nElements() ); 03224 for (int ind = 0; ind < this->nElements(); ind++) 03225 if (this->elements[ind] == 0) 03226 A[ind].clear(); 03227 } 03228 } 03229 03230 03231 /********** End of help functions for thresholding */ 03232 03233 static void gemm_upper_tr_only(const bool tA, const bool tB, 03234 const Treal alpha, 03235 const Matrix<Treal>& A, 03236 const Matrix<Treal>& B, 03237 const Treal beta, 03238 Matrix<Treal>& C); 03239 static void sytr_upper_tr_only(char const side, const Treal alpha, 03240 Matrix<Treal>& A, 03241 const Matrix<Treal>& Z); 03242 static void trmm_upper_tr_only(const char side, const char uplo, 03243 const bool tA, const Treal alpha, 03244 const Matrix<Treal>& A, 03245 Matrix<Treal>& B); 03246 static void trsytriplemm(char const side, 03247 const Matrix<Treal>& Z, 03248 Matrix<Treal>& A); 03249 03250 inline Treal frob_thresh(Treal const threshold, 03251 Matrix<Treal> * ErrorMatrix = 0) { 03252 return template_blas_sqrt 03253 (frob_squared_thresh(threshold * threshold, ErrorMatrix)); 03254 } 03255 /* Returns the Frobenius norm of the introduced error */ 03256 03257 Treal frob_squared_thresh(Treal const threshold, 03258 Matrix<Treal> * ErrorMatrix = 0); 03259 03260 03261 static void inch(const Matrix<Treal>& A, 03262 Matrix<Treal>& Z, 03263 const Treal threshold = 0, 03264 const side looking = left, 03265 const inchversion version = unstable); 03266 static void syInch(const Matrix<Treal>& A, 03267 Matrix<Treal>& Z, 03268 const Treal threshold = 0, 03269 const side looking = left, 03270 const inchversion version = unstable) { 03271 inch(A, Z, threshold, looking, version); 03272 } 03273 03274 void gersgorin(Treal& lmin, Treal& lmax) const; 03275 void sy_gersgorin(Treal& lmin, Treal& lmax) const { 03276 Matrix<Treal> tmp = (*this); 03277 tmp.symToNosym(); 03278 tmp.gersgorin(lmin, lmax); 03279 return; 03280 } 03281 void add_abs_col_sums(Treal* abscolsums) const; 03282 void get_diagonal(Treal* diag) const; /* Copy diagonal to the diag array */ 03283 03284 inline size_t memory_usage() const { /* Returns memory usage in bytes */ 03285 assert(!this->is_empty()); 03286 if (this->is_zero()) 03287 return (size_t)0; 03288 else 03289 return (size_t)this->nElements() * sizeof(Treal); 03290 } 03291 03292 inline size_t nnz() const { 03293 if (this->is_zero()) 03294 return 0; 03295 else 03296 return this->nElements(); 03297 } 03298 inline size_t sy_nnz() const { 03299 if (this->is_zero()) 03300 return 0; 03301 else 03302 return this->nElements(); 03303 } 03307 inline size_t nvalues() const { 03308 return nnz(); 03309 } 03312 size_t sy_nvalues() const { 03313 assert(this->nScalarsRows() == this->nScalarsCols()); 03314 int n = this->nrows(); 03315 return this->is_zero() ? 0 : n * (n + 1) / 2; 03316 } 03321 template<class Top> 03322 Treal syAccumulateWith(Top & op) { 03323 Treal res = 0; 03324 if (!this->is_zero()) { 03325 int rowOffset = this->rows.getOffset(); 03326 int colOffset = this->cols.getOffset(); 03327 for (int col = 0; col < this->ncols(); col++) { 03328 for (int row = 0; row < col; row++) { 03329 res += 2*op.accumulate((*this)(row, col), 03330 rowOffset + row, 03331 colOffset + col); 03332 } 03333 res += op.accumulate((*this)(col, col), 03334 rowOffset + col, 03335 colOffset + col); 03336 } 03337 } 03338 return res; 03339 } 03340 template<class Top> 03341 Treal geAccumulateWith(Top & op) { 03342 Treal res = 0; 03343 if (!this->is_zero()) { 03344 int rowOffset = this->rows.getOffset(); 03345 int colOffset = this->cols.getOffset(); 03346 for (int col = 0; col < this->ncols(); col++) 03347 for (int row = 0; row < this->nrows(); row++) 03348 res += op.accumulate((*this)(row, col), 03349 rowOffset + row, 03350 colOffset + col); 03351 } 03352 return res; 03353 } 03354 03355 static inline unsigned int level() {return 0;} 03356 03357 Treal maxAbsValue() const { 03358 if (this->is_zero()) 03359 return 0; 03360 else { 03361 Treal maxAbsGlobal = 0; 03362 Treal maxAbsLocal = 0; 03363 for (int ind = 0; ind < this->nElements(); ++ind) { 03364 maxAbsLocal = template_blas_fabs(this->elements[ind]); 03365 maxAbsGlobal = maxAbsGlobal > maxAbsLocal ? 03366 maxAbsGlobal : maxAbsLocal; 03367 } /* end for */ 03368 return maxAbsGlobal; 03369 } 03370 } 03371 03372 /* New routines above */ 03373 03374 #if 0 /* OLD ROUTINES */ 03375 03376 03377 #if 0 03378 inline Matrix<Treal>& operator=(const Matrix<Treal>& mat) { 03379 this->MatrixHierarchicBase<Treal>::operator=(mat); 03380 std::cout<<"operator="<<std::endl; 03381 } 03382 #endif 03383 03384 03385 03386 03387 03388 03389 03390 03391 03392 03393 void trim_memory_usage(); 03394 #if 1 03395 void write_to_buffer_count(int& zb_length, int& vb_length) const; 03396 void write_to_buffer(int* zerobuf, const int zb_length, 03397 Treal* valuebuf, const int vb_length, 03398 int& zb_index, int& vb_index) const; 03399 void read_from_buffer(int* zerobuf, const int zb_length, 03400 Treal* valuebuf, const int vb_length, 03401 int& zb_index, int& vb_index); 03402 #endif 03403 03404 03405 03406 03407 03408 03409 /* continue here */ 03410 03411 03412 03413 03414 03415 03416 03417 03418 inline bool lowestLevel() const {return true;} 03419 // inline unsigned int level() const {return 0;} 03420 03421 #endif /* END OF OLD ROUTINES */ 03422 protected: 03423 private: 03424 static const Treal ZERO; 03425 static const Treal ONE; 03426 }; /* end class specialization Matrix<Treal> */ 03427 03428 template<class Treal> 03429 const Treal Matrix<Treal>::ZERO = 0; 03430 template<class Treal> 03431 const Treal Matrix<Treal>::ONE = 1; 03432 03433 #if 1 03434 /* Full matrix assigns etc */ 03435 template<class Treal> 03436 void Matrix<Treal>:: 03437 assignFromFull(std::vector<Treal> const & fullMat) { 03438 int nTotalRows = this->rows.getNTotalScalars(); 03439 int nTotalCols = this->cols.getNTotalScalars(); 03440 assert((int)fullMat.size() == nTotalRows * nTotalCols); 03441 int rowOffset = this->rows.getOffset(); 03442 int colOffset = this->cols.getOffset(); 03443 if (this->is_zero()) 03444 allocate(); 03445 for (int col = 0; col < this->ncols(); col++) 03446 for (int row = 0; row < this->nrows(); row++) 03447 (*this)(row, col) = 03448 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)]; 03449 } 03450 03451 template<class Treal> 03452 void Matrix<Treal>:: 03453 fullMatrix(std::vector<Treal> & fullMat) const { 03454 int nTotalRows = this->rows.getNTotalScalars(); 03455 int nTotalCols = this->cols.getNTotalScalars(); 03456 fullMat.resize(nTotalRows * nTotalCols); 03457 int rowOffset = this->rows.getOffset(); 03458 int colOffset = this->cols.getOffset(); 03459 if (this->is_zero()) { 03460 for (int col = 0; col < this->nScalarsCols(); col++) 03461 for (int row = 0; row < this->nScalarsRows(); row++) 03462 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0; 03463 } 03464 else { 03465 for (int col = 0; col < this->ncols(); col++) 03466 for (int row = 0; row < this->nrows(); row++) 03467 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 03468 (*this)(row, col); 03469 } 03470 } 03471 03472 template<class Treal> 03473 void Matrix<Treal>:: 03474 syFullMatrix(std::vector<Treal> & fullMat) const { 03475 int nTotalRows = this->rows.getNTotalScalars(); 03476 int nTotalCols = this->cols.getNTotalScalars(); 03477 fullMat.resize(nTotalRows * nTotalCols); 03478 int rowOffset = this->rows.getOffset(); 03479 int colOffset = this->cols.getOffset(); 03480 if (this->is_zero()) { 03481 for (int col = 0; col < this->nScalarsCols(); col++) 03482 for (int row = 0; row <= col; row++) { 03483 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0; 03484 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0; 03485 } 03486 } 03487 else { 03488 for (int col = 0; col < this->ncols(); col++) 03489 for (int row = 0; row <= col; row++) { 03490 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 03491 (*this)(row, col); 03492 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 03493 (*this)(row, col); 03494 } 03495 } 03496 } 03497 03498 template<class Treal> 03499 void Matrix<Treal>:: 03500 syUpTriFullMatrix(std::vector<Treal> & fullMat) const { 03501 int nTotalRows = this->rows.getNTotalScalars(); 03502 int nTotalCols = this->cols.getNTotalScalars(); 03503 fullMat.resize(nTotalRows * nTotalCols); 03504 int rowOffset = this->rows.getOffset(); 03505 int colOffset = this->cols.getOffset(); 03506 if (this->is_zero()) { 03507 for (int col = 0; col < this->nScalarsCols(); col++) 03508 for (int row = 0; row <= this->nScalarsRows(); row++) { 03509 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0; 03510 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0; 03511 } 03512 } 03513 else { 03514 for (int col = 0; col < this->ncols(); col++) 03515 for (int row = 0; row < this->nrows(); row++) { 03516 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 03517 (*this)(row, col); 03518 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 03519 (*this)(row, col); 03520 } 03521 } 03522 } 03523 03524 #endif 03525 03526 template<class Treal> 03527 void Matrix<Treal>:: 03528 assignFromSparse(std::vector<int> const & rowind, 03529 std::vector<int> const & colind, 03530 std::vector<Treal> const & values) { 03531 assert(rowind.size() == colind.size() && 03532 rowind.size() == values.size()); 03533 std::vector<int> indexes(values.size()); 03534 for (int ind = 0; ind < values.size(); ++ind) { 03535 indexes[ind] = ind; 03536 } 03537 assignFromSparse(rowind, colind, values, indexes); 03538 } 03539 03540 template<class Treal> 03541 void Matrix<Treal>:: 03542 assignFromSparse(std::vector<int> const & rowind, 03543 std::vector<int> const & colind, 03544 std::vector<Treal> const & values, 03545 std::vector<int> const & indexes) { 03546 if (indexes.empty()) { 03547 this->clear(); 03548 return; 03549 } 03550 if (this->is_zero()) 03551 allocate(); 03552 for (int ind = 0; ind < this->nElements(); ++ind) 03553 this->elements[ind] = 0; 03554 std::vector<int>::const_iterator it; 03555 for ( it = indexes.begin(); it < indexes.end(); it++ ) { 03556 (*this)(this->rows.whichBlock( rowind[*it] ), 03557 this->cols.whichBlock( colind[*it] )) = values[*it]; 03558 } 03559 } 03560 03561 03562 template<class Treal> 03563 void Matrix<Treal>:: 03564 addValues(std::vector<int> const & rowind, 03565 std::vector<int> const & colind, 03566 std::vector<Treal> const & values) { 03567 assert(rowind.size() == colind.size() && 03568 rowind.size() == values.size()); 03569 std::vector<int> indexes(values.size()); 03570 for (int ind = 0; ind < values.size(); ++ind) { 03571 indexes[ind] = ind; 03572 } 03573 addValues(rowind, colind, values, indexes); 03574 } 03575 03576 template<class Treal> 03577 void Matrix<Treal>:: 03578 addValues(std::vector<int> const & rowind, 03579 std::vector<int> const & colind, 03580 std::vector<Treal> const & values, 03581 std::vector<int> const & indexes) { 03582 if (indexes.empty()) 03583 return; 03584 if (this->is_zero()) 03585 allocate(); 03586 std::vector<int>::const_iterator it; 03587 for ( it = indexes.begin(); it < indexes.end(); it++ ) { 03588 (*this)(this->rows.whichBlock( rowind[*it] ), 03589 this->cols.whichBlock( colind[*it] )) += values[*it]; 03590 } 03591 } 03592 03593 template<class Treal> 03594 void Matrix<Treal>:: 03595 syAssignFromSparse(std::vector<int> const & rowind, 03596 std::vector<int> const & colind, 03597 std::vector<Treal> const & values) { 03598 assert(rowind.size() == colind.size() && 03599 rowind.size() == values.size()); 03600 bool upperTriangleOnly = true; 03601 for (int ind = 0; ind < values.size(); ++ind) { 03602 upperTriangleOnly = 03603 upperTriangleOnly && colind[ind] >= rowind[ind]; 03604 } 03605 if (!upperTriangleOnly) 03606 throw Failure("Matrix<Treal>::" 03607 "syAddValues(std::vector<int> const &, " 03608 "std::vector<int> const &, " 03609 "std::vector<Treal> const &, int const): " 03610 "Only upper triangle can contain elements when assigning " 03611 "symmetric or triangular matrix "); 03612 assignFromSparse(rowind, colind, values); 03613 } 03614 03615 template<class Treal> 03616 void Matrix<Treal>:: 03617 syAddValues(std::vector<int> const & rowind, 03618 std::vector<int> const & colind, 03619 std::vector<Treal> const & values) { 03620 assert(rowind.size() == colind.size() && 03621 rowind.size() == values.size()); 03622 bool upperTriangleOnly = true; 03623 for (int ind = 0; ind < values.size(); ++ind) { 03624 upperTriangleOnly = 03625 upperTriangleOnly && colind[ind] >= rowind[ind]; 03626 } 03627 if (!upperTriangleOnly) 03628 throw Failure("Matrix<Treal>::" 03629 "syAddValues(std::vector<int> const &, " 03630 "std::vector<int> const &, " 03631 "std::vector<Treal> const &, int const): " 03632 "Only upper triangle can contain elements when assigning " 03633 "symmetric or triangular matrix "); 03634 addValues(rowind, colind, values); 03635 } 03636 03637 template<class Treal> 03638 void Matrix<Treal>:: 03639 getValues(std::vector<int> const & rowind, 03640 std::vector<int> const & colind, 03641 std::vector<Treal> & values) const { 03642 assert(rowind.size() == colind.size()); 03643 values.resize(rowind.size(),0); 03644 std::vector<int> indexes(rowind.size()); 03645 for (int ind = 0; ind < rowind.size(); ++ind) { 03646 indexes[ind] = ind; 03647 } 03648 getValues(rowind, colind, values, indexes); 03649 } 03650 03651 template<class Treal> 03652 void Matrix<Treal>:: 03653 getValues(std::vector<int> const & rowind, 03654 std::vector<int> const & colind, 03655 std::vector<Treal> & values, 03656 std::vector<int> const & indexes) const { 03657 assert(!this->is_empty()); 03658 if (indexes.empty()) 03659 return; 03660 std::vector<int>::const_iterator it; 03661 if (this->is_zero()) { 03662 for ( it = indexes.begin(); it < indexes.end(); it++ ) 03663 values[*it] = 0; 03664 return; 03665 } 03666 for ( it = indexes.begin(); it < indexes.end(); it++ ) 03667 values[*it] = (*this)(this->rows.whichBlock( rowind[*it] ), 03668 this->cols.whichBlock( colind[*it] )); 03669 } 03670 03671 03672 template<class Treal> 03673 void Matrix<Treal>:: 03674 syGetValues(std::vector<int> const & rowind, 03675 std::vector<int> const & colind, 03676 std::vector<Treal> & values) const { 03677 assert(rowind.size() == colind.size()); 03678 bool upperTriangleOnly = true; 03679 for (int ind = 0; ind < rowind.size(); ++ind) { 03680 upperTriangleOnly = 03681 upperTriangleOnly && colind[ind] >= rowind[ind]; 03682 } 03683 if (!upperTriangleOnly) 03684 throw Failure("Matrix<Treal>::" 03685 "syGetValues(std::vector<int> const &, " 03686 "std::vector<int> const &, " 03687 "std::vector<Treal> const &, int const): " 03688 "Only upper triangle when retrieving elements from " 03689 "symmetric or triangular matrix "); 03690 getValues(rowind, colind, values); 03691 } 03692 03693 template<class Treal> 03694 void Matrix<Treal>:: 03695 getAllValues(std::vector<int> & rowind, 03696 std::vector<int> & colind, 03697 std::vector<Treal> & values) const { 03698 assert(rowind.size() == colind.size() && 03699 rowind.size() == values.size()); 03700 if (!this->is_zero()) 03701 for (int col = 0; col < this->ncols(); col++) 03702 for (int row = 0; row < this->nrows(); row++) { 03703 rowind.push_back(this->rows.getOffset() + row); 03704 colind.push_back(this->cols.getOffset() + col); 03705 values.push_back((*this)(row, col)); 03706 } 03707 } 03708 03709 template<class Treal> 03710 void Matrix<Treal>:: 03711 syGetAllValues(std::vector<int> & rowind, 03712 std::vector<int> & colind, 03713 std::vector<Treal> & values) const { 03714 assert(rowind.size() == colind.size() && 03715 rowind.size() == values.size()); 03716 if (!this->is_zero()) 03717 for (int col = 0; col < this->ncols(); ++col) 03718 for (int row = 0; row <= col; ++row) { 03719 rowind.push_back(this->rows.getOffset() + row); 03720 colind.push_back(this->cols.getOffset() + col); 03721 values.push_back((*this)(row, col)); 03722 } 03723 } 03724 03725 03726 template<class Treal> 03727 void Matrix<Treal>::clear() { 03728 delete[] this->elements; 03729 this->elements = 0; 03730 } 03731 03732 template<class Treal> 03733 void Matrix<Treal>:: 03734 writeToFile(std::ofstream & file) const { 03735 int const ZERO = 0; 03736 int const ONE = 1; 03737 if (this->is_zero()) { 03738 char * tmp = (char*)&ZERO; 03739 file.write(tmp,sizeof(int)); 03740 } 03741 else { 03742 char * tmp = (char*)&ONE; 03743 file.write(tmp,sizeof(int)); 03744 char * tmpel = (char*)this->elements; 03745 file.write(tmpel,sizeof(Treal) * this->nElements()); 03746 } 03747 } 03748 03749 template<class Treal> 03750 void Matrix<Treal>:: 03751 readFromFile(std::ifstream & file) { 03752 int const ZERO = 0; 03753 int const ONE = 1; 03754 char tmp[sizeof(int)]; 03755 file.read(tmp, (std::ifstream::pos_type)sizeof(int)); 03756 switch ((int)*tmp) { 03757 case ZERO: 03758 this->clear(); 03759 break; 03760 case ONE: 03761 if (this->is_zero()) 03762 allocate(); 03763 file.read((char*)this->elements, sizeof(Treal) * this->nElements()); 03764 break; 03765 default: 03766 throw Failure("Matrix<Treal>::" 03767 "readFromFile(std::ifstream & file):" 03768 "File corruption, int value not 0 or 1"); 03769 } 03770 } 03771 03772 template<class Treal> 03773 void Matrix<Treal>::random() { 03774 if (this->is_zero()) 03775 allocate(); 03776 for (int ind = 0; ind < this->nElements(); ind++) 03777 this->elements[ind] = rand() / (Treal)RAND_MAX; 03778 } 03779 03780 template<class Treal> 03781 void Matrix<Treal>::syRandom() { 03782 if (this->is_zero()) 03783 allocate(); 03784 /* Above diagonal */ 03785 for (int col = 1; col < this->ncols(); col++) 03786 for (int row = 0; row < col; row++) 03787 (*this)(row, col) = rand() / (Treal)RAND_MAX;; 03788 /* Diagonal */ 03789 for (int rc = 0; rc < this->nrows(); rc++) 03790 (*this)(rc,rc) = rand() / (Treal)RAND_MAX;; 03791 } 03792 03793 template<class Treal> 03794 void Matrix<Treal>:: 03795 randomZeroStructure(Treal probabilityBeingZero) { 03796 if (!this->highestLevel() && 03797 probabilityBeingZero > rand() / (Treal)RAND_MAX) 03798 this->clear(); 03799 else 03800 this->random(); 03801 } 03802 03803 template<class Treal> 03804 void Matrix<Treal>:: 03805 syRandomZeroStructure(Treal probabilityBeingZero) { 03806 if (!this->highestLevel() && 03807 probabilityBeingZero > rand() / (Treal)RAND_MAX) 03808 this->clear(); 03809 else 03810 this->syRandom(); 03811 } 03812 03813 template<class Treal> 03814 template<typename TRule> 03815 void Matrix<Treal>:: 03816 setElementsByRule(TRule & rule) { 03817 if (this->is_zero()) 03818 allocate(); 03819 for (int col = 0; col < this->ncols(); col++) 03820 for (int row = 0; row < this->nrows(); row++) 03821 (*this)(row,col) = rule.set(this->rows.getOffset() + row, 03822 this->cols.getOffset() + col); 03823 } 03824 template<class Treal> 03825 template<typename TRule> 03826 void Matrix<Treal>:: 03827 sySetElementsByRule(TRule & rule) { 03828 if (this->is_zero()) 03829 allocate(); 03830 /* Upper triangle */ 03831 for (int col = 0; col < this->ncols(); col++) 03832 for (int row = 0; row <= col; row++) 03833 (*this)(row,col) = rule.set(this->rows.getOffset() + row, 03834 this->cols.getOffset() + col); 03835 } 03836 03837 03838 template<class Treal> 03839 void Matrix<Treal>:: 03840 addIdentity(Treal alpha) { 03841 if (this->is_empty()) 03842 throw Failure("Matrix<Treal>::addIdentity(Treal): " 03843 "Cannot add identity to empty matrix."); 03844 if (this->ncols() != this->nrows()) 03845 throw Failure("Matrix<Treal, Telement>::addIdentity(Treal): " 03846 "Matrix must be square to add identity"); 03847 if (this->is_zero()) 03848 allocate(); 03849 for (int ind = 0; ind < this->ncols(); ind++) 03850 (*this)(ind,ind) += alpha; 03851 } 03852 03853 template<class Treal> 03854 void Matrix<Treal>:: 03855 transpose(Matrix<Treal> const & A, Matrix<Treal> & AT) { 03856 if (A.is_zero()) { /* Condition also matches empty matrices. */ 03857 AT.rows = A.cols; 03858 AT.cols = A.rows; 03859 delete[] AT.elements; 03860 AT.elements = 0; 03861 return; 03862 } 03863 if (AT.is_zero() || (AT.nElements() != A.nElements())) { 03864 delete[] AT.elements; 03865 AT.elements = new Treal[A.nElements()]; 03866 } 03867 AT.cols = A.rows; 03868 AT.rows = A.cols; 03869 for (int row = 0; row < AT.nrows(); row++) 03870 for (int col = 0; col < AT.ncols(); col++) 03871 AT(row,col) = A(col,row); 03872 } 03873 03874 template<class Treal> 03875 void Matrix<Treal>:: 03876 symToNosym() { 03877 if (this->nrows() == this->ncols()) { 03878 if (!this->is_zero()) { 03879 /* Diagonal should be fine */ 03880 /* Fix the lower triangle */ 03881 for (int row = 1; row < this->nrows(); row++) 03882 for (int col = 0; col < row; col++) 03883 (*this)(row, col) = (*this)(col, row); 03884 } 03885 } 03886 else 03887 throw Failure("Matrix<Treal>::symToNosym(): " 03888 "Only quadratic matrices can be symmetric"); 03889 } 03890 03891 template<class Treal> 03892 void Matrix<Treal>:: 03893 nosymToSym() { 03894 if (this->nrows() == this->ncols()) { 03895 if (!this->is_zero()) { 03896 /* Diagonal should be fine */ 03897 /* Fix the lower triangle */ 03898 for (int col = 0; col < this->ncols() - 1; col++) 03899 for (int row = col + 1; row < this->nrows(); row++) 03900 (*this)(row, col) = 0; 03901 } 03902 } 03903 else 03904 throw Failure("Matrix<Treal>::nosymToSym(): " 03905 "Only quadratic matrices can be symmetric"); 03906 } 03907 03908 template<class Treal> 03909 Matrix<Treal>& 03910 Matrix<Treal>::operator=(int const k) { 03911 switch (k) { 03912 case 0: 03913 this->clear(); 03914 break; 03915 case 1: 03916 if (this->ncols() != this->nrows()) 03917 throw Failure("Matrix<Treal>::operator=(int k = 1): " 03918 "Matrix must be quadratic to " 03919 "become identity matrix."); 03920 this->clear(); 03921 this->allocate(); 03922 for (int rc = 0; rc < this->ncols(); rc++) /*Set diagonal to identity*/ 03923 (*this)(rc,rc) = 1; 03924 break; 03925 default: 03926 throw Failure 03927 ("Matrix<Treal>::operator=(int k) only implemented for k = 0, k = 1"); 03928 } 03929 return *this; 03930 } 03931 03932 template<class Treal> 03933 Matrix<Treal>& Matrix<Treal>:: 03934 operator*=(const Treal alpha) { 03935 if (!this->is_zero() && alpha != 1) { 03936 for (int ind = 0; ind < this->nElements(); ind++) 03937 this->elements[ind] *= alpha; 03938 } 03939 return *this; 03940 } 03941 03942 template<class Treal> 03943 void Matrix<Treal>:: 03944 gemm(const bool tA, const bool tB, const Treal alpha, 03945 const Matrix<Treal>& A, 03946 const Matrix<Treal>& B, 03947 const Treal beta, 03948 Matrix<Treal>& C) { 03949 assert(!A.is_empty()); 03950 assert(!B.is_empty()); 03951 if (C.is_empty()) { 03952 assert(beta == 0); 03953 if (tA) 03954 C.resetRows(A.cols); 03955 else 03956 C.resetRows(A.rows); 03957 if (tB) 03958 C.resetCols(B.rows); 03959 else 03960 C.resetCols(B.cols); 03961 } 03962 03963 if ((A.is_zero() || B.is_zero() || alpha == 0) && 03964 (C.is_zero() || beta == 0)) 03965 C = 0; 03966 else { 03967 Treal beta_tmp = beta; 03968 if (C.is_zero()) { 03969 C.allocate(); 03970 beta_tmp = 0; 03971 } 03972 03973 if (!A.is_zero() && !B.is_zero() && alpha != 0) { 03974 if (!tA && !tB) 03975 if (A.ncols() == B.nrows() && 03976 A.nrows() == C.nrows() && 03977 B.ncols() == C.ncols()) 03978 mat::gemm("N", "N", &A.nrows(), &B.ncols(), &A.ncols(), &alpha, 03979 A.elements, &A.nrows(), B.elements, &B.nrows(), 03980 &beta_tmp, C.elements, &C.nrows()); 03981 else 03982 throw Failure("Matrix<Treal>::" 03983 "gemm(bool, bool, Treal, Matrix<Treal>, " 03984 "Matrix<Treal>, Treal, " 03985 "Matrix<Treal>): " 03986 "Incorrect matrixdimensions for multiplication"); 03987 else if (tA && !tB) 03988 if (A.nrows() == B.nrows() && 03989 A.ncols() == C.nrows() && 03990 B.ncols() == C.ncols()) 03991 mat::gemm("T", "N", &A.ncols(), &B.ncols(), &A.nrows(), &alpha, 03992 A.elements, &A.nrows(), B.elements, &B.nrows(), 03993 &beta_tmp, C.elements, &C.nrows()); 03994 else 03995 throw Failure("Matrix<Treal>::" 03996 "gemm(bool, bool, Treal, Matrix<Treal>, " 03997 "Matrix<Treal>, Treal, " 03998 "Matrix<Treal>): " 03999 "Incorrect matrixdimensions for multiplication"); 04000 else if (!tA && tB) 04001 if (A.ncols() == B.ncols() && 04002 A.nrows() == C.nrows() && 04003 B.nrows() == C.ncols()) 04004 mat::gemm("N", "T", &A.nrows(), &B.nrows(), &A.ncols(), &alpha, 04005 A.elements, &A.nrows(), B.elements, &B.nrows(), 04006 &beta_tmp, C.elements, &C.nrows()); 04007 else 04008 throw Failure("Matrix<Treal>::" 04009 "gemm(bool, bool, Treal, Matrix<Treal>, " 04010 "Matrix<Treal>, Treal, " 04011 "Matrix<Treal>): " 04012 "Incorrect matrixdimensions for multiplication"); 04013 else if (tA && tB) 04014 if (A.nrows() == B.ncols() && 04015 A.ncols() == C.nrows() && 04016 B.nrows() == C.ncols()) 04017 mat::gemm("T", "T", &A.ncols(), &B.nrows(), &A.nrows(), &alpha, 04018 A.elements, &A.nrows(), B.elements, &B.nrows(), 04019 &beta_tmp, C.elements, &C.nrows()); 04020 else 04021 throw Failure("Matrix<Treal>::" 04022 "gemm(bool, bool, Treal, Matrix<Treal>, " 04023 "Matrix<Treal>, Treal, " 04024 "Matrix<Treal>): " 04025 "Incorrect matrixdimensions for multiplication"); 04026 else throw Failure("Matrix<Treal>::" 04027 "gemm(bool, bool, Treal, Matrix<Treal>, " 04028 "Matrix<Treal>, Treal, " 04029 "Matrix<Treal>):Very strange error!!"); 04030 } 04031 else 04032 C *= beta; 04033 } 04034 } 04035 04036 04037 template<class Treal> 04038 void Matrix<Treal>:: 04039 symm(const char side, const char uplo, const Treal alpha, 04040 const Matrix<Treal>& A, 04041 const Matrix<Treal>& B, 04042 const Treal beta, 04043 Matrix<Treal>& C) { 04044 assert(!A.is_empty()); 04045 assert(!B.is_empty()); 04046 assert(A.nrows() == A.ncols()); 04047 // int dimA = A.nrows(); 04048 if (C.is_empty()) { 04049 assert(beta == 0); 04050 if (side =='L') { 04051 C.resetRows(A.rows); 04052 C.resetCols(B.cols); 04053 } 04054 else { 04055 assert(side == 'R'); 04056 C.resetRows(B.rows); 04057 C.resetCols(A.cols); 04058 } 04059 } 04060 04061 if ((A.is_zero() || B.is_zero() || alpha == 0) && 04062 (C.is_zero() || beta == 0)) 04063 C = 0; 04064 else { 04065 Treal beta_tmp = beta; 04066 if (C.is_zero()) { 04067 C.allocate(); 04068 beta_tmp = 0; 04069 } 04070 if (!A.is_zero() && !B.is_zero() && alpha != 0) { 04071 if (A.nrows() == A.ncols() && C.nrows() == B.nrows() && 04072 C.ncols() == B.ncols() && 04073 ((side == 'L' && A.ncols() == B.nrows()) || 04074 (side == 'R' && B.ncols() == A.nrows()))) 04075 mat::symm(&side, &uplo, &C.nrows(), &C.ncols(), &alpha, 04076 A.elements, &A.nrows(), B.elements, &B.nrows(), 04077 &beta_tmp, C.elements, &C.nrows()); 04078 else 04079 throw Failure("Matrix<Treal>::symm: Incorrect matrix " 04080 "dimensions for symmetric multiply"); 04081 } /* end if (!A.is_zero() && !B.is_zero() && alpha != 0) */ 04082 else 04083 C *= beta; 04084 } 04085 } 04086 04087 template<class Treal> 04088 void Matrix<Treal>:: 04089 syrk(const char uplo, const bool tA, const Treal alpha, 04090 const Matrix<Treal>& A, 04091 const Treal beta, 04092 Matrix<Treal>& C) { 04093 assert(!A.is_empty()); 04094 if (C.is_empty()) { 04095 assert(beta == 0); 04096 if (tA) { 04097 C.resetRows(A.cols); 04098 C.resetCols(A.cols); 04099 } 04100 else { 04101 C.resetRows(A.rows); 04102 C.resetCols(A.rows); 04103 } 04104 } 04105 if (C.nrows() == C.ncols() && 04106 ((!tA && A.nrows() == C.nrows()) || (tA && A.ncols() == C.nrows()))) 04107 if (alpha != 0 && !A.is_zero()) { 04108 Treal beta_tmp = beta; 04109 if (C.is_zero()) { 04110 C.allocate(); 04111 beta_tmp = 0; 04112 } 04113 if (!tA) { 04114 mat::syrk(&uplo, "N", &C.nrows(), &A.ncols(), &alpha, 04115 A.elements, &A.nrows(), &beta_tmp, 04116 C.elements, &C.nrows()); 04117 } /* end if (!tA) */ 04118 else 04119 mat::syrk(&uplo, "T", &C.nrows(), &A.nrows(), &alpha, 04120 A.elements, &A.nrows(), &beta_tmp, 04121 C.elements, &C.nrows()); 04122 } 04123 else 04124 C *= beta; 04125 else 04126 throw Failure("Matrix<Treal>::syrk: Incorrect matrix " 04127 "dimensions for symmetric rank-k update"); 04128 } 04129 04130 template<class Treal> 04131 void Matrix<Treal>:: 04132 sysq(const char uplo, const Treal alpha, 04133 const Matrix<Treal>& A, 04134 const Treal beta, 04135 Matrix<Treal>& C) { 04136 assert(!A.is_empty()); 04137 if (C.is_empty()) { 04138 assert(beta == 0); 04139 C.resetRows(A.rows); 04140 C.resetCols(A.cols); 04141 } 04142 /* FIXME: slow copy */ 04143 Matrix<Treal> tmpA = A; 04144 tmpA.symToNosym(); 04145 Matrix<Treal>::syrk(uplo, false, alpha, tmpA, beta, C); 04146 } 04147 04148 /* C = alpha * A * B + beta * C where A and B are symmetric */ 04149 template<class Treal> 04150 void Matrix<Treal>:: 04151 ssmm(const Treal alpha, 04152 const Matrix<Treal>& A, 04153 const Matrix<Treal>& B, 04154 const Treal beta, 04155 Matrix<Treal>& C) { 04156 assert(!A.is_empty()); 04157 assert(!B.is_empty()); 04158 if (C.is_empty()) { 04159 assert(beta == 0); 04160 C.resetRows(A.rows); 04161 C.resetCols(B.cols); 04162 } 04163 /* FIXME: slow copy */ 04164 Matrix<Treal> tmpB = B; 04165 tmpB.symToNosym(); 04166 Matrix<Treal>::symm('L', 'U', alpha, A, tmpB, beta, C); 04167 } 04168 04169 /* C = alpha * A * B + beta * C where A and B are symmetric 04170 * and only the upper triangle of C is computed. 04171 */ 04172 template<class Treal> 04173 void Matrix<Treal>:: 04174 ssmm_upper_tr_only(const Treal alpha, 04175 const Matrix<Treal>& A, 04176 const Matrix<Treal>& B, 04177 const Treal beta, 04178 Matrix<Treal>& C) { 04179 /* FIXME: Symmetry is not utilized. */ 04180 ssmm(alpha, A, B, beta, C); 04181 C.nosymToSym(); 04182 } 04183 04184 04185 template<class Treal> 04186 void Matrix<Treal>:: 04187 trmm(const char side, const char uplo, const bool tA, 04188 const Treal alpha, 04189 const Matrix<Treal>& A, 04190 Matrix<Treal>& B) { 04191 assert(!A.is_empty()); 04192 assert(!B.is_empty()); 04193 if (alpha != 0 && !A.is_zero() && !B.is_zero()) 04194 if (((side == 'R' && B.ncols() == A.nrows()) || 04195 (side == 'L' && A.ncols() == B.nrows())) && 04196 A.nrows() == A.ncols()) 04197 if (!tA) 04198 mat::trmm(&side, &uplo, "N", "N", &B.nrows(), &B.ncols(), 04199 &alpha, A.elements, &A.nrows(), B.elements, &B.nrows()); 04200 else 04201 mat::trmm(&side, &uplo, "T", "N", &B.nrows(), &B.ncols(), 04202 &alpha, A.elements, &A.nrows(), B.elements, &B.nrows()); 04203 else 04204 throw Failure("Matrix<Treal>::trmm: " 04205 "Incorrect matrix dimensions for multiplication"); 04206 else 04207 B = 0; 04208 } 04209 04210 template<class Treal> 04211 Treal Matrix<Treal>::frobSquared() const { 04212 assert(!this->is_empty()); 04213 if (this->is_zero()) 04214 return 0; 04215 else { 04216 Treal sum(0); 04217 for (int i = 0; i < this->nElements(); i++) 04218 sum += this->elements[i] * this->elements[i]; 04219 return sum; 04220 } 04221 } 04222 04223 template<class Treal> 04224 Treal Matrix<Treal>:: 04225 syFrobSquared() const { 04226 assert(!this->is_empty()); 04227 if (this->nrows() != this->ncols()) 04228 throw Failure("Matrix<Treal>::syFrobSquared(): " 04229 "Matrix must be have equal number of rows " 04230 "and cols to be symmetric"); 04231 Treal sum(0); 04232 if (!this->is_zero()) { 04233 for (int col = 1; col < this->ncols(); col++) 04234 for (int row = 0; row < col; row++) 04235 sum += 2 * (*this)(row, col) * (*this)(row, col); 04236 for (int rc = 0; rc < this->ncols(); rc++) 04237 sum += (*this)(rc, rc) * (*this)(rc, rc); 04238 } 04239 return sum; 04240 } 04241 04242 template<class Treal> 04243 Treal Matrix<Treal>:: 04244 frobSquaredDiff 04245 (const Matrix<Treal>& A, 04246 const Matrix<Treal>& B) { 04247 assert(!A.is_empty()); 04248 assert(!B.is_empty()); 04249 if (A.nrows() != B.nrows() || A.ncols() != B.ncols()) 04250 throw Failure("Matrix<Treal>::frobSquaredDiff: " 04251 "Incorrect matrix dimensions"); 04252 Treal sum(0); 04253 if (!A.is_zero() && !B.is_zero()) { 04254 Treal diff; 04255 for (int i = 0; i < A.nElements(); i++) { 04256 diff = A.elements[i] - B.elements[i]; 04257 sum += diff * diff; 04258 } 04259 } 04260 else if (A.is_zero() && !B.is_zero()) 04261 sum = B.frobSquared(); 04262 else if (!A.is_zero() && B.is_zero()) 04263 sum = A.frobSquared(); 04264 /* sum is already zero if A.is_zero() && B.is_zero() */ 04265 return sum; 04266 } 04267 04268 04269 template<class Treal> 04270 Treal Matrix<Treal>:: 04271 syFrobSquaredDiff 04272 (const Matrix<Treal>& A, 04273 const Matrix<Treal>& B) { 04274 assert(!A.is_empty()); 04275 assert(!B.is_empty()); 04276 if (A.nrows() != B.nrows() || 04277 A.ncols() != B.ncols() || 04278 A.nrows() != A.ncols()) 04279 throw Failure("Matrix<Treal>::syFrobSquaredDiff: " 04280 "Incorrect matrix dimensions"); 04281 Treal sum(0); 04282 if (!A.is_zero() && !B.is_zero()) { 04283 Treal diff; 04284 for (int col = 1; col < A.ncols(); col++) 04285 for (int row = 0; row < col; row++) { 04286 diff = A(row, col) - B(row, col); 04287 sum += 2 * diff * diff; 04288 } 04289 for (int rc = 0; rc < A.ncols(); rc++) { 04290 diff = A(rc, rc) - B(rc, rc); 04291 sum += diff * diff; 04292 } 04293 } 04294 else if (A.is_zero() && !B.is_zero()) 04295 sum = B.syFrobSquared(); 04296 else if (!A.is_zero() && B.is_zero()) 04297 sum = A.syFrobSquared(); 04298 /* sum is already zero if A.is_zero() && B.is_zero() */ 04299 return sum; 04300 } 04301 04302 template<class Treal> 04303 Treal Matrix<Treal>:: 04304 trace() const { 04305 assert(!this->is_empty()); 04306 if (this->nrows() != this->ncols()) 04307 throw Failure("Matrix<Treal>::trace(): Matrix must be quadratic"); 04308 if (this->is_zero()) 04309 return 0; 04310 else { 04311 Treal sum = 0; 04312 for (int rc = 0; rc < this->ncols(); rc++) 04313 sum += (*this)(rc,rc); 04314 return sum; 04315 } 04316 } 04317 04318 template<class Treal> 04319 Treal Matrix<Treal>:: 04320 trace_ab(const Matrix<Treal>& A, 04321 const Matrix<Treal>& B) { 04322 assert(!A.is_empty()); 04323 assert(!B.is_empty()); 04324 if (A.nrows() != B.ncols() || A.ncols() != B.nrows()) 04325 throw Failure("Matrix<Treal>::trace_ab: " 04326 "Wrong matrix dimensions for trace(A * B)"); 04327 Treal tr = 0; 04328 if (!A.is_zero() && !B.is_zero()) 04329 for (int rc = 0; rc < A.nrows(); rc++) 04330 for (int colA = 0; colA < A.ncols(); colA++) 04331 tr += A(rc,colA) * B(colA, rc); 04332 return tr; 04333 } 04334 04335 template<class Treal> 04336 Treal Matrix<Treal>:: 04337 trace_aTb(const Matrix<Treal>& A, 04338 const Matrix<Treal>& B) { 04339 assert(!A.is_empty()); 04340 assert(!B.is_empty()); 04341 if (A.ncols() != B.ncols() || A.nrows() != B.nrows()) 04342 throw Failure("Matrix<Treal>::trace_aTb: " 04343 "Wrong matrix dimensions for trace(A' * B)"); 04344 Treal tr = 0; 04345 if (!A.is_zero() && !B.is_zero()) 04346 for (int rc = 0; rc < A.ncols(); rc++) 04347 for (int rowB = 0; rowB < B.nrows(); rowB++) 04348 tr += A(rowB,rc) * B(rowB, rc); 04349 return tr; 04350 } 04351 04352 template<class Treal> 04353 Treal Matrix<Treal>:: 04354 sy_trace_ab(const Matrix<Treal>& A, 04355 const Matrix<Treal>& B) { 04356 assert(!A.is_empty()); 04357 assert(!B.is_empty()); 04358 if (A.nrows() != B.ncols() || A.ncols() != B.nrows() || 04359 A.nrows() != A.ncols()) 04360 throw Failure("Matrix<Treal>::sy_trace_ab: " 04361 "Wrong matrix dimensions for symmetric trace(A * B)"); 04362 if (A.is_zero() || B.is_zero()) 04363 return 0; 04364 /* Now we know both A and B are nonzero */ 04365 Treal tr = 0; 04366 /* Diagonal first */ 04367 for (int rc = 0; rc < A.nrows(); rc++) 04368 tr += A(rc,rc) * B(rc, rc); 04369 /* Using that trace of transpose is equal to that without transpose: */ 04370 for (int colA = 1; colA < A.ncols(); colA++) 04371 for (int rowA = 0; rowA < colA; rowA++) 04372 tr += 2 * A(rowA, colA) * B(rowA, colA); 04373 return tr; 04374 } 04375 04376 template<class Treal> 04377 void Matrix<Treal>:: 04378 add(const Treal alpha, /* B += alpha * A */ 04379 const Matrix<Treal>& A, 04380 Matrix<Treal>& B) { 04381 assert(!A.is_empty()); 04382 assert(!B.is_empty()); 04383 if (A.nrows() != B.nrows() || A.ncols() != B.ncols()) 04384 throw Failure("Matrix<Treal>::add: " 04385 "Wrong matrix dimensions for addition"); 04386 if (!A.is_zero() && alpha != 0) { 04387 if (B.is_zero()) 04388 B.allocate(); 04389 for (int ind = 0; ind < A.nElements(); ind++) 04390 B.elements[ind] += alpha * A.elements[ind]; 04391 } 04392 } 04393 04394 template<class Treal> 04395 void Matrix<Treal>:: 04396 assign(Treal const alpha, /* *this = alpha * A */ 04397 Matrix<Treal> const & A) { 04398 assert(!A.is_empty()); 04399 if (this->is_empty()) { 04400 this->resetRows(A.rows); 04401 this->resetCols(A.cols); 04402 } 04403 Matrix<Treal>::add(alpha, A, *this); 04404 } 04405 04406 04407 /********** Help functions for thresholding */ 04408 04409 template<class Treal> 04410 void Matrix<Treal>:: 04411 getFrobSqLowestLevel(std::vector<Treal> & frobsq) const { 04412 if (!this->is_zero()) 04413 frobsq.push_back(this->frobSquared()); 04414 } 04415 04416 template<class Treal> 04417 void Matrix<Treal>:: 04418 getFrobSqElementLevel(std::vector<Treal> & frobsq) const { 04419 if (!this->is_zero()) 04420 for (int ind = 0; ind < this->nElements(); ind++) 04421 if ( this->elements[ind] != 0 ) // Add nonzero elements only 04422 frobsq.push_back(this->elements[ind] * this->elements[ind]); 04423 } 04424 04425 04426 template<class Treal> 04427 void Matrix<Treal>:: 04428 frobThreshLowestLevel 04429 (Treal const threshold, Matrix<Treal> * ErrorMatrix) { 04430 if (ErrorMatrix) { 04431 if ((!this->is_zero() && this->frobSquared() <= threshold) || 04432 (!ErrorMatrix->is_zero() && 04433 ErrorMatrix->frobSquared() > threshold)) 04434 Matrix<Treal>::swap(*this,*ErrorMatrix); 04435 } 04436 else if (!this->is_zero() && this->frobSquared() <= threshold) 04437 this->clear(); 04438 } 04439 04440 template<class Treal> 04441 void Matrix<Treal>:: 04442 frobThreshElementLevel 04443 (Treal const threshold, Matrix<Treal> * ErrorMatrix) { 04444 assert(!this->is_empty()); 04445 bool thisMatIsZero = true; 04446 if (ErrorMatrix) { 04447 assert(!ErrorMatrix->is_empty()); 04448 bool EMatIsZero = true; 04449 if (!ErrorMatrix->is_zero() || !this->is_zero()) { 04450 if (ErrorMatrix->is_zero()) 04451 ErrorMatrix->allocate(); 04452 if (this->is_zero()) 04453 this->allocate(); 04454 for (int ind = 0; ind < this->nElements(); ind++) { 04455 if ( this->elements[ind] != 0 ) { 04456 assert(ErrorMatrix->elements[ind] == 0); 04457 // ok, let's check if we want to move the element to the error matrix 04458 if ( this->elements[ind] * this->elements[ind] <= threshold ) { 04459 // we want to move the element 04460 ErrorMatrix->elements[ind] = this->elements[ind]; 04461 this->elements[ind] = 0; 04462 EMatIsZero = false; // at least one element is nonzero 04463 } 04464 else 04465 thisMatIsZero = false; // at least one element is nonzero 04466 continue; 04467 } 04468 if ( ErrorMatrix->elements[ind] != 0 ) { 04469 assert(this->elements[ind] == 0); 04470 // ok, let's check if we want to move the element from the error matrix 04471 if ( ErrorMatrix->elements[ind] * ErrorMatrix->elements[ind] > threshold ) { 04472 // we want to move the element 04473 this->elements[ind] = ErrorMatrix->elements[ind]; 04474 ErrorMatrix->elements[ind] = 0; 04475 thisMatIsZero = false; // at least one element is nonzero 04476 } 04477 else 04478 EMatIsZero = false; // at least one element is nonzero 04479 } 04480 } 04481 if (thisMatIsZero) { 04482 #if 0 04483 for (int ind = 0; ind < this->nElements(); ind++) 04484 assert( this->elements[ind] == 0); 04485 #endif 04486 this->clear(); 04487 } 04488 if (EMatIsZero) { 04489 #if 0 04490 for (int ind = 0; ind < this->nElements(); ind++) 04491 assert( ErrorMatrix->elements[ind] == 0); 04492 #endif 04493 ErrorMatrix->clear(); 04494 } 04495 } 04496 } 04497 else 04498 if (!this->is_zero()) { 04499 for (int ind = 0; ind < this->nElements(); ind++) { 04500 if ( this->elements[ind] * this->elements[ind] <= threshold ) 04501 /* FIXME BUG? EMANUEL LOOK AT THIS! */ 04502 // this->elements[ind] == 0; OLD CODE. Compiler complained that "statement has no effect". 04503 this->elements[ind] = 0; /* New code. Changed from == to =. */ 04504 else 04505 thisMatIsZero = false; 04506 } 04507 if (thisMatIsZero) 04508 this->clear(); 04509 } 04510 } 04511 04512 04513 04514 /********** End of help functions for thresholding */ 04515 04516 /* C = beta * C + alpha * A * B where only the upper triangle of C is */ 04517 /* referenced and updated */ 04518 template<class Treal> 04519 void Matrix<Treal>:: 04520 gemm_upper_tr_only(const bool tA, const bool tB, 04521 const Treal alpha, 04522 const Matrix<Treal>& A, 04523 const Matrix<Treal>& B, 04524 const Treal beta, 04525 Matrix<Treal>& C) { 04526 /* FIXME: Symmetry is not utilized. */ 04527 Matrix<Treal>::gemm(tA, tB, alpha, A, B, beta, C); 04528 C.nosymToSym(); 04529 } 04530 04531 /* A = alpha * A * Z or A = alpha * Z * A where A is symmetric, */ 04532 /* Z is upper triangular and */ 04533 /* only the upper triangle of the result is calculated */ 04534 /* side == 'R' for A = alpha * A * Z */ 04535 /* side == 'L' for A = alpha * Z * A */ 04536 template<class Treal> 04537 void Matrix<Treal>:: 04538 sytr_upper_tr_only(char const side, const Treal alpha, 04539 Matrix<Treal>& A, 04540 const Matrix<Treal>& Z) { 04541 /* FIXME: Symmetry is not utilized. */ 04542 A.symToNosym(); 04543 Matrix<Treal>::trmm(side, 'U', false, alpha, Z, A); 04544 A.nosymToSym(); 04545 } 04546 04547 /* The result B is assumed to be symmetric, i.e. only upper triangle */ 04548 /* calculated and hence only upper triangle of input B is used. */ 04549 template<class Treal> 04550 void Matrix<Treal>:: 04551 trmm_upper_tr_only(const char side, const char uplo, 04552 const bool tA, const Treal alpha, 04553 const Matrix<Treal>& A, 04554 Matrix<Treal>& B) { 04555 /* FIXME: Symmetry is not utilized. */ 04556 assert(tA); 04557 B.symToNosym(); 04558 Matrix<Treal>::trmm(side, uplo, tA, alpha, A, B); 04559 B.nosymToSym(); 04560 } 04561 04562 /* A = Z' * A * Z or A = Z * A * Z' */ 04563 /* where A is symmetric and Z is (nonunit) upper triangular */ 04564 /* side == 'R' for A = Z' * A * Z */ 04565 /* side == 'L' for A = Z * A * Z' */ 04566 template<class Treal> 04567 void Matrix<Treal>:: 04568 trsytriplemm(char const side, 04569 const Matrix<Treal>& Z, 04570 Matrix<Treal>& A) { 04571 if (side == 'R') { 04572 Matrix<Treal>:: 04573 sytr_upper_tr_only('R', 1.0, A, Z); 04574 Matrix<Treal>:: 04575 trmm_upper_tr_only('L', 'U', true, 1.0, Z, A); 04576 } 04577 else { 04578 assert(side == 'L'); 04579 Matrix<Treal>:: 04580 sytr_upper_tr_only('L', 1.0, A, Z); 04581 Matrix<Treal>:: 04582 trmm_upper_tr_only('R', 'U', true, 1.0, Z, A); 04583 } 04584 } 04585 04586 04587 template<class Treal> 04588 Treal Matrix<Treal>::frob_squared_thresh 04589 (Treal const threshold, Matrix<Treal> * ErrorMatrix) { 04590 assert(!this->is_empty()); 04591 if (ErrorMatrix && ErrorMatrix->is_empty()) { 04592 ErrorMatrix->resetRows(this->rows); 04593 ErrorMatrix->resetCols(this->cols); 04594 } 04595 Treal fs = frobSquared(); 04596 if (fs < threshold) { 04597 if (ErrorMatrix) 04598 Matrix<Treal>::swap(*this, *ErrorMatrix); 04599 return fs; 04600 } 04601 else 04602 return 0; 04603 } 04604 04605 04606 template<class Treal> 04607 void Matrix<Treal>:: 04608 inch(const Matrix<Treal>& A, 04609 Matrix<Treal>& Z, 04610 const Treal threshold, const side looking, 04611 const inchversion version) { 04612 assert(!A.is_empty()); 04613 if (A.nrows() != A.ncols()) 04614 throw Failure("Matrix<Treal>::inch: Matrix must be quadratic!"); 04615 if (A.is_zero()) 04616 throw Failure("Matrix<Treal>::inch: Matrix is zero! " 04617 "Not possible to compute inverse cholesky."); 04618 Z = A; 04619 int info; 04620 potrf("U", &A.nrows(), Z.elements, &A.nrows(), &info); 04621 if (info > 0) 04622 throw Failure("Matrix<Treal>::inch: potrf returned info > 0. The matrix is not positive definite."); 04623 if (info < 0) 04624 throw Failure("Matrix<Treal>::inch: potrf returned info < 0"); 04625 04626 trtri("U", "N", &A.nrows(), Z.elements, &A.nrows(), &info); 04627 if (info > 0) 04628 throw Failure("Matrix<Treal>::inch: trtri returned info > 0. The matrix is not positive definite."); 04629 if (info < 0) 04630 throw Failure("Matrix<Treal>::inch: trtri returned info < 0"); 04631 /* Fill lower triangle with zeroes just to be safe */ 04632 trifulltofull(Z.elements, A.nrows()); 04633 } 04634 04635 template<class Treal> 04636 void Matrix<Treal>:: 04637 gersgorin(Treal& lmin, Treal& lmax) const { 04638 assert(!this->is_empty()); 04639 if (this->nScalarsRows() == this->nScalarsCols()) { 04640 int size = this->nScalarsCols(); 04641 Treal* colsums = new Treal[size]; 04642 Treal* diag = new Treal[size]; 04643 for (int ind = 0; ind < size; ind++) 04644 colsums[ind] = 0; 04645 this->add_abs_col_sums(colsums); 04646 this->get_diagonal(diag); 04647 Treal tmp1 = colsums[0] - template_blas_fabs(diag[0]); 04648 Treal tmp2; 04649 lmin = diag[0] - tmp1; 04650 lmax = diag[0] + tmp1; 04651 for (int col = 1; col < size; col++) { 04652 tmp1 = colsums[col] - template_blas_fabs(diag[col]); 04653 tmp2 = diag[col] - tmp1; 04654 lmin = (tmp2 < lmin ? tmp2 : lmin); 04655 tmp2 = diag[col] + tmp1; 04656 lmax = (tmp2 > lmax ? tmp2 : lmax); 04657 } 04658 delete[] diag; 04659 delete[] colsums; 04660 } 04661 else 04662 throw Failure("Matrix<Treal>::gersgorin(Treal, Treal): Matrix must be quadratic"); 04663 } 04664 04665 04666 template<class Treal> 04667 void Matrix<Treal>:: 04668 add_abs_col_sums(Treal* sums) const { 04669 assert(sums); 04670 if (!this->is_zero()) 04671 for (int col = 0; col < this->ncols(); col++) 04672 for (int row = 0; row < this->nrows(); row++) 04673 sums[col] += template_blas_fabs((*this)(row,col)); 04674 } 04675 04676 template<class Treal> 04677 void Matrix<Treal>:: 04678 get_diagonal(Treal* diag) const { 04679 assert(diag); 04680 assert(this->nrows() == this->ncols()); 04681 if (this->is_zero()) 04682 for (int rc = 0; rc < this->nScalarsCols(); rc++) 04683 diag[rc] = 0; 04684 else 04685 for (int rc = 0; rc < this->ncols(); rc++) 04686 diag[rc] = (*this)(rc,rc); 04687 } 04688 04689 04690 /* New routines above */ 04691 04692 #if 0 /* OLD ROUTINES */ 04693 04694 04695 04696 04697 04698 04699 04700 04701 04702 04703 04704 04705 04706 04707 04708 04709 04710 04711 04712 04713 template<class Treal> 04714 void Matrix<Treal>::trim_memory_usage() { 04715 if (this->is_zero() && this->cap > 0) { 04716 delete[] this->elements; 04717 this->elements = NULL; 04718 this->cap = 0; 04719 this->nel = 0; 04720 } 04721 else if (this->cap > this->nel) { 04722 Treal* tmp = new Treal[this->nel]; 04723 for (int i = 0; i < this->nel; i++) 04724 tmp[i] = this->elements[i]; 04725 delete[] this->elements; 04726 this->cap = this->nel; 04727 this->elements = tmp; 04728 } 04729 assert(this->cap == this->nel); 04730 } 04731 04732 04733 04734 #if 1 04735 04736 template<class Treal> 04737 void Matrix<Treal>:: 04738 write_to_buffer_count(int& zb_length, int& vb_length) const { 04739 if (this->is_zero()) { 04740 ++zb_length; 04741 } 04742 else { 04743 ++zb_length; 04744 vb_length += this->nel; 04745 } 04746 } 04747 04748 template<class Treal> 04749 void Matrix<Treal>:: 04750 write_to_buffer(int* zerobuf, const int zb_length, 04751 Treal* valuebuf, const int vb_length, 04752 int& zb_index, int& vb_index) const { 04753 if (zb_index < zb_length) { 04754 if (this->is_zero()) { 04755 zerobuf[zb_index] = 0; 04756 ++zb_index; 04757 } 04758 else { 04759 if (vb_index + this->nel < vb_length + 1) { 04760 zerobuf[zb_index] = 1; 04761 ++zb_index; 04762 for (int i = 0; i < this->nel; i++) 04763 valuebuf[vb_index + i] = this->elements[i]; 04764 vb_index += this->nel; 04765 } 04766 else 04767 throw Failure("Matrix<Treal, Telement>::write_to_buffer: " 04768 "Insufficient space in buffers"); 04769 } 04770 } 04771 else 04772 throw Failure("Matrix<Treal, Telement>::write_to_buffer: " 04773 "Insufficient space in buffers"); 04774 } 04775 04776 template<class Treal> 04777 void Matrix<Treal>:: 04778 read_from_buffer(int* zerobuf, const int zb_length, 04779 Treal* valuebuf, const int vb_length, 04780 int& zb_index, int& vb_index) { 04781 if (zb_index < zb_length) { 04782 if (zerobuf[zb_index] == 0) { 04783 (*this) = 0; 04784 ++zb_index; 04785 } 04786 else { 04787 this->content = ful; 04788 this->nel = this->nrows() * this->ncols(); 04789 this->assert_alloc(); 04790 if (vb_index + this->nel < vb_length + 1) { 04791 assert(zerobuf[zb_index] == 1); 04792 ++zb_index; 04793 for (int i = 0; i < this->nel; i++) 04794 this->elements[i] = valuebuf[vb_index + i]; 04795 vb_index += this->nel; 04796 } 04797 else 04798 throw Failure("Matrix<Treal, Telement>::read_from_buffer: " 04799 "Mismatch, buffers does not match matrix"); 04800 } 04801 } 04802 else 04803 throw Failure("Matrix<Treal, Telement>::read_from_buffer: " 04804 "Mismatch, buffers does not match matrix"); 04805 } 04806 04807 #endif 04808 04809 04810 04811 04812 04813 04814 04815 04816 /* continue here */ 04817 04818 04819 04820 04821 04822 04823 04824 04825 04826 04827 04828 04829 04830 04831 #if 1 04832 04833 04834 04835 #endif 04836 04837 #endif /* END OF OLD ROUTINES */ 04838 04839 04840 } /* end namespace mat */ 04841 04842 #endif