Epetra Package Browser (Single Doxygen Collection) Development
example/OSKI/cxx_main.cpp
Go to the documentation of this file.
00001 /*
00002 This example compares OSKI matrix operations to native Epetra operations.  To invoke this example, use
00003 
00004   cxx_main.exe filename
00005 
00006 where filename is a Matrix Market formatted data file.  The first two lines *must* be of the form
00007 
00008 %%MatrixMarket matrix coordinate real general
00009 XX YY ZZ
00010 
00011 where the triplet XX YY ZZ contains the number of global rows, columns, and nonzeros, respectively.
00012 
00013 To compile this example, use a makefile similar to that below:
00014 
00015 ### start of makefile ####
00016 
00017 include /home/ikarlin/Trilinos/build_mpi/packages/epetraext/Makefile.export.epetraext
00018 
00019 ROOT=cxx_main
00020 
00021 FILES=/home/ikarlin/OSKI/install-debug/lib/oski/liboski.a \
00022 /home/ikarlin/OSKI/install-debug/lib/oski/liboskilt.a \
00023 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_Tid.a \
00024 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_CSR_Tid.a \
00025 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_CSC_Tid.a \
00026 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_BCSR_Tid.a \
00027 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_MBCSR_Tid.a \
00028 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_GCSR_Tid.a \
00029 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_CB_Tid.a \
00030 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_VBR_Tid.a \
00031 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_DENSE_Tid.a \
00032 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_heur_regprof_Tid.a \
00033 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_heur_symmrb_Tid.a \
00034 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_heur_mregblock_Tid.a \
00035 /home/ikarlin/OSKI/install-debug/lib/oski/liboski.a
00036 
00037 
00038 all: ${ROOT}.exe
00039 
00040 ${ROOT}.exe: ${ROOT}.o
00041         mpicxx -g -O0 -o ${ROOT}.exe ${ROOT}.o ${EPETRAEXT_INCLUDES} ${EPETRAEXT_LIBS} ${FILES} -ldl -lltdl
00042 #       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
00043 
00044 ${ROOT}.o: ${ROOT}.cpp
00045         mpicxx -g -O0 -c -I/include -DHAVE_CONFIG_H ${EPETRAEXT_INCLUDES} ${ROOT}.cpp
00046 
00047 ### end of makefile
00048 
00049 */
00050 
00051 //@HEADER
00052 // ************************************************************************
00053 // 
00054 //               Epetra: Linear Algebra Services Package 
00055 //                 Copyright 2011 Sandia Corporation
00056 // 
00057 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00058 // the U.S. Government retains certain rights in this software.
00059 //
00060 // Redistribution and use in source and binary forms, with or without
00061 // modification, are permitted provided that the following conditions are
00062 // met:
00063 //
00064 // 1. Redistributions of source code must retain the above copyright
00065 // notice, this list of conditions and the following disclaimer.
00066 //
00067 // 2. Redistributions in binary form must reproduce the above copyright
00068 // notice, this list of conditions and the following disclaimer in the
00069 // documentation and/or other materials provided with the distribution.
00070 //
00071 // 3. Neither the name of the Corporation nor the names of the
00072 // contributors may be used to endorse or promote products derived from
00073 // this software without specific prior written permission.
00074 //
00075 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00076 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00077 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00078 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00079 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00080 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00081 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00082 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00083 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00084 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00085 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00086 //
00087 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00088 // 
00089 // ************************************************************************
00090 //@HEADER
00091 
00092 
00093 #include "Epetra_Time.h"
00094 
00095 #ifdef HAVE_OSKI
00096 #ifdef HAVE_EPETRA_TEUCHOS
00097 #include "Epetra_OskiMatrix.h"
00098 #include "Epetra_OskiVector.h"
00099 #include "Epetra_OskiUtils.h"
00100 #include "Epetra_OskiPermutation.h"
00101 #include <unistd.h>
00102 
00103 #ifdef EPETRA_MPI
00104 #include "mpi.h"
00105 #include "Epetra_MpiComm.h"
00106 #else
00107 #include "Epetra_SerialComm.h"
00108 #endif
00109 #include "Epetra_Map.h"
00110 #include "Epetra_Vector.h"
00111 #include "Epetra_CrsMatrix.h"
00112 #include "Epetra_LinearProblem.h"
00113 #include "Teuchos_FileInputSource.hpp"
00114 #include "Teuchos_XMLObject.hpp"
00115 #include "Teuchos_XMLParameterListReader.hpp"
00116 #include "EpetraExt_BlockMapIn.h"
00117 #include "EpetraExt_CrsMatrixIn.h"
00118 #include "EpetraExt_VectorIn.h"
00119 #include "EpetraExt_MultiVectorIn.h"
00120 
00121 int nonzero;
00122 using namespace Teuchos;
00123 
00124 int main(int argc, char *argv[])
00125 {
00126 #ifdef HAVE_MPI
00127   MPI_Init(&argc,&argv);
00128   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00129   int mypid = Comm.MyPID();
00130 #else
00131   Epetra_SerialComm Comm;
00132   int mypid = 0;
00133 #endif
00134 
00135   ParameterList List;
00136   ParameterList List2;
00137 
00138   const char *datafile;
00139   if (argc > 1) datafile = argv[1];
00140   else          datafile = "A.dat";
00141 
00142   // ===================================================== //
00143   // READ IN MATRICES FROM FILE                            //
00144   // ===================================================== //
00145 
00146   Epetra_CrsMatrix *Amat=NULL;
00147   int errCode=0;
00148 
00149   const int lineLength = 1025;
00150   char line[lineLength];
00151   int Nrows,Ncols,NZ;
00152   FILE *handle = fopen(datafile,"r");
00153   if (handle == 0) {
00154     if (mypid==0) {
00155       printf("Cannot open file \"%s\" for reading.\n",datafile);
00156       printf("usage: cxx_main.exe <filename>, where filename defaults to \"A.dat\".\n");
00157     }
00158 #   ifdef HAVE_MPI
00159     MPI_Finalize();
00160 #   endif
00161     exit(EXIT_FAILURE);
00162   }
00163     
00164   // Strip off header lines (which start with "%")
00165   do {
00166     if(fgets(line, lineLength, handle)==0) {if (handle!=0) fclose(handle);}
00167   } while (line[0] == '%');
00168   // Get problem dimensions: #global rows, #global cols, #global nonzeros
00169   if(sscanf(line,"%d %d %d", &Nrows, &Ncols, &NZ)==0) {if (handle!=0) fclose(handle);}
00170   fclose(handle);
00171   
00172   Epetra_Map* rowmap;
00173   Epetra_Map* colmap;
00174 
00175   rowmap = new Epetra_Map (Nrows, 0, Comm);
00176   colmap = new Epetra_Map (Ncols, 0, Comm);
00177 
00178   if (mypid==0) printf("Reading matrix with %d rows, %d columns, %d nonzeros from file \"%s\".\n",Nrows,Ncols,NZ,datafile);
00179   errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, *rowmap, *colmap, Amat);
00180   Amat->OptimizeStorage();
00181 
00182   //begin epetra code
00183   Epetra_Vector x(Amat->DomainMap()); x.Random();
00184   Epetra_Vector w(Amat->DomainMap()); w.Random();
00185   Epetra_Vector z(Amat->RowMap()); z.Random();
00186   Epetra_Vector y(Amat->RowMap()); y.Random();
00187   Epetra_MultiVector X(Amat->DomainMap(), 5); X.Random();
00188   Epetra_MultiVector Y(Amat->RowMap(), 5); Y.Random();
00189 
00190   //begin example oski code
00191   
00192    
00193   Epetra_OskiUtils object; //Create an Oski utility object
00194   object.Init();  //Initialize Oski now we can make Oski matrices and vectors
00195 
00196   //Parameters to create a matrix based on.
00197   List.set("zerobased", true); //We are zero based in Epetra
00198   List.set("unique", true); //All entries are unique because optimize storage was called
00199 
00200   Epetra_OskiMatrix* OskiAmat = new Epetra_OskiMatrix(*Amat, List); //Create an OskiMatrix from an Epetra Matrix
00201   OskiAmat->Multiply(false, x, y);  //Perform y = A*x using OSKI multiply with Epetra Vectors
00202   OskiAmat->Multiply(true, y, x, 2, 1);  //Perform x = 2*A^T*y + x using OSKI multiply with Epetra Vectors
00203   
00204   //Create OskiVectors
00205   Epetra_OskiVector Oskix(x);
00206   Epetra_OskiVector Oskiy(y);
00207   Epetra_OskiVector Oskiw(w);
00208   Epetra_OskiVector Oskiz(z);
00209 
00210   OskiAmat->Multiply(false, Oskix, Oskiy);  //Perform y = A*x using OSKI multiply with Oski Vectors
00211   OskiAmat->Multiply(true, Oskiy, Oskix, 2, 1);  //Perform x = 2*A^T*y + x using OSKI multiply with Oski Vectors
00212   
00213   //Create OskiMultiVectors
00214   Epetra_OskiMultiVector OskiX(X);
00215   Epetra_OskiMultiVector OskiY(Y);
00216 
00217   OskiAmat->Multiply(false, OskiX, OskiY);  //Perform Y = A*X using OSKI multiply with Oski Vectors
00218   OskiAmat->Multiply(true, OskiY, OskiX, 2.0, 1.0);  //Perform X = 2*A^T*Y + X using OSKI multiply with Oski Vectors
00219 
00220   //Tune Multiply aggressively
00221   List2.set("singleblocksize", true); //Set machine specific hints here for blocks
00222   List2.set("row", 3); 
00223   List2.set("col", 3); 
00224   List2.set("alignedblocks", true); 
00225   OskiAmat->SetHintMultiply(false, 1.0, Oskix, 0.0, Oskiy, ALWAYS_TUNE_AGGRESSIVELY, List); //Pass routine specific hints
00226   OskiAmat->SetHint(List2); //Pass matrix specific hints
00227   OskiAmat->TuneMatrix();  //Tune matrix
00228   char* trans;  
00229   trans = OskiAmat->GetMatrixTransforms();  //Get and print out transforms performed
00230   std::cout << "Aggressive transforms performed are: " << trans << "\n";
00231   OskiAmat->Multiply(false, Oskix, Oskiy); //Perform the tuned multiply
00232 
00233   //Done for demonstration purposes
00234   delete OskiAmat;
00235   OskiAmat = new Epetra_OskiMatrix(*Amat, List); //Create an OskiMatrix from an Epetra Matrix
00236 
00237   //Tune MultiVec moderately
00238   //need tuning list params set here
00239   OskiAmat->SetHintMultiply(true, 2.0, OskiX, 1.0, OskiY, ALWAYS_TUNE, List); //Pass routine specific hints
00240   OskiAmat->SetHint(List2); //Pass matrix specific hints
00241   OskiAmat->TuneMatrix();  //Tune matrix
00242   trans = OskiAmat->GetMatrixTransforms();  //Get and print out transforms performed
00243   std::cout << "Moderate transforms performed are: " << trans << "\n";
00244   OskiAmat->Multiply(true, OskiX, OskiY, 2, 1); //Perform the tuned multiply
00245 
00246   //Done for demonstration purposes
00247   delete OskiAmat;
00248   OskiAmat = new Epetra_OskiMatrix(*Amat, List); //Create an OskiMatrix from an Epetra Matrix
00249 
00250   //Tune MultiVec based on calls
00251   OskiAmat->SetHintMultiply(true, 2.0, OskiX, 1.0, OskiY, 10, List); //Pass routine specific hints for 10 calls
00252   OskiAmat->SetHint(List2); //Pass matrix specific hints
00253   OskiAmat->TuneMatrix();  //Tune matrix
00254   trans = OskiAmat->GetMatrixTransforms();  //Get and print out transforms performed
00255   std::cout << "Moderate transforms performed are: " << trans << "\n";
00256   OskiAmat->Multiply(true, OskiX, OskiY, 2, 1); //Perform the tuned multiply
00257 
00258   //Composed multiplies
00259   //ATA
00260   OskiAmat->MatTransMatMultiply(true, Oskix, Oskiw, NULL);  //Perform the multi not saving the intermediate result.  This will not be tuned or composed by OSKI.
00261   List.set("invalidinter", true);  //Tell the tune function we're not storing the intermediate.
00262   OskiAmat->SetHintMatTransMatMultiply(true, 1.0, Oskix, 0.0, Oskiw, Oskiy, ALWAYS_TUNE, List); //Set our tuning hints.
00263   OskiAmat->TuneMatrix();  //Tune the matrix.
00264   OskiAmat->MatTransMatMultiply(true, Oskix, Oskiw, NULL); //Call the tuned matrix-vector multiply which will be composed.
00265 
00266   //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
00267   //     without communication in parallel.
00268   OskiAmat->MatTransMatMultiply(false, Oskiy, Oskiz, NULL); //Same as above
00269   OskiAmat->SetHintMatTransMatMultiply(false, 1.0, Oskiy, 0.0, Oskiz, Oskiy, ALWAYS_TUNE, List); //This time lets store the intermediate value.
00270   OskiAmat->TuneMatrix();  //Call the tune function
00271   OskiAmat->MatTransMatMultiply(false, Oskiy, Oskiz, &Oskix);  //Store the intermediate.  
00272 
00273   //2Mult
00274   OskiAmat->MultiplyAndMatTransMultiply(true, Oskix, Oskiy, Oskiz, Oskiw);  //Perform the two multiply routine.
00275   OskiAmat->SetHintMatTransMatMultiply(true, 1.0, Oskix, 0.0, Oskiy, Oskiz, ALWAYS_TUNE, List); //Tune the routine so we use the composed routine.
00276   OskiAmat->TuneMatrix();  //Tune the Matrix
00277   OskiAmat->MultiplyAndMatTransMultiply(true, Oskix, Oskiy, Oskiz, Oskiw); //Call the composed routine.
00278 
00279   OskiAmat->MultiplyAndMatTransMultiply(false, Oskix, Oskiy, Oskiw, Oskiz); //Don't transpose the second calculation.
00280 
00281   delete OskiAmat; 
00282   object.Close(); //close the OSKI object allows OSKI to do any garbage collection or freeing it needs.
00283   delete rowmap;
00284   free(trans);
00285   delete colmap;
00286   delete Amat;
00287 
00288 #ifdef HAVE_MPI
00289   MPI_Finalize();
00290 #endif
00291 
00292   return(EXIT_SUCCESS);
00293 } //main
00294 #endif
00295 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines