example/petra_transpose/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_Time.h"
00037 #include "Epetra_Map.h"
00038 #include "Epetra_BlockMap.h"
00039 #include "Epetra_MultiVector.h"
00040 #include "Epetra_Vector.h"
00041 #include "Epetra_CrsMatrix.h"
00042 #include "Epetra_CrsGraph.h"
00043 #include "Epetra_Export.h"
00044 #ifdef EPETRA_MPI
00045 #include "mpi.h"
00046 #include "Epetra_MpiComm.h"
00047 #else
00048 #include "Epetra_SerialComm.h"
00049 #endif
00050 #include "Trilinos_Util.h"
00051 #ifndef __cplusplus
00052 #define __cplusplus
00053 #endif
00054 
00055 int main(int argc, char *argv[])
00056 {
00057   
00058   int ierr = 0, i, j;
00059   bool debug = false;
00060 
00061 #ifdef EPETRA_MPI
00062 
00063   // Initialize MPI
00064 
00065   MPI_Init(&argc,&argv);
00066   int size, rank; // Number of MPI processes, My process ID
00067 
00068   MPI_Comm_size(MPI_COMM_WORLD, &size);
00069   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00070   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00071 
00072 #else
00073 
00074   int size = 1; // Serial case (not using MPI)
00075   int rank = 0;
00076   Epetra_SerialComm Comm;
00077 
00078 #endif
00079   bool verbose = true;
00080   /*
00081   bool verbose = false;
00082 
00083   // Check if we should print results to standard out
00084   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00085   */
00086 
00087 
00088  
00089   char tmp;
00090   if (rank==0) cout << "Press any key to continue..."<< endl;
00091   if (rank==0) cin >> tmp;
00092   Comm.Barrier();
00093 
00094   int MyPID = Comm.MyPID();
00095   int NumProc = Comm.NumProc();
00096   if (verbose) cout << Comm << endl;
00097 
00098   bool verbose1 = verbose;
00099 
00100   // Redefine verbose to only print on PE 0
00101   if (verbose && rank!=0) verbose = false;
00102 
00103 
00104   if(argc != 2) cout << "Error: Enter name of data file on command line." << endl; 
00105 
00106   /* Read matrix file and distribute among processors.  
00107      Returns with this processor's set of rows */ 
00108 
00109   int NumGlobalEquations, n_nonzeros, *bindx;
00110   double * val, * xguess, * b, * xexact = 0;
00111 
00112   Trilinos_Util_read_hb(argv[1], Comm.MyPID(), &NumGlobalEquations, &n_nonzeros,
00113       &val,  &bindx, &xguess, &b, &xexact);
00114 
00115   int NumMyEquations, * MyGlobalEquations;
00116 
00117   Trilinos_Util_distrib_msr_matrix(Comm, &NumGlobalEquations, &n_nonzeros, &NumMyEquations,
00118       &MyGlobalEquations, &val, &bindx, &xguess, &b, &xexact);
00119 
00120 
00121   /* Make NumNz - number of entries in each row */
00122 
00123   int * NumNz = new int[NumMyEquations];
00124   for (i=0; i<NumMyEquations; i++) NumNz[i] = bindx[i+1] - bindx[i] + 1;
00125 
00126 
00127   Epetra_Map Map(NumGlobalEquations, NumMyEquations, 
00128       MyGlobalEquations, 0, Comm);
00129  
00130   Epetra_Time timer(Comm);
00131 
00132   double start = timer.ElapsedTime();
00133   Epetra_CrsMatrix A(Copy, Map, NumNz);
00134   
00135   /* Add  rows one-at-a-time */
00136 
00137   int NumIndices;
00138   int * Indices;
00139   double * Values;
00140   for (i=0; i<NumMyEquations; i++){
00141       Values = val + bindx[i];
00142       Indices = bindx + bindx[i];
00143       NumIndices = bindx[i+1] - bindx[i];
00144      assert(A.InsertGlobalValues(MyGlobalEquations[i], NumIndices, Values, Indices)==0);
00145      assert(A.InsertGlobalValues(MyGlobalEquations[i], 1, val+i, MyGlobalEquations+i)==0);
00146     }
00147   
00148   assert(A.FillComplete()==0);
00149   
00150   if (verbose) cout << "\nTime to construct A                = " << timer.ElapsedTime() - start << endl;
00151   double * xexactt = xexact;
00152   Epetra_Vector xx(Copy, Map, xexactt);
00153 
00154   double * bt = b;
00155   Epetra_Vector bb(Copy, Map, bt);
00156   Epetra_Vector bcomp(Map);
00157 
00158   // Sanity check
00159   assert(A.Multiply(false, xx, bcomp)==0);
00160   Epetra_Vector resid(Map); 
00161  
00162   assert(resid.Update(1.0, bb, -1.0, bcomp, 0.0)==0);
00163 
00164   double residual;
00165   assert(resid.Norm2(&residual)==0);
00166   if (Comm.MyPID()==0) cout << "Sanity check: Norm of b - Ax for known x and b = " << residual << endl;
00167 
00168   // Way 1: This approach is probably most time-efficient, but is a little more complex because
00169   //        we explicitly pre-compute the local transpose.  It should not use anymore memory
00170   //        than Way 2 since we make a view of the pre-export local transpose, which we
00171   //        cannot do in Way 2.
00172 
00173   // Extract newly constructed matrix one row at a time
00174   // 1) First get the local indices to count how many nonzeros will be in the 
00175   //    transpose graph on each processor
00176 
00177   start = timer.ElapsedTime();
00178   Epetra_CrsGraph AG = A.Graph(); // Get matrix graph
00179 
00180   int NumMyCols = AG.NumMyCols();
00181   int * TransNumNz = new int[NumMyCols];
00182   for (i=0;i<NumMyCols; i++) TransNumNz[i] = 0;
00183   for (i=0; i<NumMyEquations; i++) {
00184     assert(AG.ExtractMyRowView(i, NumIndices, Indices)==0); // Get view of ith row
00185     for (j=0; j<NumIndices; j++) ++TransNumNz[Indices[j]];
00186   }
00187 
00188   int ** TransIndices = new int*[NumMyCols];
00189   double ** TransValues = new double*[NumMyCols];
00190 
00191   for(i=0; i<NumMyCols; i++) {
00192     NumIndices = TransNumNz[i];
00193     if (NumIndices>0) {
00194       TransIndices[i] = new int[NumIndices];
00195       TransValues[i] = new double[NumIndices];
00196     }
00197   }
00198 
00199   // Now copy values and global indices into newly create transpose storage
00200 
00201   for (i=0;i<NumMyCols; i++) TransNumNz[i] = 0; // Reset transpose NumNz counter
00202   for (i=0; i<NumMyEquations; i++) {
00203     assert(A.ExtractMyRowView(i, NumIndices, Values, Indices)==0);
00204     int ii = A.GRID(i);
00205     for (j=0; j<NumIndices; j++) {
00206       int TransRow = Indices[j];
00207       int loc = TransNumNz[TransRow];
00208       TransIndices[TransRow][loc] = ii;
00209       TransValues[TransRow][loc] = Values[j];
00210       ++TransNumNz[TransRow]; // increment counter into current transpose row
00211     }
00212   }
00213 
00214   //  Build Transpose matrix with some rows being shared across processors.
00215   //  We will use a view here since the matrix will not be used for anything else
00216 
00217   const Epetra_Map & TransMap = A.ImportMap();
00218 
00219   Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz);
00220   int * TransMyGlobalEquations = new int[NumMyCols];
00221   TransMap.MyGlobalElements(TransMyGlobalEquations);
00222   
00223   /* Add  rows one-at-a-time */
00224 
00225   for (i=0; i<NumMyCols; i++)
00226     {
00227      assert(TempTransA1.InsertGlobalValues(TransMyGlobalEquations[i], 
00228               TransNumNz[i], TransValues[i], TransIndices[i])==0);
00229     }
00230  
00231   // Note: The following call to FillComplete is currently necessary because
00232   //       some global constants that are needed by the Export () are computed in this routine
00233   assert(TempTransA1.FillComplete()==0);
00234 
00235   // Now that transpose matrix with shared rows is entered, create a new matrix that will
00236   // get the transpose with uniquely owned rows (using the same row distribution as A).
00237 
00238   Epetra_CrsMatrix TransA1(Copy, Map,0);
00239 
00240   // Create an Export object that will move TempTransA around
00241 
00242   Epetra_Export Export(TransMap, Map);
00243 
00244   assert(TransA1.Export(TempTransA1, Export, Add)==0);
00245   
00246   assert(TransA1.FillComplete()==0);
00247 
00248 
00249   if (verbose) cout << "\nTime to construct TransA1          = " << timer.ElapsedTime() - start << endl;
00250 
00251   // Now compute b = A' * x using the transpose option on Multiply and using 
00252   // created transpose matrix
00253 
00254   Epetra_Vector x(Map);
00255   x.Random(); // Fill with random numbers
00256 
00257   Epetra_Vector b1(Map);
00258   assert(A.Multiply(true, x, b1)==0);
00259   Epetra_Vector b2(Map);
00260   assert(TransA1.Multiply(false, x, b2)==0);
00261  
00262   assert(resid.Update(1.0, b1, -1.0, b2, 0.0)==0);
00263 
00264   assert(b1.Norm2(&residual)==0);
00265   if (verbose) cout << "Norm of RHS using Trans = true with A           = " << residual << endl;
00266   assert(b2.Norm2(&residual)==0);
00267   if (verbose) cout << "Norm of RHS using Trans = false with TransA1    = " << residual << endl;
00268   assert(resid.Norm2(&residual)==0);
00269   if (verbose) cout << "Difference between using A and TransA1          = " << residual << endl;
00270 
00271 
00272   // Way 2: This approach is probably the simplest to code, but is probably slower.
00273   //        We compute the transpose by dumping entries one-at-a-time.  
00274 
00275   // Extract newly constructed matrix one entry at a time and
00276   //  build Transpose matrix with some rows being shared across processors.
00277 
00278   // const Epetra_Map & TransMap = A.ImportMap();
00279 
00280   start = timer.ElapsedTime();
00281 
00282   Epetra_CrsMatrix TempTransA2(Copy, TransMap, 0);
00283   TransMap.MyGlobalElements(TransMyGlobalEquations);
00284 
00285   for (int LocalRow=0; LocalRow<NumMyEquations; LocalRow++) {
00286     assert(A.ExtractMyRowView(LocalRow, NumIndices, Values, Indices)==0);
00287     int TransGlobalCol = A.GRID(LocalRow);
00288     for (j=0; j<NumIndices; j++) {
00289       int TransGlobalRow = A.GCID(Indices[j]);
00290       double TransValue = Values[j];
00291       assert(TempTransA2.InsertGlobalValues(TransGlobalRow, 1, &TransValue, &TransGlobalCol)>=0);
00292     }
00293   }
00294 
00295 
00296   
00297   // Note: The following call to FillComplete is currently necessary because
00298   //       some global constants that are needed by the Export () are computed in this routine
00299   assert(TempTransA2.FillComplete()==0);
00300 
00301   if (verbose) cout << "\nTime to construct TransA2          = " << timer.ElapsedTime() - start << endl;
00302 
00303   // Now that transpose matrix with shared rows is entered, create a new matrix that will
00304   // get the transpose with uniquely owned rows (using the same row distribution as A).
00305 
00306   Epetra_CrsMatrix TransA2(Copy, Map,0);
00307 
00308   // Create an Export object that will move TempTransA around
00309 
00310   // Epetra_Export Export(TransMap, Map); // Export already built
00311 
00312   assert(TransA2.Export(TempTransA2, Export, Add)==0);
00313   
00314   assert(TransA2.FillComplete()==0);
00315 
00316   // Now compute b = A' * x using the transpose option on Multiply and using 
00317   // created transpose matrix
00318 
00319   // Epetra_Vector x(Map);
00320   // x.Random(); // Fill with random numbers
00321 
00322   // Epetra_Vector b1(Map);
00323   assert(A.Multiply(true, x, b1)==0);
00324   // Epetra_Vector b2(Map);
00325   assert(TransA2.Multiply(false, x, b2)==0);
00326  
00327   assert(resid.Update(1.0, b1, -1.0, b2, 0.0)==0);
00328 
00329   assert(b1.Norm2(&residual)==0);
00330   if (verbose) cout << "Norm of RHS using Trans = true with A           = " << residual << endl;
00331   assert(b2.Norm2(&residual)==0);
00332   if (verbose) cout << "Norm of RHS using Trans = false with TransA2    = " << residual << endl;
00333   assert(resid.Norm2(&residual)==0);
00334   if (verbose) cout << "Difference between using A and TransA2          = " << residual << endl;
00335 
00336 
00337   // The free's are needed because the HB utility routines still C-style calloc calls
00338   free ((void *) xguess);
00339   free ((void *) b);
00340   free ((void *) xexact);
00341   free ((void *) val);
00342   free ((void *) bindx);
00343   free ((void *) MyGlobalEquations);
00344 
00345   delete [] TransMyGlobalEquations;
00346 
00347   for(i=0; i<NumMyCols; i++) {
00348     NumIndices = TransNumNz[i];
00349     if (NumIndices>0) {
00350       delete [] TransIndices[i];
00351       delete [] TransValues[i];
00352     }
00353   }
00354   delete [] TransIndices;
00355   delete [] TransValues;
00356 
00357   delete [] NumNz;
00358   delete [] TransNumNz;
00359 
00360 #ifdef EPETRA_MPI
00361   MPI_Finalize() ;
00362 #endif
00363 
00364 return 0 ;
00365 }

Generated on Thu Sep 18 12:37:56 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1