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