Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
example/OSKI/cxx_main.cpp
Go to the documentation of this file.
1/*
2This example compares OSKI matrix operations to native Epetra operations. To invoke this example, use
3
4 cxx_main.exe filename
5
6where filename is a Matrix Market formatted data file. The first two lines *must* be of the form
7
8%%MatrixMarket matrix coordinate real general
9XX YY ZZ
10
11where the triplet XX YY ZZ contains the number of global rows, columns, and nonzeros, respectively.
12
13To compile this example, use a makefile similar to that below:
14
15### start of makefile ####
16
17include /home/ikarlin/Trilinos/build_mpi/packages/epetraext/Makefile.export.epetraext
18
19ROOT=cxx_main
20
21FILES=/home/ikarlin/OSKI/install-debug/lib/oski/liboski.a \
22/home/ikarlin/OSKI/install-debug/lib/oski/liboskilt.a \
23/home/ikarlin/OSKI/install-debug/lib/oski/liboski_Tid.a \
24/home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_CSR_Tid.a \
25/home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_CSC_Tid.a \
26/home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_BCSR_Tid.a \
27/home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_MBCSR_Tid.a \
28/home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_GCSR_Tid.a \
29/home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_CB_Tid.a \
30/home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_VBR_Tid.a \
31/home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_DENSE_Tid.a \
32/home/ikarlin/OSKI/install-debug/lib/oski/liboski_heur_regprof_Tid.a \
33/home/ikarlin/OSKI/install-debug/lib/oski/liboski_heur_symmrb_Tid.a \
34/home/ikarlin/OSKI/install-debug/lib/oski/liboski_heur_mregblock_Tid.a \
35/home/ikarlin/OSKI/install-debug/lib/oski/liboski.a
36
37
38all: ${ROOT}.exe
39
40${ROOT}.exe: ${ROOT}.o
41 mpicxx -g -O0 -o ${ROOT}.exe ${ROOT}.o ${EPETRAEXT_INCLUDES} ${EPETRAEXT_LIBS} ${FILES} -ldl -lltdl
42# mpicxx -o ${ROOT}.exe ${ROOT}.o ${EPETRAEXT_INCLUDES} ${EPETRAEXT_LIBS} -Wl,--whole-archive `/bin/cat /lib/oski/site-modules-static.txt` -Wl,--no-whole-archive -ldl -lltdl
43
44${ROOT}.o: ${ROOT}.cpp
45 mpicxx -g -O0 -c -I/include -DHAVE_CONFIG_H ${EPETRAEXT_INCLUDES} ${ROOT}.cpp
46
47### end of makefile
48
49*/
50
51//@HEADER
52// ************************************************************************
53//
54// Epetra: Linear Algebra Services Package
55// Copyright 2011 Sandia Corporation
56//
57// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
58// the U.S. Government retains certain rights in this software.
59//
60// Redistribution and use in source and binary forms, with or without
61// modification, are permitted provided that the following conditions are
62// met:
63//
64// 1. Redistributions of source code must retain the above copyright
65// notice, this list of conditions and the following disclaimer.
66//
67// 2. Redistributions in binary form must reproduce the above copyright
68// notice, this list of conditions and the following disclaimer in the
69// documentation and/or other materials provided with the distribution.
70//
71// 3. Neither the name of the Corporation nor the names of the
72// contributors may be used to endorse or promote products derived from
73// this software without specific prior written permission.
74//
75// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
76// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
77// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
78// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
79// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
80// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
81// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
82// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
83// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
84// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
85// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
86//
87// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
88//
89// ************************************************************************
90//@HEADER
91
92
93#include "Epetra_Time.h"
94
95#ifdef HAVE_OSKI
96#ifdef HAVE_EPETRA_TEUCHOS
97#include "Epetra_OskiMatrix.h"
98#include "Epetra_OskiVector.h"
99#include "Epetra_OskiUtils.h"
101#include <unistd.h>
102
103#ifdef EPETRA_MPI
104#include "mpi.h"
105#include "Epetra_MpiComm.h"
106#else
107#include "Epetra_SerialComm.h"
108#endif
109#include "Epetra_Map.h"
110#include "Epetra_Vector.h"
111#include "Epetra_CrsMatrix.h"
112#include "Epetra_LinearProblem.h"
113#include "Teuchos_FileInputSource.hpp"
114#include "Teuchos_XMLObject.hpp"
115#include "Teuchos_XMLParameterListReader.hpp"
116#include "EpetraExt_BlockMapIn.h"
117#include "EpetraExt_CrsMatrixIn.h"
118#include "EpetraExt_VectorIn.h"
119#include "EpetraExt_MultiVectorIn.h"
120
121int nonzero;
122using namespace Teuchos;
123
124int main(int argc, char *argv[])
125{
126#ifdef HAVE_MPI
127 MPI_Init(&argc,&argv);
128 Epetra_MpiComm Comm(MPI_COMM_WORLD);
129 int mypid = Comm.MyPID();
130#else
132 int mypid = 0;
133#endif
134
135 ParameterList List;
136 ParameterList List2;
137
138 const char *datafile;
139 if (argc > 1) datafile = argv[1];
140 else datafile = "A.dat";
141
142 // ===================================================== //
143 // READ IN MATRICES FROM FILE //
144 // ===================================================== //
145
146 Epetra_CrsMatrix *Amat=NULL;
147 int errCode=0;
148
149 const int lineLength = 1025;
150 char line[lineLength];
151 int Nrows,Ncols,NZ;
152 FILE *handle = fopen(datafile,"r");
153 if (handle == 0) {
154 if (mypid==0) {
155 printf("Cannot open file \"%s\" for reading.\n",datafile);
156 printf("usage: cxx_main.exe <filename>, where filename defaults to \"A.dat\".\n");
157 }
158# ifdef HAVE_MPI
159 MPI_Finalize();
160# endif
161 exit(EXIT_FAILURE);
162 }
163
164 // Strip off header lines (which start with "%")
165 do {
166 if(fgets(line, lineLength, handle)==0) {if (handle!=0) fclose(handle);}
167 } while (line[0] == '%');
168 // Get problem dimensions: #global rows, #global cols, #global nonzeros
169 if(sscanf(line,"%d %d %d", &Nrows, &Ncols, &NZ)==0) {if (handle!=0) fclose(handle);}
170 fclose(handle);
171
172 Epetra_Map* rowmap;
173 Epetra_Map* colmap;
174
175 rowmap = new Epetra_Map (Nrows, 0, Comm);
176 colmap = new Epetra_Map (Ncols, 0, Comm);
177
178 if (mypid==0) printf("Reading matrix with %d rows, %d columns, %d nonzeros from file \"%s\".\n",Nrows,Ncols,NZ,datafile);
179 errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, *rowmap, *colmap, Amat);
180 Amat->OptimizeStorage();
181
182 //begin epetra code
183 Epetra_Vector x(Amat->DomainMap()); x.Random();
184 Epetra_Vector w(Amat->DomainMap()); w.Random();
185 Epetra_Vector z(Amat->RowMap()); z.Random();
186 Epetra_Vector y(Amat->RowMap()); y.Random();
187 Epetra_MultiVector X(Amat->DomainMap(), 5); X.Random();
188 Epetra_MultiVector Y(Amat->RowMap(), 5); Y.Random();
189
190 //begin example oski code
191
192
193 Epetra_OskiUtils object; //Create an Oski utility object
194 object.Init(); //Initialize Oski now we can make Oski matrices and vectors
195
196 //Parameters to create a matrix based on.
197 List.set("zerobased", true); //We are zero based in Epetra
198 List.set("unique", true); //All entries are unique because optimize storage was called
199
200 Epetra_OskiMatrix* OskiAmat = new Epetra_OskiMatrix(*Amat, List); //Create an OskiMatrix from an Epetra Matrix
201 OskiAmat->Multiply(false, x, y); //Perform y = A*x using OSKI multiply with Epetra Vectors
202 OskiAmat->Multiply(true, y, x, 2, 1); //Perform x = 2*A^T*y + x using OSKI multiply with Epetra Vectors
203
204 //Create OskiVectors
205 Epetra_OskiVector Oskix(x);
206 Epetra_OskiVector Oskiy(y);
207 Epetra_OskiVector Oskiw(w);
208 Epetra_OskiVector Oskiz(z);
209
210 OskiAmat->Multiply(false, Oskix, Oskiy); //Perform y = A*x using OSKI multiply with Oski Vectors
211 OskiAmat->Multiply(true, Oskiy, Oskix, 2, 1); //Perform x = 2*A^T*y + x using OSKI multiply with Oski Vectors
212
213 //Create OskiMultiVectors
214 Epetra_OskiMultiVector OskiX(X);
215 Epetra_OskiMultiVector OskiY(Y);
216
217 OskiAmat->Multiply(false, OskiX, OskiY); //Perform Y = A*X using OSKI multiply with Oski Vectors
218 OskiAmat->Multiply(true, OskiY, OskiX, 2.0, 1.0); //Perform X = 2*A^T*Y + X using OSKI multiply with Oski Vectors
219
220 //Tune Multiply aggressively
221 List2.set("singleblocksize", true); //Set machine specific hints here for blocks
222 List2.set("row", 3);
223 List2.set("col", 3);
224 List2.set("alignedblocks", true);
225 OskiAmat->SetHintMultiply(false, 1.0, Oskix, 0.0, Oskiy, ALWAYS_TUNE_AGGRESSIVELY, List); //Pass routine specific hints
226 OskiAmat->SetHint(List2); //Pass matrix specific hints
227 OskiAmat->TuneMatrix(); //Tune matrix
228 char* trans;
229 trans = OskiAmat->GetMatrixTransforms(); //Get and print out transforms performed
230 std::cout << "Aggressive transforms performed are: " << trans << "\n";
231 OskiAmat->Multiply(false, Oskix, Oskiy); //Perform the tuned multiply
232
233 //Done for demonstration purposes
234 delete OskiAmat;
235 OskiAmat = new Epetra_OskiMatrix(*Amat, List); //Create an OskiMatrix from an Epetra Matrix
236
237 //Tune MultiVec moderately
238 //need tuning list params set here
239 OskiAmat->SetHintMultiply(true, 2.0, OskiX, 1.0, OskiY, ALWAYS_TUNE, List); //Pass routine specific hints
240 OskiAmat->SetHint(List2); //Pass matrix specific hints
241 OskiAmat->TuneMatrix(); //Tune matrix
242 trans = OskiAmat->GetMatrixTransforms(); //Get and print out transforms performed
243 std::cout << "Moderate transforms performed are: " << trans << "\n";
244 OskiAmat->Multiply(true, OskiX, OskiY, 2, 1); //Perform the tuned multiply
245
246 //Done for demonstration purposes
247 delete OskiAmat;
248 OskiAmat = new Epetra_OskiMatrix(*Amat, List); //Create an OskiMatrix from an Epetra Matrix
249
250 //Tune MultiVec based on calls
251 OskiAmat->SetHintMultiply(true, 2.0, OskiX, 1.0, OskiY, 10, List); //Pass routine specific hints for 10 calls
252 OskiAmat->SetHint(List2); //Pass matrix specific hints
253 OskiAmat->TuneMatrix(); //Tune matrix
254 trans = OskiAmat->GetMatrixTransforms(); //Get and print out transforms performed
255 std::cout << "Moderate transforms performed are: " << trans << "\n";
256 OskiAmat->Multiply(true, OskiX, OskiY, 2, 1); //Perform the tuned multiply
257
258 //Composed multiplies
259 //ATA
260 OskiAmat->MatTransMatMultiply(true, Oskix, Oskiw, NULL); //Perform the multi not saving the intermediate result. This will not be tuned or composed by OSKI.
261 List.set("invalidinter", true); //Tell the tune function we're not storing the intermediate.
262 OskiAmat->SetHintMatTransMatMultiply(true, 1.0, Oskix, 0.0, Oskiw, Oskiy, ALWAYS_TUNE, List); //Set our tuning hints.
263 OskiAmat->TuneMatrix(); //Tune the matrix.
264 OskiAmat->MatTransMatMultiply(true, Oskix, Oskiw, NULL); //Call the tuned matrix-vector multiply which will be composed.
265
266 //AAT In parallel AAT will never be composed and will instead always be two OSKI matvecs since AAT cannot be performed as an atomic operation
267 // without communication in parallel.
268 OskiAmat->MatTransMatMultiply(false, Oskiy, Oskiz, NULL); //Same as above
269 OskiAmat->SetHintMatTransMatMultiply(false, 1.0, Oskiy, 0.0, Oskiz, Oskiy, ALWAYS_TUNE, List); //This time lets store the intermediate value.
270 OskiAmat->TuneMatrix(); //Call the tune function
271 OskiAmat->MatTransMatMultiply(false, Oskiy, Oskiz, &Oskix); //Store the intermediate.
272
273 //2Mult
274 OskiAmat->MultiplyAndMatTransMultiply(true, Oskix, Oskiy, Oskiz, Oskiw); //Perform the two multiply routine.
275 OskiAmat->SetHintMatTransMatMultiply(true, 1.0, Oskix, 0.0, Oskiy, Oskiz, ALWAYS_TUNE, List); //Tune the routine so we use the composed routine.
276 OskiAmat->TuneMatrix(); //Tune the Matrix
277 OskiAmat->MultiplyAndMatTransMultiply(true, Oskix, Oskiy, Oskiz, Oskiw); //Call the composed routine.
278
279 OskiAmat->MultiplyAndMatTransMultiply(false, Oskix, Oskiy, Oskiw, Oskiz); //Don't transpose the second calculation.
280
281 delete OskiAmat;
282 object.Close(); //close the OSKI object allows OSKI to do any garbage collection or freeing it needs.
283 delete rowmap;
284 free(trans);
285 delete colmap;
286 delete Amat;
287
288#ifdef HAVE_MPI
289 MPI_Finalize();
290#endif
291
292 return(EXIT_SUCCESS);
293} //main
294#endif
295#endif
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...
const Epetra_Map & DomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_MpiComm: The Epetra MPI Communication Class.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
Epetra_OskiMatrix: A class for constructing and using OSKI Matrices within Epetra....
int SetHint(const Teuchos::ParameterList &List)
Stores the hints in List in the matrix structure.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Performs a matrix vector multiply of y = this^TransA*x.
int TuneMatrix()
Tunes the matrix multiply if its deemed profitable.
int MultiplyAndMatTransMultiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y, const Epetra_Vector &w, Epetra_Vector &z, double Alpha=1.0, double Beta=0.0, double Omega=1.0, double Zeta=0.0) const
Performs the two matrix vector multiplies of y = Alpha*this*x + Beta*y and z = Omega*this^TransA*w + ...
int SetHintMatTransMatMultiply(bool ATA, double Alpha, const Epetra_OskiMultiVector &InVec, double Beta, const Epetra_OskiMultiVector &OutVec, const Epetra_OskiMultiVector &Intermediate, int NumCalls, const Teuchos::ParameterList &List)
Workload hints for computing a two matrix-vector multiplies that are composed used by OskiTuneMat to ...
char * GetMatrixTransforms() const
Returns a string holding the transformations performed on the matrix when it was tuned.
int SetHintMultiply(bool TransA, double Alpha, const Epetra_OskiMultiVector &InVec, double Beta, const Epetra_OskiMultiVector &OutVec, int NumCalls, const Teuchos::ParameterList &List)
Workload hints for computing a matrix-vector multiply used by OskiTuneMat to optimize the data struct...
int MatTransMatMultiply(bool ATA, const Epetra_Vector &x, Epetra_Vector &y, Epetra_Vector *t, double Alpha=1.0, double Beta=0.0) const
Performs two matrix vector multiplies of y = Alpha*this^TransA*this*x + Beta*y or y = Alpha*this*this...
Epetra_OskiMultiVector: A class for constructing and using dense Oski multi-vectors on a single proce...
Epetra_OskiUtils: The Epetra OSKI Class to handle all operations that do not involve the use of a mat...
void Init()
Initializes OSKI.
Epetra_OskiVector: A class for constructing and using dense OSKI vectors on a single processor or a s...
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int main(int argc, char *argv[])