src/eigsolv/Projectors.cpp

Go to the documentation of this file.
00001 //**************************************************************************
00002 //
00003 //                                 NOTICE
00004 //
00005 // This software is a result of the research described in the report
00006 //
00007 // " A comparison of algorithms for modal analysis in the absence 
00008 //   of a sparse direct method", P. Arbenz, R. Lehoucq, and U. Hetmaniuk,
00009 //  Sandia National Laboratories, Technical report SAND2003-1028J.
00010 //
00011 // It is based on the Epetra, AztecOO, and ML packages defined in the Trilinos
00012 // framework ( http://software.sandia.gov/trilinos/ ).
00013 //
00014 // The distribution of this software follows also the rules defined in Trilinos.
00015 // This notice shall be marked on any reproduction of this software, in whole or
00016 // in part.
00017 //
00018 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00019 // license for use of this work by or on behalf of the U.S. Government.
00020 //
00021 // This program is distributed in the hope that it will be useful, but
00022 // WITHOUT ANY WARRANTY; without even the implied warranty of
00023 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00024 //
00025 // Code Authors: U. Hetmaniuk (ulhetma@sandia.gov), R. Lehoucq (rblehou@sandia.gov)
00026 //
00027 //**************************************************************************
00028 /**************************************************************************
00029  *                                                                         *
00030  *               Swiss Federal Institute of Technology (ETH)               *
00031  *                       CH-8092 Zuerich, Switzerland                      *
00032  *                                                                         *
00033  *                       (C) 1999 All Rights Reserved                      *
00034  *                                                                         *
00035  *                                NOTICE                                   *
00036  *                                                                         *
00037  *  Permission to use, copy, modify, and distribute this software and      *
00038  *  its documentation for any purpose and without fee is hereby granted    *
00039  *  provided that the above copyright notice appear in all copies and      *
00040  *  that both the copyright notice and this permission notice appear in    *
00041  *  supporting documentation.                                              *
00042  *                                                                         *
00043  *  Neither the Swiss Federal Institute of Technology nor the author make  *
00044  *  any representations about the suitability of this software for any     *
00045  *  purpose.  This software is provided ``as is'' without express or       *
00046  *  implied warranty.                                                      *
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     // This routine assumes that at least one Epetra_Operator is not empty
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     // This routine assumes that at least one Epetra_Operator is not empty
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 // Operator for the lin.eq.-solver used during the correction equation
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     }  // if (A)
00140     timeApply += MyWatch.WallTime();
00141 
00142     return 0;
00143 
00144 }
00145 
00146 
00147 /* Preconditioner for the lin.eq.-solver used during the correction equation */
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     } // if (APrec)
00169     timeApplyInverse += MyWatch.WallTime();
00170 
00171     return 0;
00172 
00173 }
00174 
00175 
00176 /* Preconditions and projects the right hand side of the correction equation */
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  * Linear Operators                                                         *
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         /* y := (I - Y'*inv(H)*Qb')*inv(K)*(A-theta*B)*x */
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         /* y := (I - Q*Qb')*(A-theta*B)*x */
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         /* y := (I - Y'*inv(H)*Q')*inv(K)*(A-theta*I)*x */
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 { /* (Pro_Bisempty == 1) && (Pro_Kisempty == 1) */
00252         /* y := (I - Q*Q')*(A-theta*I)*x */
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         /* y = (I - Q*Q')*(A - theta*I) * x */
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         /* y = (I - Qb*Q')*(A-theta*B)*(I - Q*Qb') * x */   
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         /* y = (I - Y'*inv(H)*Qb')*inv(K)*(A-theta*B) * x */   
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 /* branch_code == 11 */ {
00297         Epetra_Vector w1(View, X.Map(), work1);
00298         Epetra_Vector w2(View, X.Map(), work2);
00299         /* y = (I - Q*inv(G)*Q')* ...
00300            (I - Y'*inv(H)*Qb')*inv(K)*(A-theta*B)*(I - Q*Qb')*x */   
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         /* y := (I - Q*Q')*(A - theta*I)*x */
00319         A_theta_I(X, Y);
00320         Project1(Pro_Q, Pro_Q, Y, Y);
00321     } else {
00322         /* y := (I - Qb*Q')*(A - theta*I)*(I - Q*Qb')*x */
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     } // if (B == 0)
00335 
00336 }
00337 
00338 
00339 /****************************************************************************
00340  *                                                                          *
00341  * Preconditioners                                                          *
00342  *                                                                          *
00343  ****************************************************************************/
00344 
00345  
00346 void Projectors::Precon_UNSYM1(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {
00347 
00348     /* the Preconditioner is contained in the operator */
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     /* the Preconditioner is contained in the operator */
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         /* y := inv(K) * x */
00368         APrec->ApplyInverse(X, Y);
00369         if (B == 0) {
00370             /* y = (I-Q*G\Q') * (I-Y*H\Q') * y */
00371             Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, Y, Y);
00372         } else /* B not empty */ {
00373             /* y = (I-Q*G\Q') * (I-Y*H\Qb') * y */
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  * Right hand sides                                                         *
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         /* r := (I - Y'*inv(H)*Qb')*inv(K)*r */
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         /* r :=  (I - Q*Qb')*r */
00399         Project1(Pro_Q, Pro_Qb, r, r);
00400     else if((B == 0) && (APrec)){
00401         Epetra_Vector w1(View, r.Map(), work1);
00402         /* r := (I - Y'*inv(H)*Q')*inv(K)*r */
00403         APrec->ApplyInverse(r, w1);
00404         Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, w1, r);
00405     }
00406     else {
00407         /* r :=  (I - Q*Q')*r */
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             /* r = (I-Q*Q') * r */
00419             orthoTool.ModifiedGS(&r, Pro_Q);
00420         else {
00421             Epetra_Vector w1(View, r.Map(), work1);
00422             /* r = (I-Y*H\Q') * inv(K) * r */
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             /* r = (I - Qb*Q')*r */
00430             Project1(Pro_Qb, Pro_Q, r, r);
00431         else {
00432             /* r = (I-Q*G\Q') * (I-Y*H\Qb') * inv(K) * r */
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         /* r :=  (I - Q*Q')*r */
00446         orthoTool.ModifiedGS(&r, Pro_Q);
00447     } else { /* B not empty */
00448         /* r = r - Qb*(Q'*r); */
00449         Project1(Pro_Qb, Pro_Q, r, r); 
00450     }
00451 
00452 }
00453 
00454 
00455 /****************************************************************************
00456  *                                                                          *
00457  * Multiplication Operators                                                 *
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  * Helper routines                                                          *
00492  *                                                                          *
00493  ****************************************************************************/
00494 
00495 
00496 /* PROJECT1 - compute y := (I - A*B')*x
00497  *
00498  * sizes of matrices and vectors:
00499  *   A and B     --  r-by-c matrices, ldim = r
00500  *   x, y and w  --  r vectors
00501  *
00502  * (x, y) and also (A, B) may point to the same memory locations
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     /* if x and y do not refer to the same vector, we copy x to y */
00519     if (X.Values() != Y.Values())
00520         memcpy(Y.Values(), X.Values(), rX*cX*sizeof(double));
00521 
00522     /* W = B'*Y */
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     /* Y = Y - A*W */
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 /* PROJECT2 - compute y := (I - A*inv(B)*C')*x
00537  *
00538  * (x, y) may point to the same memory location
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     /* if x and y do not refer to the same vector, we copy x to y */
00561     if (X.Values() != Y.Values())
00562         memcpy(Y.Values(), X.Values(), rowX*colX*sizeof(double));
00563 
00564     /* W = C'*X */
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     /* W = inv(B)*W */
00570     callLAPACK.GETRS('N', rowB, colX, Blu, Pro_kmax, piv, wTmp, colC, &info);
00571     assert(info == 0);
00572 
00573     /* Y = Y - A*W = Y - A*inv(B)*C'*X */
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 

Generated on Fri Oct 26 13:35:11 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7