src/eigsolv/Projectors.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)
//
//**************************************************************************
/**************************************************************************
 *                                                                         *
 *               Swiss Federal Institute of Technology (ETH)               *
 *                       CH-8092 Zuerich, Switzerland                      *
 *                                                                         *
 *                       (C) 1999 All Rights Reserved                      *
 *                                                                         *
 *                                NOTICE                                   *
 *                                                                         *
 *  Permission to use, copy, modify, and distribute this software and      *
 *  its documentation for any purpose and without fee is hereby granted    *
 *  provided that the above copyright notice appear in all copies and      *
 *  that both the copyright notice and this permission notice appear in    *
 *  supporting documentation.                                              *
 *                                                                         *
 *  Neither the Swiss Federal Institute of Technology nor the author make  *
 *  any representations about the suitability of this software for any     *
 *  purpose.  This software is provided ``as is'' without express or       *
 *  implied warranty.                                                      *
 *                                                                         *
 ***************************************************************************/

#include "Projectors.h"


const int Projectors::OP_UNSYM1 = 1;
const int Projectors::OP_UNSYM2 = 2;
const int Projectors::OP_SYM = 3;


Projectors::Projectors(const Epetra_Comm &_Comm,
                       const Epetra_Operator *AA, const Epetra_Operator *BB, const Epetra_Operator *PP,
                       int eqType, int k, int kmax, int *Hpiv, int *Gpiv,
                       double theta, double *Hlu, double *Glu,
                       Epetra_MultiVector *Q, Epetra_MultiVector *Qb, Epetra_MultiVector *Y,
                       double *w1, double *w2)
    : callBLAS(),
      callLAPACK(),
      orthoTool(),
      MyComm(_Comm),
      MyWatch(_Comm),
      A(AA),
      B(BB),
      APrec(PP),
      Pro_k(k),
      Pro_kmax(kmax),
      Pro_EquationType(eqType),
      Pro_Hpiv(Hpiv),
      Pro_Gpiv(Gpiv),
      Pro_theta(theta),
      Pro_Hlu(Hlu),
      Pro_Glu(Glu),
      Pro_Q(Q),
      Pro_Qb(Qb),
      Pro_Y(Y),
      work1(w1),
      work2(w2),
      timeApply(0.0),
      timeApplyInverse(0.0)
{

}


const Epetra_Map& Projectors::OperatorDomainMap() const {

    // This routine assumes that at least one Epetra_Operator is not empty

    if (A)
        return A->OperatorDomainMap();
    if (B)
        return B->OperatorDomainMap();

    return APrec->OperatorDomainMap();

}


const Epetra_Map& Projectors::OperatorRangeMap() const {

    // This routine assumes that at least one Epetra_Operator is not empty

    if (A)
        return A->OperatorRangeMap();
    if (B)
        return B->OperatorRangeMap();

    return APrec->OperatorRangeMap();

}


// Operator for the lin.eq.-solver used during the correction equation
int Projectors::Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    timeApply -= MyWatch.WallTime();
    if (A) {
        switch(Pro_EquationType){
        case OP_UNSYM1:
            LinOp_UNSYM1(X, Y);
            break;
        case OP_UNSYM2:
            LinOp_UNSYM2(X, Y);
            break;
        case OP_SYM:
            LinOp_SYM(X, Y);
            break;
        }
    }
    else {
        memcpy(Y.Values(), X.Values(), X.NumVectors()*X.MyLength()*sizeof(double));
    }  // if (A)
    timeApply += MyWatch.WallTime();

    return 0;

}


/* Preconditioner for the lin.eq.-solver used during the correction equation */
int Projectors::ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    int nV = X.NumVectors();

    timeApplyInverse -= MyWatch.WallTime();
    if (APrec) {
        switch(Pro_EquationType){
        case OP_UNSYM1:
            Precon_UNSYM1(X, Y);
            break;
        case OP_UNSYM2:
            Precon_UNSYM2(X, Y);
            break;
        case OP_SYM:
            Precon_SYM(X, Y);
            break;
        }
    }
    else {
        memcpy(Y.Values(), X.Values(), nV*X.MyLength()*sizeof(double));
    } // if (APrec)
    timeApplyInverse += MyWatch.WallTime();

    return 0;

}


