ergo
MatrixGeneral.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program. If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
38 #ifndef MAT_MATRIXGENERAL
39 #define MAT_MATRIXGENERAL
40 #include "MatrixBase.h"
41 namespace mat {
58  template<typename Treal, typename Tmatrix>
59  class MatrixGeneral : public MatrixBase<Treal, Tmatrix> {
60  public:
62 
64  :MatrixBase<Treal, Tmatrix>() {}
66  :MatrixBase<Treal, Tmatrix>(matr) {}
68  :MatrixBase<Treal, Tmatrix>(symm) {
69  this->matrixPtr->symToNosym();
70  }
72  :MatrixBase<Treal, Tmatrix>(triang) {}
75 #if 0
76  template<typename Tfull>
77  inline void assign_from_full
78  (Tfull const* const fullmatrix,
79  const int totnrows, const int totncols) {
80  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
81  }
82  inline void assign_from_full
83  (Treal const* const fullmatrix,
84  const int totnrows, const int totncols) {
85  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
86  }
87 #endif
88 
89  inline void assignFromFull
90  (std::vector<Treal> const & fullMat) {
91  assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
92  this->matrixPtr->assignFromFull(fullMat);
93  }
94 
95  inline void fullMatrix(std::vector<Treal> & fullMat) const {
96  this->matrixPtr->fullMatrix(fullMat);
97  }
98 
99  inline void fullMatrix
100  (std::vector<Treal> & fullMat,
101  std::vector<int> const & rowInversePermutation,
102  std::vector<int> const & colInversePermutation) const {
103  std::vector<int> rowind;
104  std::vector<int> colind;
105  std::vector<Treal> values;
106  get_all_values(rowind, colind, values,
107  rowInversePermutation,
108  colInversePermutation);
109  fullMat.assign(this->get_nrows()*this->get_ncols(),0);
110  for (unsigned int ind = 0; ind < values.size(); ++ind)
111  fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
112  values[ind];
113  }
121  inline void assign_from_sparse
122  (std::vector<int> const & rowind,
123  std::vector<int> const & colind,
124  std::vector<Treal> const & values,
125  SizesAndBlocks const & newRows,
126  SizesAndBlocks const & newCols) {
127  this->resetSizesAndBlocks(newRows, newCols);
128  this->matrixPtr->assignFromSparse(rowind, colind, values);
129  }
137  inline void assign_from_sparse
138  (std::vector<int> const & rowind,
139  std::vector<int> const & colind,
140  std::vector<Treal> const & values,
141  std::vector<int> const & rowPermutation,
142  std::vector<int> const & colPermutation) {
143  std::vector<int> newRowind;
144  std::vector<int> newColind;
146  getPermutedIndexes(rowind, rowPermutation, newRowind);
148  getPermutedIndexes(colind, colPermutation, newColind);
149  this->matrixPtr->assignFromSparse(newRowind, newColind, values);
150  }
156  inline void assign_from_sparse
157  (std::vector<int> const & rowind,
158  std::vector<int> const & colind,
159  std::vector<Treal> const & values,
160  SizesAndBlocks const & newRows,
161  SizesAndBlocks const & newCols,
162  std::vector<int> const & rowPermutation,
163  std::vector<int> const & colPermutation) {
164  this->resetSizesAndBlocks(newRows, newCols);
165  assign_from_sparse(rowind, colind, values,
166  rowPermutation, colPermutation);
167  }
172  inline void get_values
173  (std::vector<int> const & rowind,
174  std::vector<int> const & colind,
175  std::vector<Treal> & values) const {
176  this->matrixPtr->getValues(rowind, colind, values);
177  }
185  inline void get_values
186  (std::vector<int> const & rowind,
187  std::vector<int> const & colind,
188  std::vector<Treal> & values,
189  std::vector<int> const & rowPermutation,
190  std::vector<int> const & colPermutation) const {
191  std::vector<int> newRowind;
192  std::vector<int> newColind;
194  getPermutedIndexes(rowind, rowPermutation, newRowind);
196  getPermutedIndexes(colind, colPermutation, newColind);
197  this->matrixPtr->getValues(newRowind, newColind, values);
198  }
203  inline void get_all_values
204  (std::vector<int> & rowind,
205  std::vector<int> & colind,
206  std::vector<Treal> & values) const {
207  rowind.resize(0);
208  colind.resize(0);
209  values.resize(0);
210  this->matrixPtr->getAllValues(rowind, colind, values);
211  }
221  inline void get_all_values
222  (std::vector<int> & rowind,
223  std::vector<int> & colind,
224  std::vector<Treal> & values,
225  std::vector<int> const & rowInversePermutation,
226  std::vector<int> const & colInversePermutation) const {
227  std::vector<int> tmpRowind;
228  std::vector<int> tmpColind;
229  tmpRowind.reserve(rowind.capacity());
230  tmpColind.reserve(colind.capacity());
231  values.resize(0);
232  this->matrixPtr->getAllValues(tmpRowind, tmpColind, values);
234  getPermutedIndexes(tmpRowind, rowInversePermutation, rowind);
236  getPermutedIndexes(tmpColind, colInversePermutation, colind);
237 
238  }
248 #if 0
249  inline void fullmatrix(Treal* const full,
250  const int totnrows, const int totncols) const {
251  this->matrixPtr->fullmatrix(full, totnrows, totncols);
252  }
253 #endif
257  return *this;
258  }
261  if (mt.tA)
263  else
265  return *this;
266  }
267 
271  this->matrixPtr->symToNosym();
272  return *this;
273  }
277  return *this;
278  }
279 
281  *this->matrixPtr = k;
282  return *this;
283  }
284  inline Treal frob() const {
285  return this->matrixPtr->frob();
286  }
287  static inline Treal frob_diff
290  return Tmatrix::frobDiff(*A.matrixPtr, *B.matrixPtr);
291  }
292  Treal eucl(Treal const requestedAccuracy,
293  int maxIter = -1) const;
294 
295 
296  void thresh(Treal const threshold,
297  normType const norm);
298 
299  inline void frob_thresh(Treal threshold) {
300  this->matrixPtr->frob_thresh(threshold);
301  }
302  Treal eucl_thresh(Treal const threshold);
303 
304  inline void gershgorin(Treal& lmin, Treal& lmax) {
305  this->matrixPtr->gershgorin(lmin, lmax);
306  }
307  static inline Treal trace_ab
310  return Tmatrix::trace_ab(*A.matrixPtr, *B.matrixPtr);
311  }
312  static inline Treal trace_aTb
315  return Tmatrix::trace_aTb(*A.matrixPtr, *B.matrixPtr);
316  }
317  inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
318  return this->matrixPtr->nnz();
319  }
320  inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
321  return this->matrixPtr->nvalues();
322  }
323 
324  inline void write_to_buffer(void* buffer, const int n_bytes) const {
325  this->write_to_buffer_base(buffer, n_bytes, matrix_matr);
326  }
327  inline void read_from_buffer(void* buffer, const int n_bytes) {
328  this->read_from_buffer_base(buffer, n_bytes, matrix_matr);
329  }
330 
331  /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */
334  (const XYZ<Treal,
337 
342 
344  MatrixGeneral<Treal, Tmatrix>& operator+=
345  (const XYZ<Treal,
348 
351  (const XYZpUV<Treal,
354  Treal,
356 
360  MatrixGeneral<Treal, Tmatrix> > const & mpm);
362  MatrixGeneral<Treal, Tmatrix>& operator+=
364  MatrixGeneral<Treal, Tmatrix>& operator-=
367  MatrixGeneral<Treal, Tmatrix>& operator+=
369 
370 
371  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
374  (const XYZ<Treal,
382  MatrixGeneral<Treal, Tmatrix>& operator+=
383  (const XYZ<Treal,
388  (const XYZpUV<Treal,
391  Treal,
395  (const XYZ<Treal,
403  MatrixGeneral<Treal, Tmatrix>& operator+=
404  (const XYZ<Treal,
409  (const XYZpUV<Treal,
412  Treal,
416  (const XYZ<Treal,
424  MatrixGeneral<Treal, Tmatrix>& operator+=
425  (const XYZ<Treal,
430  (const XYZpUV<Treal,
433  Treal,
435 
436  /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */
439  (const XYZ<Treal,
444  (const XYZ<Treal,
447 
448  void random() {
449  this->matrixPtr->random();
450  }
451 
452  void randomZeroStructure(Treal probabilityBeingZero) {
453  this->matrixPtr->randomZeroStructure(probabilityBeingZero);
454  }
455 
456  template<typename TRule>
457  void setElementsByRule(TRule & rule) {
458  this->matrixPtr->setElementsByRule(rule);
459  return;
460  }
461 
462  std::string obj_type_id() const {return "MatrixGeneral";}
463  protected:
464  inline void writeToFileProt(std::ofstream & file) const {
465  this->writeToFileBase(file, matrix_matr);
466  }
467  inline void readFromFileProt(std::ifstream & file) {
468  this->readFromFileBase(file, matrix_matr);
469  }
470  private:
471 
472  };
473 
474  template<typename Treal, typename Tmatrix>
476  eucl(Treal const requestedAccuracy,
477  int maxIter) const {
478  VectorType guess;
480  this->getCols(cols);
481  guess.resetSizesAndBlocks(cols);
482  guess.rand();
484  if (maxIter < 0)
485  maxIter = this->get_nrows() * 100;
488  lan(ata, guess, maxIter);
489  lan.setRelTol( 100 * mat::getMachineEpsilon<Treal>() );
490  lan.run();
491  Treal eVal = 0;
492  Treal acc = 0;
493  lan.getLargestMagnitudeEig(eVal, acc);
494  Interval<Treal> euclInt( sqrt(eVal-acc),
495  sqrt(eVal+acc) );
496  if ( euclInt.low() < 0 )
497  euclInt = Interval<Treal>( 0, sqrt(eVal+acc) );
498  if ( euclInt.length() / 2.0 > requestedAccuracy ) {
499  std::cout << "req: " << requestedAccuracy
500  << " obt: " << euclInt.length() / 2.0 << std::endl;
501  throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
502  }
503  return euclInt.midPoint();
504  }
505 
506 
507  template<typename Treal, typename Tmatrix>
509  thresh(Treal const threshold,
510  normType const norm) {
511  switch (norm) {
512  case frobNorm:
513  this->frob_thresh(threshold);
514  break;
515  default:
516  throw Failure("MatrixGeneral<Treal, Tmatrix>::"
517  "thresh(Treal const, "
518  "normType const): "
519  "Thresholding not imlpemented for selected norm");
520  }
521  }
522 
523  template<typename Treal, typename Tmatrix>
525  eucl_thresh(Treal const threshold) {
526  EuclTruncationGeneral<MatrixGeneral<Treal, Tmatrix>, Treal> TruncObj( *this );
527  return TruncObj.run( threshold );
528  }
529 
530 
531  /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */
532  /* C = alpha * op(A) * op(B) */
533  template<typename Treal, typename Tmatrix>
536  (const XYZ<Treal,
539  assert(this != &smm.B && this != &smm.C);
540  this->matrixPtr.haveDataStructureSet(true);
541  Tmatrix::gemm(smm.tB, smm.tC, smm.A,
542  *smm.B.matrixPtr, *smm.C.matrixPtr,
543  0, *this->matrixPtr);
544  return *this;
545  }
546 
547  /* C = op(A) * op(B) */
548  template<typename Treal, typename Tmatrix>
553  assert(this != &mm.A && this != &mm.B);
554  Tmatrix::gemm(mm.tA, mm.tB, 1.0,
555  *mm.A.matrixPtr, *mm.B.matrixPtr,
556  0, *this->matrixPtr);
557  return *this;
558  }
559 
560  /* C += alpha * op(A) * op(B) */
561  template<typename Treal, typename Tmatrix>
564  (const XYZ<Treal,
567  assert(this != &smm.B && this != &smm.C);
568  Tmatrix::gemm(smm.tB, smm.tC, smm.A,
569  *smm.B.matrixPtr, *smm.C.matrixPtr,
570  1, *this->matrixPtr);
571  return *this;
572  }
573 
574  /* C = alpha * op(A) * op(B) + beta * C */
575  template<typename Treal, typename Tmatrix>
578  (const XYZpUV<Treal,
581  Treal,
582  MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
583  assert(this != &smmpsm.B && this != &smmpsm.C);
584  assert(!smmpsm.tE);
585  if (this == &smmpsm.E)
586  Tmatrix::gemm(smmpsm.tB, smmpsm.tC, smmpsm.A,
587  *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
588  smmpsm.D, *this->matrixPtr);
589  else
590  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
591  "(const XYZpUV<Treal, MatrixGeneral"
592  "<Treal, Tmatrix> >&) : D = alpha "
593  "* op(A) * op(B) + beta * C not supported for C != D");
594  return *this;
595  }
596 
597 
598  /* C = A + B */
599  template<typename Treal, typename Tmatrix>
604  assert(this != &mpm.A);
605  (*this) = mpm.B;
606  Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
607  return *this;
608  }
609  /* B += A */
610  template<typename Treal, typename Tmatrix>
614  Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
615  return *this;
616  }
617  /* B -= A */
618  template<typename Treal, typename Tmatrix>
622  Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
623  return *this;
624  }
625  /* B += alpha * A */
626  template<typename Treal, typename Tmatrix>
630  assert(!sm.tB);
631  Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
632  return *this;
633  }
634 
635 
636  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
637  /* C = alpha * A * B : A is symmetric */
638  template<typename Treal, typename Tmatrix>
641  (const XYZ<Treal,
644  assert(this != &smm.C);
645  assert(!smm.tB && !smm.tC);
646  this->matrixPtr.haveDataStructureSet(true);
647  Tmatrix::symm('L', 'U', smm.A,
648  *smm.B.matrixPtr, *smm.C.matrixPtr,
649  0, *this->matrixPtr);
650  return *this;
651  }
652 
653  /* C = A * B : A is symmetric */
654  template<typename Treal, typename Tmatrix>
659  assert(this != &mm.B);
660  assert(!mm.tB);
661  Tmatrix::symm('L', 'U', 1.0,
662  *mm.A.matrixPtr, *mm.B.matrixPtr,
663  0, *this->matrixPtr);
664  return *this;
665  }
666 
667  /* C += alpha * A * B : A is symmetric */
668  template<typename Treal, typename Tmatrix>
671  (const XYZ<Treal,
674  assert(this != &smm.C);
675  assert(!smm.tB && !smm.tC);
676  Tmatrix::symm('L', 'U', smm.A,
677  *smm.B.matrixPtr, *smm.C.matrixPtr,
678  1, *this->matrixPtr);
679  return *this;
680  }
681 
682  /* C = alpha * A * B + beta * C : A is symmetric */
683  template<typename Treal, typename Tmatrix>
686  (const XYZpUV<Treal,
689  Treal,
690  MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
691  assert(this != &smmpsm.C);
692  assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
693  if (this == &smmpsm.E)
694  Tmatrix::symm('L', 'U', smmpsm.A,
695  *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
696  smmpsm.D, *this->matrixPtr);
697  else
698  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
699  "(const XYZpUV<Treal, MatrixGeneral"
700  "<Treal, Tmatrix>, MatrixSymmetric<Treal, "
701  "Tmatrix>, Treal, MatrixGeneral"
702  "<Treal, Tmatrix> >&) "
703  ": D = alpha * A * B + beta * C (with A symmetric)"
704  " not supported for C != D");
705  return *this;
706  }
707 
708  /* C = alpha * B * A : A is symmetric */
709  template<typename Treal, typename Tmatrix>
712  (const XYZ<Treal,
715  assert(this != &smm.B);
716  assert(!smm.tB && !smm.tC);
717  this->matrixPtr.haveDataStructureSet(true);
718  Tmatrix::symm('R', 'U', smm.A,
719  *smm.C.matrixPtr, *smm.B.matrixPtr,
720  0, *this->matrixPtr);
721  return *this;
722  }
723 
724  /* C = B * A : A is symmetric */
725  template<typename Treal, typename Tmatrix>
730  assert(this != &mm.A);
731  assert(!mm.tA && !mm.tB);
732  Tmatrix::symm('R', 'U', 1.0,
733  *mm.B.matrixPtr, *mm.A.matrixPtr,
734  0, *this->matrixPtr);
735  return *this;
736  }
737 
738  /* C += alpha * B * A : A is symmetric */
739  template<typename Treal, typename Tmatrix>
742  (const XYZ<Treal,
745  assert(this != &smm.B);
746  assert(!smm.tB && !smm.tC);
747  Tmatrix::symm('R', 'U', smm.A,
748  *smm.C.matrixPtr, *smm.B.matrixPtr,
749  1, *this->matrixPtr);
750  return *this;
751  }
752 
753  /* C = alpha * B * A + beta * C : A is symmetric */
754  template<typename Treal, typename Tmatrix>
757  (const XYZpUV<Treal,
760  Treal,
761  MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
762  assert(this != &smmpsm.B);
763  assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
764  if (this == &smmpsm.E)
765  Tmatrix::symm('R', 'U', smmpsm.A,
766  *smmpsm.C.matrixPtr, *smmpsm.B.matrixPtr,
767  smmpsm.D, *this->matrixPtr);
768  else
769  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
770  "(const XYZpUV<Treal, MatrixSymmetric"
771  "<Treal, Tmatrix>, MatrixGeneral<Treal, "
772  "Tmatrix>, Treal, MatrixGeneral"
773  "<Treal, Tmatrix> >&) "
774  ": D = alpha * B * A + beta * C (with A symmetric)"
775  " not supported for C != D");
776  return *this;
777  }
778 
779 
781  template<typename Treal, typename Tmatrix>
784  (const XYZ<Treal,
787  assert(!smm.tB && !smm.tC);
788  this->matrixPtr.haveDataStructureSet(true);
789  Tmatrix::ssmm(smm.A,
790  *smm.B.matrixPtr,
791  *smm.C.matrixPtr,
792  0,
793  *this->matrixPtr);
794  return *this;
795  }
796 
798  template<typename Treal, typename Tmatrix>
803  assert(!mm.tB);
804  Tmatrix::ssmm(1.0,
805  *mm.A.matrixPtr,
806  *mm.B.matrixPtr,
807  0,
808  *this->matrixPtr);
809  return *this;
810  }
811 
813  template<typename Treal, typename Tmatrix>
816  (const XYZ<Treal,
819  assert(!smm.tB && !smm.tC);
820  Tmatrix::ssmm(smm.A,
821  *smm.B.matrixPtr,
822  *smm.C.matrixPtr,
823  1,
824  *this->matrixPtr);
825  return *this;
826  }
827 
828 
830  template<typename Treal, typename Tmatrix>
833  (const XYZpUV<Treal,
836  Treal,
837  MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
838  assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
839  if (this == &smmpsm.E)
840  Tmatrix::ssmm(smmpsm.A,
841  *smmpsm.B.matrixPtr,
842  *smmpsm.C.matrixPtr,
843  smmpsm.D,
844  *this->matrixPtr);
845  else
846  throw Failure("MatrixGeneral<Treal, Tmatrix>::"
847  "operator=(const XYZpUV<"
848  "Treal, MatrixSymmetric<Treal, Tmatrix>, "
849  "MatrixSymmetric<Treal, Tmatrix>, Treal, "
850  "MatrixGeneral<Treal, Tmatrix> >&) "
851  ": D = alpha * A * B + beta * C (with A and B symmetric)"
852  " not supported for C != D");
853  return *this;
854  }
855 
856 
857 
858  /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */
859 
860  /* B = alpha * op(A) * B : A is upper triangular */
861  template<typename Treal, typename Tmatrix>
864  (const XYZ<Treal,
867  assert(!smm.tC);
868  if (this == &smm.C)
869  Tmatrix::trmm('L', 'U', smm.tB, smm.A,
870  *smm.B.matrixPtr, *this->matrixPtr);
871  else
872  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
873  "(const XYZ<Treal, MatrixTriangular"
874  "<Treal, Tmatrix>, MatrixGeneral<Treal,"
875  " Tmatrix> >& : D = alpha * op(A) * B (with"
876  " A upper triangular) not supported for B != D");
877  return *this;
878  }
879 
880 
881  /* A = alpha * A * op(B) : B is upper triangular */
882  template<typename Treal, typename Tmatrix>
885  (const XYZ<Treal,
888  assert(!smm.tB);
889  if (this == &smm.B)
890  Tmatrix::trmm('R', 'U', smm.tC, smm.A,
891  *smm.C.matrixPtr, *this->matrixPtr);
892  else
893  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
894  "(const XYZ<Treal, MatrixGeneral"
895  "<Treal, Tmatrix>, MatrixTriangular<Treal,"
896  " Tmatrix> >& : D = alpha * A * op(B) (with"
897  " B upper triangular) not supported for A != D");
898  return *this;
899  }
900 
901 
902  /******* FUNCTIONS DECLARED OUTSIDE CLASS */
903  template<typename Treal, typename Tmatrix>
904  Treal trace(const XYZ<Treal,
907  if ((!smm.tB && !smm.tC) || (smm.tB && smm.tC))
908  return smm.A * MatrixGeneral<Treal, Tmatrix>::
909  trace_ab(smm.B, smm.C);
910  else if (smm.tB)
911  return smm.A * MatrixGeneral<Treal, Tmatrix>::
912  trace_aTb(smm.B, smm.C);
913  else
914  return smm.A * MatrixGeneral<Treal, Tmatrix>::
915  trace_aTb(smm.C, smm.B);
916  }
917 
918 
919 } /* end namespace mat */
920 #endif
921 
922 
#define A
Normal matrix.
Definition: MatrixBase.h:49
normalMatrix MatrixGeneral
Definition: random_matrices.h:70
mat::SizesAndBlocks cols
Definition: test.cc:52
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:166
Treal low() const
Definition: Interval.h:144
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixGeneral.h:61
static Treal trace_ab(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition: MatrixGeneral.h:308
ValidPtr< Tmatrix > matrixPtr
Definition: MatrixBase.h:153
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixGeneral.h:476
Definition: LanczosLargestMagnitudeEig.h:46
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition: MatrixBase.h:281
MatrixGeneral(const MatrixGeneral< Treal, Tmatrix > &matr)
Copy constructor.
Definition: MatrixGeneral.h:65
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
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
MatrixGeneral< Treal, Tmatrix > & operator=(const Xtrans< MatrixGeneral< Treal, Tmatrix > > &mt)
Definition: MatrixGeneral.h:260
void assign_from_sparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Assign from sparse matrix given by three arrays.
Definition: MatrixGeneral.h:122
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixGeneral< Treal, Tmatrix > &mat)
Definition: MatrixGeneral.h:255
static Treal frob_diff(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition: MatrixGeneral.h:288
void gershgorin(Treal &lmin, Treal &lmax)
Definition: MatrixGeneral.h:304
Truncation of general matrices.
Definition: truncation.h:272
This proxy expresses the result of addition of two objects, of possibly different types...
Definition: matrix_proxy.h:246
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixSymmetric< Treal, Tmatrix > &symm)
Definition: MatrixGeneral.h:269
Definition: MatrixBase.h:55
Definition: allocate.cc:39
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
Treal frob() const
Definition: MatrixGeneral.h:284
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixTriangular< Treal, Tmatrix > &triang)
Definition: MatrixGeneral.h:275
Definition: Interval.h:46
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: MatrixGeneral.h:90
size_t nnz() const
Definition: MatrixGeneral.h:317
Definition: MatrixBase.h:56
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition: MatrixBase.h:236
MatrixGeneral(const MatrixSymmetric< Treal, Tmatrix > &symm)
Copy from symmetric matrix constructor.
Definition: MatrixGeneral.h:67
int get_ncols() const
Definition: MatrixBase.h:141
This proxy expresses the result of multiplication of three objects, of possibly different types...
Definition: matrix_proxy.h:67
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition: MatrixBase.h:254
Upper non-unit triangular matrix.
Definition: MatrixBase.h:53
This proxy expresses the result of transposition of an object of type TX.
Definition: matrix_proxy.h:118
Treal run(Treal const threshold)
Definition: truncation.h:80
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition: MatrixBase.h:221
Treal midPoint() const
Definition: Interval.h:115
void setElementsByRule(TRule &rule)
Definition: MatrixGeneral.h:457
void thresh(Treal const threshold, normType const norm)
Definition: MatrixGeneral.h:509
void rand()
Definition: VectorGeneral.h:102
size_t nvalues() const
Definition: MatrixGeneral.h:320
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:76
void randomZeroStructure(Treal probabilityBeingZero)
Definition: MatrixGeneral.h:452
void get_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Get values given by row and column index lists.
Definition: MatrixGeneral.h:173
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
std::string obj_type_id() const
Definition: MatrixGeneral.h:462
Definition: mat_utils.h:78
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
Treal eucl_thresh(Treal const threshold)
Definition: MatrixGeneral.h:525
void random()
Definition: MatrixGeneral.h:448
Base class for matrix API.
Definition: MatrixBase.h:69
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:58
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: MatrixGeneral.h:464
Base class for matrix API.
MatrixGeneral()
Default constructor.
Definition: MatrixGeneral.h:63
Treal length() const
Returns the length of the interval.
Definition: Interval.h:109
Definition: matInclude.h:139
Definition: Failure.h:57
This proxy expresses the result of multiplication of two objects, of possibly different types...
Definition: matrix_proxy.h:51
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixGeneral.h:327
#define B
MatrixGeneral(const MatrixTriangular< Treal, Tmatrix > &triang)
Copy from triangular matrix constructor.
Definition: MatrixGeneral.h:71
void frob_thresh(Treal threshold)
Definition: MatrixGeneral.h:299
int get_nrows() const
Definition: MatrixBase.h:138
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
void get_all_values(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &values) const
Get all values and corresponding row and column index lists, in matrix.
Definition: MatrixGeneral.h:204
Symmetric matrix.
Definition: MatrixBase.h:51
Treal trace(const XYZ< Treal, MatrixGeneral< Treal, Tmatrix >, MatrixGeneral< Treal, Tmatrix > > &smm)
Definition: MatrixGeneral.h:904
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition: matrix_proxy.h:88
MatrixGeneral< Treal, Tmatrix > & operator=(int const k)
Definition: MatrixGeneral.h:280
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: MatrixGeneral.h:467
static Treal trace_aTb(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition: MatrixGeneral.h:313
void fullMatrix(std::vector< Treal > &fullMat) const
Definition: MatrixGeneral.h:95
normType
Definition: matInclude.h:139
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixGeneral.h:324
static void getPermutedIndexes(std::vector< int > const &index, std::vector< int > const &permutation, std::vector< int > &newIndex)
Definition: MatrixBase.h:205