example/InverseIteration/cxx_main.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <ctype.h>
00032 #include <assert.h>
00033 #include <string.h>
00034 #include <math.h>
00035 #include "Epetra_Comm.h"
00036 #include "Epetra_Map.h"
00037 #include "Epetra_Time.h"
00038 #include "Epetra_BlockMap.h"
00039 #include "Epetra_MultiVector.h"
00040 #include "Epetra_Vector.h"
00041 #include "Epetra_CrsMatrix.h"
00042 #ifdef EPETRA_MPI
00043 #include "mpi.h"
00044 #include "Epetra_MpiComm.h"
00045 #else
00046 #include "Epetra_SerialComm.h"
00047 #endif
00048 #include "Trilinos_Util.h"
00049 #include "Ifpack_IlukGraph.h"
00050 #include "Ifpack_CrsRiluk.h"
00051 
00052 // *************************************************************
00053 // Function prototypes
00054 // *************************************************************
00055 
00056 int invIteration(Epetra_CrsMatrix& A, double &lambda, bool verbose); 
00057 int applyInverseSetup(Epetra_CrsMatrix &A, Ifpack_CrsRiluk * & M);
00058 int applyInverse(Epetra_CrsMatrix &A, Epetra_Vector & x, Epetra_Vector & b, Ifpack_CrsRiluk * M, 
00059      bool verbose);
00060 int applyInverseDestroy(Ifpack_CrsRiluk * M);
00061 void BiCGSTAB(Epetra_CrsMatrix &A, 
00062         Epetra_Vector &x, 
00063         Epetra_Vector &b, 
00064         Ifpack_CrsRiluk *M, 
00065         int Maxiter, 
00066         double Tolerance, bool verbose);
00067 
00068 // *************************************************************
00069 // main program - This benchmark code reads a Harwell-Boeing data
00070 //                set and finds the minimal eigenvalue of the matrix
00071 //                using inverse iteration.
00072 // *************************************************************
00073 int main(int argc, char *argv[]) {
00074 
00075 #ifdef EPETRA_MPI
00076   MPI_Init(&argc,&argv);
00077   Epetra_MpiComm Comm (MPI_COMM_WORLD);
00078 #else
00079   Epetra_SerialComm Comm;
00080 #endif
00081 
00082   cout << Comm << endl;
00083 
00084   int MyPID = Comm.MyPID();
00085 
00086   bool verbose = false; 
00087   if (MyPID==0) verbose = true; // Print out detailed results (turn off for best performance)
00088 
00089   if(argc != 2) {
00090     if (verbose) cerr << "Usage: " << argv[0] << " HB_data_file" << endl;
00091     exit(1); // Error
00092   }
00093 
00094   // Define pointers that will be set by HB read function
00095 
00096   Epetra_Map * readMap;
00097   Epetra_CrsMatrix * readA; 
00098   Epetra_Vector * readx; 
00099   Epetra_Vector * readb;
00100   Epetra_Vector * readxexact;
00101    
00102   // Call function to read in HB problem
00103   Trilinos_Util_ReadHb2Epetra(argv[1], Comm, readMap, readA, readx, readb, readxexact);
00104 
00105   // Not interested in x, b or xexact for an eigenvalue problem
00106   delete readx;
00107   delete readb;
00108   delete readxexact;
00109 
00110 #ifdef EPETRA_MPI // If running in parallel, we need to distribute matrix across all PEs.
00111 
00112   // Create uniform distributed map
00113   Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
00114 
00115   // Create Exporter to distribute read-in matrix and vectors
00116 
00117   Epetra_Export exporter(*readMap, map);
00118   Epetra_CrsMatrix A(Copy, map, 0);
00119 
00120   A.Export(*readA, exporter, Add);
00121   assert(A.FillComplete()==0);
00122 
00123   delete readA;
00124   delete readMap;
00125 
00126 #else // If not running in parallel, we do not need to distribute the matrix
00127   Epetra_CrsMatrix & A = *readA;
00128 #endif
00129 
00130   // Create flop counter to collect all FLOPS
00131   Epetra_Flops counter;
00132   A.SetFlopCounter(counter);
00133 
00134   double lambda = 0; // Minimal eigenvalue returned here
00135   // Call inverse iteration solver
00136   Epetra_Time timer(Comm);
00137   invIteration(A, lambda, verbose);
00138   double elapsedTime = timer.ElapsedTime();
00139   double totalFlops = counter.Flops();
00140   double MFLOPS = totalFlops/elapsedTime/1000000.0;
00141 
00142 
00143   cout << endl
00144        << "*************************************************" << endl
00145        << " Approximate smallest eigenvalue = " << lambda << endl
00146        << "    Total Time    = " << elapsedTime << endl
00147        << "    Total FLOPS   = " << totalFlops << endl
00148        << "    Total MFLOPS  = " << MFLOPS << endl
00149        << "*************************************************" << endl;
00150 
00151   // All done
00152 #ifdef EPETRA_MPI
00153   MPI_Finalize();
00154 #else
00155   delete readA;
00156   delete readMap;  
00157 #endif
00158 
00159 return (0);
00160 }
00161 
00162 // The following functions are used to solve the problem:
00163 
00164 
00165 // *************************************************************
00166 // Computes the smallest eigenvalue of the given matrix A.
00167 // Returns result in lambda
00168 // *************************************************************
00169 
00170 int invIteration(Epetra_CrsMatrix& A, double &lambda, bool verbose) {  
00171 
00172   Ifpack_CrsRiluk * M;
00173   applyInverseSetup(A, M);
00174   Epetra_Vector q(A.RowMap());
00175   Epetra_Vector z(A.RowMap());
00176   Epetra_Vector resid(A.RowMap());
00177 
00178   Epetra_Flops * counter = A.GetFlopCounter();
00179   if (counter!=0) {
00180     q.SetFlopCounter(A);
00181     z.SetFlopCounter(A);
00182     resid.SetFlopCounter(A);
00183   }
00184 
00185   // Fill z with random Numbers
00186   z.Random();
00187 
00188   // variable needed for iteration
00189   double normz, residual;
00190 
00191   int niters = 100;
00192   double tolerance = 1.0E-6;
00193   int ierr = 1;
00194 
00195   for (int iter = 0; iter < niters; iter++)
00196     {
00197       if (verbose) 
00198   cout << endl
00199        << " ***** Performing step " << iter << " of inverse iteration ***** " << endl;
00200 
00201       z.Norm2(&normz); // Compute 2-norm of z
00202       q.Scale(1.0/normz, z);
00203       applyInverse(A, z, q, M, verbose); // Compute z such that Az = q
00204       q.Dot(z, &lambda); // Approximate maximum eigenvalue
00205       if (iter%10==0 || iter+1==niters)
00206   {
00207     resid.Update(1.0, z, -lambda, q, 0.0); // Compute A(inv)*q - lambda*q
00208     resid.Norm2(&residual);
00209     cout << endl
00210          << "***** Inverse Iteration Step " << iter+1 << endl
00211          << "  Lambda = " << 1.0/lambda << endl
00212          << "  Residual of A(inv)*q - lambda*q = " 
00213          << residual << endl;
00214   } 
00215       if (residual < tolerance) {
00216   ierr = 0;
00217   break;
00218       }
00219     }
00220   // lambda is the largest eigenvalue of A(inv).  1/lambda is smallest eigenvalue of A.
00221   lambda = 1.0/lambda;
00222 
00223   // Compute A*q - lambda*q explicitly
00224   A.Multiply(false, q, z);
00225   resid.Update(1.0, z, -lambda, q, 0.0); // Compute A*q - lambda*q
00226   resid.Norm2(&residual);
00227   cout << "  Explicitly computed residual of A*q - lambda*q = " << residual << endl;
00228   applyInverseDestroy(M);
00229   return(ierr);
00230 }
00231 
00232 // Performs any setup that can be done once to reduce the cost of the applyInverse function
00233 int applyInverseSetup(Epetra_CrsMatrix &A, Ifpack_CrsRiluk * & M) {
00234   int LevelFill = 4;
00235   int Overlap = 0;
00236   Ifpack_IlukGraph * IlukGraph = new Ifpack_IlukGraph(A.Graph(), LevelFill, Overlap);
00237     assert(IlukGraph->ConstructFilledGraph()==0);
00238     M = new Ifpack_CrsRiluk(A, *IlukGraph);
00239     M->SetFlopCounter(A);
00240     assert(M->InitValues()==0);
00241     assert(M->Factor()==0);
00242   return(0);
00243 }
00244 
00245 // *************************************************************
00246 // Solves a problem Ax = b, for a given A and b.  
00247 // M is a preconditioner computed in applyInverseSetup.
00248 // *************************************************************
00249 
00250 int applyInverse(Epetra_CrsMatrix &A, 
00251      Epetra_Vector & x, Epetra_Vector & b, Ifpack_CrsRiluk * M, bool verbose) {
00252 
00253   int Maxiter = 1000;
00254   double Tolerance = 1.0E-16;
00255   BiCGSTAB(A, x, b, M, Maxiter, Tolerance, verbose);
00256 
00257   return(0);
00258 }
00259 
00260 
00261 // *************************************************************
00262 // Releases any memory associated with the preconditioner M
00263 // *************************************************************
00264 
00265 int applyInverseDestroy(Ifpack_CrsRiluk * M) {
00266   Ifpack_IlukGraph & IlukGraph = (Ifpack_IlukGraph &) M->Graph();
00267   delete M;
00268   delete &IlukGraph;
00269   return(0);
00270 }
00271 
00272 // *************************************************************
00273 // Uses the Bi-CGSTAB iterative method to solve Ax = b
00274 // *************************************************************
00275 
00276 void BiCGSTAB(Epetra_CrsMatrix &A, 
00277         Epetra_Vector &x, 
00278         Epetra_Vector &b, 
00279         Ifpack_CrsRiluk *M, 
00280         int Maxiter, 
00281         double Tolerance, bool verbose) {
00282 
00283   // Allocate vectors needed for iterations
00284   Epetra_Vector phat(x.Map()); phat.SetFlopCounter(x);
00285   Epetra_Vector p(x.Map()); p.SetFlopCounter(x);
00286   Epetra_Vector shat(x.Map()); shat.SetFlopCounter(x);
00287   Epetra_Vector s(x.Map()); s.SetFlopCounter(x);
00288   Epetra_Vector r(x.Map()); r.SetFlopCounter(x);
00289   Epetra_Vector rtilde(x.Map()); rtilde.Random(); rtilde.SetFlopCounter(x);
00290   Epetra_Vector v(x.Map()); v.SetFlopCounter(x);
00291   
00292 
00293   A.Multiply(false, x, r); // r = A*x
00294 
00295   r.Update(1.0, b, -1.0); // r = b - A*x
00296 
00297   double r_norm, b_norm, scaled_r_norm, rhon, rhonm1 = 1.0;
00298   double alpha = 1.0, omega = 1.0, sigma;
00299   double omega_num, omega_den;
00300   r.Norm2(&r_norm);
00301   b.Norm2(&b_norm);
00302   scaled_r_norm = r_norm/b_norm;
00303   r.Dot(rtilde,&rhon);
00304 
00305   if (verbose) cout << "Initial residual = " << r_norm
00306         << " Scaled residual = " << scaled_r_norm << endl;
00307 
00308 
00309   for (int i=0; i<Maxiter; i++) { // Main iteration loop   
00310 
00311     double beta = (rhon/rhonm1) * (alpha/omega);
00312     rhonm1 = rhon;
00313 
00314     /* p    = r + beta*(p - omega*v)       */
00315     /* phat = M^-1 p                       */
00316     /* v    = A phat                       */
00317 
00318     double dtemp = - beta*omega;
00319 
00320     p.Update(1.0, r, dtemp, v, beta);
00321     if (M==0) 
00322       phat.Scale(1.0, p);
00323     else
00324       M->Solve(false, p, phat);
00325     A.Multiply(false, phat, v);
00326 
00327     
00328     rtilde.Dot(v,&sigma);
00329     alpha = rhon/sigma;    
00330 
00331     /* s = r - alpha*v                     */
00332     /* shat = M^-1 s                       */
00333     /* r = A shat (r is a tmp here for t ) */
00334 
00335     s.Update(-alpha, v, 1.0, r, 0.0);
00336     if (M==0) 
00337       shat.Scale(1.0, s);
00338     else
00339       M->Solve(false, s, shat);
00340     A.Multiply(false, shat, r);
00341 
00342     r.Dot(s, &omega_num);
00343     r.Dot(r, &omega_den);
00344     omega = omega_num / omega_den;
00345 
00346     /* x = x + alpha*phat + omega*shat */
00347     /* r = s - omega*r */
00348 
00349     x.Update(alpha, phat, omega, shat, 1.0);
00350     r.Update(1.0, s, -omega); 
00351     
00352     r.Norm2(&r_norm);
00353     scaled_r_norm = r_norm/b_norm;
00354     r.Dot(rtilde,&rhon);
00355 
00356     if (verbose) cout << "Iter "<< i << " residual = " << r_norm
00357           << " Scaled residual = " << scaled_r_norm << endl;
00358 
00359     if (scaled_r_norm < Tolerance) break;
00360   }
00361   return;
00362 }

Generated on Wed May 12 21:41:04 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7