ergo
matrix_proxy.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 
39 #ifndef MAT_MATRIX_PROXY
40 #define MAT_MATRIX_PROXY
41 
42 namespace mat {
43  /*********** New code */
48  template<typename TX, typename TY>
49  struct XY {
50  TX const & A;
51  TY const & B;
52  bool const tA;
53  bool const tB;
54  XY(TX const & AA, TY const & BB,
55  bool const tAA = false, bool const tBB = false)
56  :A(AA), B(BB), tA(tAA), tB(tBB)
57  {}
58  };
59 
64  template<typename TX, typename TY, typename TZ>
65  struct XYZ {
66  TX const & A;
67  TY const & B;
68  TZ const & C;
69  bool const tA;
70  bool const tB;
71  bool const tC;
72  XYZ(TX const & AA, TY const & BB, TZ const & CC,
73  bool const tAA = false,
74  bool const tBB = false,
75  bool const tCC = false)
76  :A(AA), B(BB), C(CC), tA(tAA), tB(tBB), tC(tCC)
77  {}
78  };
79 
85  template<typename TX, typename TY, typename TZ, typename TU, typename TV>
86  struct XYZpUV {
87  TX const & A;
88  TY const & B;
89  TZ const & C;
90  TU const & D;
91  TV const & E;
92  bool const tA;
93  bool const tB;
94  bool const tC;
95  bool const tD;
96  bool const tE;
97  XYZpUV(TX const & AA, TY const & BB, TZ const & CC,
98  TU const & DD, TV const & EE,
99  bool const tAA = false,
100  bool const tBB = false,
101  bool const tCC = false,
102  bool const tDD = false,
103  bool const tEE = false)
104  :A(AA), B(BB), C(CC), D(DD), E(EE),
105  tA(tAA), tB(tBB), tC(tCC), tD(tDD), tE(tEE)
106  {}
107  };
108 
109 
115  template<typename TX>
116  struct Xtrans {
117  TX const & A;
118  bool const tA;
119  explicit Xtrans(TX const & AA, bool const tAA = false)
120  :A(AA), tA(tAA)
121  {}
122  };
123 
128  template<typename TX>
129  inline Xtrans<TX> transpose(TX const & A) {
130  return Xtrans<TX>(A,true);
131  }
139  template<typename TX>
140  inline Xtrans<TX> transpose(const Xtrans<TX>& xtrans) {
141  return Xtrans<TX>(xtrans.A, !xtrans.tA);
142  }
143 
144  /* Some operators */
154  template<typename TX, typename TY>
155  inline XY<TX, TY> operator*(Xtrans<TX> const & trAA,
156  Xtrans<TY> const & trBB) {
157  return XY<TX, TY>(trAA.A, trBB.A, trAA.tA, trBB.tA);
158  }
159 
169  template<typename TX, typename TY>
170  inline XY<TX, TY> operator*(TX const & AA,
171  Xtrans<TY> const & trBB) {
172  return XY<TX, TY>(AA, trBB.A, false, trBB.tA);
173  }
174 
184  template<typename TX, typename TY>
185  inline XY<TX, TY> operator*(Xtrans<TX> const & trAA,
186  TY const & BB) {
187  return XY<TX, TY>(trAA.A, BB, trAA.tA, false);
188  }
189 
197  template<typename TX, typename TY>
198  inline XY<TX, TY> operator*(TX const & AA,
199  TY const & BB) {
200  return XY<TX, TY>(AA, BB, false, false);
201  }
202 
210  template<typename TX, typename TY, typename TZ>
211  inline XYZ<TX, TY, TZ>
212  operator*(XY<TX, TY> const & AB, Xtrans<TZ> const & trCC) {
213  return XYZ<TX, TY, TZ>(AB.A, AB.B, trCC.A, AB.tA, AB.tB, trCC.tA);
214  }
215 
221  template<typename TX, typename TY, typename TZ>
222  inline XYZ<TX, TY, TZ>
223  operator*(XY<TX, TY> const & AB, TZ const & CC) {
224  return XYZ<TX, TY, TZ>(AB.A, AB.B, CC, AB.tA, AB.tB, false);
225  }
226 
233  template<typename TX, typename TY, typename TZ, typename TU, typename TV>
234  inline XYZpUV<TX, TY, TZ, TU, TV>
235  operator+(XYZ<TX, TY, TZ> const & ABC, XY<TU, TV> const & DE) {
236  return XYZpUV<TX, TY, TZ, TU, TV>(ABC.A, ABC.B, ABC.C, DE.A, DE.B, ABC.tA, ABC.tB, ABC.tC, DE.tA, DE.tB);
237  }
238 
243  template<typename TX, typename TY>
244  struct XpY {
245  const TX& A;
246  const TY& B;
247  XpY(const TX& AA,const TY& BB)
248  :A(AA),B(BB)
249  {}
250  };
254  template<typename TX, typename TY>
255  inline XpY<TX, TY> operator+(TX const & AA, TY const & BB) {
256  return XpY<TX, TY>(AA, BB);
257  }
258 
263  template<typename TX, typename TY>
264  struct XmY {
265  const TX& A;
266  const TY& B;
267  XmY(const TX& AA,const TY& BB)
268  :A(AA),B(BB)
269  {}
270  };
274  template<typename TX, typename TY>
275  inline XmY<TX, TY> operator-(TX const & AA, TY const & BB) {
276  return XmY<TX, TY>(AA, BB);
277  }
278 
279 
280  /************* New code ends */
281 
282 
283 #if 0
284  template<class MAT>
285  struct M2 {
286  const MAT& A;
287  M2(const MAT& AA)
288  :A(AA)
289  {}
290  };
291 
292  template<class MAT>
293  inline M2<MAT> square(const MAT& A) {
294  return M2<MAT>(A);
295  }
296 
297  template<class SCAL, class MAT>
298  struct SM2 {
299  const SCAL alpha;
300  const MAT& A;
301  SM2(const MAT& AA, const SCAL a = 1)
302  : A(AA), alpha(a)
303  {}
304  SM2(const M2<MAT>& m2)
305  :A(m2.A), alpha(1)
306  {}
307  };
308 
309  template<class SCAL, class MAT>
310  inline SM2<SCAL, MAT>
311  operator*(const SCAL s, const M2<MAT>& m2) {
312  return SM2<SCAL, MAT>(m2.A, s);
313  }
314 
315 
316 
317 
318  template<class MAT>
319  struct MT {
320  const MAT& A;
321  const bool tA;
322  MT(const MAT& AA, const bool tAA = false)
323  :A(AA), tA(tAA)
324  {}
325  };
326 
327  template<class MAT>
328  inline MT<MAT> transpose(const MAT& A) {
329  return MT<MAT>(A,true);
330  }
331  template<class MAT>
332  inline MT<MAT> transpose(const MT<MAT>& mt) {
333  return MT<MAT>(mt.A, !mt.tA);
334  }
335 
336 
337  template<class SCAL, class MAT>
338  struct SM {
339  const SCAL alpha;
340  const MAT& A;
341  const bool tA;
342  SM(const MAT& AA, const SCAL scalar = 1, const bool tAA = false)
343  :A(AA),alpha(scalar), tA(tAA)
344  {}
345  };
346 
347  template<class SCAL, class MAT>
348  inline SM<SCAL, MAT> operator*(const SCAL scalar, const MT<MAT>& mta) {
349  return SM<SCAL, MAT>(mta.A,scalar, mta.tA);
350  }
351 
352  template<class SCAL, class MAT>
353  inline SM<SCAL, MAT> operator*(const SCAL scalar, const MAT& AA) {
354  return SM<SCAL, MAT>(AA, scalar, false);
355  }
356 
357 
358 
359  template<class MAT, class MATB = MAT>
360  struct MM {
361  const MAT& A;
362  const MATB& B;
363  const bool tA;
364  const bool tB;
365  MM(const MAT& AA,const MATB& BB, const bool tAA, const bool tBB)
366  :A(AA),B(BB), tA(tAA), tB(tBB)
367  {}
368  };
369 
370  template<class MAT, class MATB = MAT>
371  struct MpM {
372  const MAT& A;
373  const MATB& B;
374  MpM(const MAT& AA,const MATB& BB)
375  :A(AA),B(BB)
376  {}
377  };
378 
379  template<class MAT, class MATB>
380  inline MpM<MAT, MATB> operator+(const MAT& AA, const MATB& BB) {
381  return MpM<MAT, MATB>(AA, BB);
382  }
383 
384  /*
385  template<class MAT, class MATB>
386  inline MM<MAT, MATB> operator*(const MT<MAT>& mta, const MT<MATB>& mtb) {
387  return MM<MAT, MATB>(mta.A, mtb.A, mta.tA, mtb.tA);
388  }
389  */
390  /*
391  template<class MAT, class MATB>
392  inline MM<MAT, MATB> operator*(const MAT& AA, const MT<MATB>& mtb) {
393  return MM<MAT, MATB>(AA, mtb.A, false, mtb.tA);
394  }
395  template<class MAT, class MATB>
396  inline MM<MAT, MATB> operator*(const MT<MAT>& mta, const MATB& BB) {
397  return MM<MAT, MATB>(mta.A, BB, mta.tA, false);
398  }
399  template<class MAT, class MATB>
400  inline MM<MAT, MATB> operator*(const MAT& AA, const MATB& BB) {
401  return MM<MAT, MATB>(AA, BB, false, false);
402  }
403  */
404 
405  template<class SCAL, class MAT, class MATB = MAT>
406  struct SMM {
407  const SCAL alpha;
408  const MAT& A;
409  const MATB& B;
410  const bool tA;
411  const bool tB;
412  SMM(const MM<MAT, MATB>& mm)
413  :A(mm.A),B(mm.B),alpha(1), tA(mm.tA), tB(mm.tB)
414  {}
415  SMM(const MAT& AA,const MATB& BB,
416  const bool tAA, const bool tBB,
417  const SCAL a = 1)
418  :A(AA), B(BB), tA(tAA), tB(tBB), alpha(a)
419  {}
420  };
421 
422  template<class SCAL, class MAT, class MATB>
423  inline SMM<SCAL, MAT, MATB>
424  operator*(const SM<SCAL, MAT>& sm,const MT<MATB>& mtb) {
425  return SMM<SCAL, MAT, MATB>(sm.A, mtb.A, sm.tA, mtb.tA, sm.alpha);
426  }
427 
428  template<class SCAL, class MAT, class MATB>
429  inline SMM<SCAL, MAT, MATB>
430  operator*(const SM<SCAL, MAT>& sm,const MATB& BB) {
431  return SMM<SCAL, MAT, MATB>(sm.A, BB, sm.tA, false, sm.alpha);
432  }
433 
434  template<class SCAL, class MATC, class MATA = MATC, class MATB = MATC>
435  struct SMMpSM {
436  const SCAL alpha;
437  const MATA& A;
438  const MATB& B;
439  const SCAL beta;
440  const MATC& C;
441  const bool tA;
442  const bool tB;
443  SMMpSM(const MATA& AA, const MATB& BB, const MATC& CC,
444  const bool tAA, const bool tBB,
445  const SCAL a=1, const SCAL b=1)
446  :A(AA), B(BB), C(CC), alpha(a), beta(b), tA(tAA), tB(tBB)
447  {}
448  };
449 
450  template<class SCAL, class MATC, class MATA, class MATB>
451  inline SMMpSM<SCAL, MATC, MATA, MATB>
452  operator+(const SMM<SCAL, MATA, MATB>& smm, const SM<SCAL, MATC>& sm) {
453  return SMMpSM<SCAL, MATC, MATA, MATB>
454  (smm.A, smm.B, sm.A, smm.tA, smm.tB, smm.alpha, sm.alpha);
455  }
456 #if 0
457 
458  template<class SCAL, class MATC, class MATA, class MATB>
459  inline SMMpSM<SCAL, MATC, MATA, MATB>
460  operator+(const SMM<SCAL, MATA, MATB>& smm, MATC& CC) {
461  return SMMpSM<SCAL, MATC, MATA, MATB>
462  (smm.A, smm.B, CC, smm.tA, smm.tB, smm.alpha, 1);
463  }
464  template<class SCAL, class MATC, class MATA, class MATB>
465  inline SMMpSM<SCAL, MATC, MATA, MATB>
466  operator+(const MM<MATA, MATB>& mm, const SM<SCAL, MATC>& sm) {
467  return SMMpSM<SCAL, MATC, MATA, MATB>
468  (mm.A, mm.B, sm.A, mm.tA, mm.tB, 1, sm.alpha);
469  }
470 #endif
471 
472 
473  template<class SCAL, class MAT>
474  struct SM2pSM {
475  const SCAL alpha;
476  const MAT& A;
477  const SCAL beta;
478  const MAT& C;
479  SM2pSM(const MAT& AA, const MAT& CC, const SCAL a = 1, const SCAL b = 0)
480  : A(AA), alpha(a), C(CC), beta(b)
481  {}
482  };
483 
484  template<class SCAL, class MAT>
485  inline SM2pSM<SCAL, MAT>
486  operator+(const SM2<SCAL, MAT>& sm2, const SM<SCAL, MAT> sm) {
487  return SM2pSM<SCAL, MAT>(sm2.A, sm.A, sm2.alpha, sm.alpha);
488  }
489 
490 
491 
492  /* Done so far with new transpose */
493 
494 template<class MAT>
495  struct MMpM {
496  const MAT& A;
497  const MAT& B;
498  const MAT& C;
499  MMpM(const MAT& AA, const MAT& BB, const MAT& CC)
500  :A(AA),B(BB),C(CC)
501  {}
502  };
503 
504 
505  template<class SCAL, class MAT>
506  struct SMpSM {
507  const SCAL alpha, beta;
508  const MAT& A, B;
509  SMpSM(const MAT& AA, const MAT& BB,
510  const SCAL scalar_a=1, const SCAL scalar_b=1)
511  :A(AA), B(BB), alpha(scalar_a), beta(scalar_b)
512  {}
513  };
514 
515  template<class SCAL, class MAT>
516  inline SMpSM<SCAL, MAT>
517  operator+(const SM<SCAL, MAT> sm1, const SM<SCAL, MAT> sm2 ) {
518  return SMpSM<SCAL, MAT>(sm1.A, sm2.A, sm1.alpha, sm2.alpha);
519  }
520 
521  /*
522  template<class MAT>
523  struct MpM {
524  const MAT& A;
525  const MAT& B;
526  MpM(const MAT& AA,const MAT& BB)
527  :A(AA),B(BB)
528  {}
529  };
530 
531  template<class MAT>
532  inline MpM<MAT> operator+(const MAT& A, const MAT& B) {
533  return MpM<MAT>(A,B);
534  }
535  */
536  template<class MAT>
537  struct MmM {
538  const MAT& A;
539  const MAT& B;
540  MmM(const MAT& AA,const MAT& BB)
541  :A(AA),B(BB)
542  {}
543  };
544 
545  template<class MAT>
546  inline MmM<MAT> operator-(const MAT& A, const MAT& B) {
547  return MmM<MAT>(A,B);
548  }
549 
550 
551 
552 
553 
554 
555  /*onodig finns redan för SMM
556  template<class MAT>
557  inline MMpM<MAT> operator+(const MM<MAT>& mm, const MAT& CC)
558  {
559  return MMpM<MAT>(mm.A,mm.B,CC);
560  }*/
561 
562 
563  /*Maste ligga i arvda klassen!!*/
564  /*
565  Matrix::Matrix(const sMMmul& mm)
566  :nrofrows(mm.A.nrofrows),nrofcols(mm.B.nrofcols)
567  {
568  this.multiply(mm.A,mm.B,*this,mm.tA,mm.tB,mm.alpha,0);
569  }
570  Matrix::Matrix(const sMMmulsMadd& mm)
571  :nrofrows(mm.A.nrofrows),nrofcols(mm.B.nrofcols)
572  {
573  this->multiply(mm.A,mm.B,mm.C,mm.tA,mm.tB,mm.alpha,mm.beta);
574  }
575 
576  */
577 #endif
578 } /* end namespace mat */
579 #endif