/* Preconditions and projects the right hand side of the correction equation */
void Projectors::MakeRHS(Epetra_MultiVector &R){

    int iV;
    int nV = R.NumVectors();

    switch(Pro_EquationType){
    case OP_UNSYM1:
        for (iV = 0; iV < nV; ++iV) {
            Epetra_Vector Rv(View, R, iV);
            Right_UNSYM1(Rv);
        }
        break;
    case OP_UNSYM2:
        for (iV = 0; iV < nV; ++iV) {
            Epetra_Vector Rv(View, R, iV);
            Right_UNSYM2(Rv);
        }
        break;
    case OP_SYM:
        for (iV = 0; iV < nV; ++iV) {
            Epetra_Vector Rv(View, R, iV);
            Right_SYM(Rv);
        }
        break;
    }

}

/****************************************************************************
 *                                                                          *
 * Linear Operators                                                         *
 *                                                                          *
 ****************************************************************************/

void Projectors::LinOp_UNSYM1(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    int branch_code;

    branch_code = 10*(B != 0) + (APrec != 0);

    Epetra_Vector w2(View, X.Map(), work2);

    int iX;
    int nX = X.NumVectors();

    if (branch_code == 11) {
        /* y := (I - Y'*inv(H)*Qb')*inv(K)*(A-theta*B)*x */
        for (iX = 0; iX < nX; ++iX) {
            Epetra_Vector Xi(View, X, iX);
            Epetra_Vector Yi(View, Y, iX);
            A_theta_B(Xi, w2, work1);
            APrec->ApplyInverse(w2, Yi);
        }
        Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, Y, Y);
    }
    else if (branch_code == 10) {
        /* y := (I - Q*Qb')*(A-theta*B)*x */
        for (iX = 0; iX < nX; ++iX) {
            Epetra_Vector Xi(View, X, iX);
            Epetra_Vector Yi(View, Y, iX);
            A_theta_B(Xi, Yi, work1);
        }
        Project1(Pro_Q, Pro_Qb, Y, Y);
    }
    else if (branch_code == 01) {
        /* y := (I - Y'*inv(H)*Q')*inv(K)*(A-theta*I)*x */
        for (iX = 0; iX < nX; ++iX) {
            Epetra_Vector Xi(View, X, iX);
            Epetra_Vector Yi(View, Y, iX);
            A_theta_I(Xi, w2);
            APrec->ApplyInverse(w2, Yi);
        }
        Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, Y, Y);
    }
    else { /* (Pro_Bisempty == 1) && (Pro_Kisempty == 1) */
        /* y := (I - Q*Q')*(A-theta*I)*x */
        A_theta_I(X, Y);
        Project1(Pro_Q, Pro_Q, Y, Y);
    }

}


void Projectors::LinOp_UNSYM2(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    int branch_code;

    branch_code = 10*(B != 0) + (APrec != 0);

    int iX;
    int nX = X.NumVectors();

    if (branch_code == 00) {
        /* y = (I - Q*Q')*(A - theta*I) * x */
        A_theta_I(X, Y);
        Project1(Pro_Q, Pro_Q, Y, Y);
    }
    else if (branch_code == 10) {
        Epetra_Vector w2(View, X.Map(), work2);
        /* y = (I - Qb*Q')*(A-theta*B)*(I - Q*Qb') * x */   
        for (iX = 0; iX < nX; ++iX) {
            Epetra_Vector Xi(View, X, iX);
            Epetra_Vector Yi(View, Y, iX);
            Project1(Pro_Q, Pro_Qb, Xi, w2);
            A_theta_B(w2, Yi, work1);
        }
        Project1(Pro_Qb, Pro_Q, Y, Y);
    }
    else if (branch_code == 01) {
        Epetra_Vector w2(View, X.Map(), work2);
        /* y = (I - Y'*inv(H)*Qb')*inv(K)*(A-theta*B) * x */   
        for (iX = 0; iX < nX; ++iX) {
            Epetra_Vector Xi(View, X, iX);
            Epetra_Vector Yi(View, Y, iX);
            A_theta_I(Xi, w2);
            APrec->ApplyInverse(w2, Yi);
        }
        Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, Y, Y);
    }
    else /* branch_code == 11 */ {
        Epetra_Vector w1(View, X.Map(), work1);
        Epetra_Vector w2(View, X.Map(), work2);
        /* y = (I - Q*inv(G)*Q')* ...
           (I - Y'*inv(H)*Qb')*inv(K)*(A-theta*B)*(I - Q*Qb')*x */   
        for (iX = 0; iX < nX; ++iX) {
            Epetra_Vector Xi(View, X, iX);
            Epetra_Vector Yi(View, Y, iX);
            Project1(Pro_Q, Pro_Qb, Xi, w2);
            A_theta_B(w2, w1, Yi.Values());
            APrec->ApplyInverse(w1, Yi);
        }
        Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, Y, Y);
        Project2(Pro_Q, Pro_k, Pro_Gpiv, Pro_Glu, Pro_Q, Y, Y);
    }

}


