src/eigsolv/CheckingTools.cpp

Go to the documentation of this file.

//**************************************************************************
//
//                                 NOTICE
//
// This software is a result of the research described in the report
//
// " A comparison of algorithms for modal analysis in the absence 
//   of a sparse direct method", P. Arbenz, R. Lehoucq, and U. Hetmaniuk,
//  Sandia National Laboratories, Technical report SAND2003-1028J.
//
// It is based on the Epetra, AztecOO, and ML packages defined in the Trilinos
// framework ( http://software.sandia.gov/trilinos/ ).
//
// The distribution of this software follows also the rules defined in Trilinos.
// This notice shall be marked on any reproduction of this software, in whole or
// in part.
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
//
// Code Authors: U. Hetmaniuk (ulhetma@sandia.gov), R. Lehoucq (rblehou@sandia.gov)
//
//**************************************************************************

#include <sstream>
#include "rlog/rlog.h"
#include "CheckingTools.h"

using namespace std;

CheckingTools::CheckingTools(const Epetra_Comm &_Comm) :
    MyComm(_Comm) 
{
}


double CheckingTools::errorOrthogonality(const Epetra_MultiVector *X, 
                                         const Epetra_MultiVector *R, 
                                         const Epetra_Operator *M) const {

    // Return the maximum value of R_i^T * M * X_j / || MR_i || || X_j ||
    // When M is not specified, the identity is used.
    double maxDot = 0.0;

    int xc = (X) ? X->NumVectors() : 0;
    int rc = (R) ? R->NumVectors() : 0;
  
    if (xc*rc == 0)
        return maxDot;

    int i, j;
    for (i = 0; i < rc; ++i) {
        Epetra_Vector MRi(Copy, (*R), i);
        if (M)
            M->Apply(*((*R)(i)), MRi);
        double normMR = 0.0;
        MRi.Norm2(&normMR);
        double dot = 0.0;
        for (j = 0; j < xc; ++j) {
            double normXj = 0.0;
            (*X)(j)->Norm2(&normXj);
            (*X)(j)->Dot(MRi, &dot);
            dot = fabs(dot)/(normMR*normXj);
            maxDot = (dot > maxDot) ? dot : maxDot;
        }
    }

    return maxDot;

}


double CheckingTools::errorOrthonormality(const Epetra_MultiVector *X, 
                                          const Epetra_Operator *M) const {

    // Return the maximum coefficient of the matrix X^T * M * X - I
    // When M is not specified, the identity is used.
    double maxDot = 0.0;

    int xc = (X) ? X->NumVectors() : 0;
    if (xc == 0)
        return maxDot;

    int i, j;
    for (i = 0; i < xc; ++i) {
        Epetra_Vector MXi(Copy, (*X), i);
        if (M)
            M->Apply(*((*X)(i)), MXi);
        double dot = 0.0;
        for (j = 0; j < xc; ++j) {
            (*X)(j)->Dot(MXi, &dot);
            dot = (i == j) ? fabs(dot - 1.0) : fabs(dot);
            maxDot = (dot > maxDot) ? dot : maxDot;
        }
    }

    return maxDot;

}


double CheckingTools::errorEquality(const Epetra_MultiVector *X,
                                    const Epetra_MultiVector *MX, const Epetra_Operator *M) const {

    // Return the maximum coefficient of the matrix M * X - MX
    // scaled by the maximum coefficient of MX.
    // When M is not specified, the identity is used.

    double maxDiff = 0.0;

    int xc = (X) ? X->NumVectors() : 0;
    int mxc = (MX) ? MX->NumVectors() : 0;

    if ((xc != mxc) || (xc*mxc == 0))
        return maxDiff;

    int i;
    double maxCoeffX = 0.0;
    for (i = 0; i < xc; ++i) {
        double tmp = 0.0;
        (*MX)(i)->NormInf(&tmp);
        maxCoeffX = (tmp > maxCoeffX) ? tmp : maxCoeffX;
    }

    for (i = 0; i < xc; ++i) {
        Epetra_Vector MtimesXi(Copy, (*X), i);
        if (M)
            M->Apply(*((*X)(i)), MtimesXi);
        MtimesXi.Update(-1.0, *((*MX)(i)), 1.0);
        double diff = 0.0;
        MtimesXi.NormInf(&diff);
        maxDiff = (diff > maxDiff) ? diff : maxDiff;
    }

    return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;

}


