src/eigsolv/LinSolvers.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 "LinSolvers.h"


const int LinSolvers::LINSOLV_CGS = 1;
const int LinSolvers::LINSOLV_MINRES = 2;
const int LinSolvers::LINSOLV_SYMMLQ = 3;
const int LinSolvers::LINSOLV_QMR = 4;


LinSolvers::LinSolvers(const Epetra_Operator *AA, const Epetra_Operator *PP, double *wTmp,
                       int _type)
    : A(AA), Prec(PP), work_s(wTmp), solverType(_type)
{

}


int LinSolvers::solve(const Epetra_Vector &b, Epetra_Vector &x, double _tol, int _itMax) {

    int info;

    switch (solverType) {
    default:
    case LINSOLV_QMR:
        info = QMRS(b, x, _tol, _itMax);
        break;
    }

    return info;

}


int LinSolvers::QMRS(const Epetra_Vector &b, Epetra_Vector &x, double _tol, int _itMax) {

    int maxIter = _itMax;
    double tol = _tol;

    // This routine uses QMRS to solve a linear system

    int nLocal = b.MyLength();
  
    Epetra_Vector wrk1(View, b.Map(), work_s);
    Epetra_Vector p(View, b.Map(), work_s + nLocal);
    Epetra_Vector d(View, b.Map(), work_s + 2*nLocal);
    Epetra_Vector v1(View, b.Map(), work_s + 3*nLocal);
    Epetra_Vector t(View, b.Map(), work_s + 4*nLocal);
    Epetra_Vector g(View, b.Map(), work_s + 5*nLocal);

    double beta;
    double c0, cc, c1;
    double delta;
    double eps0, eta0; 
    double res_init, rho0 = 0.0, rho1 = 0.0, rho1inv;
    double tau, theta0, theta;
    double xi1;

    int iter;

    memcpy(v1.Values(), b.Values(), nLocal*sizeof(double));
    v1.Norm2(&rho0);
    tau = rho0;

    v1.Scale(1.0/rho0);
    p.PutScalar(0.0);
    g.PutScalar(0.0);
    d.PutScalar(0.0);
    x.PutScalar(0.0);

    c0 = 1.0;
    eps0 = 1.0;
    xi1 = 1.0;
    theta0 = 0.0;
    eta0 = -1.0;
    res_init = rho0;

    for (iter = 0; iter < maxIter; ++iter) {

        if (eps0 == 0.0)
            return -2;

        if (Prec)
            Prec->ApplyInverse(v1, wrk1);
        else
            memcpy(wrk1.Values(), v1.Values(), nLocal*sizeof(double));

        wrk1.Dot(v1, &delta);
        if (delta == 0.0)
            return -2;

        cc = xi1*(delta/eps0);

        p.Update(1.0, v1, -cc);
        g.Update(1.0, wrk1, -cc);

        A->Apply(g, t);

        g.Dot(t, &eps0);
        beta = eps0/delta;
        v1.Update(1.0, t, -beta);

        v1.Norm2(&rho1);
        xi1 = rho1;

        theta = c0*fabs(beta);
        if (theta == 0.0)
            return -2;
        theta = rho1/theta;

        c1 = 1.0/sqrt(1.0 + theta*theta);
        eta0 = -eta0*rho0*c1*c1/(beta*c0*c0);
        tau = tau*theta*c1;

        if (rho1 == 0.0)
            return -2;

        cc = theta0*c1;
        cc = cc*cc;
        rho1inv = 1.0/rho1;

        d.Update(eta0, p, cc);
        x.Update(1.0, d, 1.0);
        v1.Scale(rho1inv);

        if (xi1 == 0)
            return -2;

        rho0 = rho1;
        c0 = c1;
        theta0 = theta;

        if (tau <= res_init*tol)
            break;

    }

    memcpy(wrk1.Values(), x.Values(), nLocal*sizeof(double));
    Prec->ApplyInverse(wrk1, x);

    if (iter < maxIter)
        iter += 1;

    if (tau > res_init*tol)
        return -1;

    return iter;

}

---