void Projectors::LinOp_SYM(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    if (B == 0) {
        /* y := (I - Q*Q')*(A - theta*I)*x */
        A_theta_I(X, Y);
        Project1(Pro_Q, Pro_Q, Y, Y);
    } else {
        /* y := (I - Qb*Q')*(A - theta*I)*(I - Q*Qb')*x */
        Project1(Pro_Q, Pro_Qb, X, Y);
        int iV;
        int nV = X.NumVectors();
        int rx = X.MyLength();
        Epetra_Vector w2(View, X.Map(), work2);
        for (iV = 0; iV < nV; ++iV) {
            Epetra_Vector Yv(View, Y, iV);
            memcpy(work2, Yv.Values(), rx*sizeof(double));
            A_theta_B(w2, Yv, work1);
        }
        Project1(Pro_Qb, Pro_Q, Y, Y);
    } // if (B == 0)

}


/****************************************************************************
 *                                                                          *
 * Preconditioners                                                          *
 *                                                                          *
 ****************************************************************************/

 
void Projectors::Precon_UNSYM1(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    /* the Preconditioner is contained in the operator */
    memcpy(Y.Values(), X.Values(), X.MyLength()*X.NumVectors()*sizeof(double));

}


void Projectors::Precon_UNSYM2(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    /* the Preconditioner is contained in the operator */
    memcpy(Y.Values(), X.Values(), X.MyLength()*X.NumVectors()*sizeof(double));

}


void Projectors::Precon_SYM(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    if (APrec == 0) {
        memcpy(Y.Values(), X.Values(), X.MyLength()*X.NumVectors()*sizeof(double));
    } else {
        /* y := inv(K) * x */
        APrec->ApplyInverse(X, Y);
        if (B == 0) {
            /* y = (I-Q*G\Q') * (I-Y*H\Q') * y */
            Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, Y, Y);
        } else /* B not empty */ {
            /* y = (I-Q*G\Q') * (I-Y*H\Qb') * y */
            Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, Y, Y);
            Project2(Pro_Q, Pro_k, Pro_Gpiv, Pro_Glu, Pro_Q, Y, Y);
        }
    }

}


/****************************************************************************
 *                                                                          *
 * Right hand sides                                                         *
 *                                                                          *
 ****************************************************************************/


void Projectors::Right_UNSYM1(Epetra_Vector &r){

    if((B) && (APrec)){
        Epetra_Vector w1(View, r.Map(), work1);
        /* r := (I - Y'*inv(H)*Qb')*inv(K)*r */
        APrec->ApplyInverse(r, w1);
        Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, w1, r);
    }
    else if((B) && (APrec == 0))
        /* r :=  (I - Q*Qb')*r */
        Project1(Pro_Q, Pro_Qb, r, r);
    else if((B == 0) && (APrec)){
        Epetra_Vector w1(View, r.Map(), work1);
        /* r := (I - Y'*inv(H)*Q')*inv(K)*r */
        APrec->ApplyInverse(r, w1);
        Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, w1, r);
    }
    else {
        /* r :=  (I - Q*Q')*r */
        orthoTool.ModifiedGS(&r, Pro_Q);
    }

}


void Projectors::Right_UNSYM2(Epetra_Vector &r) {

    if (B == 0) {
        if (APrec == 0)
            /* r = (I-Q*Q') * r */
            orthoTool.ModifiedGS(&r, Pro_Q);
        else {
            Epetra_Vector w1(View, r.Map(), work1);
            /* r = (I-Y*H\Q') * inv(K) * r */
            APrec->ApplyInverse(r, w1);
            Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Q, w1, r);
        }
    }
    else {
        if (APrec == 0)
            /* r = (I - Qb*Q')*r */
            Project1(Pro_Qb, Pro_Q, r, r);
        else {
            /* r = (I-Q*G\Q') * (I-Y*H\Qb') * inv(K) * r */
            Epetra_Vector w1(View, r.Map(), work1);
            APrec->ApplyInverse(r, w1);
            Project2(Pro_Y, Pro_k, Pro_Hpiv, Pro_Hlu, Pro_Qb, w1, r);
            Project2(Pro_Q, Pro_k, Pro_Gpiv, Pro_Glu, Pro_Q, r, r);
        }
    }
}


