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