ergo
MatrixGeneral.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
36 #ifndef MAT_MATRIXGENERAL
37 #define MAT_MATRIXGENERAL
38 #include "MatrixBase.h"
39 namespace mat {
56  template<typename Treal, typename Tmatrix>
57  class MatrixGeneral : public MatrixBase<Treal, Tmatrix> {
58  public:
60 
62  :MatrixBase<Treal, Tmatrix>() {}
64  :MatrixBase<Treal, Tmatrix>(matr) {}
66  :MatrixBase<Treal, Tmatrix>(symm) {
67  this->matrixPtr->symToNosym();
68  }
70  :MatrixBase<Treal, Tmatrix>(triang) {}
73 #if 0
74  template<typename Tfull>
75  inline void assign_from_full
76  (Tfull const* const fullmatrix,
77  const int totnrows, const int totncols) {
78  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
79  }
80  inline void assign_from_full
81  (Treal const* const fullmatrix,
82  const int totnrows, const int totncols) {
83  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
84  }
85 #endif
86 
87  inline void assignFromFull
88  (std::vector<Treal> const & fullMat) {
89  assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
90  this->matrixPtr->assignFromFull(fullMat);
91  }
92 
93  inline void fullMatrix(std::vector<Treal> & fullMat) const {
94  this->matrixPtr->fullMatrix(fullMat);
95  }
96 
97  inline void fullMatrix
98  (std::vector<Treal> & fullMat,
99  std::vector<int> const & rowInversePermutation,
100  std::vector<int> const & colInversePermutation) const {
101  std::vector<int> rowind;
102  std::vector<int> colind;
103  std::vector<Treal> values;
104  get_all_values(rowind, colind, values,
105  rowInversePermutation,
106  colInversePermutation);
107  fullMat.assign(this->get_nrows()*this->get_ncols(),0);
108  for (unsigned int ind = 0; ind < values.size(); ++ind)
109  fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
110  values[ind];
111  }
119  inline void assign_from_sparse
120  (std::vector<int> const & rowind,
121  std::vector<int> const & colind,
122  std::vector<Treal> const & values,
123  SizesAndBlocks const & newRows,
124  SizesAndBlocks const & newCols) {
125  this->resetSizesAndBlocks(newRows, newCols);
126  this->matrixPtr->assignFromSparse(rowind, colind, values);
127  }
135  inline void assign_from_sparse
136  (std::vector<int> const & rowind,
137  std::vector<int> const & colind,
138  std::vector<Treal> const & values,
139  std::vector<int> const & rowPermutation,
140  std::vector<int> const & colPermutation) {
141  std::vector<int> newRowind;
142  std::vector<int> newColind;
144  getPermutedIndexes(rowind, rowPermutation, newRowind);
146  getPermutedIndexes(colind, colPermutation, newColind);
147  this->matrixPtr->assignFromSparse(newRowind, newColind, values);
148  }
154  inline void assign_from_sparse
155  (std::vector<int> const & rowind,
156  std::vector<int> const & colind,
157  std::vector<Treal> const & values,
158  SizesAndBlocks const & newRows,
159  SizesAndBlocks const & newCols,
160  std::vector<int> const & rowPermutation,
161  std::vector<int> const & colPermutation) {
162  this->resetSizesAndBlocks(newRows, newCols);
163  assign_from_sparse(rowind, colind, values,
164  rowPermutation, colPermutation);
165  }
170  inline void get_values
171  (std::vector<int> const & rowind,
172  std::vector<int> const & colind,
173  std::vector<Treal> & values) const {
174  this->matrixPtr->getValues(rowind, colind, values);
175  }
183  inline void get_values
184  (std::vector<int> const & rowind,
185  std::vector<int> const & colind,
186  std::vector<Treal> & values,
187  std::vector<int> const & rowPermutation,
188  std::vector<int> const & colPermutation) const {
189  std::vector<int> newRowind;
190  std::vector<int> newColind;
192  getPermutedIndexes(rowind, rowPermutation, newRowind);
194  getPermutedIndexes(colind, colPermutation, newColind);
195  this->matrixPtr->getValues(newRowind, newColind, values);
196  }
201  inline void get_all_values
202  (std::vector<int> & rowind,
203  std::vector<int> & colind,
204  std::vector<Treal> & values) const {
205  rowind.resize(0);
206  colind.resize(0);
207  values.resize(0);
208  this->matrixPtr->getAllValues(rowind, colind, values);
209  }
219  inline void get_all_values
220  (std::vector<int> & rowind,
221  std::vector<int> & colind,
222  std::vector<Treal> & values,
223  std::vector<int> const & rowInversePermutation,
224  std::vector<int> const & colInversePermutation) const {
225  std::vector<int> tmpRowind;
226  std::vector<int> tmpColind;
227  tmpRowind.reserve(rowind.capacity());
228  tmpColind.reserve(colind.capacity());
229  values.resize(0);
230  this->matrixPtr->getAllValues(tmpRowind, tmpColind, values);
232  getPermutedIndexes(tmpRowind, rowInversePermutation, rowind);
234  getPermutedIndexes(tmpColind, colInversePermutation, colind);
235 
236  }
246 #if 0
247  inline void fullmatrix(Treal* const full,
248  const int totnrows, const int totncols) const {
249  this->matrixPtr->fullmatrix(full, totnrows, totncols);
250  }
251 #endif
255  return *this;
256  }
259  if (mt.tA)
261  else
263  return *this;
264  }
265 
269  this->matrixPtr->symToNosym();
270  return *this;
271  }
275  return *this;
276  }
277 
279  *this->matrixPtr = k;
280  return *this;
281  }
282  inline Treal frob() const {
283  return this->matrixPtr->frob();
284  }
285  static inline Treal frob_diff
288  return Tmatrix::frobDiff(*A.matrixPtr, *B.matrixPtr);
289  }
290  Treal eucl(Treal const requestedAccuracy,
291  int maxIter = -1) const;
292 
293 
294  void thresh(Treal const threshold,
295  normType const norm);
296 
297  inline void frob_thresh(Treal threshold) {
298  this->matrixPtr->frob_thresh(threshold);
299  }
300  Treal eucl_thresh(Treal const threshold);
301 
302  inline void gersgorin(Treal& lmin, Treal& lmax) {
303  this->matrixPtr->gersgorin(lmin, lmax);
304  }
305  static inline Treal trace_ab
308  return Tmatrix::trace_ab(*A.matrixPtr, *B.matrixPtr);
309  }
310  static inline Treal trace_aTb
313  return Tmatrix::trace_aTb(*A.matrixPtr, *B.matrixPtr);
314  }
315  inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
316  return this->matrixPtr->nnz();
317  }
318  inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
319  return this->matrixPtr->nvalues();
320  }
321 
322  inline void write_to_buffer(void* buffer, const int n_bytes) const {
323  this->write_to_buffer_base(buffer, n_bytes, matrix_matr);
324  }
325  inline void read_from_buffer(void* buffer, const int n_bytes) {
326  this->read_from_buffer_base(buffer, n_bytes, matrix_matr);
327  }
328 
329  /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */
332  (const XYZ<Treal,
334  MatrixGeneral<Treal, Tmatrix> >& smm);
335 
337  MatrixGeneral<Treal, Tmatrix>& operator=
339  MatrixGeneral<Treal, Tmatrix> >& mm);
340 
342  MatrixGeneral<Treal, Tmatrix>& operator+=
343  (const XYZ<Treal,
344  MatrixGeneral<Treal, Tmatrix>,
345  MatrixGeneral<Treal, Tmatrix> >& smm);
346 
348  MatrixGeneral<Treal, Tmatrix>& operator=
349  (const XYZpUV<Treal,
350  MatrixGeneral<Treal, Tmatrix>,
351  MatrixGeneral<Treal, Tmatrix>,
352  Treal,
353  MatrixGeneral<Treal, Tmatrix> >& smmpsm);
354 
356  MatrixGeneral<Treal, Tmatrix>& operator=
358  MatrixGeneral<Treal, Tmatrix> > const & mpm);
360  MatrixGeneral<Treal, Tmatrix>& operator+=
361  (MatrixGeneral<Treal, Tmatrix> const & A);
362  MatrixGeneral<Treal, Tmatrix>& operator-=
363  (MatrixGeneral<Treal, Tmatrix> const & A);
365  MatrixGeneral<Treal, Tmatrix>& operator+=
367 
368 
369  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
371  MatrixGeneral<Treal, Tmatrix>& operator=
372  (const XYZ<Treal,
374  MatrixGeneral<Treal, Tmatrix> >& smm);
376  MatrixGeneral<Treal, Tmatrix>& operator=
378  MatrixGeneral<Treal, Tmatrix> >& mm);
380  MatrixGeneral<Treal, Tmatrix>& operator+=
381  (const XYZ<Treal,
382  MatrixSymmetric<Treal, Tmatrix>,
383  MatrixGeneral<Treal, Tmatrix> >& smm);
385  MatrixGeneral<Treal, Tmatrix>& operator=
386  (const XYZpUV<Treal,
387  MatrixSymmetric<Treal, Tmatrix>,
388  MatrixGeneral<Treal, Tmatrix>,
389  Treal,
390  MatrixGeneral<Treal, Tmatrix> >& smmpsm);
392  MatrixGeneral<Treal, Tmatrix>& operator=
393  (const XYZ<Treal,
394  MatrixGeneral<Treal, Tmatrix>,
395  MatrixSymmetric<Treal, Tmatrix> >& smm);
397  MatrixGeneral<Treal, Tmatrix>& operator=
399  MatrixSymmetric<Treal, Tmatrix> >& mm);
401  MatrixGeneral<Treal, Tmatrix>& operator+=
402  (const XYZ<Treal,
403  MatrixGeneral<Treal, Tmatrix>,
404  MatrixSymmetric<Treal, Tmatrix> >& smm);
406  MatrixGeneral<Treal, Tmatrix>& operator=
407  (const XYZpUV<Treal,
408  MatrixGeneral<Treal, Tmatrix>,
409  MatrixSymmetric<Treal, Tmatrix>,
410  Treal,
411  MatrixGeneral<Treal, Tmatrix> >& smmpsm);
413  MatrixGeneral<Treal, Tmatrix>& operator=
414  (const XYZ<Treal,
415  MatrixSymmetric<Treal, Tmatrix>,
416  MatrixSymmetric<Treal, Tmatrix> >& smm);
418  MatrixGeneral<Treal, Tmatrix>& operator=
420  MatrixSymmetric<Treal, Tmatrix> >& mm);
422  MatrixGeneral<Treal, Tmatrix>& operator+=
423  (const XYZ<Treal,
424  MatrixSymmetric<Treal, Tmatrix>,
425  MatrixSymmetric<Treal, Tmatrix> >& smm);
427  MatrixGeneral<Treal, Tmatrix>& operator=
428  (const XYZpUV<Treal,
429  MatrixSymmetric<Treal, Tmatrix>,
430  MatrixSymmetric<Treal, Tmatrix>,
431  Treal,
432  MatrixGeneral<Treal, Tmatrix> >& smmpsm);
433 
434  /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */
436  MatrixGeneral<Treal, Tmatrix>& operator=
437  (const XYZ<Treal,
439  MatrixGeneral<Treal, Tmatrix> >& smm);
441  MatrixGeneral<Treal, Tmatrix>& operator=
442  (const XYZ<Treal,
443  MatrixGeneral<Treal, Tmatrix>,
444  MatrixTriangular<Treal, Tmatrix> >& smm);
445 
446  void random() {
447  this->matrixPtr->random();
448  }
449 
450  void randomZeroStructure(Treal probabilityBeingZero) {
451  this->matrixPtr->randomZeroStructure(probabilityBeingZero);
452  }
453 
454  template<typename TRule>
455  void setElementsByRule(TRule & rule) {
456  this->matrixPtr->setElementsByRule(rule);
457  return;
458  }
459 
460  std::string obj_type_id() const {return "MatrixGeneral";}
461  protected:
462  inline void writeToFileProt(std::ofstream & file) const {
463  this->writeToFileBase(file, matrix_matr);
464  }
465  inline void readFromFileProt(std::ifstream & file) {
466  this->readFromFileBase(file, matrix_matr);
467  }
468  private:
469 
470  };
471 
472  template<typename Treal, typename Tmatrix>
474  eucl(Treal const requestedAccuracy,
475  int maxIter) const {
476  VectorType guess;
477  SizesAndBlocks cols;
478  this->getCols(cols);
479  guess.resetSizesAndBlocks(cols);
480  guess.rand();
482  if (maxIter < 0)
483  maxIter = this->get_nrows() * 100;
485  <Treal, ATAMatrix<MatrixGeneral<Treal, Tmatrix>, Treal>, VectorType>
486  lan(ata, guess, maxIter);
487  lan.setRelTol( 100 * std::numeric_limits<Treal>::epsilon() );
488  lan.run();
489  Treal eVal = 0;
490  Treal acc = 0;
491  lan.getLargestMagnitudeEig(eVal, acc);
492  Interval<Treal> euclInt( sqrt(eVal-acc),
493  sqrt(eVal+acc) );
494  if ( euclInt.low() < 0 )
495  euclInt = Interval<Treal>( 0, sqrt(eVal+acc) );
496  if ( euclInt.length() / 2.0 > requestedAccuracy ) {
497  std::cout << "req: " << requestedAccuracy
498  << " obt: " << euclInt.length() / 2.0 << std::endl;
499  throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
500  }
501  return euclInt.midPoint();
502  }
503 
504 
505  template<typename Treal, typename Tmatrix>
507  thresh(Treal const threshold,
508  normType const norm) {
509  switch (norm) {
510  case frobNorm:
511  this->frob_thresh(threshold);
512  break;
513  default:
514  throw Failure("MatrixGeneral<Treal, Tmatrix>::"
515  "thresh(Treal const, "
516  "normType const): "
517  "Thresholding not imlpemented for selected norm");
518  }
519  }
520 
521  template<typename Treal, typename Tmatrix>
523  eucl_thresh(Treal const threshold) {
524  EuclTruncationGeneral<MatrixGeneral<Treal, Tmatrix>, Treal> TruncObj( *this );
525  return TruncObj.run( threshold );
526  }
527 
528 
529  /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */
530  /* C = alpha * op(A) * op(B) */
531  template<typename Treal, typename Tmatrix>
534  (const XYZ<Treal,
536  MatrixGeneral<Treal, Tmatrix> >& smm) {
537  assert(this != &smm.B && this != &smm.C);
538  this->matrixPtr.haveDataStructureSet(true);
539  Tmatrix::gemm(smm.tB, smm.tC, smm.A,
540  *smm.B.matrixPtr, *smm.C.matrixPtr,
541  0, *this->matrixPtr);
542  return *this;
543  }
544 
545  /* C = op(A) * op(B) */
546  template<typename Treal, typename Tmatrix>
547  MatrixGeneral<Treal, Tmatrix>&
548  MatrixGeneral<Treal, Tmatrix>::operator=
550  MatrixGeneral<Treal, Tmatrix> >& mm) {
551  assert(this != &mm.A && this != &mm.B);
552  Tmatrix::gemm(mm.tA, mm.tB, 1.0,
553  *mm.A.matrixPtr, *mm.B.matrixPtr,
554  0, *this->matrixPtr);
555  return *this;
556  }
557 
558  /* C += alpha * op(A) * op(B) */
559  template<typename Treal, typename Tmatrix>
560  MatrixGeneral<Treal, Tmatrix>&
561  MatrixGeneral<Treal, Tmatrix>::operator+=
562  (const XYZ<Treal,
563  MatrixGeneral<Treal, Tmatrix>,
564  MatrixGeneral<Treal, Tmatrix> >& smm) {
565  assert(this != &smm.B && this != &smm.C);
566  Tmatrix::gemm(smm.tB, smm.tC, smm.A,
567  *smm.B.matrixPtr, *smm.C.matrixPtr,
568  1, *this->matrixPtr);
569  return *this;
570  }
571 
572  /* C = alpha * op(A) * op(B) + beta * C */
573  template<typename Treal, typename Tmatrix>
574  MatrixGeneral<Treal, Tmatrix>&
575  MatrixGeneral<Treal, Tmatrix>::operator=
576  (const XYZpUV<Treal,
577  MatrixGeneral<Treal, Tmatrix>,
578  MatrixGeneral<Treal, Tmatrix>,
579  Treal,
580  MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
581  assert(this != &smmpsm.B && this != &smmpsm.C);
582  assert(!smmpsm.tE);
583  if (this == &smmpsm.E)
584  Tmatrix::gemm(smmpsm.tB, smmpsm.tC, smmpsm.A,
585  *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
586  smmpsm.D, *this->matrixPtr);
587  else
588  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
589  "(const XYZpUV<Treal, MatrixGeneral"
590  "<Treal, Tmatrix> >&) : D = alpha "
591  "* op(A) * op(B) + beta * C not supported for C != D");
592  return *this;
593  }
594 
595 
596  /* C = A + B */
597  template<typename Treal, typename Tmatrix>
598  inline MatrixGeneral<Treal, Tmatrix>&
599  MatrixGeneral<Treal, Tmatrix>::operator=
601  MatrixGeneral<Treal, Tmatrix> >& mpm) {
602  assert(this != &mpm.A);
603  (*this) = mpm.B;
604  Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
605  return *this;
606  }
607  /* B += A */
608  template<typename Treal, typename Tmatrix>
609  inline MatrixGeneral<Treal, Tmatrix>&
610  MatrixGeneral<Treal, Tmatrix>::operator+=
611  (MatrixGeneral<Treal, Tmatrix> const & A) {
612  Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
613  return *this;
614  }
615  /* B -= A */
616  template<typename Treal, typename Tmatrix>
617  inline MatrixGeneral<Treal, Tmatrix>&
618  MatrixGeneral<Treal, Tmatrix>::operator-=
619  (MatrixGeneral<Treal, Tmatrix> const & A) {
620  Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
621  return *this;
622  }
623  /* B += alpha * A */
624  template<typename Treal, typename Tmatrix>
625  inline MatrixGeneral<Treal, Tmatrix>&
626  MatrixGeneral<Treal, Tmatrix>::operator+=
628  assert(!sm.tB);
629  Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
630  return *this;
631  }
632 
633 
634  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
635  /* C = alpha * A * B : A is symmetric */
636  template<typename Treal, typename Tmatrix>
637  MatrixGeneral<Treal, Tmatrix>&
638  MatrixGeneral<Treal, Tmatrix>::operator=
639  (const XYZ<Treal,
641  MatrixGeneral<Treal, Tmatrix> >& smm) {
642  assert(this != &smm.C);
643  assert(!smm.tB && !smm.tC);
644  this->matrixPtr.haveDataStructureSet(true);
645  Tmatrix::symm('L', 'U', smm.A,
646  *smm.B.matrixPtr, *smm.C.matrixPtr,
647  0, *this->matrixPtr);
648  return *this;
649  }
650 
651  /* C = A * B : A is symmetric */
652  template<typename Treal, typename Tmatrix>
653  MatrixGeneral<Treal, Tmatrix>&
654  MatrixGeneral<Treal, Tmatrix>::operator=
656  MatrixGeneral<Treal, Tmatrix> >& mm) {
657  assert(this != &mm.B);
658  assert(!mm.tB);
659  Tmatrix::symm('L', 'U', 1.0,
660  *mm.A.matrixPtr, *mm.B.matrixPtr,
661  0, *this->matrixPtr);
662  return *this;
663  }
664 
665  /* C += alpha * A * B : A is symmetric */
666  template<typename Treal, typename Tmatrix>
667  MatrixGeneral<Treal, Tmatrix>&
668  MatrixGeneral<Treal, Tmatrix>::operator+=
669  (const XYZ<Treal,
670  MatrixSymmetric<Treal, Tmatrix>,
671  MatrixGeneral<Treal, Tmatrix> >& smm) {
672  assert(this != &smm.C);
673  assert(!smm.tB && !smm.tC);
674  Tmatrix::symm('L', 'U', smm.A,
675  *smm.B.matrixPtr, *smm.C.matrixPtr,
676  1, *this->matrixPtr);
677  return *this;
678  }
679 
680  /* C = alpha * A * B + beta * C : A is symmetric */
681  template<typename Treal, typename Tmatrix>
682  MatrixGeneral<Treal, Tmatrix>&
683  MatrixGeneral<Treal, Tmatrix>::operator=
684  (const XYZpUV<Treal,
685  MatrixSymmetric<Treal, Tmatrix>,
686  MatrixGeneral<Treal, Tmatrix>,
687  Treal,
688  MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
689  assert(this != &smmpsm.C);
690  assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
691  if (this == &smmpsm.E)
692  Tmatrix::symm('L', 'U', smmpsm.A,
693  *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
694  smmpsm.D, *this->matrixPtr);
695  else
696  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
697  "(const XYZpUV<Treal, MatrixGeneral"
698  "<Treal, Tmatrix>, MatrixSymmetric<Treal, "
699  "Tmatrix>, Treal, MatrixGeneral"
700  "<Treal, Tmatrix> >&) "
701  ": D = alpha * A * B + beta * C (with A symmetric)"
702  " not supported for C != D");
703  return *this;
704  }
705 
706  /* C = alpha * B * A : A is symmetric */
707  template<typename Treal, typename Tmatrix>
708  MatrixGeneral<Treal, Tmatrix>&
709  MatrixGeneral<Treal, Tmatrix>::operator=
710  (const XYZ<Treal,
711  MatrixGeneral<Treal, Tmatrix>,
712  MatrixSymmetric<Treal, Tmatrix> >& smm) {
713  assert(this != &smm.B);
714  assert(!smm.tB && !smm.tC);
715  this->matrixPtr.haveDataStructureSet(true);
716  Tmatrix::symm('R', 'U', smm.A,
717  *smm.C.matrixPtr, *smm.B.matrixPtr,
718  0, *this->matrixPtr);
719  return *this;
720  }
721 
722  /* C = B * A : A is symmetric */
723  template<typename Treal, typename Tmatrix>
724  MatrixGeneral<Treal, Tmatrix>&
725  MatrixGeneral<Treal, Tmatrix>::operator=
727  MatrixSymmetric<Treal, Tmatrix> >& mm) {
728  assert(this != &mm.A);
729  assert(!mm.tA && !mm.tB);
730  Tmatrix::symm('R', 'U', 1.0,
731  *mm.B.matrixPtr, *mm.A.matrixPtr,
732  0, *this->matrixPtr);
733  return *this;
734  }
735 
736  /* C += alpha * B * A : A is symmetric */
737  template<typename Treal, typename Tmatrix>
738  MatrixGeneral<Treal, Tmatrix>&
739  MatrixGeneral<Treal, Tmatrix>::operator+=
740  (const XYZ<Treal,
741  MatrixGeneral<Treal, Tmatrix>,
742  MatrixSymmetric<Treal, Tmatrix> >& smm) {
743  assert(this != &smm.B);
744  assert(!smm.tB && !smm.tC);
745  Tmatrix::symm('R', 'U', smm.A,
746  *smm.C.matrixPtr, *smm.B.matrixPtr,
747  1, *this->matrixPtr);
748  return *this;
749  }
750 
751  /* C = alpha * B * A + beta * C : A is symmetric */
752  template<typename Treal, typename Tmatrix>
753  MatrixGeneral<Treal, Tmatrix>&
754  MatrixGeneral<Treal, Tmatrix>::operator=
755  (const XYZpUV<Treal,
756  MatrixGeneral<Treal, Tmatrix>,
757  MatrixSymmetric<Treal, Tmatrix>,
758  Treal,
759  MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
760  assert(this != &smmpsm.B);
761  assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
762  if (this == &smmpsm.E)
763  Tmatrix::symm('R', 'U', smmpsm.A,
764  *smmpsm.C.matrixPtr, *smmpsm.B.matrixPtr,
765  smmpsm.D, *this->matrixPtr);
766  else
767  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
768  "(const XYZpUV<Treal, MatrixSymmetric"
769  "<Treal, Tmatrix>, MatrixGeneral<Treal, "
770  "Tmatrix>, Treal, MatrixGeneral"
771  "<Treal, Tmatrix> >&) "
772  ": D = alpha * B * A + beta * C (with A symmetric)"
773  " not supported for C != D");
774  return *this;
775  }
776 
777 
779  template<typename Treal, typename Tmatrix>
780  MatrixGeneral<Treal, Tmatrix>&
781  MatrixGeneral<Treal, Tmatrix>::operator=
782  (const XYZ<Treal,
783  MatrixSymmetric<Treal, Tmatrix>,
784  MatrixSymmetric<Treal, Tmatrix> >& smm) {
785  assert(!smm.tB && !smm.tC);
786  this->matrixPtr.haveDataStructureSet(true);
787  Tmatrix::ssmm(smm.A,
788  *smm.B.matrixPtr,
789  *smm.C.matrixPtr,
790  0,
791  *this->matrixPtr);
792  return *this;
793  }
794 
796  template<typename Treal, typename Tmatrix>
797  MatrixGeneral<Treal, Tmatrix>&
798  MatrixGeneral<Treal, Tmatrix>::operator=
800  MatrixSymmetric<Treal, Tmatrix> >& mm) {
801  assert(!mm.tB);
802  Tmatrix::ssmm(1.0,
803  *mm.A.matrixPtr,
804  *mm.B.matrixPtr,
805  0,
806  *this->matrixPtr);
807  return *this;
808  }
809 
811  template<typename Treal, typename Tmatrix>
812  MatrixGeneral<Treal, Tmatrix>&
813  MatrixGeneral<Treal, Tmatrix>::operator+=
814  (const XYZ<Treal,
815  MatrixSymmetric<Treal, Tmatrix>,
816  MatrixSymmetric<Treal, Tmatrix> >& smm) {
817  assert(!smm.tB && !smm.tC);
818  Tmatrix::ssmm(smm.A,
819  *smm.B.matrixPtr,
820  *smm.C.matrixPtr,
821  1,
822  *this->matrixPtr);
823  return *this;
824  }
825 
826 
828  template<typename Treal, typename Tmatrix>
829  MatrixGeneral<Treal, Tmatrix>&
830  MatrixGeneral<Treal, Tmatrix>::operator=
831  (const XYZpUV<Treal,
832  MatrixSymmetric<Treal, Tmatrix>,
833  MatrixSymmetric<Treal, Tmatrix>,
834  Treal,
835  MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
836  assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
837  if (this == &smmpsm.E)
838  Tmatrix::ssmm(smmpsm.A,
839  *smmpsm.B.matrixPtr,
840  *smmpsm.C.matrixPtr,
841  smmpsm.D,
842  *this->matrixPtr);
843  else
844  throw Failure("MatrixGeneral<Treal, Tmatrix>::"
845  "operator=(const XYZpUV<"
846  "Treal, MatrixSymmetric<Treal, Tmatrix>, "
847  "MatrixSymmetric<Treal, Tmatrix>, Treal, "
848  "MatrixGeneral<Treal, Tmatrix> >&) "
849  ": D = alpha * A * B + beta * C (with A and B symmetric)"
850  " not supported for C != D");
851  return *this;
852  }
853 
854 
855 
856  /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */
857 
858  /* B = alpha * op(A) * B : A is upper triangular */
859  template<typename Treal, typename Tmatrix>
860  MatrixGeneral<Treal, Tmatrix>&
861  MatrixGeneral<Treal, Tmatrix>::operator=
862  (const XYZ<Treal,
864  MatrixGeneral<Treal, Tmatrix> >& smm) {
865  assert(!smm.tC);
866  if (this == &smm.C)
867  Tmatrix::trmm('L', 'U', smm.tB, smm.A,
868  *smm.B.matrixPtr, *this->matrixPtr);
869  else
870  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
871  "(const XYZ<Treal, MatrixTriangular"
872  "<Treal, Tmatrix>, MatrixGeneral<Treal,"
873  " Tmatrix> >& : D = alpha * op(A) * B (with"
874  " A upper triangular) not supported for B != D");
875  return *this;
876  }
877 
878 
879  /* A = alpha * A * op(B) : B is upper triangular */
880  template<typename Treal, typename Tmatrix>
881  MatrixGeneral<Treal, Tmatrix>&
882  MatrixGeneral<Treal, Tmatrix>::operator=
883  (const XYZ<Treal,
884  MatrixGeneral<Treal, Tmatrix>,
885  MatrixTriangular<Treal, Tmatrix> >& smm) {
886  assert(!smm.tB);
887  if (this == &smm.B)
888  Tmatrix::trmm('R', 'U', smm.tC, smm.A,
889  *smm.C.matrixPtr, *this->matrixPtr);
890  else
891  throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
892  "(const XYZ<Treal, MatrixGeneral"
893  "<Treal, Tmatrix>, MatrixTriangular<Treal,"
894  " Tmatrix> >& : D = alpha * A * op(B) (with"
895  " B upper triangular) not supported for A != D");
896  return *this;
897  }
898 
899 
900  /******* FUNCTIONS DECLARED OUTSIDE CLASS */
901  template<typename Treal, typename Tmatrix>
902  Treal trace(const XYZ<Treal,
903  MatrixGeneral<Treal, Tmatrix>,
904  MatrixGeneral<Treal, Tmatrix> >& smm) {
905  if ((!smm.tB && !smm.tC) || (smm.tB && smm.tC))
906  return smm.A * MatrixGeneral<Treal, Tmatrix>::
907  trace_ab(smm.B, smm.C);
908  else if (smm.tB)
909  return smm.A * MatrixGeneral<Treal, Tmatrix>::
910  trace_aTb(smm.B, smm.C);
911  else
912  return smm.A * MatrixGeneral<Treal, Tmatrix>::
913  trace_aTb(smm.C, smm.B);
914  }
915 
916 
917 } /* end namespace mat */
918 #endif
919 
920 
#define A
void gersgorin(Treal &lmin, Treal &lmax)
Definition: MatrixGeneral.h:302
Normal matrix.
Definition: MatrixBase.h:47
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:164
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixGeneral.h:59
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixGeneral.h:474
static Treal trace_ab(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition: MatrixGeneral.h:306
ValidPtr< Tmatrix > matrixPtr
Definition: MatrixBase.h:151
Definition: LanczosLargestMagnitudeEig.h:44
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition: MatrixBase.h:279
MatrixGeneral(const MatrixGeneral< Treal, Tmatrix > &matr)
Copy constructor.
Definition: MatrixGeneral.h:63
static void symm(const char *side, const char *uplo, const int *m, const int *n, const T *alpha, const T *A, const int *lda, const T *B, const int *ldb, const T *beta, T *C, const int *ldc)
Definition: gblas.h:340
static void gemm(const char *ta, const char *tb, const int *n, const int *k, const int *l, const T *alpha, const T *A, const int *lda, const T *B, const int *ldb, const T *beta, T *C, const int *ldc)
Definition: gblas.h:230
MatrixGeneral< Treal, Tmatrix > & operator=(const Xtrans< MatrixGeneral< Treal, Tmatrix > > &mt)
Definition: MatrixGeneral.h:258
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:120
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixGeneral< Treal, Tmatrix > &mat)
Definition: MatrixGeneral.h:253
static Treal frob_diff(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition: MatrixGeneral.h:286
Truncation of general matrices.
Definition: truncation.h:270
This proxy expresses the result of addition of two objects, of possibly different types...
Definition: matrix_proxy.h:244
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixSymmetric< Treal, Tmatrix > &symm)
Definition: MatrixGeneral.h:267
Definition: MatrixBase.h:53
Definition: allocate.cc:30
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
Treal midPoint() const
Definition: Interval.h:113
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:202
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixTriangular< Treal, Tmatrix > &triang)
Definition: MatrixGeneral.h:273
friend class MatrixGeneral< Treal, Tmatrix >
Definition: MatrixBase.h:69
void fullMatrix(std::vector< Treal > &fullMat) const
Definition: MatrixGeneral.h:93
Definition: Interval.h:44
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: MatrixGeneral.h:88
Definition: MatrixBase.h:54
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition: MatrixBase.h:234
MatrixGeneral(const MatrixSymmetric< Treal, Tmatrix > &symm)
Copy from symmetric matrix constructor.
Definition: MatrixGeneral.h:65
int get_nrows() const
Definition: MatrixBase.h:136
size_t nvalues() const
Definition: MatrixGeneral.h:318
This proxy expresses the result of multiplication of three objects, of possibly different types...
Definition: matrix_proxy.h:65
Upper non-unit triangular matrix.
Definition: MatrixBase.h:51
Treal low() const
Definition: Interval.h:142
This proxy expresses the result of transposition of an object of type TX.
Definition: matrix_proxy.h:116
Treal run(Treal const threshold)
Definition: truncation.h:78
int get_ncols() const
Definition: MatrixBase.h:139
void setElementsByRule(TRule &rule)
Definition: MatrixGeneral.h:455
void thresh(Treal const threshold, normType const norm)
Definition: MatrixGeneral.h:507
void rand()
Definition: VectorGeneral.h:94
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:74
void randomZeroStructure(Treal probabilityBeingZero)
Definition: MatrixGeneral.h:450
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:49
Definition: mat_utils.h:69
static void trmm(const char *side, const char *uplo, const char *transa, const char *diag, const int *m, const int *n, const T *alpha, const T *A, const int *lda, T *B, const int *ldb)
Definition: gblas.h:279
Treal eucl_thresh(Treal const threshold)
Definition: MatrixGeneral.h:523
void random()
Definition: MatrixGeneral.h:446
Base class for matrix API.
Definition: MatrixBase.h:67
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:171
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:56
size_t nnz() const
Definition: MatrixGeneral.h:315
std::string obj_type_id() const
Definition: MatrixGeneral.h:460
Base class for matrix API.
MatrixGeneral()
Default constructor.
Definition: MatrixGeneral.h:61
Definition: matInclude.h:135
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixGeneral.h:322
Definition: Failure.h:47
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition: MatrixBase.h:219
This proxy expresses the result of multiplication of two objects, of possibly different types...
Definition: matrix_proxy.h:49
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixGeneral.h:325
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition: MatrixBase.h:252
Treal frob() const
Definition: MatrixGeneral.h:282
#define B
MatrixGeneral(const MatrixTriangular< Treal, Tmatrix > &triang)
Copy from triangular matrix constructor.
Definition: MatrixGeneral.h:69
void frob_thresh(Treal threshold)
Definition: MatrixGeneral.h:297
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:129
Symmetric matrix.
Definition: MatrixBase.h:49
Treal trace(const XYZ< Treal, MatrixGeneral< Treal, Tmatrix >, MatrixGeneral< Treal, Tmatrix > > &smm)
Definition: MatrixGeneral.h:902
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition: matrix_proxy.h:86
MatrixGeneral< Treal, Tmatrix > & operator=(int const k)
Definition: MatrixGeneral.h:278
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: MatrixGeneral.h:462
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: MatrixGeneral.h:465
static Treal trace_aTb(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition: MatrixGeneral.h:311
normType
Definition: matInclude.h:135
static void getPermutedIndexes(std::vector< int > const &index, std::vector< int > const &permutation, std::vector< int > &newIndex)
Definition: MatrixBase.h:203