00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 #include "Projectors.h"
00051
00052
00053 const int Projectors::OP_UNSYM1 = 1;
00054 const int Projectors::OP_UNSYM2 = 2;
00055 const int Projectors::OP_SYM = 3;
00056
00057
00058 Projectors::Projectors(const Epetra_Comm &_Comm,
00059 const Epetra_Operator *AA, const Epetra_Operator *BB, const Epetra_Operator *PP,
00060 int eqType, int k, int kmax, int *Hpiv, int *Gpiv,
00061 double theta, double *Hlu, double *Glu,
00062 Epetra_MultiVector *Q, Epetra_MultiVector *Qb, Epetra_MultiVector *Y,
00063 double *w1, double *w2)
00064 : callBLAS(),
00065 callLAPACK(),
00066 orthoTool(),
00067 MyComm(_Comm),
00068 MyWatch(_Comm),
00069 A(AA),
00070 B(BB),
00071 APrec(PP),
00072 Pro_k(k),
00073 Pro_kmax(kmax),
00074 Pro_EquationType(eqType),
00075 Pro_Hpiv(Hpiv),
00076 Pro_Gpiv(Gpiv),
00077 Pro_theta(theta),
00078 Pro_Hlu(Hlu),
00079 Pro_Glu(Glu),
00080 Pro_Q(Q),
00081 Pro_Qb(Qb),
00082 Pro_Y(Y),
00083 work1(w1),
00084 work2(w2),
00085 timeApply(0.0),
00086 timeApplyInverse(0.0)
00087 {
00088
00089 }
00090
00091
00092 const Epetra_Map& Projectors::OperatorDomainMap() const {
00093
00094
00095
00096 if (A)
00097 return A->OperatorDomainMap();
00098 if (B)
00099 return B->OperatorDomainMap();
00100
00101 return APrec->OperatorDomainMap();
00102
00103 }
00104
00105
00106 const Epetra_Map& Projectors::OperatorRangeMap() const {
00107
00108
00109
00110 if (A)
00111 return A->OperatorRangeMap();
00112 if (B)
00113 return B->OperatorRangeMap();
00114
00115 return APrec->OperatorRangeMap();
00116
00117 }
00118
00119
00120
00121 int Projectors::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00122
00123 timeApply -= MyWatch.WallTime();
00124 if (A) {
00125 switch(Pro_EquationType){
00126 case OP_UNSYM1:
00127 LinOp_UNSYM1(X, Y);
00128 break;
00129 case OP_UNSYM2:
00130 LinOp_UNSYM2(X, Y);
00131 break;
00132 case OP_SYM:
00133 LinOp_SYM(X, Y);
00134 break;
00135 }
00136 }
00137 else {
00138 memcpy(Y.Values(), X.Values(), X.NumVectors()*X.MyLength()*sizeof(double));
00139 }
00140 timeApply += MyWatch.WallTime();
00141
00142 return 0;
00143
00144 }
00145
00146
00147
00148 int Projectors::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00149
00150 int nV = X.NumVectors();
00151
00152 timeApplyInverse -= MyWatch.WallTime();
00153 if (APrec) {
00154 switch(Pro_EquationType){
00155 case OP_UNSYM1:
00156 Precon_UNSYM1(X, Y);
00157 break;
00158 case OP_UNSYM2:
00159 Precon_UNSYM2(X, Y);
00160 break;
00161 case OP_SYM:
00162 Precon_SYM(X, Y);
00163 break;
00164 }
00165 }
00166 else {
00167 memcpy(Y.Values(), X.Values(), nV*X.MyLength()*sizeof(double));
00168 }
00169 timeApplyInverse += MyWatch.WallTime();
00170
00171 return 0;
00172
00173 }
00174
00175
00176
00177 void Projectors::MakeRHS(Epetra_MultiVector &R){
00178
00179 int iV;
00180 int nV = R.NumVectors();
00181
00182 switch(Pro_EquationType){
00183 case OP_UNSYM1:
00184 for (iV = 0; iV < nV; ++iV) {
00185 Epetra_Vector Rv(View, R, iV);
00186 Right_UNSYM1(Rv);
00187 }
00188 break;
00189 case OP_UNSYM2:
00190 for (iV = 0; iV < nV; ++iV) {
00191 Epetra_Vector Rv(View, R, iV);
00192 Right_UNSYM2(Rv);
00193 }
00194 break;
00195 case OP_SYM:
00196 for (iV = 0; iV < nV; ++iV) {
00197 Epetra_Vector Rv(View, R, iV);
00198 Right_SYM(Rv);
00199 }
00200 break;
00201 }
00202
00203 }
00204
00205
00206
00207
00208
00209
00210
00211 void Projectors::LinOp_UNSYM1(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00212
00213 int branch_code;
00214
00215 branch_code = 10*(B != 0) + (APrec != 0);
00216
00217 Epetra_Vector w2(View, X.Map(), work2);
00218
00219 int iX;
00220 int nX = X.NumVectors();
00221
00222 if (branch_code == 11) {
00223
00224 for (iX = 0; iX < nX; ++iX) {
00225 Epetra_Vector Xi(View, X, iX);
00226 Epetra_Vector Yi(View, Y, iX);
00227 A_theta_B(Xi, w2, work1);
00228 APrec->ApplyInverse(w2, Yi);
00229 }
00230 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, Y, Y);
00231 }
00232 else if (branch_code == 10) {
00233
00234 for (iX = 0; iX < nX; ++iX) {
00235 Epetra_Vector Xi(View, X, iX);
00236 Epetra_Vector Yi(View, Y, iX);
00237 A_theta_B(Xi, Yi, work1);
00238 }
00239 Project1(Pro_Q, Pro_Qb, Y, Y);
00240 }
00241 else if (branch_code == 01) {
00242
00243 for (iX = 0; iX < nX; ++iX) {
00244 Epetra_Vector Xi(View, X, iX);
00245 Epetra_Vector Yi(View, Y, iX);
00246 A_theta_I(Xi, w2);
00247 APrec->ApplyInverse(w2, Yi);
00248 }
00249 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, Y, Y);
00250 }
00251 else {
00252
00253 A_theta_I(X, Y);
00254 Project1(Pro_Q, Pro_Q, Y, Y);
00255 }
00256
00257 }
00258
00259
00260 void Projectors::LinOp_UNSYM2(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00261
00262 int branch_code;
00263
00264 branch_code = 10*(B != 0) + (APrec != 0);
00265
00266 int iX;
00267 int nX = X.NumVectors();
00268
00269 if (branch_code == 00) {
00270
00271 A_theta_I(X, Y);
00272 Project1(Pro_Q, Pro_Q, Y, Y);
00273 }
00274 else if (branch_code == 10) {
00275 Epetra_Vector w2(View, X.Map(), work2);
00276
00277 for (iX = 0; iX < nX; ++iX) {
00278 Epetra_Vector Xi(View, X, iX);
00279 Epetra_Vector Yi(View, Y, iX);
00280 Project1(Pro_Q, Pro_Qb, Xi, w2);
00281 A_theta_B(w2, Yi, work1);
00282 }
00283 Project1(Pro_Qb, Pro_Q, Y, Y);
00284 }
00285 else if (branch_code == 01) {
00286 Epetra_Vector w2(View, X.Map(), work2);
00287
00288 for (iX = 0; iX < nX; ++iX) {
00289 Epetra_Vector Xi(View, X, iX);
00290 Epetra_Vector Yi(View, Y, iX);
00291 A_theta_I(Xi, w2);
00292 APrec->ApplyInverse(w2, Yi);
00293 }
00294 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, Y, Y);
00295 }
00296 else {
00297 Epetra_Vector w1(View, X.Map(), work1);
00298 Epetra_Vector w2(View, X.Map(), work2);
00299
00300
00301 for (iX = 0; iX < nX; ++iX) {
00302 Epetra_Vector Xi(View, X, iX);
00303 Epetra_Vector Yi(View, Y, iX);
00304 Project1(Pro_Q, Pro_Qb, Xi, w2);
00305 A_theta_B(w2, w1, Yi.Values());
00306 APrec->ApplyInverse(w1, Yi);
00307 }
00308 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, Y, Y);
00309 Project2(Pro_Q, Pro_k, Pro_Gpiv, Pro_Glu, Pro_Q, Y, Y);
00310 }
00311
00312 }
00313
00314
00315 void Projectors::LinOp_SYM(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00316
00317 if (B == 0) {
00318
00319 A_theta_I(X, Y);
00320 Project1(Pro_Q, Pro_Q, Y, Y);
00321 } else {
00322
00323 Project1(Pro_Q, Pro_Qb, X, Y);
00324 int iV;
00325 int nV = X.NumVectors();
00326 int rx = X.MyLength();
00327 Epetra_Vector w2(View, X.Map(), work2);
00328 for (iV = 0; iV < nV; ++iV) {
00329 Epetra_Vector Yv(View, Y, iV);
00330 memcpy(work2, Yv.Values(), rx*sizeof(double));
00331 A_theta_B(w2, Yv, work1);
00332 }
00333 Project1(Pro_Qb, Pro_Q, Y, Y);
00334 }
00335
00336 }
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346 void Projectors::Precon_UNSYM1(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00347
00348
00349 memcpy(Y.Values(), X.Values(), X.MyLength()*X.NumVectors()*sizeof(double));
00350
00351 }
00352
00353
00354 void Projectors::Precon_UNSYM2(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00355
00356
00357 memcpy(Y.Values(), X.Values(), X.MyLength()*X.NumVectors()*sizeof(double));
00358
00359 }
00360
00361
00362 void Projectors::Precon_SYM(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00363
00364 if (APrec == 0) {
00365 memcpy(Y.Values(), X.Values(), X.MyLength()*X.NumVectors()*sizeof(double));
00366 } else {
00367
00368 APrec->ApplyInverse(X, Y);
00369 if (B == 0) {
00370
00371 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, Y, Y);
00372 } else {
00373
00374 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, Y, Y);
00375 Project2(Pro_Q, Pro_k, Pro_Gpiv, Pro_Glu, Pro_Q, Y, Y);
00376 }
00377 }
00378
00379 }
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 void Projectors::Right_UNSYM1(Epetra_Vector &r){
00390
00391 if((B) && (APrec)){
00392 Epetra_Vector w1(View, r.Map(), work1);
00393
00394 APrec->ApplyInverse(r, w1);
00395 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, w1, r);
00396 }
00397 else if((B) && (APrec == 0))
00398
00399 Project1(Pro_Q, Pro_Qb, r, r);
00400 else if((B == 0) && (APrec)){
00401 Epetra_Vector w1(View, r.Map(), work1);
00402
00403 APrec->ApplyInverse(r, w1);
00404 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, w1, r);
00405 }
00406 else {
00407
00408 orthoTool.ModifiedGS(&r, Pro_Q);
00409 }
00410
00411 }
00412
00413
00414 void Projectors::Right_UNSYM2(Epetra_Vector &r) {
00415
00416 if (B == 0) {
00417 if (APrec == 0)
00418
00419 orthoTool.ModifiedGS(&r, Pro_Q);
00420 else {
00421 Epetra_Vector w1(View, r.Map(), work1);
00422
00423 APrec->ApplyInverse(r, w1);
00424 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, w1, r);
00425 }
00426 }
00427 else {
00428 if (APrec == 0)
00429
00430 Project1(Pro_Qb, Pro_Q, r, r);
00431 else {
00432
00433 Epetra_Vector w1(View, r.Map(), work1);
00434 APrec->ApplyInverse(r, w1);
00435 Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, w1, r);
00436 Project2(Pro_Q, Pro_k, Pro_Gpiv, Pro_Glu, Pro_Q, r, r);
00437 }
00438 }
00439 }
00440
00441
00442 void Projectors::Right_SYM(Epetra_Vector &r) {
00443
00444 if (B == 0) {
00445
00446 orthoTool.ModifiedGS(&r, Pro_Q);
00447 } else {
00448
00449 Project1(Pro_Qb, Pro_Q, r, r);
00450 }
00451
00452 }
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462 void Projectors::A_theta_B(const Epetra_MultiVector &X, Epetra_MultiVector &Y, double *wTmp)
00463 const {
00464
00465 A->Apply(X, Y);
00466
00467 int iV;
00468 int nV = X.NumVectors();
00469 int rX = X.MyLength();
00470
00471 for (iV = 0; iV < nV; ++iV) {
00472 Epetra_MultiVector Xv(View, X, iV, 1);
00473 Epetra_MultiVector w(View, X.Map(), wTmp, rX, 1);
00474 B->Apply(Xv, w);
00475 callBLAS.AXPY(rX, -Pro_theta, wTmp, Y.Values() + iV*rX);
00476 }
00477
00478 }
00479
00480
00481 void Projectors::A_theta_I(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00482
00483 A->Apply(X, Y);
00484 callBLAS.AXPY(X.MyLength()*X.NumVectors(), -Pro_theta, X.Values(), Y.Values());
00485
00486 }
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506 void Projectors::Project1(Epetra_MultiVector *AA, Epetra_MultiVector *BB,
00507 const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00508
00509 int rA = AA->MyLength();
00510 int cA = AA->NumVectors();
00511
00512 int rX = X.MyLength();
00513 int cX = X.NumVectors();
00514
00515 int size = cX*cA;
00516 double *wTmp = new double[2*size];
00517
00518
00519 if (X.Values() != Y.Values())
00520 memcpy(Y.Values(), X.Values(), rX*cX*sizeof(double));
00521
00522
00523 callBLAS.GEMM('T', 'N', cA, cX, rA, 1.0, BB->Values(), rA, Y.Values(), rX,
00524 0.0, wTmp + size, cA);
00525 MyComm.SumAll(wTmp + size, wTmp, size);
00526
00527
00528 callBLAS.GEMM('N', 'N', rA, cX, cA, -1.0, AA->Values(), rA, wTmp, cA,
00529 1.0, Y.Values(), rX);
00530
00531 delete[] wTmp;
00532
00533 }
00534
00535
00536
00537
00538
00539
00540
00541
00542 void Projectors::Project2(Epetra_MultiVector *AA, int rowB, int *piv, double *Blu,
00543 Epetra_MultiVector *CC, const Epetra_MultiVector &X,
00544 Epetra_MultiVector &Y) const {
00545
00546 int info;
00547
00548 int rowA = AA->MyLength();
00549 int colA = AA->NumVectors();
00550
00551 int rowC = CC->MyLength();
00552 int colC = CC->NumVectors();
00553
00554 int rowX = X.MyLength();
00555 int colX = X.NumVectors();
00556
00557 int size = colC*colX;
00558 double *wTmp = new double[2*size];
00559
00560
00561 if (X.Values() != Y.Values())
00562 memcpy(Y.Values(), X.Values(), rowX*colX*sizeof(double));
00563
00564
00565 callBLAS.GEMM('T', 'N', colC, colX, rowC, 1.0, CC->Values(), rowC, X.Values(), rowX,
00566 0.0, wTmp + size, colC);
00567 MyComm.SumAll(wTmp + size, wTmp, size);
00568
00569
00570 callLAPACK.GETRS('N', rowB, colX, Blu, Pro_kmax, piv, wTmp, colC, &info);
00571 assert(info == 0);
00572
00573
00574 callBLAS.GEMM('N', 'N', rowA, colX, colA, -1.0, AA->Values(), rowA, wTmp, colC,
00575 1.0, Y.Values(), rowX);
00576
00577 delete[] wTmp;
00578
00579 }
00580
00581