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 #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
00046
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
00081
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
00110
00111
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
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
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
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
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
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 }
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
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
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
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
00452
00453
00454
00455
00456
00457
00458
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