int CheckingTools::errorSubspaces(const Epetra_MultiVector &Q, const Epetra_MultiVector &Qex,
                                  const Epetra_Operator *M) const {

    int info = 0;

    int qr = Q.MyLength();
    int qc = Q.NumVectors();
    int qexc = Qex.NumVectors();

    double *mQ = new (nothrow) double[qr];
    if (mQ == 0) {
        info = -1;
        return info;
    }

    double *z = new (nothrow) double[qexc*qc];
    if (z == 0) {
        delete[] mQ;
        info = -1;
        return info;
    }

    Epetra_LocalMap lMap(qexc, 0, MyComm);
    Epetra_MultiVector QextMQ(View, lMap, z, qexc, qc);

    int j;
    for (j=0; j<qc; ++j) {
        Epetra_MultiVector Qj(View, Q, j, 1);
        Epetra_MultiVector MQ(View, Q.Map(), mQ, qr, 1);
        if (M)
            M->Apply(Qj, MQ);
        else
            memcpy(mQ, Qj.Values(), qr*sizeof(double));
        Epetra_MultiVector colJ(View, QextMQ, j, 1);
        colJ.Multiply('T', 'N', 1.0, Qex, MQ,  0.0);
    }
    delete[] mQ;

    // Compute the SVD

    int svSize = (qc > qexc) ? qexc : qc;
    double *sv = new (nothrow) double[svSize];
    if (sv == 0) {
        delete[] z;
        info  = -1;
        return info;
    }

    // lwork is larger than the value suggested by LAPACK
    int lwork = (qexc > qc) ? qexc + 5*qc : qc + 5*qexc;
    double *work = new (nothrow) double[lwork];
    if (work == 0) {
        delete[] z;
        delete[] sv;
        info = -1;
        return info;
    }

    Epetra_LAPACK call;
    call.GESVD('N','N',qexc,qc,QextMQ.Values(),qexc,sv,0,qc,0,qc,work,&lwork,&info);

    delete[] work;
    delete[] z;

    // Treat error messages

    if (info < 0) {
        rError("In DGESVD, argument %d has an illegal value", -info);
        delete[] sv;
        return info;
    }

    if (info > 0) {
        rError("In DGESVD, DBSQR did not converge (%d).", info);
        delete[] sv;
        return info;
    }

    // Output the result
    ostringstream buf;
    buf << "Principal angles between eigenspaces:" << endl;
    buf << "\tReference space built with " << Qex.NumVectors() << " vectors." << endl;
    buf << "\tDimension of computed space: " << Q.NumVectors() << endl;
    buf.setf(ios::scientific, ios::floatfield);
    buf.precision(4);
    buf << "\tSmallest singular value = " << sv[0] << endl;
    buf << "\tLargest singular value  = " << sv[svSize-1] << endl;
    buf << "\tSmallest angle between subspaces (rad) = ";
    buf << ((sv[0] > 1.0) ? 0.0 : asin(sqrt(1.0 - sv[0]*sv[0]))) << endl;
    buf << "\tLargest angle between subspaces (rad)  = ";
    buf << ((sv[0] > 1.0) ? 0.0 : asin(sqrt(1.0 - sv[svSize-1]*sv[svSize-1]))) << endl;
    rDebug(buf.str().c_str());

    delete[] sv;

    return info;
}


