main.cc

Go to the documentation of this file.
00001 //@HEADER
00002 /*
00003 ************************************************************************
00004 
00005               Epetra: Linear Algebra Services Package 
00006                 Copyright (2001) Sandia Corporation
00007 
00008 Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 license for use of this work by or on behalf of the U.S. Government.
00010 
00011 This library is free software; you can redistribute it and/or modify
00012 it under the terms of the GNU Lesser General Public License as
00013 published by the Free Software Foundation; either version 2.1 of the
00014 License, or (at your option) any later version.
00015  
00016 This library is distributed in the hope that it will be useful, but
00017 WITHOUT ANY WARRANTY; without even the implied warranty of
00018 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 Lesser General Public License for more details.
00020  
00021 You should have received a copy of the GNU Lesser General Public
00022 License along with this library; if not, write to the Free Software
00023 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 USA
00025 Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 
00027 ************************************************************************
00028 */
00029 //@HEADER
00030 
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include <ctype.h>
00034 #include <assert.h>
00035 #include <string.h>
00036 #include <math.h>
00037 #include "Petra_Comm.h"
00038 #include "Petra_Map.h"
00039 #include "Petra_RDP_MultiVector.h"
00040 #include "Petra_RDP_Vector.h"
00041 #include "Petra_RDP_DCRS_Matrix.h"
00042 #ifdef PETRA_MPI
00043 #include "mpi.h"
00044 #endif
00045 #ifndef __cplusplus
00046 #define __cplusplus
00047 #endif
00048 
00049 int main(int argc, char *argv[])
00050 {
00051   int i;
00052 
00053 #ifdef PETRA_MPI
00054   MPI_Init(&argc,&argv);
00055 #endif
00056 
00057   // get number of processors and the name of this processor
00058  
00059 #ifdef PETRA_MPI
00060   Petra_Comm& comm = *new Petra_Comm(MPI_COMM_WORLD);
00061 #else
00062   Petra_Comm& comm = *new Petra_Comm();
00063 #endif
00064 
00065   int NumProc = comm.getNumProc();
00066   int MyPID   = comm.getMyPID();
00067   cout << "Processor " << MyPID << " of " <<  NumProc << " is alive." << endl;
00068 
00069   // Get the number of local equations from the command line
00070   if (argc!=2)
00071    {
00072      if (MyPID==0) cout << "Usage: " << argv[0] << " number_of_equations" << endl;
00073     exit(1);
00074    }
00075   int numGlobalEquations = atoi(argv[1]);
00076 
00077   if (numGlobalEquations < NumProc)
00078       {
00079      if (MyPID==0)
00080        cout << "numGlobalBlocks = " << numGlobalEquations 
00081       << " cannot be < number of processors = " << NumProc << endl;
00082      exit(1);
00083       }
00084 
00085   // Construct a map that puts approximately the same number of equations on each processor
00086 
00087   Petra_Map& map = *new Petra_Map(numGlobalEquations, comm);
00088   
00089   // Get update list and number of local equations from newly created map
00090   int * UpdateList = map.getUpdateList();
00091   int numLocalEquations = map.numLocalEquations();
00092 
00093   // Create an integer vector numNz that is used to build the Petra Matrix.
00094   // numNz[i] is the number of OFF-DIAGONAL term for the ith global equation on this processor
00095 
00096   int * numNz = new int[numLocalEquations];
00097 
00098   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00099   // So we need 2 off-diagonal terms (except for the first and last equation)
00100 
00101   for (i=0; i<numLocalEquations; i++)
00102     if (UpdateList[i]==0 || UpdateList[i] == numGlobalEquations-1)
00103       numNz[i] = 1;
00104     else
00105       numNz[i] = 2;
00106 
00107   // Create a Petra_Matrix
00108 
00109   Petra_RDP_DCRS_Matrix& A = *new Petra_RDP_DCRS_Matrix(map);
00110   
00111   // Allocate space using numNz
00112   
00113   assert(A.allocate(numNz)==0);
00114 
00115   // Add  rows one-at-a-time
00116   // Need some vectors to help
00117   // Off diagonal values will always be -1
00118 
00119 
00120   double *values = new double[2];
00121   values[0] = -1.0; values[1] = -1.0;
00122   int *indices = new int[2];
00123   double two = 2.0;
00124   int numEntries;
00125   
00126   for (i=0; i<numLocalEquations; i++)
00127     {
00128     if (UpdateList[i]==0)
00129       {
00130   indices[0] = 1;
00131   numEntries = 1;
00132       }
00133     else if (UpdateList[i] == numGlobalEquations-1)
00134       {
00135   indices[0] = numGlobalEquations-2;
00136   numEntries = 1;
00137       }
00138     else
00139       {
00140   indices[0] = UpdateList[i]-1;
00141   indices[1] = UpdateList[i]+1;
00142   numEntries = 2;
00143       }
00144      assert(A.putRow(UpdateList[i], numEntries, values, indices)==0);
00145      assert(A.putRow(UpdateList[i], 1, &two, UpdateList+i)==0); // Put in the diagonal entry
00146     }
00147   
00148   // Finish up
00149   assert(A.fillComplete()==0);
00150 
00151   // Create vectors for Power method
00152 
00153   Petra_RDP_Vector& q = *new Petra_RDP_Vector(map);
00154   Petra_RDP_Vector& z = *new Petra_RDP_Vector(map);
00155   Petra_RDP_Vector& resid = *new Petra_RDP_Vector(map);
00156 
00157   // Fill z with random numbers
00158   z.random();
00159 
00160   // variable needed for iteration
00161   double normz, lambda, residual;
00162 
00163   // Iterate
00164   int niters = 500*numGlobalEquations;
00165   double tolerance = 1.0e-10;
00166   for (int iter = 0; iter < niters; iter++)
00167     {
00168       z.norm2(&normz); // Compute 2-norm of z
00169       q.scaleCopy(z, 1.0/normz);
00170       A.matvec(q, z); // Compute z = A*q
00171       q.dotProd(z, &lambda); // Approximate maximum eigenvaluE
00172       if (iter%100==0 || iter+1==niters)
00173   {
00174     resid.linComb(z, -lambda, q); // Compute A*q - lambda*q
00175     resid.norm2(&residual);
00176     if (MyPID==0) cout << "Iter = " << iter << "  Lambda = " << lambda 
00177            << "  Residual of A*q - lambda*q = " << residual << endl;
00178   } 
00179       if (residual < tolerance) break;
00180     }
00181   
00182   // Release all objects
00183 
00184   delete [] numNz;
00185   delete [] values;
00186   delete [] indices;
00187 
00188   delete &resid;
00189   delete &z;
00190   delete &q;
00191   delete &A;
00192   delete &map;
00193   delete &comm;
00194                
00195 #ifdef PETRA_MPI
00196   MPI_Finalize() ;
00197 #endif
00198 
00199 /* end main
00200 */
00201 return 0 ;
00202 }

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