Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
example/petra_transpose/cxx_main.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Epetra: Linear Algebra Services Package
5 // Copyright 2011 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 
43 #include <stdio.h>
44 #include <stdlib.h>
45 #include <ctype.h>
46 #include <assert.h>
47 #include <string.h>
48 #include <math.h>
49 #include "Epetra_Comm.h"
50 #include "Epetra_Time.h"
51 #include "Epetra_Map.h"
52 #include "Epetra_BlockMap.h"
53 #include "Epetra_MultiVector.h"
54 #include "Epetra_Vector.h"
55 #include "Epetra_CrsMatrix.h"
56 #include "Epetra_CrsGraph.h"
57 #include "Epetra_Export.h"
58 #ifdef EPETRA_MPI
59 #include "mpi.h"
60 #include "Epetra_MpiComm.h"
61 #else
62 #include "Epetra_SerialComm.h"
63 #endif
64 #include "Trilinos_Util.h"
65 #ifndef __cplusplus
66 #define __cplusplus
67 #endif
68 
69 int main(int argc, char *argv[])
70 {
71 
72  int ierr = 0, i, j;
73  bool debug = false;
74 
75 #ifdef EPETRA_MPI
76 
77  // Initialize MPI
78 
79  MPI_Init(&argc,&argv);
80  int size, rank; // Number of MPI processes, My process ID
81 
82  MPI_Comm_size(MPI_COMM_WORLD, &size);
83  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
84  Epetra_MpiComm Comm(MPI_COMM_WORLD);
85 
86 #else
87 
88  int size = 1; // Serial case (not using MPI)
89  int rank = 0;
90  Epetra_SerialComm Comm;
91 
92 #endif
93  bool verbose = true;
94  /*
95  bool verbose = false;
96 
97  // Check if we should print results to standard out
98  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
99  */
100 
101 
102 
103  char tmp;
104  if (rank==0) cout << "Press any key to continue..."<< endl;
105  if (rank==0) cin >> tmp;
106  Comm.Barrier();
107 
108  int MyPID = Comm.MyPID();
109  int NumProc = Comm.NumProc();
110  if (verbose) cout << Comm << endl;
111 
112  bool verbose1 = verbose;
113 
114  // Redefine verbose to only print on PE 0
115  if (verbose && rank!=0) verbose = false;
116 
117 
118  if(argc != 2) cout << "Error: Enter name of data file on command line." << endl;
119 
120  /* Read matrix file and distribute among processors.
121  Returns with this processor's set of rows */
122 
123  int NumGlobalEquations, n_nonzeros, *bindx;
124  double * val, * xguess, * b, * xexact = 0;
125 
126  Trilinos_Util_read_hb(argv[1], Comm.MyPID(), &NumGlobalEquations, &n_nonzeros,
127  &val, &bindx, &xguess, &b, &xexact);
128 
129  int NumMyEquations, * MyGlobalEquations;
130 
131  Trilinos_Util_distrib_msr_matrix(Comm, &NumGlobalEquations, &n_nonzeros, &NumMyEquations,
132  &MyGlobalEquations, &val, &bindx, &xguess, &b, &xexact);
133 
134 
135  /* Make NumNz - number of entries in each row */
136 
137  int * NumNz = new int[NumMyEquations];
138  for (i=0; i<NumMyEquations; i++) NumNz[i] = bindx[i+1] - bindx[i] + 1;
139 
140 
141  Epetra_Map Map(NumGlobalEquations, NumMyEquations,
142  MyGlobalEquations, 0, Comm);
143 
144  Epetra_Time timer(Comm);
145 
146  double start = timer.ElapsedTime();
147  Epetra_CrsMatrix A(Copy, Map, NumNz);
148 
149  /* Add rows one-at-a-time */
150 
151  int NumIndices;
152  int * Indices;
153  double * Values;
154  for (i=0; i<NumMyEquations; i++){
155  Values = val + bindx[i];
156  Indices = bindx + bindx[i];
157  NumIndices = bindx[i+1] - bindx[i];
158  assert(A.InsertGlobalValues(MyGlobalEquations[i], NumIndices, Values, Indices)==0);
159  assert(A.InsertGlobalValues(MyGlobalEquations[i], 1, val+i, MyGlobalEquations+i)==0);
160  }
161 
162  assert(A.FillComplete()==0);
163 
164  if (verbose) cout << "\nTime to construct A = " << timer.ElapsedTime() - start << endl;
165  double * xexactt = xexact;
166  Epetra_Vector xx(Copy, Map, xexactt);
167 
168  double * bt = b;
169  Epetra_Vector bb(Copy, Map, bt);
170  Epetra_Vector bcomp(Map);
171 
172  // Sanity check
173  assert(A.Multiply(false, xx, bcomp)==0);
174  Epetra_Vector resid(Map);
175 
176  assert(resid.Update(1.0, bb, -1.0, bcomp, 0.0)==0);
177 
178  double residual;
179  assert(resid.Norm2(&residual)==0);
180  if (Comm.MyPID()==0) cout << "Sanity check: Norm of b - Ax for known x and b = " << residual << endl;
181 
182  // Way 1: This approach is probably most time-efficient, but is a little more complex because
183  // we explicitly pre-compute the local transpose. It should not use anymore memory
184  // than Way 2 since we make a view of the pre-export local transpose, which we
185  // cannot do in Way 2.
186 
187  // Extract newly constructed matrix one row at a time
188  // 1) First get the local indices to count how many nonzeros will be in the
189  // transpose graph on each processor
190 
191  start = timer.ElapsedTime();
192  Epetra_CrsGraph AG = A.Graph(); // Get matrix graph
193 
194  int NumMyCols = AG.NumMyCols();
195  int * TransNumNz = new int[NumMyCols];
196  for (i=0;i<NumMyCols; i++) TransNumNz[i] = 0;
197  for (i=0; i<NumMyEquations; i++) {
198  assert(AG.ExtractMyRowView(i, NumIndices, Indices)==0); // Get view of ith row
199  for (j=0; j<NumIndices; j++) ++TransNumNz[Indices[j]];
200  }
201 
202  int ** TransIndices = new int*[NumMyCols];
203  double ** TransValues = new double*[NumMyCols];
204 
205  for(i=0; i<NumMyCols; i++) {
206  NumIndices = TransNumNz[i];
207  if (NumIndices>0) {
208  TransIndices[i] = new int[NumIndices];
209  TransValues[i] = new double[NumIndices];
210  }
211  }
212 
213  // Now copy values and global indices into newly create transpose storage
214 
215  for (i=0;i<NumMyCols; i++) TransNumNz[i] = 0; // Reset transpose NumNz counter
216  for (i=0; i<NumMyEquations; i++) {
217  assert(A.ExtractMyRowView(i, NumIndices, Values, Indices)==0);
218  int ii = A.GRID(i);
219  for (j=0; j<NumIndices; j++) {
220  int TransRow = Indices[j];
221  int loc = TransNumNz[TransRow];
222  TransIndices[TransRow][loc] = ii;
223  TransValues[TransRow][loc] = Values[j];
224  ++TransNumNz[TransRow]; // increment counter into current transpose row
225  }
226  }
227 
228  // Build Transpose matrix with some rows being shared across processors.
229  // We will use a view here since the matrix will not be used for anything else
230 
231  const Epetra_Map & TransMap = A.ImportMap();
232 
233  Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz);
234  int * TransMyGlobalEquations = new int[NumMyCols];
235  TransMap.MyGlobalElements(TransMyGlobalEquations);
236 
237  /* Add rows one-at-a-time */
238 
239  for (i=0; i<NumMyCols; i++)
240  {
241  assert(TempTransA1.InsertGlobalValues(TransMyGlobalEquations[i],
242  TransNumNz[i], TransValues[i], TransIndices[i])==0);
243  }
244 
245  // Note: The following call to FillComplete is currently necessary because
246  // some global constants that are needed by the Export () are computed in this routine
247  assert(TempTransA1.FillComplete()==0);
248 
249  // Now that transpose matrix with shared rows is entered, create a new matrix that will
250  // get the transpose with uniquely owned rows (using the same row distribution as A).
251 
252  Epetra_CrsMatrix TransA1(Copy, Map,0);
253 
254  // Create an Export object that will move TempTransA around
255 
256  Epetra_Export Export(TransMap, Map);
257 
258  assert(TransA1.Export(TempTransA1, Export, Add)==0);
259 
260  assert(TransA1.FillComplete()==0);
261 
262 
263  if (verbose) cout << "\nTime to construct TransA1 = " << timer.ElapsedTime() - start << endl;
264 
265  // Now compute b = A' * x using the transpose option on Multiply and using
266  // created transpose matrix
267 
268  Epetra_Vector x(Map);
269  x.Random(); // Fill with random numbers
270 
271  Epetra_Vector b1(Map);
272  assert(A.Multiply(true, x, b1)==0);
273  Epetra_Vector b2(Map);
274  assert(TransA1.Multiply(false, x, b2)==0);
275 
276  assert(resid.Update(1.0, b1, -1.0, b2, 0.0)==0);
277 
278  assert(b1.Norm2(&residual)==0);
279  if (verbose) cout << "Norm of RHS using Trans = true with A = " << residual << endl;
280  assert(b2.Norm2(&residual)==0);
281  if (verbose) cout << "Norm of RHS using Trans = false with TransA1 = " << residual << endl;
282  assert(resid.Norm2(&residual)==0);
283  if (verbose) cout << "Difference between using A and TransA1 = " << residual << endl;
284 
285 
286  // Way 2: This approach is probably the simplest to code, but is probably slower.
287  // We compute the transpose by dumping entries one-at-a-time.
288 
289  // Extract newly constructed matrix one entry at a time and
290  // build Transpose matrix with some rows being shared across processors.
291 
292  // const Epetra_Map & TransMap = A.ImportMap();
293 
294  start = timer.ElapsedTime();
295 
296  Epetra_CrsMatrix TempTransA2(Copy, TransMap, 0);
297  TransMap.MyGlobalElements(TransMyGlobalEquations);
298 
299  for (int LocalRow=0; LocalRow<NumMyEquations; LocalRow++) {
300  assert(A.ExtractMyRowView(LocalRow, NumIndices, Values, Indices)==0);
301  int TransGlobalCol = A.GRID(LocalRow);
302  for (j=0; j<NumIndices; j++) {
303  int TransGlobalRow = A.GCID(Indices[j]);
304  double TransValue = Values[j];
305  assert(TempTransA2.InsertGlobalValues(TransGlobalRow, 1, &TransValue, &TransGlobalCol)>=0);
306  }
307  }
308 
309 
310 
311  // Note: The following call to FillComplete is currently necessary because
312  // some global constants that are needed by the Export () are computed in this routine
313  assert(TempTransA2.FillComplete()==0);
314 
315  if (verbose) cout << "\nTime to construct TransA2 = " << timer.ElapsedTime() - start << endl;
316 
317  // Now that transpose matrix with shared rows is entered, create a new matrix that will
318  // get the transpose with uniquely owned rows (using the same row distribution as A).
319 
320  Epetra_CrsMatrix TransA2(Copy, Map,0);
321 
322  // Create an Export object that will move TempTransA around
323 
324  // Epetra_Export Export(TransMap, Map); // Export already built
325 
326  assert(TransA2.Export(TempTransA2, Export, Add)==0);
327 
328  assert(TransA2.FillComplete()==0);
329 
330  // Now compute b = A' * x using the transpose option on Multiply and using
331  // created transpose matrix
332 
333  // Epetra_Vector x(Map);
334  // x.Random(); // Fill with random numbers
335 
336  // Epetra_Vector b1(Map);
337  assert(A.Multiply(true, x, b1)==0);
338  // Epetra_Vector b2(Map);
339  assert(TransA2.Multiply(false, x, b2)==0);
340 
341  assert(resid.Update(1.0, b1, -1.0, b2, 0.0)==0);
342 
343  assert(b1.Norm2(&residual)==0);
344  if (verbose) cout << "Norm of RHS using Trans = true with A = " << residual << endl;
345  assert(b2.Norm2(&residual)==0);
346  if (verbose) cout << "Norm of RHS using Trans = false with TransA2 = " << residual << endl;
347  assert(resid.Norm2(&residual)==0);
348  if (verbose) cout << "Difference between using A and TransA2 = " << residual << endl;
349 
350 
351  // The free's are needed because the HB utility routines still C-style calloc calls
352  free ((void *) xguess);
353  free ((void *) b);
354  free ((void *) xexact);
355  free ((void *) val);
356  free ((void *) bindx);
357  free ((void *) MyGlobalEquations);
358 
359  delete [] TransMyGlobalEquations;
360 
361  for(i=0; i<NumMyCols; i++) {
362  NumIndices = TransNumNz[i];
363  if (NumIndices>0) {
364  delete [] TransIndices[i];
365  delete [] TransValues[i];
366  }
367  }
368  delete [] TransIndices;
369  delete [] TransValues;
370 
371  delete [] NumNz;
372  delete [] TransNumNz;
373 
374 #ifdef EPETRA_MPI
375  MPI_Finalize() ;
376 #endif
377 
378 return 0 ;
379 }
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
int Random()
Set multi-vector values to random numbers.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given global row of the matrix.
Epetra_Export: This class builds an export object for efficient exporting of off-processor elements...
Definition: Epetra_Export.h:62
const Epetra_Map & ImportMap() const
Use ColMap() instead.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int MyPID() const
Return my process ID.
Epetra_MpiComm: The Epetra MPI Communication Class.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Returns a view of the specified local row values via pointers to internal data.
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Exports an Epetra_DistObject using the Epetra_Import object.
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified local row of the graph.
Epetra_SerialComm: The Epetra Serial Communication Class.
void Barrier() const
Epetra_SerialComm Barrier function.
int GRID(int LRID_in) const
Returns the global row index for give local row index, returns IndexBase-1 if we don&#39;t have this loca...
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int main(int argc, char *argv[])
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
int GCID(int LCID_in) const
Returns the global column index for give local column index, returns IndexBase-1 if we don&#39;t have thi...