void CheckingTools::errorEigenResiduals(const Epetra_MultiVector &Q, double *lambda,
                                        const Epetra_Operator *K, const Epetra_Operator *M,
                                        double *normWeight) const {
    ostringstream buf;
    if ((K == 0) || (lambda == 0))
        return;

    int qr = Q.MyLength();
    int qc = Q.NumVectors();

    double *work = new (nothrow) double[2*qr];
    if (work == 0)
        return;

    buf << "Norms of residuals for computed eigenmodes:" << endl;
    buf << "        Eigenvalue";
    if (normWeight)
        buf << "     User Norm     Scaled User N.";
    buf << "     2-Norm     Scaled 2-Nor.\n";

    Epetra_Vector KQ(View, Q.Map(), work);
    Epetra_Vector MQ(View, Q.Map(), work + qr);
    Epetra_Vector Qj(View, Q.Map(), Q.Values());

    double maxUserNorm = 0.0;
    double minUserNorm = 1.0e+16;
    double maxL2Norm = 0.0;
    double minL2Norm = 1.0e+16;

    // Compute the residuals and norms
    int j;
    for (j=0; j<qc; ++j) {

        Qj.ResetView(Q.Values() + j*qr);
        if (M)
            M->Apply(Qj, MQ);
        else
            memcpy(MQ.Values(), Qj.Values(), qr*sizeof(double));
        K->Apply(Qj, KQ);
        KQ.Update(-lambda[j], MQ, 1.0);

        double residualL2 = 0.0;
        KQ.Norm2(&residualL2);

        double residualUser = 0.0;
        if (normWeight) {
            Epetra_Vector vectWeight(View, Q.Map(), normWeight);
            KQ.NormWeighted(vectWeight, &residualUser);
        }

        buf << " ";
        buf.width(4);
        buf.precision(8);
        buf.setf(ios::scientific, ios::floatfield);
        buf << j+1 << ". " << lambda[j] << " ";
        if (normWeight)
            buf << residualUser << " " << residualUser/lambda[j] << " ";
        buf << residualL2 << " " << residualL2/lambda[j] << " ";
        buf << endl;

        maxL2Norm = (residualL2/lambda[j] > maxL2Norm) ? residualL2/lambda[j] : maxL2Norm;
        minL2Norm = (residualL2/lambda[j] < minL2Norm) ? residualL2/lambda[j] : minL2Norm;
        if (normWeight) {
            maxUserNorm = (residualUser/lambda[j] > maxUserNorm) ? residualUser/lambda[j]
                : maxUserNorm;
            minUserNorm = (residualUser/lambda[j] < minUserNorm) ? residualUser/lambda[j]
                : minUserNorm;
        }

    } // for j=0; j<qc; ++j)

    if (normWeight) {
        buf << " >> Minimum scaled user-norm of residuals = " << minUserNorm << endl;
        buf << " >> Maximum scaled user-norm of residuals = " << maxUserNorm << endl;
    }
    buf << " >> Minimum scaled L2 norm of residuals   = " << minL2Norm << endl;
    buf << " >> Maximum scaled L2 norm of residuals   = " << maxL2Norm;

    if (work)
        delete[] work;

    rInfo(buf.str().c_str());
}


int CheckingTools::errorLambda(double *continuous, double *discrete, int numDiscrete,
                               double *lambda, int nev) const {
    ostringstream buf;
    int myPid = MyComm.MyPID();
    int nMax = 0;
    int i, j;

    // Allocate working arrays

    int *used = new (nothrow) int[numDiscrete + nev];
    if (used == 0) {
        return nMax;
    }

    int *bestMatch = used + numDiscrete;

    // Find the best match for the eigenvalues
    double eps = 0.0;
    Epetra_LAPACK call;
    call.LAMCH('E', eps);

    double gap = Epetra_MaxDouble;
    for (i=0; i<numDiscrete; ++i) {
        used[i] = -1;
        for (int j = i; j < numDiscrete; ++j)
            if (discrete[j] > (1.0 + 10.0*eps)*discrete[i]) {
                double tmp = (discrete[j] - discrete[i])/discrete[i];
                gap = (tmp < gap) ? tmp : gap;
                break;
            }
    }

    for (i=0; i<nev; ++i) 
        bestMatch[i] = -1;

    for (i=0; i<nev; ++i) {
        if (lambda[i] < continuous[0])
            continue;
        bestMatch[i] = (i == 0) ? 0 : bestMatch[i-1] + 1;
        int jStart = bestMatch[i];
        for (j = jStart; j < numDiscrete; ++j) {
            double diff = fabs(lambda[i]-discrete[j]);
            if (diff < 0.5*gap*lambda[i]) {
                bestMatch[i] = j;
                break;
            }
        }
        used[bestMatch[i]] = i;
    }

    // Print the results for the eigenvalues
    if (myPid == 0) {
        buf << endl;
        buf << " --- Relative errors on eigenvalues ---\n";
        buf << endl;
        buf << "       Exact Cont.    Exact Disc.     Computed      ";
        buf << " Alg. Err.  Mesh Err.\n";
    }

    int iCount = 0;
    for (i=0; i<nev; ++i) {
        if (bestMatch[i] == -1) {
            if (myPid == 0) {
                buf << "      ************** ************** ";
                buf.precision(8);
                buf.setf(ios::scientific, ios::floatfield);
                buf << lambda[i];
                buf << "  *********  *********" << endl;
            }
            iCount += 1;
        }
    }

    double lastDiscrete = 0.0;
    for (i=0; i<numDiscrete; ++i) {
        if ((iCount == nev) && (discrete[i] > lastDiscrete))
            break;
        if (used[i] < 0) {
            nMax += 1;
            lastDiscrete = discrete[i];
            if (myPid == 0) {
                buf << " ";
                buf.width(4);
                buf << i+1 << ". ";
                buf.precision(8);
                buf.setf(ios::scientific, ios::floatfield);
                buf << continuous[i] << " " << discrete[i] << " ";
                buf << "**************  *********  ";
                buf.precision(3);
                buf << fabs(continuous[i]-discrete[i])/continuous[i] << endl;
            }
        }
        else {
            nMax += 1;
            lastDiscrete = discrete[i];
            if (myPid == 0) {
                buf << " ";
                buf.width(4);
                buf << i+1 << ". ";
                buf.precision(8);
                buf.setf(ios::scientific, ios::floatfield);
                buf << continuous[i] << " " << discrete[i] << " " << lambda[used[i]] << "  ";
                buf.precision(3);
                buf << fabs(lambda[used[i]]-discrete[i])/discrete[i] << "  ";
                buf << fabs(continuous[i]-discrete[i])/continuous[i] << endl;
            }
            iCount += 1;
        }
    }

    delete[] used;

    if (myPid == 0)
        rDebug(buf.str().c_str());

    return nMax;
}


