src/eigsolv/OrthoPack.cpp

Go to the documentation of this file.
00001 //**************************************************************************
00002 //
00003 //                                 NOTICE
00004 //
00005 // This software is a result of the research described in the report
00006 //
00007 // " A comparison of algorithms for modal analysis in the absence 
00008 //   of a sparse direct method", P. Arbenz, R. Lehoucq, and U. Hetmaniuk,
00009 //  Sandia National Laboratories, Technical report SAND2003-1028J.
00010 //
00011 // It is based on the Epetra, AztecOO, and ML packages defined in the Trilinos
00012 // framework ( http://software.sandia.gov/trilinos/ ).
00013 //
00014 // The distribution of this software follows also the rules defined in Trilinos.
00015 // This notice shall be marked on any reproduction of this software, in whole or
00016 // in part.
00017 //
00018 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00019 // license for use of this work by or on behalf of the U.S. Government.
00020 //
00021 // This program is distributed in the hope that it will be useful, but
00022 // WITHOUT ANY WARRANTY; without even the implied warranty of
00023 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00024 //
00025 // Code Authors: U. Hetmaniuk (ulhetma@sandia.gov), R. Lehoucq (rblehou@sandia.gov)
00026 //
00027 //**************************************************************************
00028 /**************************************************************************
00029  *                                                                         *
00030  *               Swiss Federal Institute of Technology (ETH)               *
00031  *                       CH-8092 Zuerich, Switzerland                      *
00032  *                                                                         *
00033  *                       (C) 1999 All Rights Reserved                      *
00034  *                                                                         *
00035  *                                NOTICE                                   *
00036  *                                                                         *
00037  *  Permission to use, copy, modify, and distribute this software and      *
00038  *  its documentation for any purpose and without fee is hereby granted    *
00039  *  provided that the above copyright notice appear in all copies and      *
00040  *  that both the copyright notice and this permission notice appear in    *
00041  *  supporting documentation.                                              *
00042  *                                                                         *
00043  *  Neither the Swiss Federal Institute of Technology nor the author make  *
00044  *  any representations about the suitability of this software for any     *
00045  *  purpose.  This software is provided ``as is'' without express or       *
00046  *  implied warranty.                                                      *
00047  *                                                                         *
00048  ***************************************************************************/
00049 
00050 #include "OrthoPack.h"
00051 
00052 // Max. Number of Iterations in CGS
00053 const int ORTHOPACK_MAXCGSIT = 5;       
00054 
00055 //
00056 //
00057 //  Iterated Classical B-orthogonal GramSchmidt 
00058 //
00059 //  B-Orthogonalizes v with respect to A. 
00060 //  If B = 0, then use identity
00061 //
00062 //
00063 
00064 
00065 OrthoPack::OrthoPack() :
00066     callBLAS()
00067 {
00068 
00069 }
00070 
00071 
00072 void OrthoPack::IteratedClassicalGS(Epetra_Vector *v, double *vnrm, const Epetra_MultiVector *A,
00073                                     double *work1, const Epetra_Operator *B) const {
00074 
00075     // Note that A can be 0
00076     // We assume that the number of columns of A is small
00077 
00078     const double alpha = 0.5;
00079     double alpha2 = alpha*alpha;
00080 
00081     double vnrm_old = 0.0;
00082     int i = 0;
00083 
00084     int n = (A) ? A->MyLength() : v->MyLength();
00085     int m = (A) ? A->NumVectors() : 0;
00086 
00087     Epetra_Vector w1(View, v->Map(), work1);
00088     if (B)
00089         B->Apply(*v, w1);
00090     else 
00091         memcpy(work1, v->Values(), n*sizeof(double));
00092 
00093     w1.Dot(*v, &vnrm_old);
00094     (*vnrm) = vnrm_old;
00095 
00096     if (m > 0) {
00097 
00098         double *work2 = new double[2*m];
00099 
00100         for (i = 0; i < ORTHOPACK_MAXCGSIT; i ++) {
00101 
00102             callBLAS.GEMV('T', n, m, 1.0, A->Values(), n, work1, 0.0, work2 + m);
00103             (A->Comm()).SumAll(work2 + m, work2, m);
00104             callBLAS.GEMV('N', n, m, -1.0, A->Values(), n, work2, 1.0, v->Values());
00105 
00106             if (B)
00107                 B->Apply(*v, w1);
00108             else 
00109                 memcpy(work1, v->Values(), n*sizeof(double));
00110 
00111             w1.Dot((*v), vnrm);
00112 
00113             if ((*vnrm) >= alpha2*vnrm_old)
00114                 break;
00115 
00116             vnrm_old = (*vnrm);
00117 
00118         } // for (i = 0; i < ORTHOPACK_MAXCGSIT; i ++)
00119 
00120         delete[] work2;
00121 
00122     } // if (m > 0)
00123 
00124     (*vnrm) = sqrt((*vnrm));
00125 
00126     if (i >= ORTHOPACK_MAXCGSIT) {
00127         cout << "warning: loss of orthogonality. ";
00128         cout << "Iterated CGS not converged after " << ORTHOPACK_MAXCGSIT << " steps\n";
00129     }
00130 
00131 }
00132 
00133 /*
00134  *  ModifiedGramSchmidt 
00135  *
00136  *  Orthogonalizes v with respect to A
00137  */
00138 
00139 void OrthoPack::ModifiedGS(Epetra_Vector *v, const Epetra_MultiVector *A) const {
00140 
00141     int i;
00142     double s = 0.0;
00143 
00144     int n = (A) ? A->MyLength() : v->MyLength();
00145     int m = (A) ? A->NumVectors() : 0;
00146 
00147     for (i = 0; i < m; ++i) {
00148         (*A)(i)->Dot(*v, &s);
00149         callBLAS.AXPY(n, -s, A->Values() + i*n, v->Values());
00150     }
00151 
00152 }
00153 
00154 /*
00155  *  Modified B-orthogonal GramSchmidt
00156  *
00157  *  B-Orthogonalizes v with respect to A [ BA is the Product of B*A !!! ]
00158  */
00159 
00160 void OrthoPack::ModifiedGSB(Epetra_Vector *v, const Epetra_MultiVector *A, 
00161                             const Epetra_MultiVector *BA) const {
00162 
00163     int i;
00164     double s = 0.0;
00165 
00166     int n = (A) ? A->MyLength() : v->MyLength();
00167     int m = (A) ? A->NumVectors() : 0;
00168 
00169     for (i = 0; i < m; ++i) {
00170         (*BA)(i)->Dot(*v, &s);
00171         callBLAS.AXPY(n, -s, A->Values() + i*n, v->Values());
00172     }
00173 
00174 }
00175 
00176 

Generated on Fri Oct 26 13:35:11 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7