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