src/eigsolv/Projectors.h

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.                                                      *
 *                                                                         *
 ***************************************************************************/

#ifndef JDSYM_PROJECTORS_H
#define JDSYM_PROJECTORS_H

#include "Epetra_ConfigDefs.h"

#include "Epetra_BLAS.h"
#include "Epetra_Comm.h"
#include "Epetra_LAPACK.h"
#include "Epetra_MultiVector.h"
#include "Epetra_Operator.h"
#include "Epetra_Time.h"
#include "Epetra_Vector.h"

#include "OrthoPack.h"

class Projectors : public virtual Epetra_Operator {

private:

    const Epetra_BLAS callBLAS;
    const Epetra_LAPACK callLAPACK;
    OrthoPack orthoTool;

    const Epetra_Comm &MyComm;
    const Epetra_Time MyWatch;

    const Epetra_Operator *A;
    const Epetra_Operator *B;
    const Epetra_Operator *APrec;

    int Pro_k;
    int Pro_kmax;

    int Pro_EquationType;

    int *Pro_Hpiv;
    int *Pro_Gpiv;

    double Pro_theta;

    double *Pro_Hlu;
    double *Pro_Glu;

    Epetra_MultiVector *Pro_Q;
    Epetra_MultiVector *Pro_Qb;
    Epetra_MultiVector *Pro_Y;

    double *work1;
    double *work2;

    mutable double timeApply;
    mutable double timeApplyInverse;

    static const int OP_UNSYM1;
    static const int OP_UNSYM2;
    static const int OP_SYM;

    // Don't define these functions !
    Projectors(const Projectors &ref);
    Projectors& operator=(const Projectors &ref);

protected:

    void LinOp_UNSYM1(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
    void LinOp_UNSYM2(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
    void LinOp_SYM(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;

    void Precon_UNSYM1(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
    void Precon_UNSYM2(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
    void Precon_SYM(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;

    void Right_UNSYM1(Epetra_Vector &r);
    void Right_UNSYM2(Epetra_Vector &r);
    void Right_SYM(Epetra_Vector &r);

    void A_theta_B(const Epetra_MultiVector &X, Epetra_MultiVector &Y, double *wTmp) const;
    void A_theta_I(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;

    void Project1(Epetra_MultiVector *AA, Epetra_MultiVector *BB, const Epetra_MultiVector &X,
                  Epetra_MultiVector &Y) const;

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

public:

    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);

    ~Projectors() { };

    void MakeRHS(Epetra_MultiVector &R);

    void resetPro_k(int k) { Pro_k = k; };
    void resetPro_theta(double t) { Pro_theta = t; };
    void resetPro_Q(Epetra_MultiVector *dd) { Pro_Q = dd; };
    void resetPro_Qb(Epetra_MultiVector *dd) { Pro_Qb = dd; };
    void resetPro_Y(Epetra_MultiVector *dd) { Pro_Y = dd; };

    const double getTimeApply() const { return timeApply; };
    const double getTimeApplyInverse() const { return timeApplyInverse; };

    // Functions inherited from Epetra_Operator

    const char *Label() const { return "Linear operators for JDSYM"; };

    bool UseTranspose() const { return false; };
    int SetUseTranspose(bool UseTranspose) { return 0; };

    bool HasNormInf() const { return (false); };
    double NormInf() const  { return (-1.0); };

    int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
    int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;

    const Epetra_Comm& Comm() const             { return MyComm; };

    const Epetra_Map& OperatorDomainMap() const;
    const Epetra_Map& OperatorRangeMap() const;

};

#endif

---