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