Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/OSKI/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
43int main(int argc, char **argv){
44 return 0;
45}
46
47#ifdef THIS_SHOULD_NOT_BE_DEFINED
48
49#include "Epetra_Time.h"
50#include "Epetra_OskiMatrix.h"
51#include "Epetra_OskiVector.h"
52#include "Epetra_OskiUtils.h"
54
55#ifdef EPETRA_MPI
56#include "mpi.h"
57#include "Epetra_MpiComm.h"
58#else
59#include "Epetra_SerialComm.h"
60#endif
61#include "Epetra_Map.h"
62#include "Epetra_Vector.h"
63#include "Epetra_CrsMatrix.h"
65#include "Teuchos_FileInputSource.hpp"
66#include "Teuchos_XMLObject.hpp"
67#include "Teuchos_XMLParameterListReader.hpp"
68#include "EpetraExt_BlockMapIn.h"
69#include "EpetraExt_CrsMatrixIn.h"
70#include "EpetraExt_VectorIn.h"
71#include "EpetraExt_MultiVectorIn.h"
72
73using namespace Teuchos;
74
75// Turn on timing
76//#define ML_SCALING
77
78// prototypes defined after main
79void ML_Exit(int mypid, const char *msg, int code);
80void ML_Print_Help();
81void ML_Read_Matrix_Dimensions(const char *filename, int *numGlobalRows, Epetra_Comm &Comm);
82
83int main(int argc, char *argv[])
84{
85#ifdef HAVE_MPI
86 MPI_Init(&argc,&argv);
87 Epetra_MpiComm Comm(MPI_COMM_WORLD);
88 int mypid = Comm.MyPID();
89#else
91 int mypid = 0;
92#endif
93
94 // Read XML input deck
95 ParameterList masterList;
96 if (argc > 1) {
97 if (strncmp("-h",argv[1],2) == 0) {
98 cout << "help" << endl;
99 ML_Print_Help();
100 ML_Exit(mypid,0,EXIT_SUCCESS);
101 }
102 else {
103 int i=0,j;
104 FILE* fid = fopen(argv[1],"r");
105 if (fid) {
106 i++;
107 fclose(fid);
108 }
109 Comm.SumAll(&i, &j, 1);
110 if (j!=Comm.NumProc()) {
111 cout << "Could not open input file." << endl;
112 ML_Print_Help();
113 ML_Exit(mypid,0,EXIT_FAILURE);
114 }
115 FileInputSource fileSrc(argv[1]);
116 XMLObject fileXML = fileSrc.getObject();
117 XMLParameterListReader ListReader;
118 masterList = ListReader.toParameterList(fileXML);
119 }
120 } else {
121 cout << "No input file specified." << endl;
122 ML_Print_Help();
123 ML_Exit(mypid,0,EXIT_SUCCESS);
124}
125
126 ParameterList *fileList, *AztecOOList;
127 try {fileList = &(masterList.sublist("data files",true));}
128 catch(...) {ML_Exit(mypid,"Missing \"data files\" sublist.",EXIT_FAILURE);}
129 try {AztecOOList = &(masterList.sublist("AztecOO"));}
130 catch(...) {ML_Exit(mypid,"Missing \"AztecOO\" sublist.",EXIT_FAILURE);}
131
132#ifdef ML_SCALING
133 const int ntimers=4;
134 enum {total, probBuild, precBuild, solve};
135 ml_DblLoc timeVec[ntimers], maxTime[ntimers], minTime[ntimers];
136
137 for (int i=0; i<ntimers; i++) timeVec[i].rank = Comm.MyPID();
138 timeVec[total].value = MPI_Wtime();
139#endif
140 string matrixfile = fileList->get("matrix input file","A.dat");
141 const char *datafile = matrixfile.c_str();
142 int numGlobalRows;
143 ML_Read_Matrix_Dimensions(datafile, &numGlobalRows, Comm);
144
145#ifdef ML_SCALING
146 timeVec[probBuild].value = MPI_Wtime();
147#endif
148
149 // ===================================================== //
150 // READ IN MATRICES FROM FILE //
151 // ===================================================== //
152
153 if (!mypid) printf("reading %s\n",datafile); fflush(stdout);
154 Epetra_CrsMatrix *Amat=NULL;
155 //Epetra_Map *RowMap=NULL;
156 int errCode=0;
157 //if (RowMap) errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, *RowMap, Amat);
158 //else errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, Comm, Amat);
159 errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, Comm, Amat);
160 if (errCode) ML_Exit(mypid,"error reading matrix", EXIT_FAILURE);
161 Amat->OptimizeStorage();
162
163 Epetra_Vector LHS(Amat->RowMap()); LHS.Random();
164 Epetra_Vector RHS(Amat->RowMap()); RHS.PutScalar(0.0);
165 Epetra_LinearProblem Problem(Amat, &LHS, &RHS);
166
167#ifdef ML_SCALING
168 timeVec[probBuild].value = MPI_Wtime() - timeVec[probBuild].value;
169#endif
170
171 // =========================== build preconditioner ===========================
172
173#ifdef ML_SCALING
174 timeVec[precBuild].value = MPI_Wtime();
175#endif
176
177 // no preconditioner right now
178
179#ifdef ML_SCALING
180 timeVec[precBuild].value = MPI_Wtime() - timeVec[precBuild].value;
181#endif
182
183 // =========================== outer solver =============================
184
185#ifdef ML_SCALING
186 timeVec[solve].value = MPI_Wtime();
187#endif
188 solver.SetParameters(*AztecOOList);
189 int maxits = AztecOOList->get("Aztec iterations",250);
190 double tol = AztecOOList->get("Aztec tolerance",1e-10);
191#ifdef ML_SCALING
192 timeVec[solve].value = MPI_Wtime() - timeVec[solve].value;
193#endif
194
195 // compute the real residual
196 double residual;
197 LHS.Norm2(&residual);
198
199 if( Comm.MyPID()==0 ) {
200 cout << "||b-Ax||_2 = " << residual << endl;
201 }
202
203 delete Amat;
204 //delete RowMap;
205
206#ifdef ML_SCALING
207 timeVec[total].value = MPI_Wtime() - timeVec[total].value;
208
209 //avg
210 double dupTime[ntimers],avgTime[ntimers];
211 for (int i=0; i<ntimers; i++) dupTime[i] = timeVec[i].value;
212 MPI_Reduce(dupTime,avgTime,ntimers,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
213 for (int i=0; i<ntimers; i++) avgTime[i] = avgTime[i]/Comm.NumProc();
214 //min
215 MPI_Reduce(timeVec,minTime,ntimers,MPI_DOUBLE_INT,MPI_MINLOC,0,MPI_COMM_WORLD);
216 //max
217 MPI_Reduce(timeVec,maxTime,ntimers,MPI_DOUBLE_INT,MPI_MAXLOC,0,MPI_COMM_WORLD);
218
219 if (Comm.MyPID() == 0) {
220 printf("timing : max (pid) min (pid) avg\n");
221 printf("Problem build : %2.3e (%d) %2.3e (%d) %2.3e \n",
222 maxTime[probBuild].value,maxTime[probBuild].rank,
223 minTime[probBuild].value,minTime[probBuild].rank,
224 avgTime[probBuild]);
225 printf("Preconditioner build : %2.3e (%d) %2.3e (%d) %2.3e \n",
226 maxTime[precBuild].value,maxTime[precBuild].rank,
227 minTime[precBuild].value,minTime[precBuild].rank,
228 avgTime[precBuild]);
229 printf("Solve : %2.3e (%d) %2.3e (%d) %2.3e \n",
230 maxTime[solve].value,maxTime[solve].rank,
231 minTime[solve].value,minTime[solve].rank,
232 avgTime[solve]);
233 printf("Total : %2.3e (%d) %2.3e (%d) %2.3e \n",
234 maxTime[total].value,maxTime[total].rank,
235 minTime[total].value,minTime[total].rank,
236 avgTime[total]);
237 }
238#endif
239
240#ifdef HAVE_MPI
241 MPI_Finalize();
242#endif
243
244
245 return(EXIT_SUCCESS);
246} //main
247
248void ML_Exit(int mypid, const char *msg, int code)
249{
250 if (!mypid && msg != 0)
251 printf("%s\n",msg);
252#ifdef ML_MPI
253 MPI_Finalize();
254#endif
255 exit(code);
256}
257
258void ML_Print_Help()
259{
260 printf("Usage: ml_scaling.exe [XML input file]\n");
261 printf(" The XML input file must have three sublists:\n");
262 printf(" 1) \"data files\" -- input file names\n");
263 printf(" 2) \"AztecOO\" -- outer Krylov options\n");
264}
265
266void ML_Read_Matrix_Dimensions(const char *filename, int *numGlobalRows, Epetra_Comm &Comm)
267{
268 char line[35], token1[35], token2[35], token3[35], token4[35], token5[35];
269 int lineLength = 1025;
270 FILE *fid = fopen(filename,"r");
271 int N, NZ;
272 if(fgets(line, lineLength, fid)==0) {
273 if (fid!=0) fclose(fid);
274 ML_Exit(Comm.MyPID(),"error opening matrix file", EXIT_FAILURE);
275 }
276 if(sscanf(line, "%s %s %s %s %s", token1, token2, token3, token4, token5 )==0) {
277 if (fid!=0) fclose(fid);
278 ML_Exit(Comm.MyPID(),"error reading matrix file header", EXIT_FAILURE);
279 }
280 if (strcmp(token1, "%%MatrixMarket") || strcmp(token2, "matrix") ||
281 strcmp(token3, "coordinate") || strcmp(token4, "real") ||
282 strcmp(token5, "general"))
283 {
284 if (fid!=0) fclose(fid);
285 ML_Exit(Comm.MyPID(),"error reading matrix file header", EXIT_FAILURE);
286 }
287 // Next, strip off header lines (which start with "%")
288 do {
289 if(fgets(line, lineLength, fid)==0) {
290 if (fid!=0) fclose(fid);
291 ML_Exit(Comm.MyPID(),"error reading matrix file comments", EXIT_FAILURE);
292 }
293 } while (line[0] == '%');
294
295 // Next get problem dimensions: M, N, NZ
296 if(sscanf(line, "%d %d %d", numGlobalRows, &N, &NZ)==0) {
297 if (fid!=0) fclose(fid);
298 ML_Exit(Comm.MyPID(),"error reading matrix file dimensions", EXIT_FAILURE);
299 }
300} //ML_Read_Matrix_Dimensions()
301
302#endif
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
virtual int MyPID() const =0
Return my process ID.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
Epetra_LinearProblem: The Epetra Linear Problem Class.
Epetra_MpiComm: The Epetra MPI Communication Class.
Epetra_SerialComm: The Epetra Serial Communication Class.
int MyPID() const
Return my process ID.
int SumAll(double *PartialSums, double *GlobalSums, int Count) const
Epetra_SerialComm Global Sum function.
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int main(int argc, char **argv)