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