void Projectors::Right_SYM(Epetra_Vector &r) {

    if (B == 0) {
        /* r :=  (I - Q*Q')*r */
        orthoTool.ModifiedGS(&r, Pro_Q);
    } else { /* B not empty */
        /* r = r - Qb*(Q'*r); */
        Project1(Pro_Qb, Pro_Q, r, r); 
    }

}


/****************************************************************************
 *                                                                          *
 * Multiplication Operators                                                 *
 *                                                                          *
 ****************************************************************************/


void Projectors::A_theta_B(const Epetra_MultiVector &X, Epetra_MultiVector &Y, double *wTmp) 
    const {

    A->Apply(X, Y);

    int iV;
    int nV = X.NumVectors();
    int rX = X.MyLength(); 

    for (iV = 0; iV < nV; ++iV) {
        Epetra_MultiVector Xv(View, X, iV, 1);
        Epetra_MultiVector w(View, X.Map(), wTmp, rX, 1);
        B->Apply(Xv, w);
        callBLAS.AXPY(rX, -Pro_theta, wTmp, Y.Values() + iV*rX);
    }

}


void Projectors::A_theta_I(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    A->Apply(X, Y);
    callBLAS.AXPY(X.MyLength()*X.NumVectors(), -Pro_theta, X.Values(), Y.Values());

}


/****************************************************************************
 *                                                                          *
 * Helper routines                                                          *
 *                                                                          *
 ****************************************************************************/


/* PROJECT1 - compute y := (I - A*B')*x
 *
 * sizes of matrices and vectors:
 *   A and B     --  r-by-c matrices, ldim = r
 *   x, y and w  --  r vectors
 *
 * (x, y) and also (A, B) may point to the same memory locations
 */


void Projectors::Project1(Epetra_MultiVector *AA, Epetra_MultiVector *BB, 
                          const Epetra_MultiVector &X, Epetra_MultiVector &Y) const {

    int rA = AA->MyLength();
    int cA = AA->NumVectors();

    int rX = X.MyLength();
    int cX = X.NumVectors();

    int size = cX*cA;
    double *wTmp = new double[2*size];

    /* if x and y do not refer to the same vector, we copy x to y */
    if (X.Values() != Y.Values())
        memcpy(Y.Values(), X.Values(), rX*cX*sizeof(double));

    /* W = B'*Y */
    callBLAS.GEMM('T', 'N', cA, cX, rA, 1.0, BB->Values(), rA, Y.Values(), rX,
                  0.0, wTmp + size, cA);
    MyComm.SumAll(wTmp + size, wTmp, size);

    /* Y = Y - A*W */
    callBLAS.GEMM('N', 'N', rA, cX, cA, -1.0, AA->Values(), rA, wTmp, cA,
                  1.0, Y.Values(), rX);

    delete[] wTmp;

}


/* PROJECT2 - compute y := (I - A*inv(B)*C')*x
 *
 * (x, y) may point to the same memory location
 */


void Projectors::Project2(Epetra_MultiVector *AA, int rowB, int *piv, double *Blu,  
                          Epetra_MultiVector *CC, const Epetra_MultiVector &X, 
                          Epetra_MultiVector &Y) const {

    int info;

    int rowA = AA->MyLength();
    int colA = AA->NumVectors();

    int rowC = CC->MyLength();
    int colC = CC->NumVectors();

    int rowX = X.MyLength();
    int colX = X.NumVectors();

    int size = colC*colX;
    double *wTmp = new double[2*size];

    /* if x and y do not refer to the same vector, we copy x to y */
    if (X.Values() != Y.Values())
        memcpy(Y.Values(), X.Values(), rowX*colX*sizeof(double));

    /* W = C'*X */
    callBLAS.GEMM('T', 'N', colC, colX, rowC, 1.0, CC->Values(), rowC, X.Values(), rowX,
                  0.0,  wTmp + size, colC);
    MyComm.SumAll(wTmp + size, wTmp, size);

    /* W = inv(B)*W */
    callLAPACK.GETRS('N', rowB, colX, Blu, Pro_kmax, piv, wTmp, colC, &info);
    assert(info == 0);

    /* Y = Y - A*W = Y - A*inv(B)*C'*X */
    callBLAS.GEMM('N', 'N', rowA, colX, colA, -1.0, AA->Values(), rowA, wTmp, colC,
                  1.0, Y.Values(), rowX);

    delete[] wTmp;

}

---