Epetra Package Browser (Single Doxygen Collection) Development
main.cc
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //               Epetra: Linear Algebra Services Package
00005 //                 Copyright 2011 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 
00043 #include <stdio.h>
00044 #include <stdlib.h>
00045 #include <ctype.h>
00046 #include <assert.h>
00047 #include <string.h>
00048 #include <math.h>
00049 #include "Petra_Comm.h"
00050 #include "Petra_Map.h"
00051 #include "Petra_RDP_MultiVector.h"
00052 #include "Petra_RDP_Vector.h"
00053 #include "Petra_RDP_DCRS_Matrix.h"
00054 #ifdef PETRA_MPI
00055 #include "mpi.h"
00056 #endif
00057 #ifndef __cplusplus
00058 #define __cplusplus
00059 #endif
00060 
00061 int main(int argc, char *argv[])
00062 {
00063   int i;
00064 
00065 #ifdef PETRA_MPI
00066   MPI_Init(&argc,&argv);
00067 #endif
00068 
00069   // get number of processors and the name of this processor
00070 
00071 #ifdef PETRA_MPI
00072   Petra_Comm& comm = *new Petra_Comm(MPI_COMM_WORLD);
00073 #else
00074   Petra_Comm& comm = *new Petra_Comm();
00075 #endif
00076 
00077   int NumProc = comm.getNumProc();
00078   int MyPID   = comm.getMyPID();
00079   cout << "Processor " << MyPID << " of " <<  NumProc << " is alive." << endl;
00080 
00081   // Get the number of local equations from the command line
00082   if (argc!=2)
00083    {
00084      if (MyPID==0) cout << "Usage: " << argv[0] << " number_of_equations" << endl;
00085     exit(1);
00086    }
00087   int numGlobalEquations = atoi(argv[1]);
00088 
00089   if (numGlobalEquations < NumProc)
00090       {
00091      if (MyPID==0)
00092        cout << "numGlobalBlocks = " << numGlobalEquations
00093       << " cannot be < number of processors = " << NumProc << endl;
00094      exit(1);
00095       }
00096 
00097   // Construct a map that puts approximately the same number of equations on each processor
00098 
00099   Petra_Map& map = *new Petra_Map(numGlobalEquations, comm);
00100 
00101   // Get update list and number of local equations from newly created map
00102   int * UpdateList = map.getUpdateList();
00103   int numLocalEquations = map.numLocalEquations();
00104 
00105   // Create an integer vector numNz that is used to build the Petra Matrix.
00106   // numNz[i] is the number of OFF-DIAGONAL term for the ith global equation on this processor
00107 
00108   int * numNz = new int[numLocalEquations];
00109 
00110   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00111   // So we need 2 off-diagonal terms (except for the first and last equation)
00112 
00113   for (i=0; i<numLocalEquations; i++)
00114     if (UpdateList[i]==0 || UpdateList[i] == numGlobalEquations-1)
00115       numNz[i] = 1;
00116     else
00117       numNz[i] = 2;
00118 
00119   // Create a Petra_Matrix
00120 
00121   Petra_RDP_DCRS_Matrix& A = *new Petra_RDP_DCRS_Matrix(map);
00122 
00123   // Allocate space using numNz
00124 
00125   assert(A.allocate(numNz)==0);
00126 
00127   // Add  rows one-at-a-time
00128   // Need some vectors to help
00129   // Off diagonal values will always be -1
00130 
00131 
00132   double *values = new double[2];
00133   values[0] = -1.0; values[1] = -1.0;
00134   int *indices = new int[2];
00135   double two = 2.0;
00136   int numEntries;
00137 
00138   for (i=0; i<numLocalEquations; i++)
00139     {
00140     if (UpdateList[i]==0)
00141       {
00142   indices[0] = 1;
00143   numEntries = 1;
00144       }
00145     else if (UpdateList[i] == numGlobalEquations-1)
00146       {
00147   indices[0] = numGlobalEquations-2;
00148   numEntries = 1;
00149       }
00150     else
00151       {
00152   indices[0] = UpdateList[i]-1;
00153   indices[1] = UpdateList[i]+1;
00154   numEntries = 2;
00155       }
00156      assert(A.putRow(UpdateList[i], numEntries, values, indices)==0);
00157      assert(A.putRow(UpdateList[i], 1, &two, UpdateList+i)==0); // Put in the diagonal entry
00158     }
00159 
00160   // Finish up
00161   assert(A.fillComplete()==0);
00162 
00163   // Create vectors for Power method
00164 
00165   Petra_RDP_Vector& q = *new Petra_RDP_Vector(map);
00166   Petra_RDP_Vector& z = *new Petra_RDP_Vector(map);
00167   Petra_RDP_Vector& resid = *new Petra_RDP_Vector(map);
00168 
00169   // Fill z with random numbers
00170   z.random();
00171 
00172   // variable needed for iteration
00173   double normz, lambda, residual;
00174 
00175   // Iterate
00176   int niters = 500*numGlobalEquations;
00177   double tolerance = 1.0e-10;
00178   for (int iter = 0; iter < niters; iter++)
00179     {
00180       z.norm2(&normz); // Compute 2-norm of z
00181       q.scaleCopy(z, 1.0/normz);
00182       A.matvec(q, z); // Compute z = A*q
00183       q.dotProd(z, &lambda); // Approximate maximum eigenvaluE
00184       if (iter%100==0 || iter+1==niters)
00185   {
00186     resid.linComb(z, -lambda, q); // Compute A*q - lambda*q
00187     resid.norm2(&residual);
00188     if (MyPID==0) cout << "Iter = " << iter << "  Lambda = " << lambda
00189            << "  Residual of A*q - lambda*q = " << residual << endl;
00190   }
00191       if (residual < tolerance) break;
00192     }
00193 
00194   // Release all objects
00195 
00196   delete [] numNz;
00197   delete [] values;
00198   delete [] indices;
00199 
00200   delete &resid;
00201   delete &z;
00202   delete &q;
00203   delete &A;
00204   delete &map;
00205   delete &comm;
00206         
00207 #ifdef PETRA_MPI
00208   MPI_Finalize() ;
00209 #endif
00210 
00211 /* end main
00212 */
00213 return 0 ;
00214 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines