ergo
MatrixGeneral.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
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,
335 
340 
342  MatrixGeneral<Treal, Tmatrix>& operator+=
343  (const XYZ<Treal,
346 
349  (const XYZpUV<Treal,
352  Treal,
354 
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 */
372  (const XYZ<Treal,
380  MatrixGeneral<Treal, Tmatrix>& operator+=
381  (const XYZ<Treal,
386  (const XYZpUV<Treal,
389  Treal,
393  (const XYZ<Treal,
401  MatrixGeneral<Treal, Tmatrix>& operator+=
402  (const XYZ<Treal,
407  (const XYZpUV<Treal,
410  Treal,
414  (const XYZ<Treal,
422  MatrixGeneral<Treal, Tmatrix>& operator+=
423  (const XYZ<Treal,
428  (const XYZpUV<Treal,
431  Treal,
433 
434  /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */
437  (const XYZ<Treal,
442  (const XYZ<Treal,
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;
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,
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>
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>
562  (const XYZ<Treal,
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>
576  (const XYZpUV<Treal,
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>
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>
612  Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
613  return *this;
614  }
615  /* B -= A */
616  template<typename Treal, typename Tmatrix>
620  Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
621  return *this;
622  }
623  /* B += alpha * A */
624  template<typename Treal, typename Tmatrix>
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>
639  (const XYZ<Treal,
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>
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>
669  (const XYZ<Treal,
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>
684  (const XYZpUV<Treal,
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>
710  (const XYZ<Treal,
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>
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>
740  (const XYZ<Treal,
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>
755  (const XYZpUV<Treal,
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>
782  (const XYZ<Treal,
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>
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>
814  (const XYZ<Treal,
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>
831  (const XYZpUV<Treal,
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>
862  (const XYZ<Treal,
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>
883  (const XYZ<Treal,
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,
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