ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * Kohn−Sham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00039 #ifndef MAT_MATRIX_PROXY 00040 #define MAT_MATRIX_PROXY 00041 00042 namespace mat { 00043 /*********** New code */ 00048 template<typename TX, typename TY> 00049 struct XY { 00050 TX const & A; 00051 TY const & B; 00052 bool const tA; 00053 bool const tB; 00054 XY(TX const & AA, TY const & BB, 00055 bool const tAA = false, bool const tBB = false) 00056 :A(AA), B(BB), tA(tAA), tB(tBB) 00057 {} 00058 }; 00059 00064 template<typename TX, typename TY, typename TZ> 00065 struct XYZ { 00066 TX const & A; 00067 TY const & B; 00068 TZ const & C; 00069 bool const tA; 00070 bool const tB; 00071 bool const tC; 00072 XYZ(TX const & AA, TY const & BB, TZ const & CC, 00073 bool const tAA = false, 00074 bool const tBB = false, 00075 bool const tCC = false) 00076 :A(AA), B(BB), C(CC), tA(tAA), tB(tBB), tC(tCC) 00077 {} 00078 }; 00079 00085 template<typename TX, typename TY, typename TZ, typename TU, typename TV> 00086 struct XYZpUV { 00087 TX const & A; 00088 TY const & B; 00089 TZ const & C; 00090 TU const & D; 00091 TV const & E; 00092 bool const tA; 00093 bool const tB; 00094 bool const tC; 00095 bool const tD; 00096 bool const tE; 00097 XYZpUV(TX const & AA, TY const & BB, TZ const & CC, 00098 TU const & DD, TV const & EE, 00099 bool const tAA = false, 00100 bool const tBB = false, 00101 bool const tCC = false, 00102 bool const tDD = false, 00103 bool const tEE = false) 00104 :A(AA), B(BB), C(CC), D(DD), E(EE), 00105 tA(tAA), tB(tBB), tC(tCC), tD(tDD), tE(tEE) 00106 {} 00107 }; 00108 00109 00115 template<typename TX> 00116 struct Xtrans { 00117 TX const & A; 00118 bool const tA; 00119 explicit Xtrans(TX const & AA, bool const tAA = false) 00120 :A(AA), tA(tAA) 00121 {} 00122 }; 00123 00128 template<typename TX> 00129 inline Xtrans<TX> transpose(TX const & A) { 00130 return Xtrans<TX>(A,true); 00131 } 00139 template<typename TX> 00140 inline Xtrans<TX> transpose(const Xtrans<TX>& xtrans) { 00141 return Xtrans<TX>(xtrans.A, !xtrans.tA); 00142 } 00143 00144 /* Some operators */ 00154 template<typename TX, typename TY> 00155 inline XY<TX, TY> operator*(Xtrans<TX> const & trAA, 00156 Xtrans<TY> const & trBB) { 00157 return XY<TX, TY>(trAA.A, trBB.A, trAA.tA, trBB.tA); 00158 } 00159 00169 template<typename TX, typename TY> 00170 inline XY<TX, TY> operator*(TX const & AA, 00171 Xtrans<TY> const & trBB) { 00172 return XY<TX, TY>(AA, trBB.A, false, trBB.tA); 00173 } 00174 00184 template<typename TX, typename TY> 00185 inline XY<TX, TY> operator*(Xtrans<TX> const & trAA, 00186 TY const & BB) { 00187 return XY<TX, TY>(trAA.A, BB, trAA.tA, false); 00188 } 00189 00197 template<typename TX, typename TY> 00198 inline XY<TX, TY> operator*(TX const & AA, 00199 TY const & BB) { 00200 return XY<TX, TY>(AA, BB, false, false); 00201 } 00202 00210 template<typename TX, typename TY, typename TZ> 00211 inline XYZ<TX, TY, TZ> 00212 operator*(XY<TX, TY> const & AB, Xtrans<TZ> const & trCC) { 00213 return XYZ<TX, TY, TZ>(AB.A, AB.B, trCC.A, AB.tA, AB.tB, trCC.tA); 00214 } 00215 00221 template<typename TX, typename TY, typename TZ> 00222 inline XYZ<TX, TY, TZ> 00223 operator*(XY<TX, TY> const & AB, TZ const & CC) { 00224 return XYZ<TX, TY, TZ>(AB.A, AB.B, CC, AB.tA, AB.tB, false); 00225 } 00226 00233 template<typename TX, typename TY, typename TZ, typename TU, typename TV> 00234 inline XYZpUV<TX, TY, TZ, TU, TV> 00235 operator+(XYZ<TX, TY, TZ> const & ABC, XY<TU, TV> const & DE) { 00236 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); 00237 } 00238 00243 template<typename TX, typename TY> 00244 struct XpY { 00245 const TX& A; 00246 const TY& B; 00247 XpY(const TX& AA,const TY& BB) 00248 :A(AA),B(BB) 00249 {} 00250 }; 00254 template<typename TX, typename TY> 00255 inline XpY<TX, TY> operator+(TX const & AA, TY const & BB) { 00256 return XpY<TX, TY>(AA, BB); 00257 } 00258 00263 template<typename TX, typename TY> 00264 struct XmY { 00265 const TX& A; 00266 const TY& B; 00267 XmY(const TX& AA,const TY& BB) 00268 :A(AA),B(BB) 00269 {} 00270 }; 00274 template<typename TX, typename TY> 00275 inline XmY<TX, TY> operator-(TX const & AA, TY const & BB) { 00276 return XmY<TX, TY>(AA, BB); 00277 } 00278 00279 00280 /************* New code ends */ 00281 00282 00283 #if 0 00284 template<class MAT> 00285 struct M2 { 00286 const MAT& A; 00287 M2(const MAT& AA) 00288 :A(AA) 00289 {} 00290 }; 00291 00292 template<class MAT> 00293 inline M2<MAT> square(const MAT& A) { 00294 return M2<MAT>(A); 00295 } 00296 00297 template<class SCAL, class MAT> 00298 struct SM2 { 00299 const SCAL alpha; 00300 const MAT& A; 00301 SM2(const MAT& AA, const SCAL a = 1) 00302 : A(AA), alpha(a) 00303 {} 00304 SM2(const M2<MAT>& m2) 00305 :A(m2.A), alpha(1) 00306 {} 00307 }; 00308 00309 template<class SCAL, class MAT> 00310 inline SM2<SCAL, MAT> 00311 operator*(const SCAL s, const M2<MAT>& m2) { 00312 return SM2<SCAL, MAT>(m2.A, s); 00313 } 00314 00315 00316 00317 00318 template<class MAT> 00319 struct MT { 00320 const MAT& A; 00321 const bool tA; 00322 MT(const MAT& AA, const bool tAA = false) 00323 :A(AA), tA(tAA) 00324 {} 00325 }; 00326 00327 template<class MAT> 00328 inline MT<MAT> transpose(const MAT& A) { 00329 return MT<MAT>(A,true); 00330 } 00331 template<class MAT> 00332 inline MT<MAT> transpose(const MT<MAT>& mt) { 00333 return MT<MAT>(mt.A, !mt.tA); 00334 } 00335 00336 00337 template<class SCAL, class MAT> 00338 struct SM { 00339 const SCAL alpha; 00340 const MAT& A; 00341 const bool tA; 00342 SM(const MAT& AA, const SCAL scalar = 1, const bool tAA = false) 00343 :A(AA),alpha(scalar), tA(tAA) 00344 {} 00345 }; 00346 00347 template<class SCAL, class MAT> 00348 inline SM<SCAL, MAT> operator*(const SCAL scalar, const MT<MAT>& mta) { 00349 return SM<SCAL, MAT>(mta.A,scalar, mta.tA); 00350 } 00351 00352 template<class SCAL, class MAT> 00353 inline SM<SCAL, MAT> operator*(const SCAL scalar, const MAT& AA) { 00354 return SM<SCAL, MAT>(AA, scalar, false); 00355 } 00356 00357 00358 00359 template<class MAT, class MATB = MAT> 00360 struct MM { 00361 const MAT& A; 00362 const MATB& B; 00363 const bool tA; 00364 const bool tB; 00365 MM(const MAT& AA,const MATB& BB, const bool tAA, const bool tBB) 00366 :A(AA),B(BB), tA(tAA), tB(tBB) 00367 {} 00368 }; 00369 00370 template<class MAT, class MATB = MAT> 00371 struct MpM { 00372 const MAT& A; 00373 const MATB& B; 00374 MpM(const MAT& AA,const MATB& BB) 00375 :A(AA),B(BB) 00376 {} 00377 }; 00378 00379 template<class MAT, class MATB> 00380 inline MpM<MAT, MATB> operator+(const MAT& AA, const MATB& BB) { 00381 return MpM<MAT, MATB>(AA, BB); 00382 } 00383 00384 /* 00385 template<class MAT, class MATB> 00386 inline MM<MAT, MATB> operator*(const MT<MAT>& mta, const MT<MATB>& mtb) { 00387 return MM<MAT, MATB>(mta.A, mtb.A, mta.tA, mtb.tA); 00388 } 00389 */ 00390 /* 00391 template<class MAT, class MATB> 00392 inline MM<MAT, MATB> operator*(const MAT& AA, const MT<MATB>& mtb) { 00393 return MM<MAT, MATB>(AA, mtb.A, false, mtb.tA); 00394 } 00395 template<class MAT, class MATB> 00396 inline MM<MAT, MATB> operator*(const MT<MAT>& mta, const MATB& BB) { 00397 return MM<MAT, MATB>(mta.A, BB, mta.tA, false); 00398 } 00399 template<class MAT, class MATB> 00400 inline MM<MAT, MATB> operator*(const MAT& AA, const MATB& BB) { 00401 return MM<MAT, MATB>(AA, BB, false, false); 00402 } 00403 */ 00404 00405 template<class SCAL, class MAT, class MATB = MAT> 00406 struct SMM { 00407 const SCAL alpha; 00408 const MAT& A; 00409 const MATB& B; 00410 const bool tA; 00411 const bool tB; 00412 SMM(const MM<MAT, MATB>& mm) 00413 :A(mm.A),B(mm.B),alpha(1), tA(mm.tA), tB(mm.tB) 00414 {} 00415 SMM(const MAT& AA,const MATB& BB, 00416 const bool tAA, const bool tBB, 00417 const SCAL a = 1) 00418 :A(AA), B(BB), tA(tAA), tB(tBB), alpha(a) 00419 {} 00420 }; 00421 00422 template<class SCAL, class MAT, class MATB> 00423 inline SMM<SCAL, MAT, MATB> 00424 operator*(const SM<SCAL, MAT>& sm,const MT<MATB>& mtb) { 00425 return SMM<SCAL, MAT, MATB>(sm.A, mtb.A, sm.tA, mtb.tA, sm.alpha); 00426 } 00427 00428 template<class SCAL, class MAT, class MATB> 00429 inline SMM<SCAL, MAT, MATB> 00430 operator*(const SM<SCAL, MAT>& sm,const MATB& BB) { 00431 return SMM<SCAL, MAT, MATB>(sm.A, BB, sm.tA, false, sm.alpha); 00432 } 00433 00434 template<class SCAL, class MATC, class MATA = MATC, class MATB = MATC> 00435 struct SMMpSM { 00436 const SCAL alpha; 00437 const MATA& A; 00438 const MATB& B; 00439 const SCAL beta; 00440 const MATC& C; 00441 const bool tA; 00442 const bool tB; 00443 SMMpSM(const MATA& AA, const MATB& BB, const MATC& CC, 00444 const bool tAA, const bool tBB, 00445 const SCAL a=1, const SCAL b=1) 00446 :A(AA), B(BB), C(CC), alpha(a), beta(b), tA(tAA), tB(tBB) 00447 {} 00448 }; 00449 00450 template<class SCAL, class MATC, class MATA, class MATB> 00451 inline SMMpSM<SCAL, MATC, MATA, MATB> 00452 operator+(const SMM<SCAL, MATA, MATB>& smm, const SM<SCAL, MATC>& sm) { 00453 return SMMpSM<SCAL, MATC, MATA, MATB> 00454 (smm.A, smm.B, sm.A, smm.tA, smm.tB, smm.alpha, sm.alpha); 00455 } 00456 #if 0 00457 00458 template<class SCAL, class MATC, class MATA, class MATB> 00459 inline SMMpSM<SCAL, MATC, MATA, MATB> 00460 operator+(const SMM<SCAL, MATA, MATB>& smm, MATC& CC) { 00461 return SMMpSM<SCAL, MATC, MATA, MATB> 00462 (smm.A, smm.B, CC, smm.tA, smm.tB, smm.alpha, 1); 00463 } 00464 template<class SCAL, class MATC, class MATA, class MATB> 00465 inline SMMpSM<SCAL, MATC, MATA, MATB> 00466 operator+(const MM<MATA, MATB>& mm, const SM<SCAL, MATC>& sm) { 00467 return SMMpSM<SCAL, MATC, MATA, MATB> 00468 (mm.A, mm.B, sm.A, mm.tA, mm.tB, 1, sm.alpha); 00469 } 00470 #endif 00471 00472 00473 template<class SCAL, class MAT> 00474 struct SM2pSM { 00475 const SCAL alpha; 00476 const MAT& A; 00477 const SCAL beta; 00478 const MAT& C; 00479 SM2pSM(const MAT& AA, const MAT& CC, const SCAL a = 1, const SCAL b = 0) 00480 : A(AA), alpha(a), C(CC), beta(b) 00481 {} 00482 }; 00483 00484 template<class SCAL, class MAT> 00485 inline SM2pSM<SCAL, MAT> 00486 operator+(const SM2<SCAL, MAT>& sm2, const SM<SCAL, MAT> sm) { 00487 return SM2pSM<SCAL, MAT>(sm2.A, sm.A, sm2.alpha, sm.alpha); 00488 } 00489 00490 00491 00492 /* Done so far with new transpose */ 00493 00494 template<class MAT> 00495 struct MMpM { 00496 const MAT& A; 00497 const MAT& B; 00498 const MAT& C; 00499 MMpM(const MAT& AA, const MAT& BB, const MAT& CC) 00500 :A(AA),B(BB),C(CC) 00501 {} 00502 }; 00503 00504 00505 template<class SCAL, class MAT> 00506 struct SMpSM { 00507 const SCAL alpha, beta; 00508 const MAT& A, B; 00509 SMpSM(const MAT& AA, const MAT& BB, 00510 const SCAL scalar_a=1, const SCAL scalar_b=1) 00511 :A(AA), B(BB), alpha(scalar_a), beta(scalar_b) 00512 {} 00513 }; 00514 00515 template<class SCAL, class MAT> 00516 inline SMpSM<SCAL, MAT> 00517 operator+(const SM<SCAL, MAT> sm1, const SM<SCAL, MAT> sm2 ) { 00518 return SMpSM<SCAL, MAT>(sm1.A, sm2.A, sm1.alpha, sm2.alpha); 00519 } 00520 00521 /* 00522 template<class MAT> 00523 struct MpM { 00524 const MAT& A; 00525 const MAT& B; 00526 MpM(const MAT& AA,const MAT& BB) 00527 :A(AA),B(BB) 00528 {} 00529 }; 00530 00531 template<class MAT> 00532 inline MpM<MAT> operator+(const MAT& A, const MAT& B) { 00533 return MpM<MAT>(A,B); 00534 } 00535 */ 00536 template<class MAT> 00537 struct MmM { 00538 const MAT& A; 00539 const MAT& B; 00540 MmM(const MAT& AA,const MAT& BB) 00541 :A(AA),B(BB) 00542 {} 00543 }; 00544 00545 template<class MAT> 00546 inline MmM<MAT> operator-(const MAT& A, const MAT& B) { 00547 return MmM<MAT>(A,B); 00548 } 00549 00550 00551 00552 00553 00554 00555 /*onodig finns redan för SMM 00556 template<class MAT> 00557 inline MMpM<MAT> operator+(const MM<MAT>& mm, const MAT& CC) 00558 { 00559 return MMpM<MAT>(mm.A,mm.B,CC); 00560 }*/ 00561 00562 00563 /*Maste ligga i arvda klassen!!*/ 00564 /* 00565 Matrix::Matrix(const sMMmul& mm) 00566 :nrofrows(mm.A.nrofrows),nrofcols(mm.B.nrofcols) 00567 { 00568 this.multiply(mm.A,mm.B,*this,mm.tA,mm.tB,mm.alpha,0); 00569 } 00570 Matrix::Matrix(const sMMmulsMadd& mm) 00571 :nrofrows(mm.A.nrofrows),nrofcols(mm.B.nrofcols) 00572 { 00573 this->multiply(mm.A,mm.B,mm.C,mm.tA,mm.tB,mm.alpha,mm.beta); 00574 } 00575 00576 */ 00577 #endif 00578 } /* end namespace mat */ 00579 #endif