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