int CheckingTools::inputArguments(const int &numEigen, const Epetra_Operator *K, 
                                  const Epetra_Operator *M, const Epetra_Operator *P,
                                  const Epetra_MultiVector &Q, const int &minSize) const {

    // Routine to check some arguments
    // 
    // info = -  1 >> The stiffness matrix K has not been specified.
    // info = -  2 >> The maps for the matrix K and the matrix M differ.
    // info = -  3 >> The maps for the matrix K and the preconditioner P differ.
    // info = -  4 >> The maps for the vectors and the matrix K differ.
    // info = -  5 >> Q is too small for the number of eigenvalues requested.
    // info = -  6 >> Q is too small for the computation parameters.
    //

    int myPid = MyComm.MyPID();

    if (K == 0) {
        if (myPid == 0)
            cerr << "\n !!! The matrix K to analyze has not been specified !!! \n\n";
        return -1;
    }

    if (M) {
        int mGlobal = (M->OperatorDomainMap()).NumGlobalElements();
        int mLocal = (M->OperatorDomainMap()).NumMyElements();
        int kGlobal = (K->OperatorDomainMap()).NumGlobalElements();
        int kLocal = (K->OperatorDomainMap()).NumMyElements();
        if ((mGlobal != kGlobal) || (mLocal != kLocal)) {
            if (myPid == 0) {
                cerr << endl;
                cerr << " !!! The maps for the matrice K and the mass M are different !!!\n";
                cerr << endl;
            }
            return -2;
        }
    }

    if (P) {
        int pGlobal = (P->OperatorDomainMap()).NumGlobalElements();
        int pLocal = (P->OperatorDomainMap()).NumMyElements();
        int kGlobal = (K->OperatorDomainMap()).NumGlobalElements();
        int kLocal = (K->OperatorDomainMap()).NumMyElements();
        if ((pGlobal != kGlobal) || (pLocal != kLocal)) {
            if (myPid == 0) {
                cerr << endl;
                cerr << " !!! The maps for the matrice K and the preconditioner P are different !!!\n";
                cerr << endl;
            }
            return -3;
        }
    }

    if ((Q.MyLength() != (K->OperatorDomainMap()).NumMyElements()) || 
        (Q.GlobalLength() != (K->OperatorDomainMap()).NumGlobalElements())) {
        if (myPid == 0) {
            cerr << "\n !!! The maps for the vectors and the matrix are different !!! \n\n";
        }
        return -4;
    }

    if (Q.NumVectors() < numEigen) {
        if (myPid == 0) {
            cerr << endl;
            cerr << " !!! The number of eigenvalues is too large for the space allocated !!! \n\n";
            cerr << " The recommended size for " << numEigen << " eigenvalues is ";
            cerr << minSize << endl;
            cerr << endl;
        }
        return -5;
    }

    if (Q.NumVectors() < minSize) {
        if (myPid == 0) {
            cerr << endl;
            cerr << " !!! The space allocated is too small for the number of eigenvalues";
            cerr << " and the size of blocks specified !!! \n";
            cerr << " The recommended size for " << numEigen << " eigenvalues is ";
            cerr << minSize << endl;
            cerr << endl;
        }
        return -6;
    }

    return 0;

}

---