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

// Max. Number of Iterations in CGS
const int ORTHOPACK_MAXCGSIT = 5;       

//
//
//  Iterated Classical B-orthogonal GramSchmidt 
//
//  B-Orthogonalizes v with respect to A. 
//  If B = 0, then use identity
//
//


OrthoPack::OrthoPack() :
    callBLAS()
{

}


void OrthoPack::IteratedClassicalGS(Epetra_Vector *v, double *vnrm, const Epetra_MultiVector *A,
                                    double *work1, const Epetra_Operator *B) const {

    // Note that A can be 0
    // We assume that the number of columns of A is small

    const double alpha = 0.5;
    double alpha2 = alpha*alpha;

    double vnrm_old = 0.0;
    int i = 0;

    int n = (A) ? A->MyLength() : v->MyLength();
    int m = (A) ? A->NumVectors() : 0;

    Epetra_Vector w1(View, v->Map(), work1);
    if (B)
        B->Apply(*v, w1);
    else 
        memcpy(work1, v->Values(), n*sizeof(double));

    w1.Dot(*v, &vnrm_old);
    (*vnrm) = vnrm_old;

    if (m > 0) {

        double *work2 = new double[2*m];

        for (i = 0; i < ORTHOPACK_MAXCGSIT; i ++) {

            callBLAS.GEMV('T', n, m, 1.0, A->Values(), n, work1, 0.0, work2 + m);
            (A->Comm()).SumAll(work2 + m, work2, m);
            callBLAS.GEMV('N', n, m, -1.0, A->Values(), n, work2, 1.0, v->Values());

            if (B)
                B->Apply(*v, w1);
            else 
                memcpy(work1, v->Values(), n*sizeof(double));

            w1.Dot((*v), vnrm);

            if ((*vnrm) >= alpha2*vnrm_old)
                break;

            vnrm_old = (*vnrm);

        } // for (i = 0; i < ORTHOPACK_MAXCGSIT; i ++)

        delete[] work2;

    } // if (m > 0)

    (*vnrm) = sqrt((*vnrm));

    if (i >= ORTHOPACK_MAXCGSIT) {
        cout << "warning: loss of orthogonality. ";
        cout << "Iterated CGS not converged after " << ORTHOPACK_MAXCGSIT << " steps\n";
    }

}

/*
 *  ModifiedGramSchmidt 
 *
 *  Orthogonalizes v with respect to A
 */

void OrthoPack::ModifiedGS(Epetra_Vector *v, const Epetra_MultiVector *A) const {

    int i;
    double s = 0.0;

    int n = (A) ? A->MyLength() : v->MyLength();
    int m = (A) ? A->NumVectors() : 0;

    for (i = 0; i < m; ++i) {
        (*A)(i)->Dot(*v, &s);
        callBLAS.AXPY(n, -s, A->Values() + i*n, v->Values());
    }

}

/*
 *  Modified B-orthogonal GramSchmidt
 *
 *  B-Orthogonalizes v with respect to A [ BA is the Product of B*A !!! ]
 */

void OrthoPack::ModifiedGSB(Epetra_Vector *v, const Epetra_MultiVector *A, 
                            const Epetra_MultiVector *BA) const {

    int i;
    double s = 0.0;

    int n = (A) ? A->MyLength() : v->MyLength();
    int m = (A) ? A->NumVectors() : 0;

    for (i = 0; i < m; ++i) {
        (*BA)(i)->Dot(*v, &s);
        callBLAS.AXPY(n, -s, A->Values() + i*n, v->Values());
    }

}

---