src/eigsolv/CheckingTools.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 #include <sstream>
00030 #include "rlog/rlog.h"
00031 #include "CheckingTools.h"
00032 
00033 using namespace std;
00034 
00035 CheckingTools::CheckingTools(const Epetra_Comm &_Comm) :
00036     MyComm(_Comm) 
00037 {
00038 }
00039 
00040 
00041 double CheckingTools::errorOrthogonality(const Epetra_MultiVector *X, 
00042                                          const Epetra_MultiVector *R, 
00043                                          const Epetra_Operator *M) const {
00044 
00045     // Return the maximum value of R_i^T * M * X_j / || MR_i || || X_j ||
00046     // When M is not specified, the identity is used.
00047     double maxDot = 0.0;
00048 
00049     int xc = (X) ? X->NumVectors() : 0;
00050     int rc = (R) ? R->NumVectors() : 0;
00051   
00052     if (xc*rc == 0)
00053         return maxDot;
00054 
00055     int i, j;
00056     for (i = 0; i < rc; ++i) {
00057         Epetra_Vector MRi(Copy, (*R), i);
00058         if (M)
00059             M->Apply(*((*R)(i)), MRi);
00060         double normMR = 0.0;
00061         MRi.Norm2(&normMR);
00062         double dot = 0.0;
00063         for (j = 0; j < xc; ++j) {
00064             double normXj = 0.0;
00065             (*X)(j)->Norm2(&normXj);
00066             (*X)(j)->Dot(MRi, &dot);
00067             dot = fabs(dot)/(normMR*normXj);
00068             maxDot = (dot > maxDot) ? dot : maxDot;
00069         }
00070     }
00071 
00072     return maxDot;
00073 
00074 }
00075 
00076 
00077 double CheckingTools::errorOrthonormality(const Epetra_MultiVector *X, 
00078                                           const Epetra_Operator *M) const {
00079 
00080     // Return the maximum coefficient of the matrix X^T * M * X - I
00081     // When M is not specified, the identity is used.
00082     double maxDot = 0.0;
00083 
00084     int xc = (X) ? X->NumVectors() : 0;
00085     if (xc == 0)
00086         return maxDot;
00087 
00088     int i, j;
00089     for (i = 0; i < xc; ++i) {
00090         Epetra_Vector MXi(Copy, (*X), i);
00091         if (M)
00092             M->Apply(*((*X)(i)), MXi);
00093         double dot = 0.0;
00094         for (j = 0; j < xc; ++j) {
00095             (*X)(j)->Dot(MXi, &dot);
00096             dot = (i == j) ? fabs(dot - 1.0) : fabs(dot);
00097             maxDot = (dot > maxDot) ? dot : maxDot;
00098         }
00099     }
00100 
00101     return maxDot;
00102 
00103 }
00104 
00105 
00106 double CheckingTools::errorEquality(const Epetra_MultiVector *X,
00107                                     const Epetra_MultiVector *MX, const Epetra_Operator *M) const {
00108 
00109     // Return the maximum coefficient of the matrix M * X - MX
00110     // scaled by the maximum coefficient of MX.
00111     // When M is not specified, the identity is used.
00112 
00113     double maxDiff = 0.0;
00114 
00115     int xc = (X) ? X->NumVectors() : 0;
00116     int mxc = (MX) ? MX->NumVectors() : 0;
00117 
00118     if ((xc != mxc) || (xc*mxc == 0))
00119         return maxDiff;
00120 
00121     int i;
00122     double maxCoeffX = 0.0;
00123     for (i = 0; i < xc; ++i) {
00124         double tmp = 0.0;
00125         (*MX)(i)->NormInf(&tmp);
00126         maxCoeffX = (tmp > maxCoeffX) ? tmp : maxCoeffX;
00127     }
00128 
00129     for (i = 0; i < xc; ++i) {
00130         Epetra_Vector MtimesXi(Copy, (*X), i);
00131         if (M)
00132             M->Apply(*((*X)(i)), MtimesXi);
00133         MtimesXi.Update(-1.0, *((*MX)(i)), 1.0);
00134         double diff = 0.0;
00135         MtimesXi.NormInf(&diff);
00136         maxDiff = (diff > maxDiff) ? diff : maxDiff;
00137     }
00138 
00139     return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;
00140 
00141 }
00142 
00143 
00144 int CheckingTools::errorSubspaces(const Epetra_MultiVector &Q, const Epetra_MultiVector &Qex,
00145                                   const Epetra_Operator *M) const {
00146 
00147     int info = 0;
00148 
00149     int qr = Q.MyLength();
00150     int qc = Q.NumVectors();
00151     int qexc = Qex.NumVectors();
00152 
00153     double *mQ = new (nothrow) double[qr];
00154     if (mQ == 0) {
00155         info = -1;
00156         return info;
00157     }
00158 
00159     double *z = new (nothrow) double[qexc*qc];
00160     if (z == 0) {
00161         delete[] mQ;
00162         info = -1;
00163         return info;
00164     }
00165 
00166     Epetra_LocalMap lMap(qexc, 0, MyComm);
00167     Epetra_MultiVector QextMQ(View, lMap, z, qexc, qc);
00168 
00169     int j;
00170     for (j=0; j<qc; ++j) {
00171         Epetra_MultiVector Qj(View, Q, j, 1);
00172         Epetra_MultiVector MQ(View, Q.Map(), mQ, qr, 1);
00173         if (M)
00174             M->Apply(Qj, MQ);
00175         else
00176             memcpy(mQ, Qj.Values(), qr*sizeof(double));
00177         Epetra_MultiVector colJ(View, QextMQ, j, 1);
00178         colJ.Multiply('T', 'N', 1.0, Qex, MQ,  0.0);
00179     }
00180     delete[] mQ;
00181 
00182     // Compute the SVD
00183 
00184     int svSize = (qc > qexc) ? qexc : qc;
00185     double *sv = new (nothrow) double[svSize];
00186     if (sv == 0) {
00187         delete[] z;
00188         info  = -1;
00189         return info;
00190     }
00191 
00192     // lwork is larger than the value suggested by LAPACK
00193     int lwork = (qexc > qc) ? qexc + 5*qc : qc + 5*qexc;
00194     double *work = new (nothrow) double[lwork];
00195     if (work == 0) {
00196         delete[] z;
00197         delete[] sv;
00198         info = -1;
00199         return info;
00200     }
00201 
00202     Epetra_LAPACK call;
00203     call.GESVD('N','N',qexc,qc,QextMQ.Values(),qexc,sv,0,qc,0,qc,work,&lwork,&info);
00204 
00205     delete[] work;
00206     delete[] z;
00207 
00208     // Treat error messages
00209 
00210     if (info < 0) {
00211         rError("In DGESVD, argument %d has an illegal value", -info);
00212         delete[] sv;
00213         return info;
00214     }
00215 
00216     if (info > 0) {
00217         rError("In DGESVD, DBSQR did not converge (%d).", info);
00218         delete[] sv;
00219         return info;
00220     }
00221 
00222     // Output the result
00223     ostringstream buf;
00224     buf << "Principal angles between eigenspaces:" << endl;
00225     buf << "\tReference space built with " << Qex.NumVectors() << " vectors." << endl;
00226     buf << "\tDimension of computed space: " << Q.NumVectors() << endl;
00227     buf.setf(ios::scientific, ios::floatfield);
00228     buf.precision(4);
00229     buf << "\tSmallest singular value = " << sv[0] << endl;
00230     buf << "\tLargest singular value  = " << sv[svSize-1] << endl;
00231     buf << "\tSmallest angle between subspaces (rad) = ";
00232     buf << ((sv[0] > 1.0) ? 0.0 : asin(sqrt(1.0 - sv[0]*sv[0]))) << endl;
00233     buf << "\tLargest angle between subspaces (rad)  = ";
00234     buf << ((sv[0] > 1.0) ? 0.0 : asin(sqrt(1.0 - sv[svSize-1]*sv[svSize-1]))) << endl;
00235     rDebug(buf.str().c_str());
00236 
00237     delete[] sv;
00238 
00239     return info;
00240 }
00241 
00242 
00243 void CheckingTools::errorEigenResiduals(const Epetra_MultiVector &Q, double *lambda,
00244                                         const Epetra_Operator *K, const Epetra_Operator *M,
00245                                         double *normWeight) const {
00246     ostringstream buf;
00247     if ((K == 0) || (lambda == 0))
00248         return;
00249 
00250     int qr = Q.MyLength();
00251     int qc = Q.NumVectors();
00252 
00253     double *work = new (nothrow) double[2*qr];
00254     if (work == 0)
00255         return;
00256 
00257     buf << "Norms of residuals for computed eigenmodes:" << endl;
00258     buf << "        Eigenvalue";
00259     if (normWeight)
00260         buf << "     User Norm     Scaled User N.";
00261     buf << "     2-Norm     Scaled 2-Nor.\n";
00262 
00263     Epetra_Vector KQ(View, Q.Map(), work);
00264     Epetra_Vector MQ(View, Q.Map(), work + qr);
00265     Epetra_Vector Qj(View, Q.Map(), Q.Values());
00266 
00267     double maxUserNorm = 0.0;
00268     double minUserNorm = 1.0e+16;
00269     double maxL2Norm = 0.0;
00270     double minL2Norm = 1.0e+16;
00271 
00272     // Compute the residuals and norms
00273     int j;
00274     for (j=0; j<qc; ++j) {
00275 
00276         Qj.ResetView(Q.Values() + j*qr);
00277         if (M)
00278             M->Apply(Qj, MQ);
00279         else
00280             memcpy(MQ.Values(), Qj.Values(), qr*sizeof(double));
00281         K->Apply(Qj, KQ);
00282         KQ.Update(-lambda[j], MQ, 1.0);
00283 
00284         double residualL2 = 0.0;
00285         KQ.Norm2(&residualL2);
00286 
00287         double residualUser = 0.0;
00288         if (normWeight) {
00289             Epetra_Vector vectWeight(View, Q.Map(), normWeight);
00290             KQ.NormWeighted(vectWeight, &residualUser);
00291         }
00292 
00293         buf << " ";
00294         buf.width(4);
00295         buf.precision(8);
00296         buf.setf(ios::scientific, ios::floatfield);
00297         buf << j+1 << ". " << lambda[j] << " ";
00298         if (normWeight)
00299             buf << residualUser << " " << residualUser/lambda[j] << " ";
00300         buf << residualL2 << " " << residualL2/lambda[j] << " ";
00301         buf << endl;
00302 
00303         maxL2Norm = (residualL2/lambda[j] > maxL2Norm) ? residualL2/lambda[j] : maxL2Norm;
00304         minL2Norm = (residualL2/lambda[j] < minL2Norm) ? residualL2/lambda[j] : minL2Norm;
00305         if (normWeight) {
00306             maxUserNorm = (residualUser/lambda[j] > maxUserNorm) ? residualUser/lambda[j]
00307                 : maxUserNorm;
00308             minUserNorm = (residualUser/lambda[j] < minUserNorm) ? residualUser/lambda[j]
00309                 : minUserNorm;
00310         }
00311 
00312     } // for j=0; j<qc; ++j)
00313 
00314     if (normWeight) {
00315         buf << " >> Minimum scaled user-norm of residuals = " << minUserNorm << endl;
00316         buf << " >> Maximum scaled user-norm of residuals = " << maxUserNorm << endl;
00317     }
00318     buf << " >> Minimum scaled L2 norm of residuals   = " << minL2Norm << endl;
00319     buf << " >> Maximum scaled L2 norm of residuals   = " << maxL2Norm;
00320 
00321     if (work)
00322         delete[] work;
00323 
00324     rInfo(buf.str().c_str());
00325 }
00326 
00327 
00328 int CheckingTools::errorLambda(double *continuous, double *discrete, int numDiscrete,
00329                                double *lambda, int nev) const {
00330     ostringstream buf;
00331     int myPid = MyComm.MyPID();
00332     int nMax = 0;
00333     int i, j;
00334 
00335     // Allocate working arrays
00336 
00337     int *used = new (nothrow) int[numDiscrete + nev];
00338     if (used == 0) {
00339         return nMax;
00340     }
00341 
00342     int *bestMatch = used + numDiscrete;
00343 
00344     // Find the best match for the eigenvalues
00345     double eps = 0.0;
00346     Epetra_LAPACK call;
00347     call.LAMCH('E', eps);
00348 
00349     double gap = Epetra_MaxDouble;
00350     for (i=0; i<numDiscrete; ++i) {
00351         used[i] = -1;
00352         for (int j = i; j < numDiscrete; ++j)
00353             if (discrete[j] > (1.0 + 10.0*eps)*discrete[i]) {
00354                 double tmp = (discrete[j] - discrete[i])/discrete[i];
00355                 gap = (tmp < gap) ? tmp : gap;
00356                 break;
00357             }
00358     }
00359 
00360     for (i=0; i<nev; ++i) 
00361         bestMatch[i] = -1;
00362 
00363     for (i=0; i<nev; ++i) {
00364         if (lambda[i] < continuous[0])
00365             continue;
00366         bestMatch[i] = (i == 0) ? 0 : bestMatch[i-1] + 1;
00367         int jStart = bestMatch[i];
00368         for (j = jStart; j < numDiscrete; ++j) {
00369             double diff = fabs(lambda[i]-discrete[j]);
00370             if (diff < 0.5*gap*lambda[i]) {
00371                 bestMatch[i] = j;
00372                 break;
00373             }
00374         }
00375         used[bestMatch[i]] = i;
00376     }
00377 
00378     // Print the results for the eigenvalues
00379     if (myPid == 0) {
00380         buf << endl;
00381         buf << " --- Relative errors on eigenvalues ---\n";
00382         buf << endl;
00383         buf << "       Exact Cont.    Exact Disc.     Computed      ";
00384         buf << " Alg. Err.  Mesh Err.\n";
00385     }
00386 
00387     int iCount = 0;
00388     for (i=0; i<nev; ++i) {
00389         if (bestMatch[i] == -1) {
00390             if (myPid == 0) {
00391                 buf << "      ************** ************** ";
00392                 buf.precision(8);
00393                 buf.setf(ios::scientific, ios::floatfield);
00394                 buf << lambda[i];
00395                 buf << "  *********  *********" << endl;
00396             }
00397             iCount += 1;
00398         }
00399     }
00400 
00401     double lastDiscrete = 0.0;
00402     for (i=0; i<numDiscrete; ++i) {
00403         if ((iCount == nev) && (discrete[i] > lastDiscrete))
00404             break;
00405         if (used[i] < 0) {
00406             nMax += 1;
00407             lastDiscrete = discrete[i];
00408             if (myPid == 0) {
00409                 buf << " ";
00410                 buf.width(4);
00411                 buf << i+1 << ". ";
00412                 buf.precision(8);
00413                 buf.setf(ios::scientific, ios::floatfield);
00414                 buf << continuous[i] << " " << discrete[i] << " ";
00415                 buf << "**************  *********  ";
00416                 buf.precision(3);
00417                 buf << fabs(continuous[i]-discrete[i])/continuous[i] << endl;
00418             }
00419         }
00420         else {
00421             nMax += 1;
00422             lastDiscrete = discrete[i];
00423             if (myPid == 0) {
00424                 buf << " ";
00425                 buf.width(4);
00426                 buf << i+1 << ". ";
00427                 buf.precision(8);
00428                 buf.setf(ios::scientific, ios::floatfield);
00429                 buf << continuous[i] << " " << discrete[i] << " " << lambda[used[i]] << "  ";
00430                 buf.precision(3);
00431                 buf << fabs(lambda[used[i]]-discrete[i])/discrete[i] << "  ";
00432                 buf << fabs(continuous[i]-discrete[i])/continuous[i] << endl;
00433             }
00434             iCount += 1;
00435         }
00436     }
00437 
00438     delete[] used;
00439 
00440     if (myPid == 0)
00441         rDebug(buf.str().c_str());
00442 
00443     return nMax;
00444 }
00445 
00446 
00447 int CheckingTools::inputArguments(const int &numEigen, const Epetra_Operator *K, 
00448                                   const Epetra_Operator *M, const Epetra_Operator *P,
00449                                   const Epetra_MultiVector &Q, const int &minSize) const {
00450 
00451     // Routine to check some arguments
00452     // 
00453     // info = -  1 >> The stiffness matrix K has not been specified.
00454     // info = -  2 >> The maps for the matrix K and the matrix M differ.
00455     // info = -  3 >> The maps for the matrix K and the preconditioner P differ.
00456     // info = -  4 >> The maps for the vectors and the matrix K differ.
00457     // info = -  5 >> Q is too small for the number of eigenvalues requested.
00458     // info = -  6 >> Q is too small for the computation parameters.
00459     //
00460 
00461     int myPid = MyComm.MyPID();
00462 
00463     if (K == 0) {
00464         if (myPid == 0)
00465             cerr << "\n !!! The matrix K to analyze has not been specified !!! \n\n";
00466         return -1;
00467     }
00468 
00469     if (M) {
00470         int mGlobal = (M->OperatorDomainMap()).NumGlobalElements();
00471         int mLocal = (M->OperatorDomainMap()).NumMyElements();
00472         int kGlobal = (K->OperatorDomainMap()).NumGlobalElements();
00473         int kLocal = (K->OperatorDomainMap()).NumMyElements();
00474         if ((mGlobal != kGlobal) || (mLocal != kLocal)) {
00475             if (myPid == 0) {
00476                 cerr << endl;
00477                 cerr << " !!! The maps for the matrice K and the mass M are different !!!\n";
00478                 cerr << endl;
00479             }
00480             return -2;
00481         }
00482     }
00483 
00484     if (P) {
00485         int pGlobal = (P->OperatorDomainMap()).NumGlobalElements();
00486         int pLocal = (P->OperatorDomainMap()).NumMyElements();
00487         int kGlobal = (K->OperatorDomainMap()).NumGlobalElements();
00488         int kLocal = (K->OperatorDomainMap()).NumMyElements();
00489         if ((pGlobal != kGlobal) || (pLocal != kLocal)) {
00490             if (myPid == 0) {
00491                 cerr << endl;
00492                 cerr << " !!! The maps for the matrice K and the preconditioner P are different !!!\n";
00493                 cerr << endl;
00494             }
00495             return -3;
00496         }
00497     }
00498 
00499     if ((Q.MyLength() != (K->OperatorDomainMap()).NumMyElements()) || 
00500         (Q.GlobalLength() != (K->OperatorDomainMap()).NumGlobalElements())) {
00501         if (myPid == 0) {
00502             cerr << "\n !!! The maps for the vectors and the matrix are different !!! \n\n";
00503         }
00504         return -4;
00505     }
00506 
00507     if (Q.NumVectors() < numEigen) {
00508         if (myPid == 0) {
00509             cerr << endl;
00510             cerr << " !!! The number of eigenvalues is too large for the space allocated !!! \n\n";
00511             cerr << " The recommended size for " << numEigen << " eigenvalues is ";
00512             cerr << minSize << endl;
00513             cerr << endl;
00514         }
00515         return -5;
00516     }
00517 
00518     if (Q.NumVectors() < minSize) {
00519         if (myPid == 0) {
00520             cerr << endl;
00521             cerr << " !!! The space allocated is too small for the number of eigenvalues";
00522             cerr << " and the size of blocks specified !!! \n";
00523             cerr << " The recommended size for " << numEigen << " eigenvalues is ";
00524             cerr << minSize << endl;
00525             cerr << endl;
00526         }
00527         return -6;
00528     }
00529 
00530     return 0;
00531 
00532 }
00533 
00534 

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