ergo
Matrix.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
58 #ifndef MAT_MATRIX
59 #define MAT_MATRIX
60 
61 #include <math.h>
62 #include <cstdlib>
63 #include <algorithm>
64 
65 #include "MatrixHierarchicBase.h"
66 #include "matrix_proxy.h"
67 #include "gblas.h"
68 #include "sort.h"
69 #include "allocate.h"
70 //#define MAT_NAIVE
71 
72 namespace mat{
73  enum side {left, right};
75  template<class Treal, class Telement>
76  class Vector;
77  template<typename Tperm>
78  struct AccessMap;
79 
81  private:
83  public:
84  void reset() { accumulatedTime = 0; }
85  double getAccumulatedTime() { return accumulatedTime; }
86  void addTime(double timeToAdd) {
87 #ifdef _OPENMP
88  #pragma omp critical
89 #endif
90  {
91  accumulatedTime += timeToAdd;
92  }
93  }
95  static SingletonForTimings theInstance;
96  return theInstance;
97  }
98  private:
101  };
102 
103 
112  template<class Treal, class Telement = Treal>
113  class Matrix: public MatrixHierarchicBase<Treal, Telement> {
114  public:
115  typedef Telement ElementType;
117 
118  friend class Vector<Treal, Telement>;
119  Matrix():MatrixHierarchicBase<Treal, Telement>(){}
120 
121 
122  void allocate() {
123  assert(!this->is_empty());
124  assert(this->is_zero());
125  //#define MAT_USE_ALLOC_TIMER
126 #ifdef MAT_USE_ALLOC_TIMER
127  mat::Time theTimer; theTimer.tic();
128 #endif
129  this->elements = allocateElements<Telement>(this->nElements());
130 #ifdef MAT_USE_ALLOC_TIMER
132 #endif
133  SizesAndBlocks colSAB;
134  SizesAndBlocks rowSAB;
135  for (int col = 0; col < this->cols.getNBlocks(); col++) {
136  colSAB = this->cols.getSizesAndBlocksForLowerLevel(col);
137  for (int row = 0; row < this->rows.getNBlocks(); row++) {
138  /* This could be improved - precompute rowSAB as well as colSAB */
139  rowSAB = this->rows.getSizesAndBlocksForLowerLevel(row);
140  (*this)(row,col).resetRows(rowSAB);
141  (*this)(row,col).resetCols(colSAB);
142  }
143  }
144  }
145 
146  /* Full matrix assigns etc */
147  void assignFromFull(std::vector<Treal> const & fullMat);
148  void fullMatrix(std::vector<Treal> & fullMat) const;
149  void syFullMatrix(std::vector<Treal> & fullMat) const;
150  void syUpTriFullMatrix(std::vector<Treal> & fullMat) const;
151 
152  /* Sparse matrix assigns etc */
153  void assignFromSparse(std::vector<int> const & rowind,
154  std::vector<int> const & colind,
155  std::vector<Treal> const & values);
156  void assignFromSparse(std::vector<int> const & rowind,
157  std::vector<int> const & colind,
158  std::vector<Treal> const & values,
159  std::vector<int> const & indexes);
160  /* Adds values (+=) to elements */
161  void addValues(std::vector<int> const & rowind,
162  std::vector<int> const & colind,
163  std::vector<Treal> const & values);
164  void addValues(std::vector<int> const & rowind,
165  std::vector<int> const & colind,
166  std::vector<Treal> const & values,
167  std::vector<int> const & indexes);
168 
169  void syAssignFromSparse(std::vector<int> const & rowind,
170  std::vector<int> const & colind,
171  std::vector<Treal> const & values);
172 
173  void syAddValues(std::vector<int> const & rowind,
174  std::vector<int> const & colind,
175  std::vector<Treal> const & values);
176 
177  void getValues(std::vector<int> const & rowind,
178  std::vector<int> const & colind,
179  std::vector<Treal> & values) const;
180  void getValues(std::vector<int> const & rowind,
181  std::vector<int> const & colind,
182  std::vector<Treal> &,
183  std::vector<int> const & indexes) const;
184  void syGetValues(std::vector<int> const & rowind,
185  std::vector<int> const & colind,
186  std::vector<Treal> & values) const;
187  void getAllValues(std::vector<int> & rowind,
188  std::vector<int> & colind,
189  std::vector<Treal> &) const;
190  void syGetAllValues(std::vector<int> & rowind,
191  std::vector<int> & colind,
192  std::vector<Treal> &) const;
193 
194 
198  return *this;
199  }
200 
201 
202  void clear();
204  clear();
205  }
206  void writeToFile(std::ofstream & file) const;
207  void readFromFile(std::ifstream & file);
208 
209  void random();
210  void syRandom();
214  void randomZeroStructure(Treal probabilityBeingZero);
215  void syRandomZeroStructure(Treal probabilityBeingZero);
216 
217  template<typename TRule>
218  void setElementsByRule(TRule & rule);
219  template<typename TRule>
220  void sySetElementsByRule(TRule & rule);
221  template<typename TRule>
222  void trSetElementsByRule(TRule & rule) {
223  // Setting elements for triangular is the same as for symmetric
224  sySetElementsByRule(rule);
225  }
226 
227  void addIdentity(Treal alpha); /* C += alpha * I */
228 
229  static void transpose(Matrix<Treal, Telement> const & A,
231 
232  void symToNosym();
233  void nosymToSym();
234 
235  /* Basic linear algebra routines */
236 
237  /* Set matrix to zero (k = 0) or identity (k = 1) */
238  Matrix<Treal, Telement>& operator=(int const k);
239 
240  Matrix<Treal, Telement>& operator*=(const Treal alpha);
241 
242  static void gemm(const bool tA, const bool tB, const Treal alpha,
243  const Matrix<Treal, Telement>& A,
244  const Matrix<Treal, Telement>& B,
245  const Treal beta,
247  static void symm(const char side, const char uplo, const Treal alpha,
248  const Matrix<Treal, Telement>& A,
249  const Matrix<Treal, Telement>& B,
250  const Treal beta,
252  static void syrk(const char uplo, const bool tA, const Treal alpha,
253  const Matrix<Treal, Telement>& A,
254  const Treal beta,
256  /* C = alpha * A * A + beta * C where A and C are symmetric */
257  static void sysq(const char uplo, const Treal alpha,
258  const Matrix<Treal, Telement>& A,
259  const Treal beta,
261  /* C = alpha * A * B + beta * C where A and B are symmetric */
262  static void ssmm(const Treal alpha,
263  const Matrix<Treal, Telement>& A,
264  const Matrix<Treal, Telement>& B,
265  const Treal beta,
267  /* C = alpha * A * B + beta * C where A and B are symmetric
268  * and only the upper triangle of C is computed.
269  */
270  static void ssmm_upper_tr_only(const Treal alpha,
271  const Matrix<Treal, Telement>& A,
272  const Matrix<Treal, Telement>& B,
273  const Treal beta,
275 
276  static void trmm(const char side, const char uplo, const bool tA,
277  const Treal alpha,
278  const Matrix<Treal, Telement>& A,
280 
281  /* Frobenius norms */
282 
283  /* Returns the Frobenius norm of the matrix. */
284  inline Treal frob() const {
285  return template_blas_sqrt(this->frobSquared());
286  }
287  /* Returns the squared Frobenius norm */
288  Treal frobSquared() const;
289  inline Treal syFrob() const {
290  return template_blas_sqrt(this->syFrobSquared());
291  }
292  Treal syFrobSquared() const;
293 
294  inline static Treal frobDiff
296  const Matrix<Treal, Telement>& B) {
297  return template_blas_sqrt(frobSquaredDiff(A, B));
298  }
299  static Treal frobSquaredDiff
300  (const Matrix<Treal, Telement>& A,
301  const Matrix<Treal, Telement>& B);
302 
303  inline static Treal syFrobDiff
305  const Matrix<Treal, Telement>& B) {
306  return template_blas_sqrt(syFrobSquaredDiff(A, B));
307  }
308  static Treal syFrobSquaredDiff
309  (const Matrix<Treal, Telement>& A,
310  const Matrix<Treal, Telement>& B);
311 
312  /* Traces */
313  Treal trace() const;
314  static Treal trace_ab(const Matrix<Treal, Telement>& A,
315  const Matrix<Treal, Telement>& B);
316  static Treal trace_aTb(const Matrix<Treal, Telement>& A,
317  const Matrix<Treal, Telement>& B);
318  static Treal sy_trace_ab(const Matrix<Treal, Telement>& A,
319  const Matrix<Treal, Telement>& B);
320 
321  static void add(const Treal alpha, /* B += alpha * A */
322  const Matrix<Treal, Telement>& A,
324  void assign(Treal const alpha, /* *this = alpha * A */
325  Matrix<Treal, Telement> const & A);
326 
327 
328  /********** Help functions for thresholding */
329  // int getnnzBlocksLowestLevel() const;
330  void getFrobSqLowestLevel(std::vector<Treal> & frobsq) const;
332  (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix);
333 
334  void getFrobSqElementLevel(std::vector<Treal> & frobsq) const;
336  (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix);
337 
338 
339 #if 0
340  inline void frobThreshLowestLevel
341  (Treal const threshold,
342  Matrix<Treal, Telement> * ErrorMatrix) {
343  bool a,b;
344  frobThreshLowestLevel(threshold, ErrorMatrix, a, b);
345  }
346 #endif
347 
351  ( Matrix<Treal, Matrix<Treal, Telement> > const & A );
354  ( Matrix<Treal, Matrix<Treal, Telement> > const & A );
355 
360  ( Matrix<Treal, Matrix<Treal, Telement> > const & A,
361  Matrix<Treal, Matrix<Treal, Telement> > const & B );
366  ( Matrix<Treal, Matrix<Treal, Telement> > const & A,
367  Matrix<Treal, Matrix<Treal, Telement> > const & B );
368 
371  void truncateAccordingToSparsityPattern( Matrix<Treal, Matrix<Treal, Telement> > & A ) const;
372 
373 
374  /********** End of help functions for thresholding */
375 
376  static void gemm_upper_tr_only(const bool tA, const bool tB,
377  const Treal alpha,
378  const Matrix<Treal, Telement>& A,
379  const Matrix<Treal, Telement>& B,
380  const Treal beta,
381  Matrix<Treal, Telement>& C);
382  static void sytr_upper_tr_only(char const side, const Treal alpha,
383  Matrix<Treal, Telement>& A,
384  const Matrix<Treal, Telement>& Z);
385  static void trmm_upper_tr_only(const char side, const char uplo,
386  const bool tA, const Treal alpha,
387  const Matrix<Treal, Telement>& A,
388  Matrix<Treal, Telement>& B);
389  static void trsytriplemm(char const side,
390  const Matrix<Treal, Telement>& Z,
391  Matrix<Treal, Telement>& A);
392 
393  inline Treal frob_thresh
394  (Treal const threshold,
395  Matrix<Treal, Telement> * ErrorMatrix = 0) {
396  return template_blas_sqrt(frob_squared_thresh(threshold * threshold,
397  ErrorMatrix));
398  }
404  Treal frob_squared_thresh
405  (Treal const threshold,
406  Matrix<Treal, Telement> * ErrorMatrix = 0);
412  static void syInch(const Matrix<Treal, Telement>& A,
414  const Treal threshold = 0,
415  const side looking = left,
416  const inchversion version = unstable);
417 
418  void gersgorin(Treal& lmin, Treal& lmax) const; /* Computes bounds for*/
419  /* real part of eigenvalues. The matrix must be quadratic (of course) */
420  void sy_gersgorin(Treal& lmin, Treal& lmax) const {
421  Matrix<Treal, Telement> tmp = (*this);
422  tmp.symToNosym();
423  tmp.gersgorin(lmin, lmax);
424  return;
425  }
426 
427  void add_abs_col_sums(Treal* abscolsums) const; /* Adds the absolute */
428  /* column sums to the abscolsums array. */
429  /* abscolsums(i) += sum(abs(C(:,i))) for all i (C == *this) */
430  /* Used by e.g. gersgorin eigenvalue inclusion */
431  void get_diagonal(Treal* diag) const; /*Copy diagonal to the diag array*/
432 
433  size_t memory_usage() const; /* Returns memory usage in bytes */
434 
435  /* Note: we use size_t instead of int for nnz and nvalues to avoid
436  integer overflow. */
437  size_t nnz() const;
438  size_t sy_nnz() const;
442  inline size_t nvalues() const {
443  return nnz();
444  }
447  size_t sy_nvalues() const;
454  template<typename Top>
455  Treal syAccumulateWith(Top & op) {
456  Treal res = 0;
457  if (!this->is_zero()) {
458  for (int col = 0; col < this->ncols(); col++) {
459  for (int row = 0; row < col; row++)
460  res += 2 * (*this)(row, col).geAccumulateWith(op);
461  res += (*this)(col, col).syAccumulateWith(op);
462  }
463  }
464  return res;
465  }
467  template<typename Top>
468  Treal geAccumulateWith(Top & op) {
469  Treal res = 0;
470  if (!this->is_zero()) {
471  for (int col = 0; col < this->ncols(); col++)
472  for (int row = 0; row < this->nrows(); row++)
473  res += (*this)(row, col).geAccumulateWith(op);
474  }
475  return res;
476  }
477 
478  static inline unsigned int level() {return Telement::level() + 1;}
479 
480  Treal maxAbsValue() const {
481  if (this->is_zero())
482  return 0;
483  else {
484  Treal maxAbsGlobal = 0;
485  Treal maxAbsLocal = 0;
486  for (int ind = 0; ind < this->nElements(); ++ind) {
487  maxAbsLocal = this->elements[ind].maxAbsValue();
488  maxAbsGlobal = maxAbsGlobal > maxAbsLocal ?
489  maxAbsGlobal : maxAbsLocal;
490  } /* end for */
491  return maxAbsGlobal;
492  }
493  }
494 
495  protected:
496  private:
497  }; /* end class */
498 
499 
500 #if 1
501  /* Full matrix assigns etc */
502  template<class Treal, class Telement>
504  assignFromFull(std::vector<Treal> const & fullMat) {
505  int nTotalRows = this->rows.getNTotalScalars();
506  int nTotalCols = this->cols.getNTotalScalars();
507  assert((int)fullMat.size() == nTotalRows * nTotalCols);
508  if (this->is_zero())
509  allocate();
510  for (int col = 0; col < this->ncols(); col++)
511  for (int row = 0; row < this->nrows(); row++)
512  (*this)(row, col).assignFromFull(fullMat);
513  }
514 
515  template<class Treal, class Telement>
517  fullMatrix(std::vector<Treal> & fullMat) const {
518  int nTotalRows = this->rows.getNTotalScalars();
519  int nTotalCols = this->cols.getNTotalScalars();
520  fullMat.resize(nTotalRows * nTotalCols);
521  if (this->is_zero()) {
522  int rowOffset = this->rows.getOffset();
523  int colOffset = this->cols.getOffset();
524  for (int col = 0; col < this->nScalarsCols(); col++)
525  for (int row = 0; row < this->nScalarsRows(); row++)
526  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
527  }
528  else {
529  for (int col = 0; col < this->ncols(); col++)
530  for (int row = 0; row < this->nrows(); row++)
531  (*this)(row, col).fullMatrix(fullMat);
532  }
533  }
534 
535  template<class Treal, class Telement>
537  syFullMatrix(std::vector<Treal> & fullMat) const {
538  int nTotalRows = this->rows.getNTotalScalars();
539  int nTotalCols = this->cols.getNTotalScalars();
540  fullMat.resize(nTotalRows * nTotalCols);
541  if (this->is_zero()) {
542  int rowOffset = this->rows.getOffset();
543  int colOffset = this->cols.getOffset();
544  for (int col = 0; col < this->nScalarsCols(); col++)
545  for (int row = 0; row <= col; row++) {
546  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
547  fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
548  }
549  }
550  else {
551  for (int col = 0; col < this->ncols(); col++) {
552  for (int row = 0; row < col; row++)
553  (*this)(row, col).syUpTriFullMatrix(fullMat);
554  (*this)(col, col).syFullMatrix(fullMat);
555  }
556  }
557  }
558 
559  template<class Treal, class Telement>
561  syUpTriFullMatrix(std::vector<Treal> & fullMat) const {
562  int nTotalRows = this->rows.getNTotalScalars();
563  int nTotalCols = this->cols.getNTotalScalars();
564  fullMat.resize(nTotalRows * nTotalCols);
565  if (this->is_zero()) {
566  int rowOffset = this->rows.getOffset();
567  int colOffset = this->cols.getOffset();
568  for (int col = 0; col < this->nScalarsCols(); col++)
569  for (int row = 0; row <= this->nScalarsRows(); row++) {
570  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
571  fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
572  }
573  }
574  else {
575  for (int col = 0; col < this->ncols(); col++)
576  for (int row = 0; row < this->nrows(); row++)
577  (*this)(row, col).syUpTriFullMatrix(fullMat);
578  }
579  }
580 
581 #endif
582 
583 
584  template<class Treal, class Telement>
586  assignFromSparse(std::vector<int> const & rowind,
587  std::vector<int> const & colind,
588  std::vector<Treal> const & values) {
589  assert(rowind.size() == colind.size() &&
590  rowind.size() == values.size());
591  std::vector<int> indexes(values.size());
592  for (unsigned int ind = 0; ind < values.size(); ++ind)
593  indexes[ind] = ind;
594  assignFromSparse(rowind, colind, values, indexes);
595  }
596 
597  template<class Treal, class Telement>
599  assignFromSparse(std::vector<int> const & rowind,
600  std::vector<int> const & colind,
601  std::vector<Treal> const & values,
602  std::vector<int> const & indexes) {
603  if (indexes.empty()) {
604  this->clear();
605  return;
606  }
607  if (this->is_zero())
608  allocate();
609 
610  std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
611  std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());
612  int currentInd = 0;
613 
614 
615  std::vector<int>::const_iterator it;
616  for ( it = indexes.begin(); it < indexes.end(); it++ )
617  columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
618 
619  /* Go through all column buckets. */
620  for (int col = 0; col < this->ncols(); col++) {
621  /* Go through current column bucket and distribute to row buckets. */
622  while (!columnBuckets[col].empty()) {
623  currentInd = columnBuckets[col].back();
624  columnBuckets[col].pop_back();
625  rowBuckets[ this->rows.whichBlock
626  ( rowind[currentInd] ) ].push_back (currentInd);
627  }
628  /* Make calls to lower level for every row bucket. */
629  for (int row = 0; row < this->nrows(); row++) {
630  (*this)(row,col).assignFromSparse(rowind, colind, values, rowBuckets[row]);
631  rowBuckets[row].clear();
632  } /* end row loop */
633  } /* end column loop */
634  }
635 
636  template<class Treal, class Telement>
638  addValues(std::vector<int> const & rowind,
639  std::vector<int> const & colind,
640  std::vector<Treal> const & values) {
641  assert(rowind.size() == colind.size() &&
642  rowind.size() == values.size());
643  std::vector<int> indexes(values.size());
644  for (unsigned int ind = 0; ind < values.size(); ++ind)
645  indexes[ind] = ind;
646  addValues(rowind, colind, values, indexes);
647  }
648 
649  template<class Treal, class Telement>
651  addValues(std::vector<int> const & rowind,
652  std::vector<int> const & colind,
653  std::vector<Treal> const & values,
654  std::vector<int> const & indexes) {
655  if (indexes.empty())
656  return;
657  if (this->is_zero())
658  allocate();
659 
660  std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
661  std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());
662  int currentInd = 0;
663 
664  std::vector<int>::const_iterator it;
665  for ( it = indexes.begin(); it < indexes.end(); it++ )
666  columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
667 
668  /* Go through all column buckets. */
669  for (int col = 0; col < this->ncols(); col++) {
670  /* Go through current column bucket and distribute to row buckets. */
671  while (!columnBuckets[col].empty()) {
672  currentInd = columnBuckets[col].back();
673  columnBuckets[col].pop_back();
674  rowBuckets[ this->rows.whichBlock
675  ( rowind[currentInd] ) ].push_back (currentInd);
676  }
677  /* Make calls to lower level for every row bucket. */
678  for (int row = 0; row < this->nrows(); row++) {
679  (*this)(row,col).addValues(rowind, colind, values, rowBuckets[row]);
680  rowBuckets[row].clear();
681  } /* end row loop */
682  } /* end column loop */
683  }
684 
685  template<class Treal, class Telement>
687  syAssignFromSparse(std::vector<int> const & rowind,
688  std::vector<int> const & colind,
689  std::vector<Treal> const & values) {
690  assert(rowind.size() == colind.size() &&
691  rowind.size() == values.size());
692  bool upperTriangleOnly = true;
693  for (unsigned int ind = 0; ind < values.size(); ++ind) {
694  upperTriangleOnly =
695  upperTriangleOnly && colind[ind] >= rowind[ind];
696  }
697  if (!upperTriangleOnly)
698  throw Failure("Matrix<Treal, Telement>::"
699  "syAddValues(std::vector<int> const &, "
700  "std::vector<int> const &, "
701  "std::vector<Treal> const &, int const): "
702  "Only upper triangle can contain elements when assigning "
703  "symmetric or triangular matrix ");
704  assignFromSparse(rowind, colind, values);
705  }
706 
707  template<class Treal, class Telement>
709  syAddValues(std::vector<int> const & rowind,
710  std::vector<int> const & colind,
711  std::vector<Treal> const & values) {
712  assert(rowind.size() == colind.size() &&
713  rowind.size() == values.size());
714  bool upperTriangleOnly = true;
715  for (unsigned int ind = 0; ind < values.size(); ++ind) {
716  upperTriangleOnly =
717  upperTriangleOnly && colind[ind] >= rowind[ind];
718  }
719  if (!upperTriangleOnly)
720  throw Failure("Matrix<Treal, Telement>::"
721  "syAddValues(std::vector<int> const &, "
722  "std::vector<int> const &, "
723  "std::vector<Treal> const &, int const): "
724  "Only upper triangle can contain elements when assigning "
725  "symmetric or triangular matrix ");
726  addValues(rowind, colind, values);
727  }
728 
729  template<class Treal, class Telement>
731  getValues(std::vector<int> const & rowind,
732  std::vector<int> const & colind,
733  std::vector<Treal> & values) const {
734  assert(rowind.size() == colind.size());
735  values.resize(rowind.size(),0);
736  std::vector<int> indexes(rowind.size());
737  for (unsigned int ind = 0; ind < rowind.size(); ++ind)
738  indexes[ind] = ind;
739  getValues(rowind, colind, values, indexes);
740  }
741 
742  template<class Treal, class Telement>
744  getValues(std::vector<int> const & rowind,
745  std::vector<int> const & colind,
746  std::vector<Treal> & values,
747  std::vector<int> const & indexes) const {
748  assert(!this->is_empty());
749  if (indexes.empty())
750  return;
751  std::vector<int>::const_iterator it;
752  if (this->is_zero()) {
753  for ( it = indexes.begin(); it < indexes.end(); it++ )
754  values[*it] = 0;
755  return;
756  }
757 
758  std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
759  std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());
760  int currentInd = 0;
761 
762  for ( it = indexes.begin(); it < indexes.end(); it++ )
763  columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
764 
765  /* Go through all column buckets. */
766  for (int col = 0; col < this->ncols(); col++) {
767  /* Go through current column bucket and distribute to row buckets. */
768  while (!columnBuckets[col].empty()) {
769  currentInd = columnBuckets[col].back();
770  columnBuckets[col].pop_back();
771  rowBuckets[ this->rows.whichBlock
772  ( rowind[currentInd] ) ].push_back (currentInd);
773  }
774  /* Make calls to lower level for every row bucket. */
775  for (int row = 0; row < this->nrows(); row++) {
776  (*this)(row,col).getValues(rowind, colind, values, rowBuckets[row]);
777  rowBuckets[row].clear();
778  } /* end row loop */
779  } /* end column loop */
780  }
781 
782  template<class Treal, class Telement>
784  syGetValues(std::vector<int> const & rowind,
785  std::vector<int> const & colind,
786  std::vector<Treal> & values) const {
787  assert(rowind.size() == colind.size());
788  bool upperTriangleOnly = true;
789  for (int unsigned ind = 0; ind < rowind.size(); ++ind) {
790  upperTriangleOnly =
791  upperTriangleOnly && colind[ind] >= rowind[ind];
792  }
793  if (!upperTriangleOnly)
794  throw Failure("Matrix<Treal, Telement>::"
795  "syGetValues(std::vector<int> const &, "
796  "std::vector<int> const &, "
797  "std::vector<Treal> const &, int const): "
798  "Only upper triangle when retrieving elements from "
799  "symmetric or triangular matrix ");
800  getValues(rowind, colind, values);
801  }
802 
803 
804  template<class Treal, class Telement>
806  getAllValues(std::vector<int> & rowind,
807  std::vector<int> & colind,
808  std::vector<Treal> & values) const {
809  assert(rowind.size() == colind.size() &&
810  rowind.size() == values.size());
811  if (!this->is_zero())
812  for (int ind = 0; ind < this->nElements(); ++ind)
813  this->elements[ind].getAllValues(rowind, colind, values);
814  }
815 
816  template<class Treal, class Telement>
818  syGetAllValues(std::vector<int> & rowind,
819  std::vector<int> & colind,
820  std::vector<Treal> & values) const {
821  assert(rowind.size() == colind.size() &&
822  rowind.size() == values.size());
823  if (!this->is_zero())
824  for (int col = 0; col < this->ncols(); ++col) {
825  for (int row = 0; row < col; ++row)
826  (*this)(row, col).getAllValues(rowind, colind, values);
827  (*this)(col, col).syGetAllValues(rowind, colind, values);
828  }
829  }
830 
831 
832  template<class Treal, class Telement>
834  if (!this->is_zero())
835  for (int i = 0; i < this->nElements(); i++)
836  this->elements[i].clear();
837  freeElements(this->elements);
838  this->elements = 0;
839  }
840 
841  template<class Treal, class Telement>
843  writeToFile(std::ofstream & file) const {
844  int const ZERO = 0;
845  int const ONE = 1;
846  if (this->is_zero()) {
847  char * tmp = (char*)&ZERO;
848  file.write(tmp,sizeof(int));
849  }
850  else {
851  char * tmp = (char*)&ONE;
852  file.write(tmp,sizeof(int));
853  for (int i = 0; i < this->nElements(); i++)
854  this->elements[i].writeToFile(file);
855  }
856  }
857  template<class Treal, class Telement>
859  readFromFile(std::ifstream & file) {
860  int const ZERO = 0;
861  int const ONE = 1;
862  char tmp[sizeof(int)];
863  file.read(tmp, (std::ifstream::pos_type)sizeof(int));
864  switch ((int)*tmp) {
865  case ZERO:
866  this->clear();
867  break;
868  case ONE:
869  if (this->is_zero())
870  allocate();
871  for (int i = 0; i < this->nElements(); i++)
872  this->elements[i].readFromFile(file);
873  break;
874  default:
875  throw Failure("Matrix<Treal, Telement>::"
876  "readFromFile(std::ifstream & file):"
877  "File corruption int value not 0 or 1");
878  }
879  }
880 
881  template<class Treal, class Telement>
883  if (this->is_zero())
884  allocate();
885  for (int ind = 0; ind < this->nElements(); ind++)
886  this->elements[ind].random();
887  }
888 
889  template<class Treal, class Telement>
891  if (this->is_zero())
892  allocate();
893  /* Above diagonal */
894  for (int col = 1; col < this->ncols(); col++)
895  for (int row = 0; row < col; row++)
896  (*this)(row, col).random();
897  /* Diagonal */
898  for (int rc = 0; rc < this->nrows(); rc++)
899  (*this)(rc,rc).syRandom();
900  }
901 
902  template<class Treal, class Telement>
904  randomZeroStructure(Treal probabilityBeingZero) {
905  if (!this->highestLevel() &&
906  probabilityBeingZero > rand() / (Treal)RAND_MAX)
907  this->clear();
908  else {
909  if (this->is_zero())
910  allocate();
911  for (int ind = 0; ind < this->nElements(); ind++)
912  this->elements[ind].randomZeroStructure(probabilityBeingZero);
913  }
914  }
915 
916  template<class Treal, class Telement>
918  syRandomZeroStructure(Treal probabilityBeingZero) {
919  if (!this->highestLevel() &&
920  probabilityBeingZero > rand() / (Treal)RAND_MAX)
921  this->clear();
922  else {
923  if (this->is_zero())
924  allocate();
925  /* Above diagonal */
926  for (int col = 1; col < this->ncols(); col++)
927  for (int row = 0; row < col; row++)
928  (*this)(row, col).randomZeroStructure(probabilityBeingZero);
929  /* Diagonal */
930  for (int rc = 0; rc < this->nrows(); rc++)
931  (*this)(rc,rc).syRandomZeroStructure(probabilityBeingZero);
932  }
933  }
934 
935  template<class Treal, class Telement>
936  template<typename TRule>
938  setElementsByRule(TRule & rule) {
939  if (this->is_zero())
940  allocate();
941  for (int ind = 0; ind < this->nElements(); ind++)
942  this->elements[ind].setElementsByRule(rule);
943  }
944  template<class Treal, class Telement>
945  template<typename TRule>
947  sySetElementsByRule(TRule & rule) {
948  if (this->is_zero())
949  allocate();
950  /* Above diagonal */
951  for (int col = 1; col < this->ncols(); col++)
952  for (int row = 0; row < col; row++)
953  (*this)(row, col).setElementsByRule(rule);
954  /* Diagonal */
955  for (int rc = 0; rc < this->nrows(); rc++)
956  (*this)(rc,rc).sySetElementsByRule(rule);
957  }
958 
959 
960  template<class Treal, class Telement>
962  addIdentity(Treal alpha) {
963  if (this->is_empty())
964  throw Failure("Matrix<Treal, Telement>::addIdentity(Treal): "
965  "Cannot add identity to empty matrix.");
966  if (this->ncols() != this->nrows())
967  throw Failure("Matrix<Treal, Telement>::addIdentity(Treal): "
968  "Matrix must be square to add identity");
969  if (this->is_zero())
970  allocate();
971  for (int ind = 0; ind < this->ncols(); ind++)
972  (*this)(ind,ind).addIdentity(alpha);
973  }
974 
975  template<class Treal, class Telement>
979  if (A.is_zero()) { /* Condition also matches empty matrices. */
980  AT.rows = A.cols;
981  AT.cols = A.rows;
983  AT.elements = 0;
984  return;
985  }
986  if (AT.is_zero() || (AT.nElements() != A.nElements())) {
988  AT.elements = allocateElements<Telement>(A.nElements());
989  }
990  AT.cols = A.rows;
991  AT.rows = A.cols;
992  for (int row = 0; row < AT.nrows(); row++)
993  for (int col = 0; col < AT.ncols(); col++)
994  Telement::transpose(A(col,row), AT(row,col));
995  }
996 
997  template<class Treal, class Telement>
1000  try {
1001  if (this->nrows() == this->ncols()) {
1002  if (!this->is_zero()) {
1003  /* Fix the diagonal: */
1004  for (int rc = 0; rc < this->ncols(); rc++)
1005  (*this)(rc, rc).symToNosym();
1006  /* Fix the lower triangle */
1007  for (int row = 1; row < this->nrows(); row++)
1008  for (int col = 0; col < row; col++)
1009  Telement::transpose((*this)(col, row), (*this)(row, col));
1010  }
1011  }
1012  else
1013  throw Failure("Matrix<Treal, Telement>::symToNosym(): "
1014  "Only quadratic matrices can be symmetric");
1015  }
1016  catch(Failure& f) {
1017  std::cout<<"Failure caught:Matrix<Treal, Telement>::symToNosym()"
1018  <<std::endl;
1019  throw f;
1020  }
1021  }
1022 
1023  template<class Treal, class Telement>
1026  if (this->nrows() == this->ncols()) {
1027  if (!this->is_zero()) {
1028  /* Fix the diagonal: */
1029  for (int rc = 0; rc < this->ncols(); rc++)
1030  (*this)(rc, rc).nosymToSym();
1031  /* Fix the lower triangle */
1032  for (int col = 0; col < this->ncols() - 1; col++)
1033  for (int row = col + 1; row < this->nrows(); row++)
1034  (*this)(row, col).clear();
1035  }
1036  }
1037  else
1038  throw Failure("Matrix<Treal, Telement>::nosymToSym(): "
1039  "Only quadratic matrices can be symmetric");
1040  }
1041 
1042  /* Basic linear algebra routines */
1043 
1044  template<class Treal, class Telement>
1047  switch (k) {
1048  case 0:
1049  this->clear();
1050  break;
1051  case 1:
1052  if (this->ncols() != this->nrows())
1053  throw Failure
1054  ("Matrix::operator=(int k = 1): "
1055  "Matrix must be quadratic to become identity matrix.");
1056  this->clear();
1057  this->allocate();
1058  for (int rc = 0; rc < this->ncols(); rc++) /*Set diagonal to identity*/
1059  (*this)(rc,rc) = 1;
1060  break;
1061  default:
1062  throw Failure("Matrix::operator=(int k) only "
1063  "implemented for k = 0, k = 1");
1064  }
1065  return *this;
1066  }
1067 
1068 
1069  template<class Treal, class Telement>
1071  operator*=(const Treal alpha) {
1072  if (!this->is_zero() && alpha != 1) {
1073  for (int ind = 0; ind < this->nElements(); ind++)
1074  this->elements[ind] *= alpha;
1075  }
1076  return *this;
1077  }
1078 
1079  /* C = beta * C + alpha * A * B */
1080  template<class Treal, class Telement>
1082  gemm(const bool tA, const bool tB, const Treal alpha,
1083  const Matrix<Treal, Telement>& A,
1084  const Matrix<Treal, Telement>& B,
1085  const Treal beta,
1087  assert(!A.is_empty());
1088  assert(!B.is_empty());
1089  if (C.is_empty()) {
1090  assert(beta == 0);
1091  if (tA)
1092  C.resetRows(A.cols);
1093  else
1094  C.resetRows(A.rows);
1095  if (tB)
1096  C.resetCols(B.rows);
1097  else
1098  C.resetCols(B.cols);
1099  }
1100 
1101  if ((A.is_zero() || B.is_zero() || alpha == 0) &&
1102  (C.is_zero() || beta == 0))
1103  C = 0;
1104  else {
1105  Treal beta_tmp = beta;
1106  if (C.is_zero()) {
1107  C.allocate();
1108  beta_tmp = 0;
1109  }
1110  if (!A.is_zero() && !B.is_zero() && alpha != 0) {
1111  MAT_OMP_INIT;
1112  if (!tA && !tB)
1113  if (A.ncols() == B.nrows() &&
1114  A.nrows() == C.nrows() &&
1115  B.ncols() == C.ncols()) {
1116 #ifdef _OPENMP
1117 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1118 #endif
1119  for (int col = 0; col < C.ncols(); col++) {
1120  MAT_OMP_START;
1121  for (int row = 0; row < C.nrows(); row++) {
1122  Telement::gemm(tA, tB, alpha,
1123  A(row, 0), B(0, col),
1124  beta_tmp,
1125  C(row, col));
1126  for(int cola = 1; cola < A.ncols(); cola++)
1127  Telement::gemm(tA, tB, alpha,
1128  A(row, cola), B(cola, col),
1129  1.0,
1130  C(row, col));
1131  }
1132  MAT_OMP_END;
1133  } /* end omp for */
1134  }
1135  else
1136  throw Failure("Matrix<class Treal, class Telement>::"
1137  "gemm(bool, bool, Treal, "
1138  "Matrix<Treal, Telement>, "
1139  "Matrix<Treal, Telement>, Treal, "
1140  "Matrix<Treal, Telement>): "
1141  "Incorrect matrixdimensions for multiplication");
1142  else if (tA && !tB)
1143  if (A.nrows() == B.nrows() &&
1144  A.ncols() == C.nrows() &&
1145  B.ncols() == C.ncols()) {
1146 #ifdef _OPENMP
1147 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1148 #endif
1149  for (int col = 0; col < C.ncols(); col++) {
1150  MAT_OMP_START;
1151  for (int row = 0; row < C.nrows(); row++) {
1152  Telement::gemm(tA, tB, alpha,
1153  A(0,row), B(0,col),
1154  beta_tmp,
1155  C(row,col));
1156  for(int cola = 1; cola < A.nrows(); cola++)
1157  Telement::gemm(tA, tB, alpha,
1158  A(cola, row), B(cola, col),
1159  1.0,
1160  C(row,col));
1161  }
1162  MAT_OMP_END;
1163  } /* end omp for */
1164  }
1165  else
1166  throw Failure("Matrix<class Treal, class Telement>::"
1167  "gemm(bool, bool, Treal, "
1168  "Matrix<Treal, Telement>, "
1169  "Matrix<Treal, Telement>, Treal, "
1170  "Matrix<Treal, Telement>): "
1171  "Incorrect matrixdimensions for multiplication");
1172  else if (!tA && tB)
1173  if (A.ncols() == B.ncols() &&
1174  A.nrows() == C.nrows() &&
1175  B.nrows() == C.ncols()) {
1176 #ifdef _OPENMP
1177 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1178 #endif
1179  for (int col = 0; col < C.ncols(); col++) {
1180  MAT_OMP_START;
1181  for (int row = 0; row < C.nrows(); row++) {
1182  Telement::gemm(tA, tB, alpha,
1183  A(row, 0), B(col, 0),
1184  beta_tmp,
1185  C(row,col));
1186  for(int cola = 1; cola < A.ncols(); cola++)
1187  Telement::gemm(tA, tB, alpha,
1188  A(row, cola), B(col, cola),
1189  1.0,
1190  C(row,col));
1191  }
1192  MAT_OMP_END;
1193  } /* end omp for */
1194  }
1195  else
1196  throw Failure("Matrix<class Treal, class Telement>::"
1197  "gemm(bool, bool, Treal, "
1198  "Matrix<Treal, Telement>, "
1199  "Matrix<Treal, Telement>, Treal, "
1200  "Matrix<Treal, Telement>): "
1201  "Incorrect matrixdimensions for multiplication");
1202  else if (tA && tB)
1203  if (A.nrows() == B.ncols() &&
1204  A.ncols() == C.nrows() &&
1205  B.nrows() == C.ncols()) {
1206 #ifdef _OPENMP
1207 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1208 #endif
1209  for (int col = 0; col < C.ncols(); col++) {
1210  MAT_OMP_START;
1211  for (int row = 0; row < C.nrows(); row++) {
1212  Telement::gemm(tA, tB, alpha,
1213  A(0, row), B(col, 0),
1214  beta_tmp,
1215  C(row,col));
1216  for(int cola = 1; cola < A.nrows(); cola++)
1217  Telement::gemm(tA, tB, alpha,
1218  A(cola, row), B(col, cola),
1219  1.0,
1220  C(row,col));
1221  }
1222  MAT_OMP_END;
1223  } /* end omp for */
1224  }
1225  else
1226  throw Failure("Matrix<class Treal, class Telement>::"
1227  "gemm(bool, bool, Treal, "
1228  "Matrix<Treal, Telement>, "
1229  "Matrix<Treal, Telement>, "
1230  "Treal, Matrix<Treal, Telement>): "
1231  "Incorrect matrixdimensions for multiplication");
1232  else throw Failure("Matrix<class Treal, class Telement>::"
1233  "gemm(bool, bool, Treal, "
1234  "Matrix<Treal, Telement>, "
1235  "Matrix<Treal, Telement>, Treal, "
1236  "Matrix<Treal, Telement>):"
1237  "Very strange error!!");
1239  }
1240  else
1241  C *= beta;
1242  }
1243  }
1244 
1245  template<class Treal, class Telement>
1247  symm(const char side, const char uplo, const Treal alpha,
1248  const Matrix<Treal, Telement>& A,
1249  const Matrix<Treal, Telement>& B,
1250  const Treal beta,
1252  assert(!A.is_empty());
1253  assert(!B.is_empty());
1254  assert(A.nrows() == A.ncols());
1255  int dimA = A.nrows();
1256  if (C.is_empty()) {
1257  assert(beta == 0);
1258  if (side =='L') {
1259  C.resetRows(A.rows);
1260  C.resetCols(B.cols);
1261  }
1262  else {
1263  assert(side == 'R');
1264  C.resetRows(B.rows);
1265  C.resetCols(A.cols);
1266  }
1267  }
1268  if ((A.is_zero() || B.is_zero() || alpha == 0) &&
1269  (C.is_zero() || beta == 0))
1270  C = 0;
1271  else {
1272  if (uplo == 'U') {
1273  Treal beta_tmp = beta;
1274  if (C.is_zero()) {
1275  C.allocate();
1276  beta_tmp = 0;
1277  }
1278  if (!A.is_zero() && !B.is_zero() && alpha != 0) {
1279  MAT_OMP_INIT;
1280  if (side =='L')
1281  if (dimA == B.nrows() &&
1282  dimA == C.nrows() &&
1283  B.ncols() == C.ncols()) {
1284 #ifdef _OPENMP
1285 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1286 #endif
1287  for (int col = 0; col < C.ncols(); col++) {
1288  MAT_OMP_START;
1289  for (int row = 0; row < C.nrows(); row++) {
1290  /* Diagonal element in matrix A */
1291  Telement::symm(side, uplo, alpha, A(row, row),
1292  B(row, col), beta_tmp, C(row, col));
1293  /* Elements below diagonal in A */
1294  for(int ind = 0; ind < row; ind++)
1295  Telement::gemm(true, false, alpha,
1296  A(ind, row), B(ind, col),
1297  1.0, C(row,col));
1298  /* Elements above diagonal in A */
1299  for(int ind = row + 1; ind < dimA; ind++)
1300  Telement::gemm(false, false, alpha,
1301  A(row, ind), B(ind, col),
1302  1.0, C(row,col));
1303  }
1304  MAT_OMP_END;
1305  } /* end omp for */
1306  }
1307  else
1308  throw Failure("Matrix<class Treal, class Telement>"
1309  "::symm(bool, bool, Treal, Matrix<Treal, "
1310  "Telement>, Matrix<Treal, Telement>, "
1311  "Treal, Matrix<Treal, Telement>): "
1312  "Incorrect matrixdimensions for multiplication");
1313  else { /* side == 'R' */
1314  assert(side == 'R');
1315  if (B.ncols() == dimA &&
1316  B.nrows() == C.nrows() &&
1317  dimA == C.ncols()) {
1318 #ifdef _OPENMP
1319 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1320 #endif
1321  for (int col = 0; col < C.ncols(); col++) {
1322  MAT_OMP_START;
1323  for (int row = 0; row < C.nrows(); row++) {
1324  /* Diagonal element in matrix A */
1325  Telement::symm(side, uplo, alpha, A(col, col),
1326  B(row, col), beta_tmp, C(row, col));
1327  /* Elements below diagonal in A */
1328  for(int ind = col + 1; ind < dimA; ind++)
1329  Telement::gemm(false, true, alpha,
1330  B(row, ind), A(col, ind),
1331  1.0, C(row,col));
1332  /* Elements above diagonal in A */
1333  for(int ind = 0; ind < col; ind++)
1334  Telement::gemm(false, false, alpha,
1335  B(row, ind), A(ind, col),
1336  1.0, C(row,col));
1337  }
1338  MAT_OMP_END;
1339  } /* end omp for */
1340  }
1341  else
1342  throw Failure("Matrix<class Treal, class Telement>"
1343  "::symm(bool, bool, Treal, Matrix<Treal, "
1344  "Telement>, Matrix<Treal, Telement>, "
1345  "Treal, Matrix<Treal, Telement>): "
1346  "Incorrect matrixdimensions for multiplication");
1347  }
1349  }
1350  else
1351  C *= beta;
1352  }
1353  else
1354  throw Failure("Matrix<class Treal, class Telement>::"
1355  "symm only implemented for symmetric matrices in "
1356  "upper triangular storage");
1357  }
1358  }
1359 
1360  template<class Treal, class Telement>
1362  syrk(const char uplo, const bool tA, const Treal alpha,
1363  const Matrix<Treal, Telement>& A,
1364  const Treal beta,
1366  assert(!A.is_empty());
1367  if (C.is_empty()) {
1368  assert(beta == 0);
1369  if (tA) {
1370  C.resetRows(A.cols);
1371  C.resetCols(A.cols);
1372  }
1373  else {
1374  C.resetRows(A.rows);
1375  C.resetCols(A.rows);
1376  }
1377  }
1378 
1379  if (C.nrows() == C.ncols() &&
1380  ((!tA && A.nrows() == C.nrows()) || (tA && A.ncols() == C.nrows())))
1381  if (alpha != 0 && !A.is_zero()) {
1382  Treal beta_tmp = beta;
1383  if (C.is_zero()) {
1384  C.allocate();
1385  beta_tmp = 0;
1386  }
1387  MAT_OMP_INIT;
1388  if (!tA && uplo == 'U') { /* C = alpha * A * A' + beta * C */
1389 #ifdef _OPENMP
1390 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
1391 #endif
1392  {
1393 #ifdef _OPENMP
1394 #pragma omp for schedule(dynamic) nowait
1395 #endif
1396  for (int rc = 0; rc < C.ncols(); rc++) {
1397  MAT_OMP_START;
1398  Telement::syrk(uplo, tA, alpha, A(rc, 0), beta_tmp, C(rc, rc));
1399  for (int cola = 1; cola < A.ncols(); cola++)
1400  Telement::syrk(uplo, tA, alpha, A(rc, cola), 1.0, C(rc, rc));
1401  MAT_OMP_END;
1402  }
1403 #ifdef _OPENMP
1404 #pragma omp for schedule(dynamic) nowait
1405 #endif
1406  for (int row = 0; row < C.nrows() - 1; row++) {
1407  MAT_OMP_START;
1408  for (int col = row + 1; col < C.ncols(); col++) {
1409  Telement::gemm(tA, !tA, alpha,
1410  A(row, 0), A(col,0), beta_tmp, C(row,col));
1411  for (int ind = 1; ind < A.ncols(); ind++)
1412  Telement::gemm(tA, !tA, alpha,
1413  A(row, ind), A(col,ind), 1.0, C(row,col));
1414  }
1415  MAT_OMP_END;
1416  }
1417  } /* end omp parallel */
1418  } /* end if (!tA) */
1419  else if (tA && uplo == 'U') { /* C = alpha * A' * A + beta * C */
1420 #ifdef _OPENMP
1421 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
1422 #endif
1423  {
1424 #ifdef _OPENMP
1425 #pragma omp for schedule(dynamic) nowait
1426 #endif
1427  for (int rc = 0; rc < C.ncols(); rc++) {
1428  MAT_OMP_START;
1429  Telement::syrk(uplo, tA, alpha, A(0, rc), beta_tmp, C(rc, rc));
1430  for (int rowa = 1; rowa < A.nrows(); rowa++)
1431  Telement::syrk(uplo, tA, alpha, A(rowa, rc), 1.0, C(rc, rc));
1432  MAT_OMP_END;
1433  }
1434 #ifdef _OPENMP
1435 #pragma omp for schedule(dynamic) nowait
1436 #endif
1437  for (int row = 0; row < C.nrows() - 1; row++) {
1438  MAT_OMP_START;
1439  for (int col = row + 1; col < C.ncols(); col++) {
1440  Telement::gemm(tA, !tA, alpha,
1441  A(0, row), A(0, col), beta_tmp, C(row,col));
1442  for (int ind = 1; ind < A.nrows(); ind++)
1443  Telement::gemm(tA, !tA, alpha,
1444  A(ind, row), A(ind, col), 1.0, C(row,col));
1445  }
1446  MAT_OMP_END;
1447  }
1448  } /* end omp parallel */
1449  } /* end if (tA) */
1450  else
1451  throw Failure("Matrix<class Treal, class Telement>::"
1452  "syrk not implemented for lower triangular storage");
1454  }
1455  else {
1456  C *= beta;
1457  }
1458  else
1459  throw Failure("Matrix<class Treal, class Telement>::syrk: "
1460  "Incorrect matrix dimensions for symmetric rank-k update");
1461  }
1462 
1463  template<class Treal, class Telement>
1465  sysq(const char uplo, const Treal alpha,
1466  const Matrix<Treal, Telement>& A,
1467  const Treal beta,
1469  assert(!A.is_empty());
1470  if (C.is_empty()) {
1471  assert(beta == 0);
1472  C.resetRows(A.rows);
1473  C.resetCols(A.cols);
1474  }
1475  if (C.nrows() == C.ncols() &&
1476  A.nrows() == C.nrows() && A.nrows() == A.ncols())
1477  if (alpha != 0 && !A.is_zero()) {
1478  if (uplo == 'U') {
1479  Treal beta_tmp = beta;
1480  if (C.is_zero()) {
1481  C.allocate();
1482  beta_tmp = 0;
1483  }
1484  MAT_OMP_INIT;
1485 #ifdef _OPENMP
1486 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
1487 #endif
1488  {
1489 #ifdef _OPENMP
1490 #pragma omp for schedule(dynamic) nowait
1491 #endif
1492  for (int rc = 0; rc < C.ncols(); rc++) {
1493  MAT_OMP_START;
1494  Telement::sysq(uplo, alpha, A(rc, rc), beta_tmp, C(rc, rc));
1495  for (int cola = 0; cola < rc; cola++)
1496  Telement::syrk(uplo, true, alpha, A(cola, rc), 1.0, C(rc, rc));
1497  for (int cola = rc + 1; cola < A.ncols(); cola++)
1498  Telement::syrk(uplo, false, alpha, A(rc, cola), 1.0, C(rc, rc));
1499  MAT_OMP_END;
1500  }
1501  /* Maste anvanda symm? */
1502 #ifdef _OPENMP
1503 #pragma omp for schedule(dynamic) nowait
1504 #endif
1505  for (int row = 0; row < C.nrows() - 1; row++) {
1506  MAT_OMP_START;
1507  for (int col = row + 1; col < C.ncols(); col++) {
1508  /* First the two operations involving diagonal elements */
1509  Telement::symm('L', 'U', alpha, A(row, row), A(row, col),
1510  beta_tmp, C(row, col));
1511  Telement::symm('R', 'U', alpha, A(col, col), A(row, col),
1512  1.0, C(row, col));
1513  /* First element in product is below the diagonal */
1514  for (int ind = 0; ind < row; ind++)
1515  Telement::gemm(true, false, alpha,
1516  A(ind, row), A(ind, col), 1.0, C(row, col));
1517  /* None of the elements are below the diagonal */
1518  for (int ind = row + 1; ind < col; ind++)
1519  Telement::gemm(false, false, alpha,
1520  A(row, ind), A(ind, col), 1.0, C(row, col));
1521  /* Second element is below the diagonal */
1522  for (int ind = col + 1; ind < A.ncols(); ind++)
1523  Telement::gemm(false, true, alpha,
1524  A(row, ind), A(col, ind), 1.0, C(row, col));
1525  }
1526  MAT_OMP_END;
1527  } /* end omp for */
1528  } /* end omp parallel */
1530  }
1531  else
1532  throw Failure("Matrix<class Treal, class Telement>::"
1533  "sysq only implemented for symmetric matrices in "
1534  "upper triangular storage");;
1535  }
1536  else {
1537  C *= beta;
1538  }
1539  else
1540  throw Failure("Matrix<class Treal, class Telement>::sysq: "
1541  "Incorrect matrix dimensions for symmetric square "
1542  "operation");
1543  }
1544 
1545  /* C = alpha * A * B + beta * C where A and B are symmetric */
1546  template<class Treal, class Telement>
1548  ssmm(const Treal alpha,
1549  const Matrix<Treal, Telement>& A,
1550  const Matrix<Treal, Telement>& B,
1551  const Treal beta,
1553  assert(!A.is_empty());
1554  assert(!B.is_empty());
1555  if (C.is_empty()) {
1556  assert(beta == 0);
1557  C.resetRows(A.rows);
1558  C.resetCols(B.cols);
1559  }
1560  if (A.ncols() != B.nrows() ||
1561  A.nrows() != C.nrows() ||
1562  B.ncols() != C.ncols() ||
1563  A.nrows() != A.ncols() ||
1564  B.nrows() != B.ncols()) {
1565  throw Failure("Matrix<class Treal, class Telement>::"
1566  "ssmm(Treal, "
1567  "Matrix<Treal, Telement>, "
1568  "Matrix<Treal, Telement>, Treal, "
1569  "Matrix<Treal, Telement>): "
1570  "Incorrect matrixdimensions for ssmm multiplication");
1571  }
1572  if ((A.is_zero() || B.is_zero() || alpha == 0) &&
1573  (C.is_zero() || beta == 0)) {
1574  C = 0;
1575  return;
1576  }
1577  Treal beta_tmp = beta;
1578  if (C.is_zero()) {
1579  C.allocate();
1580  beta_tmp = 0;
1581  }
1582  if (A.is_zero() || B.is_zero() || alpha == 0) {
1583  C *= beta;
1584  return;
1585  }
1586  MAT_OMP_INIT;
1587 #ifdef _OPENMP
1588 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1589 #endif
1590  for (int col = 0; col < C.ncols(); col++) {
1591  MAT_OMP_START;
1592  /* Diagonal */
1593  /* C(col, col) = alpha * A(col, col) * B(col, col) + beta * C(col, col)*/
1594  Telement::ssmm(alpha, A(col,col), B(col, col),
1595  beta_tmp, C(col,col));
1596  for (int ind = 0; ind < col; ind++)
1597  /* C(col, col) += alpha * A(ind, col)' * B(ind, col) */
1598  Telement::gemm(true, false,
1599  alpha, A(ind,col), B(ind, col),
1600  1.0, C(col,col));
1601  for (int ind = col + 1; ind < A.ncols(); ind++)
1602  /* C(col, col) += alpha * A(col, ind) * B(col, ind)' */
1603  Telement::gemm(false, true,
1604  alpha, A(col, ind), B(col, ind),
1605  1.0, C(col,col));
1606  /* Upper half */
1607  for (int row = 0; row < col; row++) {
1608  /* C(row, col) = alpha * A(row, row) * B(row, col) +
1609  * beta * C(row, col)
1610  */
1611  Telement::symm('L', 'U',
1612  alpha, A(row, row), B(row, col),
1613  beta_tmp, C(row, col));
1614  /* C(row, col) += alpha * A(row, col) * B(col, col) */
1615  Telement::symm('R', 'U',
1616  alpha, B(col, col), A(row, col),
1617  1.0, C(row, col));
1618  for (int ind = 0; ind < row; ind++)
1619  /* C(row, col) += alpha * A(ind, row)' * B(ind, col) */
1620  Telement::gemm(true, false,
1621  alpha, A(ind, row), B(ind, col),
1622  1.0, C(row,col));
1623  for (int ind = row + 1; ind < col; ind++)
1624  /* C(row, col) += alpha * A(row, ind) * B(ind, col) */
1625  Telement::gemm(false, false,
1626  alpha, A(row, ind), B(ind, col),
1627  1.0, C(row,col));
1628  for (int ind = col + 1; ind < A.ncols(); ind++)
1629  /* C(row, col) += alpha * A(row, ind) * B(col, ind)' */
1630  Telement::gemm(false, true,
1631  alpha, A(row, ind), B(col, ind),
1632  1.0, C(row,col));
1633  }
1634  /* Lower half */
1635  Telement tmpSubMat;
1636  for (int row = col + 1; row < C.nrows(); row++) {
1637  Telement::transpose(C(row, col), tmpSubMat);
1638  /* tmpSubMat = alpha * B(col, col) * A(col, row) +
1639  * beta * tmpSubMat
1640  */
1641  Telement::symm('L', 'U',
1642  alpha, B(col, col), A(col, row),
1643  beta_tmp, tmpSubMat);
1644  /* tmpSubMat += alpha * B(col, row) * A(row, row) */
1645  Telement::symm('R', 'U',
1646  alpha, A(row, row), B(col, row),
1647  1.0, tmpSubMat);
1648  for (int ind = 0; ind < col; ind++)
1649  /* tmpSubMat += alpha * B(ind, col)' * A(ind, row) */
1650  Telement::gemm(true, false,
1651  alpha, B(ind, col), A(ind, row),
1652  1.0, tmpSubMat);
1653  for (int ind = col + 1; ind < row; ind++)
1654  /* tmpSubMat += alpha * B(col, ind) * A(ind, row) */
1655  Telement::gemm(false, false,
1656  alpha, B(col, ind), A(ind, row),
1657  1.0, tmpSubMat);
1658  for (int ind = row + 1; ind < B.nrows(); ind++)
1659  /* tmpSubMat += alpha * B(col, ind) * A(row, ind)' */
1660  Telement::gemm(false, true,
1661  alpha, B(col, ind), A(row, ind),
1662  1.0, tmpSubMat);
1663  Telement::transpose(tmpSubMat, C(row, col));
1664  }
1665  MAT_OMP_END;
1666  }
1668  } /* end ssmm */
1669 
1670 
1671  /* C = alpha * A * B + beta * C where A and B are symmetric
1672  * and only the upper triangle of C is computed.
1673  */
1674  template<class Treal, class Telement>
1676  ssmm_upper_tr_only(const Treal alpha,
1677  const Matrix<Treal, Telement>& A,
1678  const Matrix<Treal, Telement>& B,
1679  const Treal beta,
1681  assert(!A.is_empty());
1682  assert(!B.is_empty());
1683  if (C.is_empty()) {
1684  assert(beta == 0);
1685  C.resetRows(A.rows);
1686  C.resetCols(B.cols);
1687  }
1688  if (A.ncols() != B.nrows() ||
1689  A.nrows() != C.nrows() ||
1690  B.ncols() != C.ncols() ||
1691  A.nrows() != A.ncols() ||
1692  B.nrows() != B.ncols()) {
1693  throw Failure("Matrix<class Treal, class Telement>::"
1694  "ssmm_upper_tr_only(Treal, "
1695  "Matrix<Treal, Telement>, "
1696  "Matrix<Treal, Telement>, Treal, "
1697  "Matrix<Treal, Telement>): "
1698  "Incorrect matrixdimensions for ssmm_upper_tr_only "
1699  "multiplication");
1700  }
1701  if ((A.is_zero() || B.is_zero() || alpha == 0) &&
1702  (C.is_zero() || beta == 0)) {
1703  C = 0;
1704  return;
1705  }
1706  Treal beta_tmp = beta;
1707  if (C.is_zero()) {
1708  C.allocate();
1709  beta_tmp = 0;
1710  }
1711  if (A.is_zero() || B.is_zero() || alpha == 0) {
1712  C *= beta;
1713  return;
1714  }
1715  MAT_OMP_INIT;
1716 #ifdef _OPENMP
1717 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1718 #endif
1719  for (int col = 0; col < C.ncols(); col++) {
1720  MAT_OMP_START;
1721  /* Diagonal */
1722  /* C(col, col) = alpha * A(col, col) * B(col, col) + beta * C(col, col)*/
1723  Telement::ssmm_upper_tr_only(alpha, A(col,col), B(col, col),
1724  beta_tmp, C(col,col));
1725  for (int ind = 0; ind < col; ind++)
1726  /* C(col, col) += alpha * A(ind, col)' * B(ind, col) */
1727  Telement::gemm_upper_tr_only(true, false,
1728  alpha, A(ind,col), B(ind, col),
1729  1.0, C(col,col));
1730  for (int ind = col + 1; ind < A.ncols(); ind++)
1731  /* C(col, col) += alpha * A(col, ind) * B(col, ind)' */
1732  Telement::gemm_upper_tr_only(false, true,
1733  alpha, A(col, ind), B(col, ind),
1734  1.0, C(col,col));
1735  /* Upper half */
1736  for (int row = 0; row < col; row++) {
1737  /* C(row, col) = alpha * A(row, row) * B(row, col) +
1738  * beta * C(row, col)
1739  */
1740  Telement::symm('L', 'U',
1741  alpha, A(row, row), B(row, col),
1742  beta_tmp, C(row, col));
1743  /* C(row, col) += alpha * A(row, col) * B(col, col) */
1744  Telement::symm('R', 'U',
1745  alpha, B(col, col), A(row, col),
1746  1.0, C(row, col));
1747  for (int ind = 0; ind < row; ind++)
1748  /* C(row, col) += alpha * A(ind, row)' * B(ind, col) */
1749  Telement::gemm(true, false,
1750  alpha, A(ind, row), B(ind, col),
1751  1.0, C(row,col));
1752  for (int ind = row + 1; ind < col; ind++)
1753  /* C(row, col) += alpha * A(row, ind) * B(ind, col) */
1754  Telement::gemm(false, false,
1755  alpha, A(row, ind), B(ind, col),
1756  1.0, C(row,col));
1757  for (int ind = col + 1; ind < A.ncols(); ind++)
1758  /* C(row, col) += alpha * A(row, ind) * B(col, ind)' */
1759  Telement::gemm(false, true,
1760  alpha, A(row, ind), B(col, ind),
1761  1.0, C(row,col));
1762  }
1763  MAT_OMP_END;
1764  }
1766  } /* end ssmm_upper_tr_only */
1767 
1768 
1769 
1770  template<class Treal, class Telement>
1772  trmm(const char side, const char uplo, const bool tA, const Treal alpha,
1773  const Matrix<Treal, Telement>& A,
1775  assert(!A.is_empty());
1776  assert(!B.is_empty());
1777  if (alpha != 0 && !A.is_zero() && !B.is_zero())
1778  if (((side == 'R' && B.ncols() == A.nrows()) ||
1779  (side == 'L' && A.ncols() == B.nrows())) &&
1780  A.nrows() == A.ncols())
1781  if (uplo == 'U')
1782  if (!tA) {
1783  if (side == 'R') {
1784  /* Last column must be calculated first */
1785  for (int col = B.ncols() - 1; col >= 0; col--)
1786  for (int row = 0; row < B.nrows(); row++) {
1787  /* Use first the diagonal element in A */
1788  /* Otherwise the faster call to trmm cannot be utilized */
1789  Telement::trmm(side, uplo, tA, alpha,
1790  A(col, col), B(row,col));
1791  /* And the rest: */
1792  for (int ind = 0; ind < col; ind++)
1793  Telement::gemm(false, tA, alpha,
1794  B(row,ind), A(ind, col),
1795  1.0, B(row,col));
1796  }
1797  } /* end if (side == 'R')*/
1798  else {
1799  assert(side == 'L');
1800  /* First row must be calculated first */
1801  for (int row = 0; row < B.nrows(); row++ )
1802  for (int col = 0; col < B.ncols(); col++) {
1803  Telement::trmm(side, uplo, tA, alpha,
1804  A(row,row), B(row,col));
1805  for (int ind = row + 1 ; ind < B.nrows(); ind++)
1806  Telement::gemm(tA, false, alpha,
1807  A(row,ind), B(ind,col),
1808  1.0, B(row,col));
1809  }
1810  } /* end else (side == 'L')*/
1811  } /* end if (tA == false) */
1812  else {
1813  assert(tA == true);
1814  if (side == 'R') {
1815  /* First column must be calculated first */
1816  for (int col = 0; col < B.ncols(); col++)
1817  for (int row = 0; row < B.nrows(); row++) {
1818  Telement::trmm(side, uplo, tA, alpha,
1819  A(col,col), B(row,col));
1820  for (int ind = col + 1; ind < A.ncols(); ind++)
1821  Telement::gemm(false, tA, alpha,
1822  B(row,ind), A(col,ind),
1823  1.0, B(row,col));
1824  }
1825  } /* end if (side == 'R')*/
1826  else {
1827  assert(side == 'L');
1828  /* Last row must be calculated first */
1829  for (int row = B.nrows() - 1; row >= 0; row--)
1830  for (int col = 0; col < B.ncols(); col++) {
1831  Telement::trmm(side, uplo, tA, alpha,
1832  A(row,row), B(row,col));
1833  for (int ind = 0; ind < row; ind++)
1834  Telement::gemm(tA, false, alpha,
1835  A(ind,row), B(ind,col),
1836  1.0, B(row,col));
1837  }
1838  } /* end else (side == 'L')*/
1839  } /* end else (tA == true)*/
1840  else
1841  throw Failure("Matrix<class Treal, class Telement>::"
1842  "trmm not implemented for lower triangular matrices");
1843  else
1844  throw Failure("Matrix<class Treal, class Telement>::trmm"
1845  ": Incorrect matrix dimensions for multiplication");
1846  else
1847  B = 0;
1848  }
1849 
1850 
1851  template<class Treal, class Telement>
1853  assert(!this->is_empty());
1854  if (this->is_zero())
1855  return 0;
1856  else {
1857  Treal sum(0);
1858  for (int i = 0; i < this->nElements(); i++)
1859  sum += this->elements[i].frobSquared();
1860  return sum;
1861  }
1862  }
1863 
1864  template<class Treal, class Telement>
1867  assert(!this->is_empty());
1868  if (this->nrows() != this->ncols())
1869  throw Failure("Matrix<Treal, Telement>::syFrobSquared(): "
1870  "Matrix must be have equal number of rows "
1871  "and cols to be symmetric");
1872  Treal sum(0);
1873  if (!this->is_zero()) {
1874  for (int col = 1; col < this->ncols(); col++)
1875  for (int row = 0; row < col; row++)
1876  sum += 2 * (*this)(row, col).frobSquared();
1877  for (int rc = 0; rc < this->ncols(); rc++)
1878  sum += (*this)(rc, rc).syFrobSquared();
1879  }
1880  return sum;
1881  }
1882 
1883  template<class Treal, class Telement>
1887  const Matrix<Treal, Telement>& B) {
1888  assert(!A.is_empty());
1889  assert(!B.is_empty());
1890  if (A.nrows() != B.nrows() || A.ncols() != B.ncols())
1891  throw Failure("Matrix<Treal, Telement>::"
1892  "frobSquaredDiff: Incorrect matrix dimensions");
1893  Treal sum(0);
1894  if (!A.is_zero() && !B.is_zero())
1895  for (int i = 0; i < A.nElements(); i++)
1896  sum += Telement::frobSquaredDiff(A.elements[i],B.elements[i]);
1897  else if (A.is_zero() && !B.is_zero())
1898  sum = B.frobSquared();
1899  else if (!A.is_zero() && B.is_zero())
1900  sum = A.frobSquared();
1901  /* sum is already zero if A.is_zero() && B.is_zero() */
1902  return sum;
1903  }
1904 
1905  template<class Treal, class Telement>
1909  const Matrix<Treal, Telement>& B) {
1910  assert(!A.is_empty());
1911  assert(!B.is_empty());
1912  if (A.nrows() != B.nrows() ||
1913  A.ncols() != B.ncols() ||
1914  A.nrows() != A.ncols())
1915  throw Failure("Matrix<Treal, Telement>::syFrobSquaredDiff: "
1916  "Incorrect matrix dimensions");
1917  Treal sum(0);
1918  if (!A.is_zero() && !B.is_zero()) {
1919  for (int col = 1; col < A.ncols(); col++)
1920  for (int row = 0; row < col; row++)
1921  sum += 2 * Telement::frobSquaredDiff(A(row, col), B(row, col));
1922  for (int rc = 0; rc < A.ncols(); rc++)
1923  sum += Telement::syFrobSquaredDiff(A(rc, rc),B(rc, rc));
1924  }
1925  else if (A.is_zero() && !B.is_zero())
1926  sum = B.syFrobSquared();
1927  else if (!A.is_zero() && B.is_zero())
1928  sum = A.syFrobSquared();
1929  /* sum is already zero if A.is_zero() && B.is_zero() */
1930  return sum;
1931  }
1932 
1933  template<class Treal, class Telement>
1935  trace() const {
1936  assert(!this->is_empty());
1937  if (this->nrows() != this->ncols())
1938  throw Failure("Matrix<Treal, Telement>::trace(): "
1939  "Matrix must be quadratic");
1940  if (this->is_zero())
1941  return 0;
1942  else {
1943  Treal sum = 0;
1944  for (int rc = 0; rc < this->ncols(); rc++)
1945  sum += (*this)(rc,rc).trace();
1946  return sum;
1947  }
1948  }
1949 
1950  template<class Treal, class Telement>
1953  const Matrix<Treal, Telement>& B) {
1954  assert(!A.is_empty());
1955  assert(!B.is_empty());
1956  if (A.nrows() != B.ncols() || A.ncols() != B.nrows())
1957  throw Failure("Matrix<Treal, Telement>::trace_ab: "
1958  "Wrong matrix dimensions for trace(A * B)");
1959  Treal tr = 0;
1960  if (!A.is_zero() && !B.is_zero())
1961  for (int rc = 0; rc < A.nrows(); rc++)
1962  for (int colA = 0; colA < A.ncols(); colA++)
1963  tr += Telement::trace_ab(A(rc,colA), B(colA, rc));
1964  return tr;
1965  }
1966 
1967  template<class Treal, class Telement>
1970  const Matrix<Treal, Telement>& B) {
1971  assert(!A.is_empty());
1972  assert(!B.is_empty());
1973  if (A.ncols() != B.ncols() || A.nrows() != B.nrows())
1974  throw Failure("Matrix<Treal, Telement>::trace_aTb: "
1975  "Wrong matrix dimensions for trace(A' * B)");
1976  Treal tr = 0;
1977  if (!A.is_zero() && !B.is_zero())
1978  for (int rc = 0; rc < A.ncols(); rc++)
1979  for (int rowB = 0; rowB < B.nrows(); rowB++)
1980  tr += Telement::trace_aTb(A(rowB,rc), B(rowB, rc));
1981  return tr;
1982  }
1983 
1984  template<class Treal, class Telement>
1987  const Matrix<Treal, Telement>& B) {
1988  assert(!A.is_empty());
1989  assert(!B.is_empty());
1990  if (A.nrows() != B.ncols() || A.ncols() != B.nrows() ||
1991  A.nrows() != A.ncols())
1992  throw Failure("Matrix<Treal, Telement>::sy_trace_ab: "
1993  "Wrong matrix dimensions for symmetric trace(A * B)");
1994  Treal tr = 0;
1995  if (!A.is_zero() && !B.is_zero()) {
1996  /* Diagonal first */
1997  for (int rc = 0; rc < A.nrows(); rc++)
1998  tr += Telement::sy_trace_ab(A(rc,rc), B(rc, rc));
1999  /* Using that trace of transpose is equal to that without transpose: */
2000  for (int colA = 1; colA < A.ncols(); colA++)
2001  for (int rowA = 0; rowA < colA; rowA++)
2002  tr += 2 * Telement::trace_aTb(A(rowA, colA), B(rowA, colA));
2003  }
2004  return tr;
2005  }
2006 
2007  template<class Treal, class Telement>
2009  add(const Treal alpha, /* B += alpha * A */
2010  const Matrix<Treal, Telement>& A,
2012  assert(!A.is_empty());
2013  assert(!B.is_empty());
2014  if (A.nrows() != B.nrows() || A.ncols() != B.ncols())
2015  throw Failure("Matrix<Treal, Telement>::add: "
2016  "Wrong matrix dimensions for addition");
2017  if (!A.is_zero() && alpha != 0) {
2018  if (B.is_zero())
2019  B.allocate();
2020  for (int ind = 0; ind < A.nElements(); ind++)
2021  Telement::add(alpha, A.elements[ind], B.elements[ind]);
2022  }
2023  }
2024 
2025  template<class Treal, class Telement>
2027  assign(Treal const alpha, /* *this = alpha * A */
2028  Matrix<Treal, Telement> const & A) {
2029  assert(!A.is_empty());
2030  if (this->is_empty()) {
2031  this->resetRows(A.rows);
2032  this->resetCols(A.cols);
2033  }
2034  *this = 0;
2036  add(alpha, A, *this);
2037  }
2038 
2039 
2040  /********** Help functions for thresholding */
2041 
2042 
2043  template<class Treal, class Telement>
2045  getFrobSqLowestLevel(std::vector<Treal> & frobsq) const {
2046  if (!this->is_zero())
2047  for (int ind = 0; ind < this->nElements(); ind++)
2048  this->elements[ind].getFrobSqLowestLevel(frobsq);
2049  }
2050 
2051  template<class Treal, class Telement>
2054  (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix) {
2055  assert(!this->is_empty());
2056  bool thisMatIsZero = true;
2057  if (ErrorMatrix) {
2058  assert(!ErrorMatrix->is_empty());
2059  bool EMatIsZero = true;
2060  if (!ErrorMatrix->is_zero() || !this->is_zero()) {
2061  if (ErrorMatrix->is_zero())
2062  ErrorMatrix->allocate();
2063  if (this->is_zero())
2064  this->allocate();
2065  for (int ind = 0; ind < this->nElements(); ind++) {
2066  this->elements[ind].
2067  frobThreshLowestLevel(threshold, &ErrorMatrix->elements[ind]);
2068  thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2069  EMatIsZero = EMatIsZero && ErrorMatrix->elements[ind].is_zero();
2070  }
2071  if (thisMatIsZero)
2072  this->clear();
2073  if (EMatIsZero)
2074  ErrorMatrix->clear();
2075  }
2076  }
2077  else
2078  if (!this->is_zero()) {
2079  for (int ind = 0; ind < this->nElements(); ind++) {
2080  this->elements[ind].frobThreshLowestLevel(threshold, 0);
2081  thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2082  }
2083  if (thisMatIsZero)
2084  this->clear();
2085  }
2086  }
2087 
2088  template<class Treal, class Telement>
2090  getFrobSqElementLevel(std::vector<Treal> & frobsq) const {
2091  if (!this->is_zero())
2092  for (int ind = 0; ind < this->nElements(); ind++)
2093  this->elements[ind].getFrobSqElementLevel(frobsq);
2094  }
2095 
2096  template<class Treal, class Telement>
2099  (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix) {
2100  assert(!this->is_empty());
2101  bool thisMatIsZero = true;
2102  if (ErrorMatrix) {
2103  assert(!ErrorMatrix->is_empty());
2104  bool EMatIsZero = true;
2105  if (!ErrorMatrix->is_zero() || !this->is_zero()) {
2106  if (ErrorMatrix->is_zero())
2107  ErrorMatrix->allocate();
2108  if (this->is_zero())
2109  this->allocate();
2110  for (int ind = 0; ind < this->nElements(); ind++) {
2111  this->elements[ind].
2112  frobThreshElementLevel(threshold, &ErrorMatrix->elements[ind]);
2113  thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2114  EMatIsZero = EMatIsZero && ErrorMatrix->elements[ind].is_zero();
2115  }
2116  if (thisMatIsZero)
2117  this->clear();
2118  if (EMatIsZero)
2119  ErrorMatrix->clear();
2120  }
2121  }
2122  else
2123  if (!this->is_zero()) {
2124  for (int ind = 0; ind < this->nElements(); ind++) {
2125  this->elements[ind].frobThreshElementLevel(threshold, 0);
2126  thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2127  }
2128  if (thisMatIsZero)
2129  this->clear();
2130  }
2131  }
2132 
2133 
2134 
2135  template<class Treal, class Telement>
2137  ( Matrix<Treal, Matrix<Treal, Telement> > const & A ) {
2138  if (!A.is_zero()) {
2139  if ( this->is_zero() )
2140  this->allocate();
2141  assert( this->nElements() == A.nElements() );
2142  for (int ind = 0; ind < this->nElements(); ind++)
2143  this->elements[ind].assignFrobNormsLowestLevel(A[ind]);
2144  }
2145  else
2146  this->clear();
2147  }
2148 
2149  template<class Treal, class Telement>
2151  ( Matrix<Treal, Matrix<Treal, Telement> > const & A ) {
2152  assert(!A.is_empty());
2153  if (A.nrows() != A.ncols())
2154  throw Failure("Matrix<Treal, Telement>::syAssignFrobNormsLowestLevel(...): "
2155  "Matrix must be have equal number of rows "
2156  "and cols to be symmetric");
2157  if (!A.is_zero()) {
2158  if ( this->is_zero() )
2159  this->allocate();
2160  assert( this->nElements() == A.nElements() );
2161  for (int col = 1; col < this->ncols(); col++)
2162  for (int row = 0; row < col; row++)
2163  (*this)(row, col).assignFrobNormsLowestLevel( A(row,col) );
2164  for (int rc = 0; rc < this->ncols(); rc++)
2165  (*this)(rc, rc).syAssignFrobNormsLowestLevel( A(rc,rc) );
2166  }
2167  else
2168  this->clear();
2169  }
2170 
2171  template<class Treal, class Telement>
2173  ( Matrix<Treal, Matrix<Treal, Telement> > const & A,
2174  Matrix<Treal, Matrix<Treal, Telement> > const & B ) {
2175  if ( A.is_zero() && B.is_zero() ) {
2176  // Both A and B are zero
2177  this->clear();
2178  return;
2179  }
2180  // At least one of A and B is nonzero
2181  if ( this->is_zero() )
2182  this->allocate();
2183  if ( !A.is_zero() && !B.is_zero() ) {
2184  // Both are nonzero
2185  assert( this->nElements() == A.nElements() );
2186  assert( this->nElements() == B.nElements() );
2187  for (int ind = 0; ind < this->nElements(); ind++)
2188  this->elements[ind].assignDiffFrobNormsLowestLevel( A[ind], B[ind] );
2189  return;
2190  }
2191  if ( !A.is_zero() ) {
2192  // A is nonzero
2193  this->assignFrobNormsLowestLevel( A );
2194  return;
2195  }
2196  if ( !B.is_zero() ) {
2197  // B is nonzero
2198  this->assignFrobNormsLowestLevel( B );
2199  return;
2200  }
2201  }
2202  template<class Treal, class Telement>
2204  ( Matrix<Treal, Matrix<Treal, Telement> > const & A,
2205  Matrix<Treal, Matrix<Treal, Telement> > const & B ) {
2206  if ( A.is_zero() && B.is_zero() ) {
2207  // Both A and B are zero
2208  this->clear();
2209  return;
2210  }
2211  // At least one of A and B is nonzero
2212  if ( this->is_zero() )
2213  this->allocate();
2214  if ( !A.is_zero() && !B.is_zero() ) {
2215  // Both are nonzero
2216  assert( this->nElements() == A.nElements() );
2217  assert( this->nElements() == B.nElements() );
2218  for (int col = 1; col < this->ncols(); col++)
2219  for (int row = 0; row < col; row++)
2220  (*this)(row, col).assignDiffFrobNormsLowestLevel( A(row,col), B(row,col) );
2221  for (int rc = 0; rc < this->ncols(); rc++)
2222  (*this)(rc, rc).syAssignDiffFrobNormsLowestLevel( A(rc,rc), B(rc,rc) );
2223  return;
2224  }
2225  if ( !A.is_zero() ) {
2226  // A is nonzero
2227  this->syAssignFrobNormsLowestLevel( A );
2228  return;
2229  }
2230  if ( !B.is_zero() ) {
2231  // B is nonzero
2232  this->syAssignFrobNormsLowestLevel( B );
2233  return;
2234  }
2235  }
2236 
2237 
2238 
2239  template<class Treal, class Telement>
2242  if ( this->is_zero() ) {
2243  A.clear();
2244  }
2245  else {
2246  assert( !A.is_zero() );
2247  assert( this->nElements() == A.nElements() );
2248  for (int ind = 0; ind < this->nElements(); ind++)
2249  this->elements[ind].truncateAccordingToSparsityPattern( A[ind] );
2250  }
2251  }
2252 
2253 
2254 
2255  /********** End of help functions for thresholding */
2256 
2257  /* C = beta * C + alpha * A * B where only the upper triangle of C is */
2258  /* referenced and updated */
2259  template<class Treal, class Telement>
2261  gemm_upper_tr_only(const bool tA, const bool tB, const Treal alpha,
2262  const Matrix<Treal, Telement>& A,
2263  const Matrix<Treal, Telement>& B,
2264  const Treal beta,
2266  assert(!A.is_empty());
2267  assert(!B.is_empty());
2268  if (C.is_empty()) {
2269  assert(beta == 0);
2270  if (tA)
2271  C.resetRows(A.cols);
2272  else
2273  C.resetRows(A.rows);
2274  if (tB)
2275  C.resetCols(B.rows);
2276  else
2277  C.resetCols(B.cols);
2278  }
2279  if ((A.is_zero() || B.is_zero() || alpha == 0) &&
2280  (C.is_zero() || beta == 0))
2281  C = 0;
2282  else {
2283  Treal beta_tmp = beta;
2284  if (C.is_zero()) {
2285  C.allocate();
2286  beta_tmp = 0;
2287  }
2288  if (!A.is_zero() && !B.is_zero() && alpha != 0) {
2289  if (!tA && !tB)
2290  if (A.ncols() == B.nrows() &&
2291  A.nrows() == C.nrows() &&
2292  B.ncols() == C.ncols()) {
2293  for (int col = 0; col < C.ncols(); col++) {
2294  Telement::gemm_upper_tr_only(tA, tB, alpha,
2295  A(col, 0), B(0, col),
2296  beta_tmp,
2297  C(col, col));
2298  for(int cola = 1; cola < A.ncols(); cola++)
2299  Telement::gemm_upper_tr_only(tA, tB, alpha,
2300  A(col, cola), B(cola, col),
2301  1.0,
2302  C(col,col));
2303  for (int row = 0; row < col; row++) {
2304  Telement::gemm(tA, tB, alpha,
2305  A(row, 0), B(0, col),
2306  beta_tmp,
2307  C(row,col));
2308  for(int cola = 1; cola < A.ncols(); cola++)
2309  Telement::gemm(tA, tB, alpha,
2310  A(row, cola), B(cola, col),
2311  1.0,
2312  C(row,col));
2313  }
2314  }
2315  }
2316  else
2317  throw Failure("Matrix<class Treal, class Telement>::"
2318  "gemm_upper_tr_only(bool, bool, Treal, "
2319  "Matrix<Treal, Telement>, "
2320  "Matrix<Treal, Telement>, "
2321  "Treal, Matrix<Treal, Telement>): "
2322  "Incorrect matrixdimensions for multiplication");
2323  else if (tA && !tB)
2324  if (A.nrows() == B.nrows() &&
2325  A.ncols() == C.nrows() &&
2326  B.ncols() == C.ncols()) {
2327  for (int col = 0; col < C.ncols(); col++) {
2328  Telement::gemm_upper_tr_only(tA, tB, alpha,
2329  A(0,col), B(0,col),
2330  beta_tmp,
2331  C(col,col));
2332  for(int cola = 1; cola < A.nrows(); cola++)
2333  Telement::gemm_upper_tr_only(tA, tB, alpha,
2334  A(cola, col), B(cola, col),
2335  1.0,
2336  C(col, col));
2337  for (int row = 0; row < col; row++) {
2338  Telement::gemm(tA, tB, alpha,
2339  A(0,row), B(0,col),
2340  beta_tmp,
2341  C(row,col));
2342  for(int cola = 1; cola < A.nrows(); cola++)
2343  Telement::gemm(tA, tB, alpha,
2344  A(cola, row), B(cola, col),
2345  1.0,
2346  C(row,col));
2347  }
2348  }
2349  }
2350  else
2351  throw Failure("Matrix<class Treal, class Telement>::"
2352  "gemm_upper_tr_only(bool, bool, Treal, "
2353  "Matrix<Treal, Telement>, "
2354  "Matrix<Treal, Telement>, "
2355  "Treal, Matrix<Treal, Telement>): "
2356  "Incorrect matrixdimensions for multiplication");
2357  else if (!tA && tB)
2358  if (A.ncols() == B.ncols() &&
2359  A.nrows() == C.nrows() &&
2360  B.nrows() == C.ncols()) {
2361  for (int col = 0; col < C.ncols(); col++) {
2362  Telement::gemm_upper_tr_only(tA, tB, alpha,
2363  A(col, 0), B(col, 0),
2364  beta_tmp,
2365  C(col,col));
2366  for(int cola = 1; cola < A.ncols(); cola++)
2367  Telement::gemm_upper_tr_only(tA, tB, alpha,
2368  A(col, cola), B(col, cola),
2369  1.0,
2370  C(col,col));
2371  for (int row = 0; row < col; row++) {
2372  Telement::gemm(tA, tB, alpha,
2373  A(row, 0), B(col, 0),
2374  beta_tmp,
2375  C(row,col));
2376  for(int cola = 1; cola < A.ncols(); cola++)
2377  Telement::gemm(tA, tB, alpha,
2378  A(row, cola), B(col, cola),
2379  1.0,
2380  C(row,col));
2381  }
2382  }
2383  }
2384  else
2385  throw Failure("Matrix<class Treal, class Telement>::"
2386  "gemm_upper_tr_only(bool, bool, Treal, "
2387  "Matrix<Treal, Telement>, "
2388  "Matrix<Treal, Telement>, "
2389  "Treal, Matrix<Treal, Telement>): "
2390  "Incorrect matrixdimensions for multiplication");
2391  else if (tA && tB)
2392  if (A.nrows() == B.ncols() &&
2393  A.ncols() == C.nrows() &&
2394  B.nrows() == C.ncols()) {
2395  for (int col = 0; col < C.ncols(); col++) {
2396  Telement::gemm_upper_tr_only(tA, tB, alpha,
2397  A(0, col), B(col, 0),
2398  beta_tmp,
2399  C(col,col));
2400  for(int cola = 1; cola < A.nrows(); cola++)
2401  Telement::gemm_upper_tr_only(tA, tB, alpha,
2402  A(cola, col), B(col, cola),
2403  1.0,
2404  C(col,col));
2405  for (int row = 0; row < col; row++) {
2406  Telement::gemm(tA, tB, alpha,
2407  A(0, row), B(col, 0),
2408  beta_tmp,
2409  C(row,col));
2410  for(int cola = 1; cola < A.nrows(); cola++)
2411  Telement::gemm(tA, tB, alpha,
2412  A(cola, row), B(col, cola),
2413  1.0,
2414  C(row,col));
2415  }
2416  }
2417  }
2418  else
2419  throw Failure("Matrix<class Treal, class Telement>::"
2420  "gemm_upper_tr_only(bool, bool, Treal, "
2421  "Matrix<Treal, Telement>, "
2422  "Matrix<Treal, Telement>, Treal, "
2423  "Matrix<Treal, Telement>): "
2424  "Incorrect matrixdimensions for multiplication");
2425  else throw Failure("Matrix<class Treal, class Telement>::"
2426  "gemm_upper_tr_only(bool, bool, Treal, "
2427  "Matrix<Treal, Telement>, "
2428  "Matrix<Treal, Telement>, Treal, "
2429  "Matrix<Treal, Telement>):"
2430  "Very strange error!!");
2431  }
2432  else
2433  C *= beta;
2434  }
2435  }
2436 
2437  /* A = alpha * A * Z or A = alpha * Z * A where A is symmetric, */
2438  /* Z is upper triangular and */
2439  /* only the upper triangle of the result is calculated */
2440  /* side == 'R' for A = alpha * A * Z */
2441  /* side == 'L' for A = alpha * Z * A */
2442  template<class Treal, class Telement>
2444  sytr_upper_tr_only(char const side, const Treal alpha,
2446  const Matrix<Treal, Telement>& Z) {
2447  assert(!A.is_empty());
2448  assert(!Z.is_empty());
2449  if (alpha != 0 && !A.is_zero() && !Z.is_zero()) {
2450  if (A.nrows() == A.ncols() &&
2451  Z.nrows() == Z.ncols() &&
2452  A.nrows() == Z.nrows()) {
2453  if (side == 'R') {
2454  /* Last column must be calculated first */
2455  for (int col = A.ncols() - 1; col >= 0; col--) {
2456  // A(col, col) = alpha * A(col, col) * Z(col, col)
2457  Telement::sytr_upper_tr_only(side, alpha,
2458  A(col, col), Z(col, col));
2459  for (int ind = 0; ind < col; ind++) {
2460  // A(col, col) += alpha * A(ind, col)' * Z(ind, col)
2461  Telement::gemm_upper_tr_only(true, false, alpha, A(ind, col),
2462  Z(ind, col), 1.0, A(col, col));
2463  }
2464  /* Last row must be calculated first */
2465  for (int row = col - 1; row >= 0; row--) {
2466  // A(row, col) = alpha * A(row, col) * Z(col, col);
2467  Telement::trmm(side, 'U', false,
2468  alpha, Z(col, col), A(row, col));
2469  // A(row, col) += alpha * A(row, row) * Z(row, col);
2470  Telement::symm('L', 'U', alpha, A(row, row), Z(row, col),
2471  1.0, A(row, col));
2472  for (int ind = 0; ind < row; ind++) {
2473  // A(row, col) += alpha * A(ind, row)' * Z(ind, col);
2474  Telement::gemm(true, false, alpha, A(ind, row), Z(ind, col),
2475  1.0, A(row, col));
2476  }
2477  for (int ind = row + 1; ind < col; ind++) {
2478  // A(row, col) += alpha * A(row, ind) * Z(ind, col);
2479  Telement::gemm(false, false, alpha, A(row, ind), Z(ind, col),
2480  1.0, A(row, col));
2481  }
2482  }
2483  }
2484  }
2485  else { /* side == 'L' */
2486  assert(side == 'L');
2487  /* First column must be calculated first */
2488  for (int col = 0; col < A.ncols(); col++) {
2489  /* First row must be calculated first */
2490  for (int row = 0; row < col; row++) {
2491  // A(row, col) = alpha * Z(row, row) * A(row, col)
2492  Telement::trmm(side, 'U', false, alpha,
2493  Z(row, row), A(row, col));
2494  // A(row, col) += alpha * Z(row, col) * A(col, col)
2495  Telement::symm('R', 'U', alpha, A(col, col), Z(row, col),
2496  1.0, A(row, col));
2497  for (int ind = row + 1; ind < col; ind++)
2498  // A(row, col) += alpha * Z(row, ind) * A(ind, col)
2499  Telement::gemm(false, false, alpha, Z(row, ind), A(ind, col),
2500  1.0, A(row, col));
2501  for (int ind = col + 1; ind < A.nrows(); ind++)
2502  // A(row, col) += alpha * Z(row, ind) * A(col, ind)'
2503  Telement::gemm(false, true, alpha, Z(row, ind), A(col, ind),
2504  1.0, A(row, col));
2505  }
2506  Telement::sytr_upper_tr_only(side, alpha,
2507  A(col, col), Z(col, col));
2508  for (int ind = col + 1; ind < A.ncols(); ind++)
2509  Telement::gemm_upper_tr_only(false, true, alpha, Z(col, ind),
2510  A(col, ind), 1.0, A(col, col));
2511  }
2512  }
2513  }
2514  else
2515  throw Failure("Matrix<class Treal, class Telement>::"
2516  "sytr_upper_tr_only: Incorrect matrix dimensions "
2517  "for symmetric-triangular multiplication");
2518  }
2519  else
2520  A = 0;
2521  }
2522 
2523  /* The result B is assumed to be symmetric, i.e. only upper triangle */
2524  /* calculated and hence only upper triangle of input B is used. */
2525  template<class Treal, class Telement>
2527  trmm_upper_tr_only(const char side, const char uplo,
2528  const bool tA, const Treal alpha,
2529  const Matrix<Treal, Telement>& A,
2531  assert(!A.is_empty());
2532  assert(!B.is_empty());
2533  if (alpha != 0 && !A.is_zero() && !B.is_zero())
2534  if (((side == 'R' && B.ncols() == A.nrows()) ||
2535  (side == 'L' && A.ncols() == B.nrows())) &&
2536  A.nrows() == A.ncols())
2537  if (uplo == 'U')
2538  if (!tA) {
2539  throw Failure("Matrix<Treal, class Telement>::"
2540  "trmm_upper_tr_only : "
2541  "not possible for upper triangular not transposed "
2542  "matrices A since lower triangle of B is needed");
2543  } /* end if (tA == false) */
2544  else {
2545  assert(tA == true);
2546  if (side == 'R') {
2547  /* First column must be calculated first */
2548  for (int col = 0; col < B.ncols(); col++) {
2549  Telement::trmm_upper_tr_only(side, uplo, tA, alpha,
2550  A(col,col), B(col,col));
2551  for (int ind = col + 1; ind < A.ncols(); ind++)
2552  Telement::gemm_upper_tr_only(false, tA, alpha,
2553  B(col,ind), A(col,ind),
2554  1.0, B(col,col));
2555  for (int row = 0; row < col; row++) {
2556  Telement::trmm(side, uplo, tA, alpha,
2557  A(col,col), B(row,col));
2558  for (int ind = col + 1; ind < A.ncols(); ind++)
2559  Telement::gemm(false, tA, alpha,
2560  B(row,ind), A(col,ind),
2561  1.0, B(row,col));
2562  }
2563  }
2564  } /* end if (side == 'R')*/
2565  else {
2566  assert(side == 'L');
2567  /* Last row must be calculated first */
2568  for (int row = B.nrows() - 1; row >= 0; row--) {
2569  Telement::trmm_upper_tr_only(side, uplo, tA, alpha,
2570  A(row,row), B(row,row));
2571  for (int ind = 0; ind < row; ind++)
2572  Telement::gemm_upper_tr_only(tA, false, alpha,
2573  A(ind,row), B(ind,row),
2574  1.0, B(row,row));
2575  for (int col = row + 1; col < B.ncols(); col++) {
2576  Telement::trmm(side, uplo, tA, alpha,
2577  A(row,row), B(row,col));
2578  for (int ind = 0; ind < row; ind++)
2579  Telement::gemm(tA, false, alpha,
2580  A(ind,row), B(ind,col),
2581  1.0, B(row,col));
2582  }
2583  }
2584  } /* end else (side == 'L')*/
2585  } /* end else (tA == true)*/
2586  else
2587  throw Failure("Matrix<class Treal, class Telement>::"
2588  "trmm_upper_tr_only not implemented for lower "
2589  "triangular matrices");
2590  else
2591  throw Failure("Matrix<class Treal, class Telement>::"
2592  "trmm_upper_tr_only: Incorrect matrix dimensions "
2593  "for multiplication");
2594  else
2595  B = 0;
2596  }
2597 
2598  /* A = Z' * A * Z or A = Z * A * Z' */
2599  /* where A is symmetric and Z is (nonunit) upper triangular */
2600  /* side == 'R' for A = Z' * A * Z */
2601  /* side == 'L' for A = Z * A * Z' */
2602  template<class Treal, class Telement>
2604  trsytriplemm(char const side,
2605  const Matrix<Treal, Telement>& Z,
2607  if (side == 'R') {
2609  sytr_upper_tr_only('R', 1.0, A, Z);
2611  trmm_upper_tr_only('L', 'U', true, 1.0, Z, A);
2612  }
2613  else {
2614  assert(side == 'L');
2616  sytr_upper_tr_only('L', 1.0, A, Z);
2618  trmm_upper_tr_only('R', 'U', true, 1.0, Z, A);
2619  }
2620  }
2621 
2622  template<class Treal, class Telement>
2624  frob_squared_thresh(Treal const threshold,
2625  Matrix<Treal, Telement> * ErrorMatrix) {
2626  assert(!this->is_empty());
2627  if (ErrorMatrix && ErrorMatrix->is_empty()) {
2628  ErrorMatrix->resetRows(this->rows);
2629  ErrorMatrix->resetCols(this->cols);
2630  }
2631  assert(threshold >= (Treal)0.0);
2632  if (threshold == (Treal)0.0)
2633  return 0;
2634 
2635  std::vector<Treal> frobsq_norms;
2636  getFrobSqLowestLevel(frobsq_norms);
2637  /* Sort in ascending order */
2638  std::sort(frobsq_norms.begin(), frobsq_norms.end());
2639  int lastRemoved = -1;
2640  Treal frobsqSum = 0;
2641  int nnzBlocks = frobsq_norms.size();
2642  while(lastRemoved < nnzBlocks - 1 && frobsqSum < threshold) {
2643  ++lastRemoved;
2644  frobsqSum += frobsq_norms[lastRemoved];
2645  }
2646 
2647  /* Check if entire matrix will be removed */
2648  if (lastRemoved == nnzBlocks - 1 && frobsqSum < threshold) {
2649  if (ErrorMatrix)
2650  Matrix<Treal, Telement>::swap(*this, *ErrorMatrix);
2651  else
2652  *this = 0;
2653  }
2654  else {
2655  frobsqSum -= frobsq_norms[lastRemoved];
2656  --lastRemoved;
2657  while(lastRemoved > -1 && frobsq_norms[lastRemoved] ==
2658  frobsq_norms[lastRemoved + 1]) {
2659  frobsqSum -= frobsq_norms[lastRemoved];
2660  --lastRemoved;
2661  }
2662  if (lastRemoved > -1) {
2663  Treal threshLowestLevel =
2664  (frobsq_norms[lastRemoved + 1] +
2665  frobsq_norms[lastRemoved]) / 2;
2666  this->frobThreshLowestLevel(threshLowestLevel, ErrorMatrix);
2667  }
2668  }
2669  return frobsqSum;
2670  }
2671 
2672  template<class Treal, class Telement>
2676  const Treal threshold, const side looking,
2677  const inchversion version) {
2678  assert(!A.is_empty());
2679  if (A.nrows() != A.ncols())
2680  throw Failure("Matrix<Treal, Telement>::sy_inch: "
2681  "Matrix must be quadratic!");
2682  if (A.is_zero())
2683  throw Failure("Matrix<Treal>::sy_inch: Matrix is zero! "
2684  "Not possible to compute inverse cholesky.");
2685  if (version == stable) /* STABILIZED */
2686  throw Failure("Matrix<Treal>::sy_inch: Only unstable "
2687  "version of sy_inch implemented.");
2688  Treal myThresh = 0;
2689  if (threshold != 0)
2690  myThresh = threshold / (Z.ncols() * Z.nrows());
2691  Z.resetRows(A.rows);
2692  Z.resetCols(A.cols);
2693  Z.allocate();
2694 
2695  if (looking == left) { /* LEFT-LOOKING INCH */
2696  if (threshold != 0)
2697  throw Failure("Matrix<Treal, Telement>::syInch: "
2698  "Thresholding not implemented for left-looking inch.");
2699  /* Left and unstable */
2700  Telement::syInch(A(0,0), Z(0,0), threshold, looking, version);
2701  Telement Ptmp;//, tmp;
2702  for (int i = 1; i < Z.ncols(); i++) {
2703  for (int j = 0; j < i; j++) {
2704  /* Z(k,i) is nonzero for k = 0, 1, ...,j - 2, j - 1, i */
2705  /* and Z(i,i) = I (yes it is i ^) */
2706  Ptmp = A(j,i); /* (Z(i,i) == I) */
2707  for (int k = 0; k < j; k++) /* Ptmp = A(j,:) * Z(:,i) */
2708  Telement::gemm(true, false, 1.0, /* SYMMETRY USED */
2709  A(k,j), Z(k,i), 1.0, Ptmp);
2710  Telement::trmm('L', 'U', true, 1.0, Z(j,j), Ptmp);
2711 
2712  for (int k = 0; k < j; k++) /* Z(:,i) -= Z(:,j) * Ptmp */
2713  Telement::gemm(false, false, -1.0,
2714  Z(k,j), Ptmp, 1.0, Z(k,i));
2715  /* Z(j,j) is triangular: */
2716  Telement::trmm('L', 'U', false, -1.0, Z(j,j), Ptmp);
2717  Telement::add(1.0, Ptmp, Z(j,i));
2718  }
2719  Ptmp = A(i,i); /* Z(i,i) == I */
2720  for (int k = 0; k < i; k++) /* SYMMETRY USED */
2721  Telement::gemm_upper_tr_only(true, false, 1.0,
2722  A(k,i), Z(k,i), 1.0, Ptmp);
2723  /* Z(i,i) == I !! */
2724  /* Z(:,i) *= INCH(Ptmp) */
2725  Telement::syInch(Ptmp, Z(i,i), threshold, looking, version);
2726  for (int k = 0; k < i; k++) {
2727  Telement::trmm('R', 'U', false, 1.0, Z(i,i), Z(k,i));
2728  /* INCH(Ptmp) is upper triangular*/
2729  }
2730  }
2731  } /* end if left-looking inch */
2732  else { /* RIGHT-LOOKING INCH */
2733  assert(looking == right); /* right and unstable */
2734  Telement Ptmp;
2735  Treal newThresh = 0;
2736  if (myThresh != 0 && Z.ncols() > 1)
2737  newThresh = myThresh * 0.0001;
2738  else
2739  newThresh = myThresh;
2740 
2741  for (int i = 0; i < Z.ncols(); i++) {
2742  /* Ptmp = A(i,:) * Z(:,i) */
2743  Ptmp = A(i,i); /* Z(i,i) == I */
2744  for (int k = 0; k < i; k++)
2745  Telement::gemm_upper_tr_only(true, false, 1.0, /* SYMMETRY USED */
2746  A(k,i), Z(k,i), 1.0, Ptmp);
2747 
2748  /* Z(:,i) *= INCH(Ptmp) */
2749  Telement::syInch(Ptmp, Z(i,i), newThresh, looking, version);
2750  for (int k = 0; k < i; k++) {
2751  Telement::trmm('R', 'U', false, 1.0, Z(i,i), Z(k,i));
2752  /* INCH(Ptmp) is upper triangular */
2753  }
2754 
2755  for (int j = i + 1; j < Z.ncols(); j++) {
2756  /* Compute Ptmp = Z(i,i)^T * A(i,:) * Z(:,j) */
2757  /* Z(k,j) is nonzero for k = 0, 1, ...,i - 2, i - 1, j */
2758  /* and Z(j,j) = I */
2759  Ptmp = A(i,j); /* (Z(j,j) == I) */
2760  for (int k = 0; k < i; k++) /* Ptmp = A(i,:) * Z(:,j) */
2761  Telement::gemm(true, false, 1.0, /* SYMMETRY USED */
2762  A(k,i), Z(k,j), 1.0, Ptmp);
2763  Telement::trmm('L', 'U', true, 1.0, Z(i,i), Ptmp);
2764 
2765  for (int k = 0; k < i; k++) /* Z(:,j) -= Z(:,i) * Ptmp */
2766  Telement::gemm(false, false, -1.0,
2767  Z(k,i), Ptmp, 1.0, Z(k,j));
2768  /* Z(i,i) is triangular: */
2769  Telement::trmm('L', 'U', false, -1.0, Z(i,i), Ptmp);
2770  Telement::add(1.0, Ptmp, Z(i,j));
2771  } /* end for j */
2772 
2773  /* Truncation starts here */
2774  if (threshold != 0) {
2775  for (int k = 0; k < i; k++)
2776  Z(k,i).frob_thresh(myThresh);
2777  }
2778  } /* end for i */
2779  } /* end else right-looking inch */
2780  }
2781 
2782  template<class Treal, class Telement>
2784  gersgorin(Treal& lmin, Treal& lmax) const {
2785  assert(!this->is_empty());
2786  if (this->nScalarsRows() == this->nScalarsCols()) {
2787  int size = this->nScalarsCols();
2788  Treal* colsums = new Treal[size];
2789  Treal* diag = new Treal[size];
2790  for (int ind = 0; ind < size; ind++)
2791  colsums[ind] = 0;
2792  this->add_abs_col_sums(colsums);
2793  this->get_diagonal(diag);
2794  Treal tmp1 = colsums[0] - template_blas_fabs(diag[0]);
2795  Treal tmp2;
2796  lmin = diag[0] - tmp1;
2797  lmax = diag[0] + tmp1;
2798  for (int col = 1; col < size; col++) {
2799  tmp1 = colsums[col] - template_blas_fabs(diag[col]);
2800  tmp2 = diag[col] - tmp1;
2801  lmin = (tmp2 < lmin ? tmp2 : lmin);
2802  tmp2 = diag[col] + tmp1;
2803  lmax = (tmp2 > lmax ? tmp2 : lmax);
2804  }
2805  delete[] diag;
2806  delete[] colsums;
2807  }
2808  else
2809  throw Failure("Matrix<Treal, Telement, int>::gersgorin(Treal, Treal): "
2810  "Matrix must be quadratic");
2811  }
2812 
2813 
2814  template<class Treal, class Telement>
2816  add_abs_col_sums(Treal* sums) const {
2817  assert(sums);
2818  if (!this->is_zero()) {
2819  int offset = 0;
2820  for (int col = 0; col < this->ncols(); col++) {
2821  for (int row = 0; row < this->nrows(); row++) {
2822  (*this)(row,col).add_abs_col_sums(&sums[offset]);
2823  }
2824  offset += (*this)(0,col).nScalarsCols();
2825  }
2826  }
2827  }
2828 
2829  template<class Treal, class Telement>
2831  get_diagonal(Treal* diag) const {
2832  assert(diag);
2833  assert(this->nrows() == this->ncols());
2834  if (this->is_zero())
2835  for (int rc = 0; rc < this->nScalarsCols(); rc++)
2836  diag[rc] = 0;
2837  else {
2838  int offset = 0;
2839  for (int rc = 0; rc < this->ncols(); rc++) {
2840  (*this)(rc,rc).get_diagonal(&diag[offset]);
2841  offset += (*this)(rc,rc).nScalarsCols();
2842  }
2843  }
2844  }
2845 
2846  template<class Treal, class Telement>
2848  assert(!this->is_empty());
2849  size_t sum = 0;
2850  if (this->is_zero())
2851  return (size_t)0;
2852  for (int ind = 0; ind < this->nElements(); ind++)
2853  sum += this->elements[ind].memory_usage();
2854  return sum;
2855  }
2856 
2857  template<class Treal, class Telement>
2859  size_t sum = 0;
2860  if (!this->is_zero()) {
2861  for (int ind = 0; ind < this->nElements(); ind++)
2862  sum += this->elements[ind].nnz();
2863  }
2864  return sum;
2865  }
2866  template<class Treal, class Telement>
2868  size_t sum = 0;
2869  if (!this->is_zero()) {
2870  /* Above diagonal */
2871  for (int col = 1; col < this->ncols(); col++)
2872  for (int row = 0; row < col; row++)
2873  sum += (*this)(row, col).nnz();
2874  /* Below diagonal */
2875  sum *= 2;
2876  /* Diagonal */
2877  for (int rc = 0; rc < this->nrows(); rc++)
2878  sum += (*this)(rc,rc).sy_nnz();
2879  }
2880  return sum;
2881  }
2882 
2883  template<class Treal, class Telement>
2885  size_t sum = 0;
2886  if (!this->is_zero()) {
2887  /* Above diagonal */
2888  for (int col = 1; col < this->ncols(); col++)
2889  for (int row = 0; row < col; row++)
2890  sum += (*this)(row, col).nvalues();
2891  /* Diagonal */
2892  for (int rc = 0; rc < this->nrows(); rc++)
2893  sum += (*this)(rc,rc).sy_nvalues();
2894  }
2895  return sum;
2896  }
2897 
2898 
2899 
2900 
2901 
2902  /***************************************************************************/
2903  /***************************************************************************/
2904  /* Specialization for Telement = Treal */
2905  /***************************************************************************/
2906  /***************************************************************************/
2907 
2908  template<class Treal>
2909  class Matrix<Treal>: public MatrixHierarchicBase<Treal> {
2910 
2911  public:
2913  friend class Vector<Treal, Treal>;
2914 
2916  :MatrixHierarchicBase<Treal>() {
2917  }
2918  /* Matrix(const Atomblock<Treal>& row_atoms,
2919  const Atomblock<Treal>& col_atoms,
2920  bool z = true, int nr = 0, int nc = 0, char tr = 'N')
2921  :MatrixHierarchicBase<Treal>(row_atoms, col_atoms, z, nr, nc,tr) {}
2922  */
2923 
2924  void allocate() {
2925  assert(!this->is_empty());
2926  assert(this->is_zero());
2927  this->elements = allocateElements<Treal>(this->nElements());
2928  for (int ind = 0; ind < this->nElements(); ++ind)
2929  this->elements[ind] = 0;
2930  }
2931 
2932  /* Full matrix assigns etc */
2933  void assignFromFull(std::vector<Treal> const & fullMat);
2934  void fullMatrix(std::vector<Treal> & fullMat) const;
2935  void syFullMatrix(std::vector<Treal> & fullMat) const;
2936  void syUpTriFullMatrix(std::vector<Treal> & fullMat) const;
2937 
2938  /* Sparse matrix assigns etc */
2939  void assignFromSparse(std::vector<int> const & rowind,
2940  std::vector<int> const & colind,
2941  std::vector<Treal> const & values);
2942  void assignFromSparse(std::vector<int> const & rowind,
2943  std::vector<int> const & colind,
2944  std::vector<Treal> const & values,
2945  std::vector<int> const & indexes);
2946 
2947  /* Adds values (+=) to elements */
2948  void addValues(std::vector<int> const & rowind,
2949  std::vector<int> const & colind,
2950  std::vector<Treal> const & values);
2951  void addValues(std::vector<int> const & rowind,
2952  std::vector<int> const & colind,
2953  std::vector<Treal> const & values,
2954  std::vector<int> const & indexes);
2955 
2956  void syAssignFromSparse(std::vector<int> const & rowind,
2957  std::vector<int> const & colind,
2958  std::vector<Treal> const & values);
2959 
2960  void syAddValues(std::vector<int> const & rowind,
2961  std::vector<int> const & colind,
2962  std::vector<Treal> const & values);
2963 
2964  void getValues(std::vector<int> const & rowind,
2965  std::vector<int> const & colind,
2966  std::vector<Treal> & values) const;
2967  void getValues(std::vector<int> const & rowind,
2968  std::vector<int> const & colind,
2969  std::vector<Treal> & values,
2970  std::vector<int> const & indexes) const;
2971  void syGetValues(std::vector<int> const & rowind,
2972  std::vector<int> const & colind,
2973  std::vector<Treal> & values) const;
2974 
2975  void getAllValues(std::vector<int> & rowind,
2976  std::vector<int> & colind,
2977  std::vector<Treal> & values) const;
2978  void syGetAllValues(std::vector<int> & rowind,
2979  std::vector<int> & colind,
2980  std::vector<Treal> & values) const;
2981 
2982  Matrix<Treal>&
2983  operator=(const Matrix<Treal>& mat) {
2985  return *this;
2986  }
2987 
2988  void clear();
2990  clear();
2991  }
2992 
2993  void writeToFile(std::ofstream & file) const;
2994  void readFromFile(std::ifstream & file);
2995 
2996  void random();
2997  void syRandom();
2998  void randomZeroStructure(Treal probabilityBeingZero);
2999  void syRandomZeroStructure(Treal probabilityBeingZero);
3000  template<typename TRule>
3001  void setElementsByRule(TRule & rule);
3002  template<typename TRule>
3003  void sySetElementsByRule(TRule & rule);
3004 
3005 
3006  void addIdentity(Treal alpha); /* C += alpha * I */
3007 
3008  static void transpose(Matrix<Treal> const & A,
3009  Matrix<Treal> & AT);
3010 
3011  void symToNosym();
3012  void nosymToSym();
3013 
3014  /* Set matrix to zero (k = 0) or identity (k = 1) */
3015  Matrix<Treal>& operator=(int const k);
3016 
3017  Matrix<Treal>& operator*=(const Treal alpha);
3018 
3019  static void gemm(const bool tA, const bool tB, const Treal alpha,
3020  const Matrix<Treal>& A,
3021  const Matrix<Treal>& B,
3022  const Treal beta,
3023  Matrix<Treal>& C);
3024  static void symm(const char side, const char uplo, const Treal alpha,
3025  const Matrix<Treal>& A,
3026  const Matrix<Treal>& B,
3027  const Treal beta,
3028  Matrix<Treal>& C);
3029  static void syrk(const char uplo, const bool tA, const Treal alpha,
3030  const Matrix<Treal>& A,
3031  const Treal beta,
3032  Matrix<Treal>& C);
3033  /* C = beta * C + alpha * A * A where A and C are symmetric */
3034  static void sysq(const char uplo, const Treal alpha,
3035  const Matrix<Treal>& A,
3036  const Treal beta,
3037  Matrix<Treal>& C);
3038  /* C = alpha * A * B + beta * C where A and B are symmetric */
3039  static void ssmm(const Treal alpha,
3040  const Matrix<Treal>& A,
3041  const Matrix<Treal>& B,
3042  const Treal beta,
3043  Matrix<Treal>& C);
3044  /* C = alpha * A * B + beta * C where A and B are symmetric
3045  * and only the upper triangle of C is computed.
3046  */
3047  static void ssmm_upper_tr_only(const Treal alpha,
3048  const Matrix<Treal>& A,
3049  const Matrix<Treal>& B,
3050  const Treal beta,
3051  Matrix<Treal>& C);
3052 
3053  static void trmm(const char side, const char uplo, const bool tA,
3054  const Treal alpha,
3055  const Matrix<Treal>& A,
3056  Matrix<Treal>& B);
3057 
3058  /* Frobenius norms */
3059  inline Treal frob() const {return template_blas_sqrt(frobSquared());}
3060  Treal frobSquared() const;
3061  inline Treal syFrob() const {
3062  return template_blas_sqrt(this->syFrobSquared());
3063  }
3064  Treal syFrobSquared() const;
3065 
3066  inline static Treal frobDiff
3067  (const Matrix<Treal>& A,
3068  const Matrix<Treal>& B) {
3069  return template_blas_sqrt(frobSquaredDiff(A, B));
3070  }
3071  static Treal frobSquaredDiff
3072  (const Matrix<Treal>& A,
3073  const Matrix<Treal>& B);
3074 
3075  inline static Treal syFrobDiff
3076  (const Matrix<Treal>& A,
3077  const Matrix<Treal>& B) {
3078  return template_blas_sqrt(syFrobSquaredDiff(A, B));
3079  }
3080  static Treal syFrobSquaredDiff
3081  (const Matrix<Treal>& A,
3082  const Matrix<Treal>& B);
3083 
3084  Treal trace() const;
3085  static Treal trace_ab(const Matrix<Treal>& A,
3086  const Matrix<Treal>& B);
3087  static Treal trace_aTb(const Matrix<Treal>& A,
3088  const Matrix<Treal>& B);
3089  static Treal sy_trace_ab(const Matrix<Treal>& A,
3090  const Matrix<Treal>& B);
3091 
3092  static void add(const Treal alpha, /* B += alpha * A */
3093  const Matrix<Treal>& A,
3094  Matrix<Treal>& B);
3095  void assign(Treal const alpha, /* *this = alpha * A */
3096  Matrix<Treal> const & A);
3097 
3098 
3099  /********** Help functions for thresholding */
3100  // int getnnzBlocksLowestLevel() const;
3101  void getFrobSqLowestLevel(std::vector<Treal> & frobsq) const;
3103  (Treal const threshold, Matrix<Treal> * ErrorMatrix);
3104 
3105  void getFrobSqElementLevel(std::vector<Treal> & frobsq) const;
3107  (Treal const threshold, Matrix<Treal> * ErrorMatrix);
3108 
3109 
3110 #if 0
3111  inline void frobThreshLowestLevel
3112  (Treal const threshold,
3113  Matrix<Treal> * ErrorMatrix) {
3114  bool a,b;
3115  frobThreshLowestLevel(threshold, ErrorMatrix, a, b);
3116  }
3117 #endif
3118 
3120  ( Matrix<Treal, Matrix<Treal> > const & A ) {
3121  if (!A.is_zero()) {
3122  if ( this->is_zero() )
3123  this->allocate();
3124  assert( this->nElements() == A.nElements() );
3125  for (int ind = 0; ind < this->nElements(); ind++)
3126  this->elements[ind] = A[ind].frob();
3127  }
3128  else
3129  this->clear();
3130  }
3131 
3133  if (!A.is_zero()) {
3134  if ( this->is_zero() )
3135  this->allocate();
3136  assert( this->nElements() == A.nElements() );
3137  for (int col = 1; col < this->ncols(); col++)
3138  for (int row = 0; row < col; row++)
3139  (*this)(row,col) = A(row,col).frob();
3140  for (int rc = 0; rc < this->nrows(); rc++)
3141  (*this)(rc,rc) = A(rc,rc).syFrob();
3142  }
3143  else
3144  this->clear();
3145  }
3146 
3148  Matrix<Treal, Matrix<Treal> > const & B ) {
3149  if ( A.is_zero() && B.is_zero() ) {
3150  // Both A and B are zero
3151  this->clear();
3152  return;
3153  }
3154  // At least one of A and B is nonzero
3155  if ( this->is_zero() )
3156  this->allocate();
3157  if ( !A.is_zero() && !B.is_zero() ) {
3158  // Both are nonzero
3159  assert( this->nElements() == A.nElements() );
3160  assert( this->nElements() == B.nElements() );
3161  for (int ind = 0; ind < this->nElements(); ind++) {
3162  Matrix<Treal> Diff(A[ind]);
3163  Matrix<Treal>::add( -1.0, B[ind], Diff );
3164  this->elements[ind] = Diff.frob();
3165  }
3166  return;
3167  }
3168  if ( !A.is_zero() ) {
3169  // A is nonzero
3170  this->assignFrobNormsLowestLevel( A );
3171  return;
3172  }
3173  if ( !B.is_zero() ) {
3174  // B is nonzero
3175  this->assignFrobNormsLowestLevel( B );
3176  return;
3177  }
3178  }
3180  Matrix<Treal, Matrix<Treal> > const & B ) {
3181  if ( A.is_zero() && B.is_zero() ) {
3182  // Both A and B are zero
3183  this->clear();
3184  return;
3185  }
3186  // At least one of A and B is nonzero
3187  if ( this->is_zero() )
3188  this->allocate();
3189  if ( !A.is_zero() && !B.is_zero() ) {
3190  // Both are nonzero
3191  assert( this->nElements() == A.nElements() );
3192  assert( this->nElements() == B.nElements() );
3193  for (int col = 1; col < this->ncols(); col++)
3194  for (int row = 0; row < col; row++) {
3195  Matrix<Treal> Diff(A(row,col));
3196  Matrix<Treal>::add( -1.0, B(row,col), Diff );
3197  (*this)(row, col) = Diff.frob();
3198  }
3199  for (int rc = 0; rc < this->ncols(); rc++) {
3200  Matrix<Treal> Diff( A(rc,rc) );
3201  Matrix<Treal>::add( -1.0, B(rc,rc), Diff );
3202  (*this)(rc, rc) = Diff.syFrob();
3203  }
3204  return;
3205  }
3206  if ( !A.is_zero() ) {
3207  // A is nonzero
3208  this->syAssignFrobNormsLowestLevel( A );
3209  return;
3210  }
3211  if ( !B.is_zero() ) {
3212  // B is nonzero
3213  this->syAssignFrobNormsLowestLevel( B );
3214  return;
3215  }
3216  }
3217 
3218 
3220  if ( this->is_zero() )
3221  A.clear();
3222  else {
3223  assert( !A.is_zero() );
3224  assert( this->nElements() == A.nElements() );
3225  for (int ind = 0; ind < this->nElements(); ind++)
3226  if (this->elements[ind] == 0)
3227  A[ind].clear();
3228  }
3229  }
3230 
3231 
3232  /********** End of help functions for thresholding */
3233 
3234  static void gemm_upper_tr_only(const bool tA, const bool tB,
3235  const Treal alpha,
3236  const Matrix<Treal>& A,
3237  const Matrix<Treal>& B,
3238  const Treal beta,
3239  Matrix<Treal>& C);
3240  static void sytr_upper_tr_only(char const side, const Treal alpha,
3241  Matrix<Treal>& A,
3242  const Matrix<Treal>& Z);
3243  static void trmm_upper_tr_only(const char side, const char uplo,
3244  const bool tA, const Treal alpha,
3245  const Matrix<Treal>& A,
3246  Matrix<Treal>& B);
3247  static void trsytriplemm(char const side,
3248  const Matrix<Treal>& Z,
3249  Matrix<Treal>& A);
3250 
3251  inline Treal frob_thresh(Treal const threshold,
3252  Matrix<Treal> * ErrorMatrix = 0) {
3253  return template_blas_sqrt
3254  (frob_squared_thresh(threshold * threshold, ErrorMatrix));
3255  }
3256  /* Returns the Frobenius norm of the introduced error */
3257 
3258  Treal frob_squared_thresh(Treal const threshold,
3259  Matrix<Treal> * ErrorMatrix = 0);
3260 
3261 
3262  static void inch(const Matrix<Treal>& A,
3263  Matrix<Treal>& Z,
3264  const Treal threshold = 0,
3265  const side looking = left,
3266  const inchversion version = unstable);
3267  static void syInch(const Matrix<Treal>& A,
3268  Matrix<Treal>& Z,
3269  const Treal threshold = 0,
3270  const side looking = left,
3271  const inchversion version = unstable) {
3272  inch(A, Z, threshold, looking, version);
3273  }
3274 
3275  void gersgorin(Treal& lmin, Treal& lmax) const;
3276  void sy_gersgorin(Treal& lmin, Treal& lmax) const {
3277  Matrix<Treal> tmp = (*this);
3278  tmp.symToNosym();
3279  tmp.gersgorin(lmin, lmax);
3280  return;
3281  }
3282  void add_abs_col_sums(Treal* abscolsums) const;
3283  void get_diagonal(Treal* diag) const; /* Copy diagonal to the diag array */
3284 
3285  inline size_t memory_usage() const { /* Returns memory usage in bytes */
3286  assert(!this->is_empty());
3287  if (this->is_zero())
3288  return (size_t)0;
3289  else
3290  return (size_t)this->nElements() * sizeof(Treal);
3291  }
3292 
3293  inline size_t nnz() const {
3294  if (this->is_zero())
3295  return 0;
3296  else
3297  return this->nElements();
3298  }
3299  inline size_t sy_nnz() const {
3300  if (this->is_zero())
3301  return 0;
3302  else
3303  return this->nElements();
3304  }
3308  inline size_t nvalues() const {
3309  return nnz();
3310  }
3313  size_t sy_nvalues() const {
3314  assert(this->nScalarsRows() == this->nScalarsCols());
3315  int n = this->nrows();
3316  return this->is_zero() ? 0 : n * (n + 1) / 2;
3317  }
3322  template<class Top>
3323  Treal syAccumulateWith(Top & op) {
3324  Treal res = 0;
3325  if (!this->is_zero()) {
3326  int rowOffset = this->rows.getOffset();
3327  int colOffset = this->cols.getOffset();
3328  for (int col = 0; col < this->ncols(); col++) {
3329  for (int row = 0; row < col; row++) {
3330  res += 2*op.accumulate((*this)(row, col),
3331  rowOffset + row,
3332  colOffset + col);
3333  }
3334  res += op.accumulate((*this)(col, col),
3335  rowOffset + col,
3336  colOffset + col);
3337  }
3338  }
3339  return res;
3340  }
3341  template<class Top>
3342  Treal geAccumulateWith(Top & op) {
3343  Treal res = 0;
3344  if (!this->is_zero()) {
3345  int rowOffset = this->rows.getOffset();
3346  int colOffset = this->cols.getOffset();
3347  for (int col = 0; col < this->ncols(); col++)
3348  for (int row = 0; row < this->nrows(); row++)
3349  res += op.accumulate((*this)(row, col),
3350  rowOffset + row,
3351  colOffset + col);
3352  }
3353  return res;
3354  }
3355 
3356  static inline unsigned int level() {return 0;}
3357 
3358  Treal maxAbsValue() const {
3359  if (this->is_zero())
3360  return 0;
3361  else {
3362  Treal maxAbsGlobal = 0;
3363  Treal maxAbsLocal = 0;
3364  for (int ind = 0; ind < this->nElements(); ++ind) {
3365  maxAbsLocal = template_blas_fabs(this->elements[ind]);
3366  maxAbsGlobal = maxAbsGlobal > maxAbsLocal ?
3367  maxAbsGlobal : maxAbsLocal;
3368  } /* end for */
3369  return maxAbsGlobal;
3370  }
3371  }
3372 
3373  /* New routines above */
3374 
3375 #if 0 /* OLD ROUTINES */
3376 
3377 
3378 #if 0
3379  inline Matrix<Treal>& operator=(const Matrix<Treal>& mat) {
3381  std::cout<<"operator="<<std::endl;
3382  }
3383 #endif
3384 
3385 
3386 
3387 
3388 
3389 
3390 
3391 
3392 
3393 
3394  void trim_memory_usage();
3395 #if 1
3396  void write_to_buffer_count(int& zb_length, int& vb_length) const;
3397  void write_to_buffer(int* zerobuf, const int zb_length,
3398  Treal* valuebuf, const int vb_length,
3399  int& zb_index, int& vb_index) const;
3400  void read_from_buffer(int* zerobuf, const int zb_length,
3401  Treal* valuebuf, const int vb_length,
3402  int& zb_index, int& vb_index);
3403 #endif
3404 
3405 
3406 
3407 
3408 
3409 
3410  /* continue here */
3411 
3412 
3413 
3414 
3415 
3416 
3417 
3418 
3419  inline bool lowestLevel() const {return true;}
3420  // inline unsigned int level() const {return 0;}
3421 
3422 #endif /* END OF OLD ROUTINES */
3423  protected:
3424  private:
3425  static const Treal ZERO;
3426  static const Treal ONE;
3427  }; /* end class specialization Matrix<Treal> */
3428 
3429  template<class Treal>
3430  const Treal Matrix<Treal>::ZERO = 0;
3431  template<class Treal>
3432  const Treal Matrix<Treal>::ONE = 1;
3433 
3434 #if 1
3435  /* Full matrix assigns etc */
3436  template<class Treal>
3437  void Matrix<Treal>::
3438  assignFromFull(std::vector<Treal> const & fullMat) {
3439  int nTotalRows = this->rows.getNTotalScalars();
3440  int nTotalCols = this->cols.getNTotalScalars();
3441  assert((int)fullMat.size() == nTotalRows * nTotalCols);
3442  int rowOffset = this->rows.getOffset();
3443  int colOffset = this->cols.getOffset();
3444  if (this->is_zero())
3445  allocate();
3446  for (int col = 0; col < this->ncols(); col++)
3447  for (int row = 0; row < this->nrows(); row++)
3448  (*this)(row, col) =
3449  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)];
3450  }
3451 
3452  template<class Treal>
3453  void Matrix<Treal>::
3454  fullMatrix(std::vector<Treal> & fullMat) const {
3455  int nTotalRows = this->rows.getNTotalScalars();
3456  int nTotalCols = this->cols.getNTotalScalars();
3457  fullMat.resize(nTotalRows * nTotalCols);
3458  int rowOffset = this->rows.getOffset();
3459  int colOffset = this->cols.getOffset();
3460  if (this->is_zero()) {
3461  for (int col = 0; col < this->nScalarsCols(); col++)
3462  for (int row = 0; row < this->nScalarsRows(); row++)
3463  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3464  }
3465  else {
3466  for (int col = 0; col < this->ncols(); col++)
3467  for (int row = 0; row < this->nrows(); row++)
3468  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] =
3469  (*this)(row, col);
3470  }
3471  }
3472 
3473  template<class Treal>
3474  void Matrix<Treal>::
3475  syFullMatrix(std::vector<Treal> & fullMat) const {
3476  int nTotalRows = this->rows.getNTotalScalars();
3477  int nTotalCols = this->cols.getNTotalScalars();
3478  fullMat.resize(nTotalRows * nTotalCols);
3479  int rowOffset = this->rows.getOffset();
3480  int colOffset = this->cols.getOffset();
3481  if (this->is_zero()) {
3482  for (int col = 0; col < this->nScalarsCols(); col++)
3483  for (int row = 0; row <= col; row++) {
3484  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3485  fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
3486  }
3487  }
3488  else {
3489  for (int col = 0; col < this->ncols(); col++)
3490  for (int row = 0; row <= col; row++) {
3491  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] =
3492  (*this)(row, col);
3493  fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] =
3494  (*this)(row, col);
3495  }
3496  }
3497  }
3498 
3499  template<class Treal>
3500  void Matrix<Treal>::
3501  syUpTriFullMatrix(std::vector<Treal> & fullMat) const {
3502  int nTotalRows = this->rows.getNTotalScalars();
3503  int nTotalCols = this->cols.getNTotalScalars();
3504  fullMat.resize(nTotalRows * nTotalCols);
3505  int rowOffset = this->rows.getOffset();
3506  int colOffset = this->cols.getOffset();
3507  if (this->is_zero()) {
3508  for (int col = 0; col < this->nScalarsCols(); col++)
3509  for (int row = 0; row <= this->nScalarsRows(); row++) {
3510  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3511  fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
3512  }
3513  }
3514  else {
3515  for (int col = 0; col < this->ncols(); col++)
3516  for (int row = 0; row < this->nrows(); row++) {
3517  fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] =
3518  (*this)(row, col);
3519  fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] =
3520  (*this)(row, col);
3521  }
3522  }
3523  }
3524 
3525 #endif
3526 
3527  template<class Treal>
3528  void Matrix<Treal>::
3529  assignFromSparse(std::vector<int> const & rowind,
3530  std::vector<int> const & colind,
3531  std::vector<Treal> const & values) {
3532  assert(rowind.size() == colind.size() &&
3533  rowind.size() == values.size());
3534  std::vector<int> indexes(values.size());
3535  for (int ind = 0; ind < values.size(); ++ind) {
3536  indexes[ind] = ind;
3537  }
3538  assignFromSparse(rowind, colind, values, indexes);
3539  }
3540 
3541  template<class Treal>
3542  void Matrix<Treal>::
3543  assignFromSparse(std::vector<int> const & rowind,
3544  std::vector<int> const & colind,
3545  std::vector<Treal> const & values,
3546  std::vector<int> const & indexes) {
3547  if (indexes.empty()) {
3548  this->clear();
3549  return;
3550  }
3551  if (this->is_zero())
3552  allocate();
3553  for (int ind = 0; ind < this->nElements(); ++ind)
3554  this->elements[ind] = 0;
3555  std::vector<int>::const_iterator it;
3556  for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3557  (*this)(this->rows.whichBlock( rowind[*it] ),
3558  this->cols.whichBlock( colind[*it] )) = values[*it];
3559  }
3560  }
3561 
3562 
3563  template<class Treal>
3564  void Matrix<Treal>::
3565  addValues(std::vector<int> const & rowind,
3566  std::vector<int> const & colind,
3567  std::vector<Treal> const & values) {
3568  assert(rowind.size() == colind.size() &&
3569  rowind.size() == values.size());
3570  std::vector<int> indexes(values.size());
3571  for (int ind = 0; ind < values.size(); ++ind) {
3572  indexes[ind] = ind;
3573  }
3574  addValues(rowind, colind, values, indexes);
3575  }
3576 
3577  template<class Treal>
3578  void Matrix<Treal>::
3579  addValues(std::vector<int> const & rowind,
3580  std::vector<int> const & colind,
3581  std::vector<Treal> const & values,
3582  std::vector<int> const & indexes) {
3583  if (indexes.empty())
3584  return;
3585  if (this->is_zero())
3586  allocate();
3587  std::vector<int>::const_iterator it;
3588  for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3589  (*this)(this->rows.whichBlock( rowind[*it] ),
3590  this->cols.whichBlock( colind[*it] )) += values[*it];
3591  }
3592  }
3593 
3594  template<class Treal>
3595  void Matrix<Treal>::
3596  syAssignFromSparse(std::vector<int> const & rowind,
3597  std::vector<int> const & colind,
3598  std::vector<Treal> const & values) {
3599  assert(rowind.size() == colind.size() &&
3600  rowind.size() == values.size());
3601  bool upperTriangleOnly = true;
3602  for (int ind = 0; ind < values.size(); ++ind) {
3603  upperTriangleOnly =
3604  upperTriangleOnly && colind[ind] >= rowind[ind];
3605  }
3606  if (!upperTriangleOnly)
3607  throw Failure("Matrix<Treal>::"
3608  "syAddValues(std::vector<int> const &, "
3609  "std::vector<int> const &, "
3610  "std::vector<Treal> const &, int const): "
3611  "Only upper triangle can contain elements when assigning "
3612  "symmetric or triangular matrix ");
3613  assignFromSparse(rowind, colind, values);
3614  }
3615 
3616  template<class Treal>
3617  void Matrix<Treal>::
3618  syAddValues(std::vector<int> const & rowind,
3619  std::vector<int> const & colind,
3620  std::vector<Treal> const & values) {
3621  assert(rowind.size() == colind.size() &&
3622  rowind.size() == values.size());
3623  bool upperTriangleOnly = true;
3624  for (int ind = 0; ind < values.size(); ++ind) {
3625  upperTriangleOnly =
3626  upperTriangleOnly && colind[ind] >= rowind[ind];
3627  }
3628  if (!upperTriangleOnly)
3629  throw Failure("Matrix<Treal>::"
3630  "syAddValues(std::vector<int> const &, "
3631  "std::vector<int> const &, "
3632  "std::vector<Treal> const &, int const): "
3633  "Only upper triangle can contain elements when assigning "
3634  "symmetric or triangular matrix ");
3635  addValues(rowind, colind, values);
3636  }
3637 
3638  template<class Treal>
3639  void Matrix<Treal>::
3640  getValues(std::vector<int> const & rowind,
3641  std::vector<int> const & colind,
3642  std::vector<Treal> & values) const {
3643  assert(rowind.size() == colind.size());
3644  values.resize(rowind.size(),0);
3645  std::vector<int> indexes(rowind.size());
3646  for (int ind = 0; ind < rowind.size(); ++ind) {
3647  indexes[ind] = ind;
3648  }
3649  getValues(rowind, colind, values, indexes);
3650  }
3651 
3652  template<class Treal>
3653  void Matrix<Treal>::
3654  getValues(std::vector<int> const & rowind,
3655  std::vector<int> const & colind,
3656  std::vector<Treal> & values,
3657  std::vector<int> const & indexes) const {
3658  assert(!this->is_empty());
3659  if (indexes.empty())
3660  return;
3661  std::vector<int>::const_iterator it;
3662  if (this->is_zero()) {
3663  for ( it = indexes.begin(); it < indexes.end(); it++ )
3664  values[*it] = 0;
3665  return;
3666  }
3667  for ( it = indexes.begin(); it < indexes.end(); it++ )
3668  values[*it] = (*this)(this->rows.whichBlock( rowind[*it] ),
3669  this->cols.whichBlock( colind[*it] ));
3670  }
3671 
3672 
3673  template<class Treal>
3674  void Matrix<Treal>::
3675  syGetValues(std::vector<int> const & rowind,
3676  std::vector<int> const & colind,
3677  std::vector<Treal> & values) const {
3678  assert(rowind.size() == colind.size());
3679  bool upperTriangleOnly = true;
3680  for (int ind = 0; ind < rowind.size(); ++ind) {
3681  upperTriangleOnly =
3682  upperTriangleOnly && colind[ind] >= rowind[ind];
3683  }
3684  if (!upperTriangleOnly)
3685  throw Failure("Matrix<Treal>::"
3686  "syGetValues(std::vector<int> const &, "
3687  "std::vector<int> const &, "
3688  "std::vector<Treal> const &, int const): "
3689  "Only upper triangle when retrieving elements from "
3690  "symmetric or triangular matrix ");
3691  getValues(rowind, colind, values);
3692  }
3693 
3694  template<class Treal>
3695  void Matrix<Treal>::
3696  getAllValues(std::vector<int> & rowind,
3697  std::vector<int> & colind,
3698  std::vector<Treal> & values) const {
3699  assert(rowind.size() == colind.size() &&
3700  rowind.size() == values.size());
3701  if (!this->is_zero())
3702  for (int col = 0; col < this->ncols(); col++)
3703  for (int row = 0; row < this->nrows(); row++) {
3704  rowind.push_back(this->rows.getOffset() + row);
3705  colind.push_back(this->cols.getOffset() + col);
3706  values.push_back((*this)(row, col));
3707  }
3708  }
3709 
3710  template<class Treal>
3711  void Matrix<Treal>::
3712  syGetAllValues(std::vector<int> & rowind,
3713  std::vector<int> & colind,
3714  std::vector<Treal> & values) const {
3715  assert(rowind.size() == colind.size() &&
3716  rowind.size() == values.size());
3717  if (!this->is_zero())
3718  for (int col = 0; col < this->ncols(); ++col)
3719  for (int row = 0; row <= col; ++row) {
3720  rowind.push_back(this->rows.getOffset() + row);
3721  colind.push_back(this->cols.getOffset() + col);
3722  values.push_back((*this)(row, col));
3723  }
3724  }
3725 
3726 
3727  template<class Treal>
3729  freeElements(this->elements);
3730  this->elements = 0;
3731  }
3732 
3733  template<class Treal>
3734  void Matrix<Treal>::
3735  writeToFile(std::ofstream & file) const {
3736  int const ZERO = 0;
3737  int const ONE = 1;
3738  if (this->is_zero()) {
3739  char * tmp = (char*)&ZERO;
3740  file.write(tmp,sizeof(int));
3741  }
3742  else {
3743  char * tmp = (char*)&ONE;
3744  file.write(tmp,sizeof(int));
3745  char * tmpel = (char*)this->elements;
3746  file.write(tmpel,sizeof(Treal) * this->nElements());
3747  }
3748  }
3749 
3750  template<class Treal>
3751  void Matrix<Treal>::
3752  readFromFile(std::ifstream & file) {
3753  int const ZERO = 0;
3754  int const ONE = 1;
3755  char tmp[sizeof(int)];
3756  file.read(tmp, (std::ifstream::pos_type)sizeof(int));
3757  switch ((int)*tmp) {
3758  case ZERO:
3759  this->clear();
3760  break;
3761  case ONE:
3762  if (this->is_zero())
3763  allocate();
3764  file.read((char*)this->elements, sizeof(Treal) * this->nElements());
3765  break;
3766  default:
3767  throw Failure("Matrix<Treal>::"
3768  "readFromFile(std::ifstream & file):"
3769  "File corruption, int value not 0 or 1");
3770  }
3771  }
3772 
3773  template<class Treal>
3775  if (this->is_zero())
3776  allocate();
3777  for (int ind = 0; ind < this->nElements(); ind++)
3778  this->elements[ind] = rand() / (Treal)RAND_MAX;
3779  }
3780 
3781  template<class Treal>
3783  if (this->is_zero())
3784  allocate();
3785  /* Above diagonal */
3786  for (int col = 1; col < this->ncols(); col++)
3787  for (int row = 0; row < col; row++)
3788  (*this)(row, col) = rand() / (Treal)RAND_MAX;;
3789  /* Diagonal */
3790  for (int rc = 0; rc < this->nrows(); rc++)
3791  (*this)(rc,rc) = rand() / (Treal)RAND_MAX;;
3792  }
3793 
3794  template<class Treal>
3795  void Matrix<Treal>::
3796  randomZeroStructure(Treal probabilityBeingZero) {
3797  if (!this->highestLevel() &&
3798  probabilityBeingZero > rand() / (Treal)RAND_MAX)
3799  this->clear();
3800  else
3801  this->random();
3802  }
3803 
3804  template<class Treal>
3805  void Matrix<Treal>::
3806  syRandomZeroStructure(Treal probabilityBeingZero) {
3807  if (!this->highestLevel() &&
3808  probabilityBeingZero > rand() / (Treal)RAND_MAX)
3809  this->clear();
3810  else
3811  this->syRandom();
3812  }
3813 
3814  template<class Treal>
3815  template<typename TRule>
3816  void Matrix<Treal>::
3817  setElementsByRule(TRule & rule) {
3818  if (this->is_zero())
3819  allocate();
3820  for (int col = 0; col < this->ncols(); col++)
3821  for (int row = 0; row < this->nrows(); row++)
3822  (*this)(row,col) = rule.set(this->rows.getOffset() + row,
3823  this->cols.getOffset() + col);
3824  }
3825  template<class Treal>
3826  template<typename TRule>
3827  void Matrix<Treal>::
3828  sySetElementsByRule(TRule & rule) {
3829  if (this->is_zero())
3830  allocate();
3831  /* Upper triangle */
3832  for (int col = 0; col < this->ncols(); col++)
3833  for (int row = 0; row <= col; row++)
3834  (*this)(row,col) = rule.set(this->rows.getOffset() + row,
3835  this->cols.getOffset() + col);
3836  }
3837 
3838 
3839  template<class Treal>
3840  void Matrix<Treal>::
3841  addIdentity(Treal alpha) {
3842  if (this->is_empty())
3843  throw Failure("Matrix<Treal>::addIdentity(Treal): "
3844  "Cannot add identity to empty matrix.");
3845  if (this->ncols() != this->nrows())
3846  throw Failure("Matrix<Treal, Telement>::addIdentity(Treal): "
3847  "Matrix must be square to add identity");
3848  if (this->is_zero())
3849  allocate();
3850  for (int ind = 0; ind < this->ncols(); ind++)
3851  (*this)(ind,ind) += alpha;
3852  }
3853 
3854  template<class Treal>
3855  void Matrix<Treal>::
3857  if (A.is_zero()) { /* Condition also matches empty matrices. */
3858  AT.rows = A.cols;
3859  AT.cols = A.rows;
3860  freeElements(AT.elements);
3861  AT.elements = 0;
3862  return;
3863  }
3864  if (AT.is_zero() || (AT.nElements() != A.nElements())) {
3865  freeElements(AT.elements);
3866  AT.elements = allocateElements<Treal>(A.nElements());
3867  }
3868  AT.cols = A.rows;
3869  AT.rows = A.cols;
3870  for (int row = 0; row < AT.nrows(); row++)
3871  for (int col = 0; col < AT.ncols(); col++)
3872  AT(row,col) = A(col,row);
3873  }
3874 
3875  template<class Treal>
3876  void Matrix<Treal>::
3878  if (this->nrows() == this->ncols()) {
3879  if (!this->is_zero()) {
3880  /* Diagonal should be fine */
3881  /* Fix the lower triangle */
3882  for (int row = 1; row < this->nrows(); row++)
3883  for (int col = 0; col < row; col++)
3884  (*this)(row, col) = (*this)(col, row);
3885  }
3886  }
3887  else
3888  throw Failure("Matrix<Treal>::symToNosym(): "
3889  "Only quadratic matrices can be symmetric");
3890  }
3891 
3892  template<class Treal>
3893  void Matrix<Treal>::
3895  if (this->nrows() == this->ncols()) {
3896  if (!this->is_zero()) {
3897  /* Diagonal should be fine */
3898  /* Fix the lower triangle */
3899  for (int col = 0; col < this->ncols() - 1; col++)
3900  for (int row = col + 1; row < this->nrows(); row++)
3901  (*this)(row, col) = 0;
3902  }
3903  }
3904  else
3905  throw Failure("Matrix<Treal>::nosymToSym(): "
3906  "Only quadratic matrices can be symmetric");
3907  }
3908 
3909  template<class Treal>
3910  Matrix<Treal>&
3912  switch (k) {
3913  case 0:
3914  this->clear();
3915  break;
3916  case 1:
3917  if (this->ncols() != this->nrows())
3918  throw Failure("Matrix<Treal>::operator=(int k = 1): "
3919  "Matrix must be quadratic to "
3920  "become identity matrix.");
3921  this->clear();
3922  this->allocate();
3923  for (int rc = 0; rc < this->ncols(); rc++) /*Set diagonal to identity*/
3924  (*this)(rc,rc) = 1;
3925  break;
3926  default:
3927  throw Failure
3928  ("Matrix<Treal>::operator=(int k) only implemented for k = 0, k = 1");
3929  }
3930  return *this;
3931  }
3932 
3933  template<class Treal>
3935  operator*=(const Treal alpha) {
3936  if (!this->is_zero() && alpha != 1) {
3937  for (int ind = 0; ind < this->nElements(); ind++)
3938  this->elements[ind] *= alpha;
3939  }
3940  return *this;
3941  }
3942 
3943  template<class Treal>
3944  void Matrix<Treal>::
3945  gemm(const bool tA, const bool tB, const Treal alpha,
3946  const Matrix<Treal>& A,
3947  const Matrix<Treal>& B,
3948  const Treal beta,
3949  Matrix<Treal>& C) {
3950  assert(!A.is_empty());
3951  assert(!B.is_empty());
3952  if (C.is_empty()) {
3953  assert(beta == 0);
3954  if (tA)
3955  C.resetRows(A.cols);
3956  else
3957  C.resetRows(A.rows);
3958  if (tB)
3959  C.resetCols(B.rows);
3960  else
3961  C.resetCols(B.cols);
3962  }
3963 
3964  if ((A.is_zero() || B.is_zero() || alpha == 0) &&
3965  (C.is_zero() || beta == 0))
3966  C = 0;
3967  else {
3968  Treal beta_tmp = beta;
3969  if (C.is_zero()) {
3970  C.allocate();
3971  beta_tmp = 0;
3972  }
3973 
3974  if (!A.is_zero() && !B.is_zero() && alpha != 0) {
3975  if (!tA && !tB)
3976  if (A.ncols() == B.nrows() &&
3977  A.nrows() == C.nrows() &&
3978  B.ncols() == C.ncols())
3979  mat::gemm("N", "N", &A.nrows(), &B.ncols(), &A.ncols(), &alpha,
3980  A.elements, &A.nrows(), B.elements, &B.nrows(),
3981  &beta_tmp, C.elements, &C.nrows());
3982  else
3983  throw Failure("Matrix<Treal>::"
3984  "gemm(bool, bool, Treal, Matrix<Treal>, "
3985  "Matrix<Treal>, Treal, "
3986  "Matrix<Treal>): "
3987  "Incorrect matrixdimensions for multiplication");
3988  else if (tA && !tB)
3989  if (A.nrows() == B.nrows() &&
3990  A.ncols() == C.nrows() &&
3991  B.ncols() == C.ncols())
3992  mat::gemm("T", "N", &A.ncols(), &B.ncols(), &A.nrows(), &alpha,
3993  A.elements, &A.nrows(), B.elements, &B.nrows(),
3994  &beta_tmp, C.elements, &C.nrows());
3995  else
3996  throw Failure("Matrix<Treal>::"
3997  "gemm(bool, bool, Treal, Matrix<Treal>, "
3998  "Matrix<Treal>, Treal, "
3999  "Matrix<Treal>): "
4000  "Incorrect matrixdimensions for multiplication");
4001  else if (!tA && tB)
4002  if (A.ncols() == B.ncols() &&
4003  A.nrows() == C.nrows() &&
4004  B.nrows() == C.ncols())
4005  mat::gemm("N", "T", &A.nrows(), &B.nrows(), &A.ncols(), &alpha,
4006  A.elements, &A.nrows(), B.elements, &B.nrows(),
4007  &beta_tmp, C.elements, &C.nrows());
4008  else
4009  throw Failure("Matrix<Treal>::"
4010  "gemm(bool, bool, Treal, Matrix<Treal>, "
4011  "Matrix<Treal>, Treal, "
4012  "Matrix<Treal>): "
4013  "Incorrect matrixdimensions for multiplication");
4014  else if (tA && tB)
4015  if (A.nrows() == B.ncols() &&
4016  A.ncols() == C.nrows() &&
4017  B.nrows() == C.ncols())
4018  mat::gemm("T", "T", &A.ncols(), &B.nrows(), &A.nrows(), &alpha,
4019  A.elements, &A.nrows(), B.elements, &B.nrows(),
4020  &beta_tmp, C.elements, &C.nrows());
4021  else
4022  throw Failure("Matrix<Treal>::"
4023  "gemm(bool, bool, Treal, Matrix<Treal>, "
4024  "Matrix<Treal>, Treal, "
4025  "Matrix<Treal>): "
4026  "Incorrect matrixdimensions for multiplication");
4027  else throw Failure("Matrix<Treal>::"
4028  "gemm(bool, bool, Treal, Matrix<Treal>, "
4029  "Matrix<Treal>, Treal, "
4030  "Matrix<Treal>):Very strange error!!");
4031  }
4032  else
4033  C *= beta;
4034  }
4035  }
4036 
4037 
4038  template<class Treal>
4039  void Matrix<Treal>::
4040  symm(const char side, const char uplo, const Treal alpha,
4041  const Matrix<Treal>& A,
4042  const Matrix<Treal>& B,
4043  const Treal beta,
4044  Matrix<Treal>& C) {
4045  assert(!A.is_empty());
4046  assert(!B.is_empty());
4047  assert(A.nrows() == A.ncols());
4048  // int dimA = A.nrows();
4049  if (C.is_empty()) {
4050  assert(beta == 0);
4051  if (side =='L') {
4052  C.resetRows(A.rows);
4053  C.resetCols(B.cols);
4054  }
4055  else {
4056  assert(side == 'R');
4057  C.resetRows(B.rows);
4058  C.resetCols(A.cols);
4059  }
4060  }
4061 
4062  if ((A.is_zero() || B.is_zero() || alpha == 0) &&
4063  (C.is_zero() || beta == 0))
4064  C = 0;
4065  else {
4066  Treal beta_tmp = beta;
4067  if (C.is_zero()) {
4068  C.allocate();
4069  beta_tmp = 0;
4070  }
4071  if (!A.is_zero() && !B.is_zero() && alpha != 0) {
4072  if (A.nrows() == A.ncols() && C.nrows() == B.nrows() &&
4073  C.ncols() == B.ncols() &&
4074  ((side == 'L' && A.ncols() == B.nrows()) ||
4075  (side == 'R' && B.ncols() == A.nrows())))
4076  mat::symm(&side, &uplo, &C.nrows(), &C.ncols(), &alpha,
4077  A.elements, &A.nrows(), B.elements, &B.nrows(),
4078  &beta_tmp, C.elements, &C.nrows());
4079  else
4080  throw Failure("Matrix<Treal>::symm: Incorrect matrix "
4081  "dimensions for symmetric multiply");
4082  } /* end if (!A.is_zero() && !B.is_zero() && alpha != 0) */
4083  else
4084  C *= beta;
4085  }
4086  }
4087 
4088  template<class Treal>
4089  void Matrix<Treal>::
4090  syrk(const char uplo, const bool tA, const Treal alpha,
4091  const Matrix<Treal>& A,
4092  const Treal beta,
4093  Matrix<Treal>& C) {
4094  assert(!A.is_empty());
4095  if (C.is_empty()) {
4096  assert(beta == 0);
4097  if (tA) {
4098  C.resetRows(A.cols);
4099  C.resetCols(A.cols);
4100  }
4101  else {
4102  C.resetRows(A.rows);
4103  C.resetCols(A.rows);
4104  }
4105  }
4106  if (C.nrows() == C.ncols() &&
4107  ((!tA && A.nrows() == C.nrows()) || (tA && A.ncols() == C.nrows())))
4108  if (alpha != 0 && !A.is_zero()) {
4109  Treal beta_tmp = beta;
4110  if (C.is_zero()) {
4111  C.allocate();
4112  beta_tmp = 0;
4113  }
4114  if (!tA) {
4115  mat::syrk(&uplo, "N", &C.nrows(), &A.ncols(), &alpha,
4116  A.elements, &A.nrows(), &beta_tmp,
4117  C.elements, &C.nrows());
4118  } /* end if (!tA) */
4119  else
4120  mat::syrk(&uplo, "T", &C.nrows(), &A.nrows(), &alpha,
4121  A.elements, &A.nrows(), &beta_tmp,
4122  C.elements, &C.nrows());
4123  }
4124  else
4125  C *= beta;
4126  else
4127  throw Failure("Matrix<Treal>::syrk: Incorrect matrix "
4128  "dimensions for symmetric rank-k update");
4129  }
4130 
4131  template<class Treal>
4132  void Matrix<Treal>::
4133  sysq(const char uplo, const Treal alpha,
4134  const Matrix<Treal>& A,
4135  const Treal beta,
4136  Matrix<Treal>& C) {
4137  assert(!A.is_empty());
4138  if (C.is_empty()) {
4139  assert(beta == 0);
4140  C.resetRows(A.rows);
4141  C.resetCols(A.cols);
4142  }
4143  /* FIXME: slow copy */
4144  Matrix<Treal> tmpA = A;
4145  tmpA.symToNosym();
4146  Matrix<Treal>::syrk(uplo, false, alpha, tmpA, beta, C);
4147  }
4148 
4149  /* C = alpha * A * B + beta * C where A and B are symmetric */
4150  template<class Treal>
4151  void Matrix<Treal>::
4152  ssmm(const Treal alpha,
4153  const Matrix<Treal>& A,
4154  const Matrix<Treal>& B,
4155  const Treal beta,
4156  Matrix<Treal>& C) {
4157  assert(!A.is_empty());
4158  assert(!B.is_empty());
4159  if (C.is_empty()) {
4160  assert(beta == 0);
4161  C.resetRows(A.rows);
4162  C.resetCols(B.cols);
4163  }
4164  /* FIXME: slow copy */
4165  Matrix<Treal> tmpB = B;
4166  tmpB.symToNosym();
4167  Matrix<Treal>::symm('L', 'U', alpha, A, tmpB, beta, C);
4168  }
4169 
4170  /* C = alpha * A * B + beta * C where A and B are symmetric
4171  * and only the upper triangle of C is computed.
4172  */
4173  template<class Treal>
4174  void Matrix<Treal>::
4175  ssmm_upper_tr_only(const Treal alpha,
4176  const Matrix<Treal>& A,
4177  const Matrix<Treal>& B,
4178  const Treal beta,
4179  Matrix<Treal>& C) {
4180  /* FIXME: Symmetry is not utilized. */
4181  ssmm(alpha, A, B, beta, C);
4182  C.nosymToSym();
4183  }
4184 
4185 
4186  template<class Treal>
4187  void Matrix<Treal>::
4188  trmm(const char side, const char uplo, const bool tA,
4189  const Treal alpha,
4190  const Matrix<Treal>& A,
4191  Matrix<Treal>& B) {
4192  assert(!A.is_empty());
4193  assert(!B.is_empty());
4194  if (alpha != 0 && !A.is_zero() && !B.is_zero())
4195  if (((side == 'R' && B.ncols() == A.nrows()) ||
4196  (side == 'L' && A.ncols() == B.nrows())) &&
4197  A.nrows() == A.ncols())
4198  if (!tA)
4199  mat::trmm(&side, &uplo, "N", "N", &B.nrows(), &B.ncols(),
4200  &alpha, A.elements, &A.nrows(), B.elements, &B.nrows());
4201  else
4202  mat::trmm(&side, &uplo, "T", "N", &B.nrows(), &B.ncols(),
4203  &alpha, A.elements, &A.nrows(), B.elements, &B.nrows());
4204  else
4205  throw Failure("Matrix<Treal>::trmm: "
4206  "Incorrect matrix dimensions for multiplication");
4207  else
4208  B = 0;
4209  }
4210 
4211  template<class Treal>
4213  assert(!this->is_empty());
4214  if (this->is_zero())
4215  return 0;
4216  else {
4217  Treal sum(0);
4218  for (int i = 0; i < this->nElements(); i++)
4219  sum += this->elements[i] * this->elements[i];
4220  return sum;
4221  }
4222  }
4223 
4224  template<class Treal>
4225  Treal Matrix<Treal>::
4227  assert(!this->is_empty());
4228  if (this->nrows() != this->ncols())
4229  throw Failure("Matrix<Treal>::syFrobSquared(): "
4230  "Matrix must be have equal number of rows "
4231  "and cols to be symmetric");
4232  Treal sum(0);
4233  if (!this->is_zero()) {
4234  for (int col = 1; col < this->ncols(); col++)
4235  for (int row = 0; row < col; row++)
4236  sum += 2 * (*this)(row, col) * (*this)(row, col);
4237  for (int rc = 0; rc < this->ncols(); rc++)
4238  sum += (*this)(rc, rc) * (*this)(rc, rc);
4239  }
4240  return sum;
4241  }
4242 
4243  template<class Treal>
4244  Treal Matrix<Treal>::
4246  (const Matrix<Treal>& A,
4247  const Matrix<Treal>& B) {
4248  assert(!A.is_empty());
4249  assert(!B.is_empty());
4250  if (A.nrows() != B.nrows() || A.ncols() != B.ncols())
4251  throw Failure("Matrix<Treal>::frobSquaredDiff: "
4252  "Incorrect matrix dimensions");
4253  Treal sum(0);
4254  if (!A.is_zero() && !B.is_zero()) {
4255  Treal diff;
4256  for (int i = 0; i < A.nElements(); i++) {
4257  diff = A.elements[i] - B.elements[i];
4258  sum += diff * diff;
4259  }
4260  }
4261  else if (A.is_zero() && !B.is_zero())
4262  sum = B.frobSquared();
4263  else if (!A.is_zero() && B.is_zero())
4264  sum = A.frobSquared();
4265  /* sum is already zero if A.is_zero() && B.is_zero() */
4266  return sum;
4267  }
4268 
4269 
4270  template<class Treal>
4271  Treal Matrix<Treal>::
4273  (const Matrix<Treal>& A,
4274  const Matrix<Treal>& B) {
4275  assert(!A.is_empty());
4276  assert(!B.is_empty());
4277  if (A.nrows() != B.nrows() ||
4278  A.ncols() != B.ncols() ||
4279  A.nrows() != A.ncols())
4280  throw Failure("Matrix<Treal>::syFrobSquaredDiff: "
4281  "Incorrect matrix dimensions");
4282  Treal sum(0);
4283  if (!A.is_zero() && !B.is_zero()) {
4284  Treal diff;
4285  for (int col = 1; col < A.ncols(); col++)
4286  for (int row = 0; row < col; row++) {
4287  diff = A(row, col) - B(row, col);
4288  sum += 2 * diff * diff;
4289  }
4290  for (int rc = 0; rc < A.ncols(); rc++) {
4291  diff = A(rc, rc) - B(rc, rc);
4292  sum += diff * diff;
4293  }
4294  }
4295  else if (A.is_zero() && !B.is_zero())
4296  sum = B.syFrobSquared();
4297  else if (!A.is_zero() && B.is_zero())
4298  sum = A.syFrobSquared();
4299  /* sum is already zero if A.is_zero() && B.is_zero() */
4300  return sum;
4301  }
4302 
4303  template<class Treal>
4304  Treal Matrix<Treal>::
4305  trace() const {
4306  assert(!this->is_empty());
4307  if (this->nrows() != this->ncols())
4308  throw Failure("Matrix<Treal>::trace(): Matrix must be quadratic");
4309  if (this->is_zero())
4310  return 0;
4311  else {
4312  Treal sum = 0;
4313  for (int rc = 0; rc < this->ncols(); rc++)
4314  sum += (*this)(rc,rc);
4315  return sum;
4316  }
4317  }
4318 
4319  template<class Treal>
4320  Treal Matrix<Treal>::
4322  const Matrix<Treal>& B) {
4323  assert(!A.is_empty());
4324  assert(!B.is_empty());
4325  if (A.nrows() != B.ncols() || A.ncols() != B.nrows())
4326  throw Failure("Matrix<Treal>::trace_ab: "
4327  "Wrong matrix dimensions for trace(A * B)");
4328  Treal tr = 0;
4329  if (!A.is_zero() && !B.is_zero())
4330  for (int rc = 0; rc < A.nrows(); rc++)
4331  for (int colA = 0; colA < A.ncols(); colA++)
4332  tr += A(rc,colA) * B(colA, rc);
4333  return tr;
4334  }
4335 
4336  template<class Treal>
4337  Treal Matrix<Treal>::
4339  const Matrix<Treal>& B) {
4340  assert(!A.is_empty());
4341  assert(!B.is_empty());
4342  if (A.ncols() != B.ncols() || A.nrows() != B.nrows())
4343  throw Failure("Matrix<Treal>::trace_aTb: "
4344  "Wrong matrix dimensions for trace(A' * B)");
4345  Treal tr = 0;
4346  if (!A.is_zero() && !B.is_zero())
4347  for (int rc = 0; rc < A.ncols(); rc++)
4348  for (int rowB = 0; rowB < B.nrows(); rowB++)
4349  tr += A(rowB,rc) * B(rowB, rc);
4350  return tr;
4351  }
4352 
4353  template<class Treal>
4354  Treal Matrix<Treal>::
4356  const Matrix<Treal>& B) {
4357  assert(!A.is_empty());
4358  assert(!B.is_empty());
4359  if (A.nrows() != B.ncols() || A.ncols() != B.nrows() ||
4360  A.nrows() != A.ncols())
4361  throw Failure("Matrix<Treal>::sy_trace_ab: "
4362  "Wrong matrix dimensions for symmetric trace(A * B)");
4363  if (A.is_zero() || B.is_zero())
4364  return 0;
4365  /* Now we know both A and B are nonzero */
4366  Treal tr = 0;
4367  /* Diagonal first */
4368  for (int rc = 0; rc < A.nrows(); rc++)
4369  tr += A(rc,rc) * B(rc, rc);
4370  /* Using that trace of transpose is equal to that without transpose: */
4371  for (int colA = 1; colA < A.ncols(); colA++)
4372  for (int rowA = 0; rowA < colA; rowA++)
4373  tr += 2 * A(rowA, colA) * B(rowA, colA);
4374  return tr;
4375  }
4376 
4377  template<class Treal>
4378  void Matrix<Treal>::
4379  add(const Treal alpha, /* B += alpha * A */
4380  const Matrix<Treal>& A,
4381  Matrix<Treal>& B) {
4382  assert(!A.is_empty());
4383  assert(!B.is_empty());
4384  if (A.nrows() != B.nrows() || A.ncols() != B.ncols())
4385  throw Failure("Matrix<Treal>::add: "
4386  "Wrong matrix dimensions for addition");
4387  if (!A.is_zero() && alpha != 0) {
4388  if (B.is_zero())
4389  B.allocate();
4390  for (int ind = 0; ind < A.nElements(); ind++)
4391  B.elements[ind] += alpha * A.elements[ind];
4392  }
4393  }
4394 
4395  template<class Treal>
4396  void Matrix<Treal>::
4397  assign(Treal const alpha, /* *this = alpha * A */
4398  Matrix<Treal> const & A) {
4399  assert(!A.is_empty());
4400  if (this->is_empty()) {
4401  this->resetRows(A.rows);
4402  this->resetCols(A.cols);
4403  }
4404  Matrix<Treal>::add(alpha, A, *this);
4405  }
4406 
4407 
4408  /********** Help functions for thresholding */
4409 
4410  template<class Treal>
4411  void Matrix<Treal>::
4412  getFrobSqLowestLevel(std::vector<Treal> & frobsq) const {
4413  if (!this->is_zero())
4414  frobsq.push_back(this->frobSquared());
4415  }
4416 
4417  template<class Treal>
4418  void Matrix<Treal>::
4419  getFrobSqElementLevel(std::vector<Treal> & frobsq) const {
4420  if (!this->is_zero())
4421  for (int ind = 0; ind < this->nElements(); ind++)
4422  if ( this->elements[ind] != 0 ) // Add nonzero elements only
4423  frobsq.push_back(this->elements[ind] * this->elements[ind]);
4424  }
4425 
4426 
4427  template<class Treal>
4428  void Matrix<Treal>::
4430  (Treal const threshold, Matrix<Treal> * ErrorMatrix) {
4431  if (ErrorMatrix) {
4432  if ((!this->is_zero() && this->frobSquared() <= threshold) ||
4433  (!ErrorMatrix->is_zero() &&
4434  ErrorMatrix->frobSquared() > threshold))
4435  Matrix<Treal>::swap(*this,*ErrorMatrix);
4436  }
4437  else if (!this->is_zero() && this->frobSquared() <= threshold)
4438  this->clear();
4439  }
4440 
4441  template<class Treal>
4442  void Matrix<Treal>::
4444  (Treal const threshold, Matrix<Treal> * ErrorMatrix) {
4445  assert(!this->is_empty());
4446  bool thisMatIsZero = true;
4447  if (ErrorMatrix) {
4448  assert(!ErrorMatrix->is_empty());
4449  bool EMatIsZero = true;
4450  if (!ErrorMatrix->is_zero() || !this->is_zero()) {
4451  if (ErrorMatrix->is_zero())
4452  ErrorMatrix->allocate();
4453  if (this->is_zero())
4454  this->allocate();
4455  for (int ind = 0; ind < this->nElements(); ind++) {
4456  if ( this->elements[ind] != 0 ) {
4457  assert(ErrorMatrix->elements[ind] == 0);
4458  // ok, let's check if we want to move the element to the error matrix
4459  if ( this->elements[ind] * this->elements[ind] <= threshold ) {
4460  // we want to move the element
4461  ErrorMatrix->elements[ind] = this->elements[ind];
4462  this->elements[ind] = 0;
4463  EMatIsZero = false; // at least one element is nonzero
4464  }
4465  else
4466  thisMatIsZero = false; // at least one element is nonzero
4467  continue;
4468  }
4469  if ( ErrorMatrix->elements[ind] != 0 ) {
4470  assert(this->elements[ind] == 0);
4471  // ok, let's check if we want to move the element from the error matrix
4472  if ( ErrorMatrix->elements[ind] * ErrorMatrix->elements[ind] > threshold ) {
4473  // we want to move the element
4474  this->elements[ind] = ErrorMatrix->elements[ind];
4475  ErrorMatrix->elements[ind] = 0;
4476  thisMatIsZero = false; // at least one element is nonzero
4477  }
4478  else
4479  EMatIsZero = false; // at least one element is nonzero
4480  }
4481  }
4482  if (thisMatIsZero) {
4483 #if 0
4484  for (int ind = 0; ind < this->nElements(); ind++)
4485  assert( this->elements[ind] == 0);
4486 #endif
4487  this->clear();
4488  }
4489  if (EMatIsZero) {
4490 #if 0
4491  for (int ind = 0; ind < this->nElements(); ind++)
4492  assert( ErrorMatrix->elements[ind] == 0);
4493 #endif
4494  ErrorMatrix->clear();
4495  }
4496  }
4497  }
4498  else
4499  if (!this->is_zero()) {
4500  for (int ind = 0; ind < this->nElements(); ind++) {
4501  if ( this->elements[ind] * this->elements[ind] <= threshold )
4502  /* FIXME BUG? EMANUEL LOOK AT THIS! */
4503  // this->elements[ind] == 0; OLD CODE. Compiler complained that "statement has no effect".
4504  this->elements[ind] = 0; /* New code. Changed from == to =. */
4505  else
4506  thisMatIsZero = false;
4507  }
4508  if (thisMatIsZero)
4509  this->clear();
4510  }
4511  }
4512 
4513 
4514 
4515  /********** End of help functions for thresholding */
4516 
4517  /* C = beta * C + alpha * A * B where only the upper triangle of C is */
4518  /* referenced and updated */
4519  template<class Treal>
4520  void Matrix<Treal>::
4521  gemm_upper_tr_only(const bool tA, const bool tB,
4522  const Treal alpha,
4523  const Matrix<Treal>& A,
4524  const Matrix<Treal>& B,
4525  const Treal beta,
4526  Matrix<Treal>& C) {
4527  /* FIXME: Symmetry is not utilized. */
4528  Matrix<Treal>::gemm(tA, tB, alpha, A, B, beta, C);
4529  C.nosymToSym();
4530  }
4531 
4532  /* A = alpha * A * Z or A = alpha * Z * A where A is symmetric, */
4533  /* Z is upper triangular and */
4534  /* only the upper triangle of the result is calculated */
4535  /* side == 'R' for A = alpha * A * Z */
4536  /* side == 'L' for A = alpha * Z * A */
4537  template<class Treal>
4538  void Matrix<Treal>::
4539  sytr_upper_tr_only(char const side, const Treal alpha,
4540  Matrix<Treal>& A,
4541  const Matrix<Treal>& Z) {
4542  /* FIXME: Symmetry is not utilized. */
4543  A.symToNosym();
4544  Matrix<Treal>::trmm(side, 'U', false, alpha, Z, A);
4545  A.nosymToSym();
4546  }
4547 
4548  /* The result B is assumed to be symmetric, i.e. only upper triangle */
4549  /* calculated and hence only upper triangle of input B is used. */
4550  template<class Treal>
4551  void Matrix<Treal>::
4552  trmm_upper_tr_only(const char side, const char uplo,
4553  const bool tA, const Treal alpha,
4554  const Matrix<Treal>& A,
4555  Matrix<Treal>& B) {
4556  /* FIXME: Symmetry is not utilized. */
4557  assert(tA);
4558  B.symToNosym();
4559  Matrix<Treal>::trmm(side, uplo, tA, alpha, A, B);
4560  B.nosymToSym();
4561  }
4562 
4563  /* A = Z' * A * Z or A = Z * A * Z' */
4564  /* where A is symmetric and Z is (nonunit) upper triangular */
4565  /* side == 'R' for A = Z' * A * Z */
4566  /* side == 'L' for A = Z * A * Z' */
4567  template<class Treal>
4568  void Matrix<Treal>::
4569  trsytriplemm(char const side,
4570  const Matrix<Treal>& Z,
4571  Matrix<Treal>& A) {
4572  if (side == 'R') {
4574  sytr_upper_tr_only('R', 1.0, A, Z);
4576  trmm_upper_tr_only('L', 'U', true, 1.0, Z, A);
4577  }
4578  else {
4579  assert(side == 'L');
4581  sytr_upper_tr_only('L', 1.0, A, Z);
4583  trmm_upper_tr_only('R', 'U', true, 1.0, Z, A);
4584  }
4585  }
4586 
4587 
4588  template<class Treal>
4590  (Treal const threshold, Matrix<Treal> * ErrorMatrix) {
4591  assert(!this->is_empty());
4592  if (ErrorMatrix && ErrorMatrix->is_empty()) {
4593  ErrorMatrix->resetRows(this->rows);
4594  ErrorMatrix->resetCols(this->cols);
4595  }
4596  Treal fs = frobSquared();
4597  if (fs < threshold) {
4598  if (ErrorMatrix)
4599  Matrix<Treal>::swap(*this, *ErrorMatrix);
4600  return fs;
4601  }
4602  else
4603  return 0;
4604  }
4605 
4606 
4607  template<class Treal>
4608  void Matrix<Treal>::
4610  Matrix<Treal>& Z,
4611  const Treal threshold, const side looking,
4612  const inchversion version) {
4613  assert(!A.is_empty());
4614  if (A.nrows() != A.ncols())
4615  throw Failure("Matrix<Treal>::inch: Matrix must be quadratic!");
4616  if (A.is_zero())
4617  throw Failure("Matrix<Treal>::inch: Matrix is zero! "
4618  "Not possible to compute inverse cholesky.");
4619  Z = A;
4620  int info;
4621  potrf("U", &A.nrows(), Z.elements, &A.nrows(), &info);
4622  if (info > 0)
4623  throw Failure("Matrix<Treal>::inch: potrf returned info > 0. The matrix is not positive definite.");
4624  if (info < 0)
4625  throw Failure("Matrix<Treal>::inch: potrf returned info < 0");
4626 
4627  trtri("U", "N", &A.nrows(), Z.elements, &A.nrows(), &info);
4628  if (info > 0)
4629  throw Failure("Matrix<Treal>::inch: trtri returned info > 0. The matrix is not positive definite.");
4630  if (info < 0)
4631  throw Failure("Matrix<Treal>::inch: trtri returned info < 0");
4632  /* Fill lower triangle with zeroes just to be safe */
4633  trifulltofull(Z.elements, A.nrows());
4634  }
4635 
4636  template<class Treal>
4637  void Matrix<Treal>::
4638  gersgorin(Treal& lmin, Treal& lmax) const {
4639  assert(!this->is_empty());
4640  if (this->nScalarsRows() == this->nScalarsCols()) {
4641  int size = this->nScalarsCols();
4642  Treal* colsums = new Treal[size];
4643  Treal* diag = new Treal[size];
4644  for (int ind = 0; ind < size; ind++)
4645  colsums[ind] = 0;
4646  this->add_abs_col_sums(colsums);
4647  this->get_diagonal(diag);
4648  Treal tmp1 = colsums[0] - template_blas_fabs(diag[0]);
4649  Treal tmp2;
4650  lmin = diag[0] - tmp1;
4651  lmax = diag[0] + tmp1;
4652  for (int col = 1; col < size; col++) {
4653  tmp1 = colsums[col] - template_blas_fabs(diag[col]);
4654  tmp2 = diag[col] - tmp1;
4655  lmin = (tmp2 < lmin ? tmp2 : lmin);
4656  tmp2 = diag[col] + tmp1;
4657  lmax = (tmp2 > lmax ? tmp2 : lmax);
4658  }
4659  delete[] diag;
4660  delete[] colsums;
4661  }
4662  else
4663  throw Failure("Matrix<Treal>::gersgorin(Treal, Treal): Matrix must be quadratic");
4664  }
4665 
4666 
4667  template<class Treal>
4668  void Matrix<Treal>::
4669  add_abs_col_sums(Treal* sums) const {
4670  assert(sums);
4671  if (!this->is_zero())
4672  for (int col = 0; col < this->ncols(); col++)
4673  for (int row = 0; row < this->nrows(); row++)
4674  sums[col] += template_blas_fabs((*this)(row,col));
4675  }
4676 
4677  template<class Treal>
4678  void Matrix<Treal>::
4679  get_diagonal(Treal* diag) const {
4680  assert(diag);
4681  assert(this->nrows() == this->ncols());
4682  if (this->is_zero())
4683  for (int rc = 0; rc < this->nScalarsCols(); rc++)
4684  diag[rc] = 0;
4685  else
4686  for (int rc = 0; rc < this->ncols(); rc++)
4687  diag[rc] = (*this)(rc,rc);
4688  }
4689 
4690 
4691  /* New routines above */
4692 
4693 #if 0 /* OLD ROUTINES */
4694 
4695 
4696 
4697 
4698 
4699 
4700 
4701 
4702 
4703 
4704 
4705 
4706 
4707 
4708 
4709 
4710 
4711 
4712 
4713 
4714  template<class Treal>
4716  if (this->is_zero() && this->cap > 0) {
4717  freeElements(this->elements);
4718  this->elements = NULL;
4719  this->cap = 0;
4720  this->nel = 0;
4721  }
4722  else if (this->cap > this->nel) {
4723  Treal* tmp = new Treal[this->nel];
4724  for (int i = 0; i < this->nel; i++)
4725  tmp[i] = this->elements[i];
4726  freeElements(this->elements);
4727  this->cap = this->nel;
4728  this->elements = tmp;
4729  }
4730  assert(this->cap == this->nel);
4731  }
4732 
4733 
4734 
4735 #if 1
4736 
4737  template<class Treal>
4738  void Matrix<Treal>::
4739  write_to_buffer_count(int& zb_length, int& vb_length) const {
4740  if (this->is_zero()) {
4741  ++zb_length;
4742  }
4743  else {
4744  ++zb_length;
4745  vb_length += this->nel;
4746  }
4747  }
4748 
4749  template<class Treal>
4750  void Matrix<Treal>::
4751  write_to_buffer(int* zerobuf, const int zb_length,
4752  Treal* valuebuf, const int vb_length,
4753  int& zb_index, int& vb_index) const {
4754  if (zb_index < zb_length) {
4755  if (this->is_zero()) {
4756  zerobuf[zb_index] = 0;
4757  ++zb_index;
4758  }
4759  else {
4760  if (vb_index + this->nel < vb_length + 1) {
4761  zerobuf[zb_index] = 1;
4762  ++zb_index;
4763  for (int i = 0; i < this->nel; i++)
4764  valuebuf[vb_index + i] = this->elements[i];
4765  vb_index += this->nel;
4766  }
4767  else
4768  throw Failure("Matrix<Treal, Telement>::write_to_buffer: "
4769  "Insufficient space in buffers");
4770  }
4771  }
4772  else
4773  throw Failure("Matrix<Treal, Telement>::write_to_buffer: "
4774  "Insufficient space in buffers");
4775  }
4776 
4777  template<class Treal>
4778  void Matrix<Treal>::
4779  read_from_buffer(int* zerobuf, const int zb_length,
4780  Treal* valuebuf, const int vb_length,
4781  int& zb_index, int& vb_index) {
4782  if (zb_index < zb_length) {
4783  if (zerobuf[zb_index] == 0) {
4784  (*this) = 0;
4785  ++zb_index;
4786  }
4787  else {
4788  this->content = ful;
4789  this->nel = this->nrows() * this->ncols();
4790  this->assert_alloc();
4791  if (vb_index + this->nel < vb_length + 1) {
4792  assert(zerobuf[zb_index] == 1);
4793  ++zb_index;
4794  for (int i = 0; i < this->nel; i++)
4795  this->elements[i] = valuebuf[vb_index + i];
4796  vb_index += this->nel;
4797  }
4798  else
4799  throw Failure("Matrix<Treal, Telement>::read_from_buffer: "
4800  "Mismatch, buffers does not match matrix");
4801  }
4802  }
4803  else
4804  throw Failure("Matrix<Treal, Telement>::read_from_buffer: "
4805  "Mismatch, buffers does not match matrix");
4806  }
4807 
4808 #endif
4809 
4810 
4811 
4812 
4813 
4814 
4815 
4816 
4817  /* continue here */
4818 
4819 
4820 
4821 
4822 
4823 
4824 
4825 
4826 
4827 
4828 
4829 
4830 
4831 
4832 #if 1
4833 
4834 
4835 
4836 #endif
4837 
4838 #endif /* END OF OLD ROUTINES */
4839 
4840 
4841 } /* end namespace mat */
4842 
4843 #endif
#define A
void syFullMatrix(std::vector< Treal > &fullMat) const
Definition: Matrix.h:537
Vector< Treal, typename ElementType::VectorType > VectorType
Definition: Matrix.h:116
static Treal sy_trace_ab(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1986
static void trifulltofull(Treal *trifull, const int size)
Definition: gblas.h:1191
void syAssignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A)
Definition: Matrix.h:3132
size_t nnz() const
Returns number of nonzeros in matrix.
Definition: Matrix.h:2858
Definition: Matrix.h:2909
int whichBlock(int const globalIndex) const
Returns the blocknumber (between 0 and nBlocks-1) that contains elements with the given global index...
Definition: SizesAndBlocks.h:71
void assign(Treal const alpha, Matrix< Treal, Telement > const &A)
Definition: Matrix.h:2027
void frobThreshElementLevel(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix)
Definition: Matrix.h:2099
SizesAndBlocks getSizesAndBlocksForLowerLevel(int const blockNumber) const
Definition: SizesAndBlocks.cc:63
Treal frob_squared_thresh(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix=0)
Removes small elements so that the introduced error is smaller than threshold in the squared Frobeniu...
Definition: Matrix.h:2624
void getValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Definition: Matrix.h:731
size_t nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:442
Definition: Matrix.h:73
Vector< Treal, Treal > VectorType
Definition: Matrix.h:2912
SizesAndBlocks cols
Definition: MatrixHierarchicBase.h:163
Treal frob() const
Definition: Matrix.h:284
double accumulatedTime
Definition: Matrix.h:82
void assignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A, Matrix< Treal, Matrix< Treal, Telement > > const &B)
Same as assignFrobNormsLowestLevel except that the Frobenius norms of the differences between submatr...
Definition: Matrix.h:2173
Treal syFrobSquared() const
Definition: Matrix.h:1866
#define MAT_OMP_FINALIZE
Definition: matInclude.h:88
inchversion
Definition: Matrix.h:74
void writeToFile(std::ofstream &file) const
Definition: Matrix.h:843
static void syInch(const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &Z, const Treal threshold=0, const side looking=left, const inchversion version=unstable)
Definition: Matrix.h:2674
void symToNosym()
Definition: Matrix.h:3877
void tic()
Definition: matInclude.cc:77
void getAllValues(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &) const
Definition: Matrix.h:806
static void add(const Treal alpha, const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &B)
Definition: Matrix.h:2009
void getFrobSqLowestLevel(std::vector< Treal > &frobsq) const
Definition: Matrix.h:2045
bool highestLevel() const
Definition: MatrixHierarchicBase.h:131
int getOffset() const
Definition: SizesAndBlocks.h:75
void getFrobSqElementLevel(std::vector< Treal > &frobsq) const
Definition: Matrix.h:2090
Base class for Matrix.
int getNTotalScalars() const
Definition: SizesAndBlocks.h:76
Proxy structs used by the matrix API.
Treal syFrobSquared() const
Definition: Matrix.h:4226
void fullMatrix(std::vector< Treal > &fullMat) const
Definition: Matrix.h:517
size_t sy_nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:2884
static void symm(const char *side, const char *uplo, const int *m, const int *n, const T *alpha, const T *A, const int *lda, const T *B, const int *ldb, const T *beta, T *C, const int *ldc)
Definition: gblas.h:340
static void transpose(Matrix< Treal, Telement > const &A, Matrix< Treal, Telement > &AT)
Definition: Matrix.h:977
static void gemm(const char *ta, const char *tb, const int *n, const int *k, const int *l, const T *alpha, const T *A, const int *lda, const T *B, const int *ldb, const T *beta, T *C, const int *ldc)
Definition: gblas.h:230
Matrix()
Definition: Matrix.h:2915
void setElementsByRule(TRule &rule)
Definition: Matrix.h:938
Vector class.
Definition: Matrix.h:76
MatrixHierarchicBase< Treal, Telement > & operator=(const MatrixHierarchicBase< Treal, Telement > &mat)
Definition: MatrixHierarchicBase.h:200
return elements[row+col *nrows()]
Definition: MatrixHierarchicBase.h:79
void syRandomZeroStructure(Treal probabilityBeingZero)
Definition: Matrix.h:918
size_t sy_nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:3313
static void ssmm(const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1548
void reset()
Definition: Matrix.h:84
Treal frob_thresh(Treal const threshold, Matrix< Treal > *ErrorMatrix=0)
Definition: Matrix.h:3251
void truncateAccordingToSparsityPattern(Matrix< Treal, Matrix< Treal, Telement > > &A) const
Truncate matrix A according to the sparsity pattern of the this matrix (frobNormMat).
Definition: Matrix.h:2241
Matrix< Treal, Telement > & operator=(const Matrix< Treal, Telement > &mat)
Definition: Matrix.h:196
size_t memory_usage() const
Definition: Matrix.h:3285
void syAssignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A, Matrix< Treal, Matrix< Treal > > const &B)
Definition: Matrix.h:3179
const int & ncols() const
Definition: MatrixHierarchicBase.h:69
const int & nScalarsRows() const
Definition: MatrixHierarchicBase.h:61
static void syrk(const char *uplo, const char *trans, const int *n, const int *k, const T *alpha, const T *A, const int *lda, const T *beta, T *C, const int *ldc)
Definition: gblas.h:332
Treal frob() const
Definition: Matrix.h:3059
const int & nrows() const
Definition: MatrixHierarchicBase.h:67
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
void get_diagonal(Treal *diag) const
Definition: Matrix.h:2831
void addValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:638
static Treal syFrobSquaredDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1908
void syAddValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:709
void clear()
Set matrix to zero and delete all arrays.
Definition: Matrix.h:3728
Definition: matInclude.h:158
Treal template_blas_fabs(Treal x)
float toc()
Definition: matInclude.cc:81
Treal syFrob() const
Definition: Matrix.h:289
Definition: Matrix.h:80
Definition: matInclude.h:134
void allocate()
Definition: Matrix.h:122
void symToNosym()
Definition: Matrix.h:999
static void trmm_upper_tr_only(const char side, const char uplo, const bool tA, const Treal alpha, const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &B)
Definition: Matrix.h:2527
Matrix< Treal, Telement > & operator*=(const Treal alpha)
Definition: Matrix.h:1071
static Treal trace_ab(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1952
C++ interface to a subset of BLAS and LAPACK.
static void gemm_upper_tr_only(const bool tA, const bool tB, const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:2261
Base class for Matrix and Matrix specialization.
Definition: MatrixHierarchicBase.h:50
void sy_gersgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:420
static void trtri(const char *uplo, const char *diag, const int *n, T *A, const int *lda, int *info)
Definition: gblas.h:319
void addTime(double timeToAdd)
Definition: Matrix.h:86
static Treal trace_aTb(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1969
void sy_gersgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:3276
void syAssignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A)
Version of assignFrobNormsLowestLevelToMatrix for symmetric matrices.
Definition: Matrix.h:2151
void random()
Definition: Matrix.h:882
static unsigned int level()
Definition: Matrix.h:3356
Treal trace() const
Definition: Matrix.h:1935
static Treal syFrobDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:304
void syRandom()
Definition: Matrix.h:890
void truncateAccordingToSparsityPattern(Matrix< Treal, Matrix< Treal > > &A) const
Definition: Matrix.h:3219
Treal geAccumulateWith(Top &op)
Definition: Matrix.h:3342
static Treal frobDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:295
static void potrf(const char *uplo, const int *n, T *A, const int *lda, int *info)
Definition: gblas.h:312
void allocate()
Definition: Matrix.h:2924
static const Treal ONE
Definition: Matrix.h:3426
void freeElements(float *ptr)
Definition: allocate.cc:40
void gersgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:4638
int nElements() const
Definition: MatrixHierarchicBase.h:108
static void ssmm_upper_tr_only(const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1676
size_t sy_nnz() const
Returns number of nonzeros in matrix including lower triangle elements.
Definition: Matrix.h:2867
#define MAT_OMP_START
Definition: matInclude.h:86
static unsigned int level()
Definition: Matrix.h:478
void syGetAllValues(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &) const
Definition: Matrix.h:818
Telement ElementType
Definition: Matrix.h:115
Matrix class and heart of the matrix library.
Definition: Matrix.h:113
Matrix< Treal > & operator=(const Matrix< Treal > &mat)
Definition: Matrix.h:2983
void randomZeroStructure(Treal probabilityBeingZero)
Get a random zero structure with a specified probability that each submatrix is zero.
Definition: Matrix.h:904
double getAccumulatedTime()
Definition: Matrix.h:85
static Treal frobSquaredDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1886
void frobThreshLowestLevel(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix)
Definition: Matrix.h:2054
void assignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:586
#define MAT_OMP_END
Definition: matInclude.h:87
static void syInch(const Matrix< Treal > &A, Matrix< Treal > &Z, const Treal threshold=0, const side looking=left, const inchversion version=unstable)
Definition: Matrix.h:3267
~Matrix()
Definition: Matrix.h:203
size_t nnz() const
Returns number of nonzeros in matrix.
Definition: Matrix.h:3293
static void trmm(const char *side, const char *uplo, const char *transa, const char *diag, const int *m, const int *n, const T *alpha, const T *A, const int *lda, T *B, const int *ldb)
Definition: gblas.h:279
size_t nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:3308
Treal frobSquared() const
Definition: Matrix.h:1852
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: Matrix.h:504
SingletonForTimings()
Definition: Matrix.h:100
size_t sy_nnz() const
Returns number of nonzeros in matrix including lower triangle elements.
Definition: Matrix.h:3299
bool is_empty() const
Check if matrix is empty Empty is different from zero, a zero matrix contains information about block...
Definition: MatrixHierarchicBase.h:141
void syAssignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:687
void syGetValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Definition: Matrix.h:784
void sySetElementsByRule(TRule &rule)
Definition: Matrix.h:947
static void trsytriplemm(char const side, const Matrix< Treal, Telement > &Z, Matrix< Treal, Telement > &A)
Definition: Matrix.h:2604
static void syrk(const char uplo, const bool tA, const Treal alpha, const Matrix< Treal, Telement > &A, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1362
side
Definition: Matrix.h:73
static SingletonForTimings & instance()
Definition: Matrix.h:94
Definition: Matrix.h:78
void syUpTriFullMatrix(std::vector< Treal > &fullMat) const
Definition: Matrix.h:561
void resetCols(SizesAndBlocks const &newCols)
Definition: MatrixHierarchicBase.h:117
Definition: Failure.h:47
int const & getNBlocks() const
Definition: SizesAndBlocks.h:64
static void trmm(const char side, const char uplo, const bool tA, const Treal alpha, const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &B)
Definition: Matrix.h:1772
size_t memory_usage() const
Definition: Matrix.h:2847
void assignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A, Matrix< Treal, Matrix< Treal > > const &B)
Definition: Matrix.h:3147
void syAssignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A, Matrix< Treal, Matrix< Treal, Telement > > const &B)
Same as syAssignFrobNormsLowestLevel except that the Frobenius norms of the differences between subma...
Definition: Matrix.h:2204
void add_abs_col_sums(Treal *abscolsums) const
Definition: Matrix.h:2816
~Matrix()
Definition: Matrix.h:2989
void assignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A)
Build a matrix with single matrix elements at the lowest level containing the Frobenius norms of the ...
Definition: Matrix.h:2137
void trSetElementsByRule(TRule &rule)
Definition: Matrix.h:222
Treal geAccumulateWith(Top &op)
Accumulation algorithm for general matrices.
Definition: Matrix.h:468
#define MAT_OMP_INIT
Definition: matInclude.h:85
Treal frobSquared() const
Definition: Matrix.h:4212
static void symm(const char side, const char uplo, const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1247
void readFromFile(std::ifstream &file)
Definition: Matrix.h:859
void nosymToSym()
Definition: Matrix.h:3894
Matrix()
Definition: Matrix.h:119
#define B
bool is_zero() const
Definition: MatrixHierarchicBase.h:106
static void sysq(const char uplo, const Treal alpha, const Matrix< Treal, Telement > &A, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1465
Definition: Matrix.h:73
static const Treal ZERO
Definition: Matrix.h:3425
SizesAndBlocks rows
Definition: MatrixHierarchicBase.h:162
Treal syAccumulateWith(Top &op)
Definition: Matrix.h:3323
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:129
const int & nScalarsCols() const
Definition: MatrixHierarchicBase.h:63
Treal syFrob() const
Definition: Matrix.h:3061
Treal trace(const XYZ< Treal, MatrixGeneral< Treal, Tmatrix >, MatrixGeneral< Treal, Tmatrix > > &smm)
Definition: MatrixGeneral.h:902
void resetRows(SizesAndBlocks const &newRows)
Definition: MatrixHierarchicBase.h:112
void addIdentity(Treal alpha)
Definition: Matrix.h:962
Definition: Matrix.h:74
static void sytr_upper_tr_only(char const side, const Treal alpha, Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &Z)
Definition: Matrix.h:2444
Treal template_blas_sqrt(Treal x)
Treal maxAbsValue() const
Definition: Matrix.h:480
void nosymToSym()
Definition: Matrix.h:1025
void clear()
Definition: Matrix.h:833
Treal maxAbsValue() const
Definition: Matrix.h:3358
static void gemm(const bool tA, const bool tB, const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1082
Treal frob_thresh(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix=0)
Removes small elements so that the introduced error is smaller than the threshold in the Frobenius no...
Definition: Matrix.h:394
void gersgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:2784
Definition: Matrix.h:74
Treal syAccumulateWith(Top &op)
Definition: Matrix.h:455
static void swap(MatrixHierarchicBase< Treal, Telement > &A, MatrixHierarchicBase< Treal, Telement > &B)
Definition: MatrixHierarchicBase.h:221