Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
69int 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;
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
378return 0 ;
379}
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs.
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified local row of the graph.
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int GCID(int LCID_in) const
Returns the global column index for give local column index, returns IndexBase-1 if we don't have thi...
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 FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
const Epetra_Map & ImportMap() const
Use ColMap() instead.
int GRID(int LRID_in) const
Returns the global row index for give local row index, returns IndexBase-1 if we don't have this loca...
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
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.
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.
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.
Epetra_Export: This class builds an export object for efficient exporting of off-processor elements.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
Epetra_MpiComm: The Epetra MPI Communication Class.
void Barrier() const
Epetra_MpiComm Barrier function.
int NumProc() const
Returns total number of processes.
int MyPID() const
Return my process ID.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int Random()
Set multi-vector values to random numbers.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Time: The Epetra Timing Class.
Definition Epetra_Time.h:75
double ElapsedTime(void) const
Epetra_Time elapsed time function.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int main(int argc, char *argv[])