example/OSKI/cxx_main.cpp

Go to the documentation of this file.
00001 /*
00002 This comment includes my makefile I used to compile this example.  It is ugly but this is what I had to do to get things to work.  Good luck if you try.  IK 08-06-2008
00003 
00004 include /home/ikarlin/Trilinos/build_mpi/packages/epetraext/Makefile.export.epetraext
00005 
00006 ROOT=cxx_main
00007 
00008 FILES=/home/ikarlin/OSKI/install-debug/lib/oski/liboski.a \
00009 /home/ikarlin/OSKI/install-debug/lib/oski/liboskilt.a \
00010 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_Tid.a \
00011 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_CSR_Tid.a \
00012 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_CSC_Tid.a \
00013 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_BCSR_Tid.a \
00014 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_MBCSR_Tid.a \
00015 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_GCSR_Tid.a \
00016 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_CB_Tid.a \
00017 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_VBR_Tid.a \
00018 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_mat_DENSE_Tid.a \
00019 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_heur_regprof_Tid.a \
00020 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_heur_symmrb_Tid.a \
00021 /home/ikarlin/OSKI/install-debug/lib/oski/liboski_heur_mregblock_Tid.a \
00022 /home/ikarlin/OSKI/install-debug/lib/oski/liboski.a
00023 
00024 
00025 all: ${ROOT}.exe
00026 
00027 ${ROOT}.exe: ${ROOT}.o
00028         mpicxx -g -O0 -o ${ROOT}.exe ${ROOT}.o ${EPETRAEXT_INCLUDES} ${EPETRAEXT_LIBS} ${FILES} -ldl -lltdl
00029 #       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
00030 
00031 ${ROOT}.o: ${ROOT}.cpp
00032         mpicxx -g -O0 -c -I/include -DHAVE_CONFIG_H ${EPETRAEXT_INCLUDES} ${ROOT}.cpp
00033 */
00034 
00035 //@HEADER
00036 // ************************************************************************
00037 // 
00038 //               ML: A Multilevel Preconditioner Package
00039 //                 Copyright (2002) Sandia Corporation
00040 // 
00041 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00042 // license for use of this work by or on behalf of the U.S. Government.
00043 // 
00044 // This library is free software; you can redistribute it and/or modify
00045 // it under the terms of the GNU Lesser General Public License as
00046 // published by the Free Software Foundation; either version 2.1 of the
00047 // License, or (at your option) any later version.
00048 //  
00049 // This library is distributed in the hope that it will be useful, but
00050 // WITHOUT ANY WARRANTY; without even the implied warranty of
00051 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00052 // Lesser General Public License for more details.
00053 //  
00054 // You should have received a copy of the GNU Lesser General Public
00055 // License along with this library; if not, write to the Free Software
00056 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00057 // USA
00058 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00059 // 
00060 // ************************************************************************
00061 //@HEADER
00062 
00063 #include "Epetra_Time.h"
00064 
00065 #ifdef HAVE_OSKI
00066 #ifdef HAVE_EPETRA_TEUCHOS
00067 #include "Epetra_OskiMatrix.h"
00068 #include "Epetra_OskiVector.h"
00069 #include "Epetra_OskiUtils.h"
00070 #include "Epetra_OskiPermutation.h"
00071 #include <unistd.h>
00072 
00073 #ifdef EPETRA_MPI
00074 #include "mpi.h"
00075 #include "Epetra_MpiComm.h"
00076 #else
00077 #include "Epetra_SerialComm.h"
00078 #endif
00079 #include "Epetra_Map.h"
00080 #include "Epetra_Vector.h"
00081 #include "Epetra_CrsMatrix.h"
00082 #include "Epetra_LinearProblem.h"
00083 #include "Teuchos_FileInputSource.hpp"
00084 #include "Teuchos_XMLObject.hpp"
00085 #include "Teuchos_XMLParameterListReader.hpp"
00086 #include "EpetraExt_BlockMapIn.h"
00087 #include "EpetraExt_CrsMatrixIn.h"
00088 #include "EpetraExt_VectorIn.h"
00089 #include "EpetraExt_MultiVectorIn.h"
00090 
00091 int nonzero;
00092 using namespace Teuchos;
00093 
00094 int main(int argc, char *argv[])
00095 {
00096 #ifdef HAVE_MPI
00097   MPI_Init(&argc,&argv);
00098   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00099   int mypid = Comm.MyPID();
00100 #else
00101   Epetra_SerialComm Comm;
00102   int mypid = 0;
00103 #endif
00104 
00105   ParameterList List;
00106   ParameterList List2;
00107   string matrixfile = "A.dat";
00108   const char *datafile = matrixfile.c_str();
00109 
00110   // ===================================================== //
00111   // READ IN MATRICES FROM FILE                            //
00112   // ===================================================== //
00113 
00114   if (!mypid) printf("reading %s\n",datafile); fflush(stdout);
00115   Epetra_CrsMatrix *Amat=NULL;
00116   int errCode=0;
00117   
00118   int N[1];
00119   int NZ[1];
00120 
00121   N[0] = 90;
00122   NZ[0] = 90;
00123 
00124   std::ifstream inFile("A.dat", std::ios::in);
00125 
00126   Epetra_Map* rowmap;
00127   Epetra_Map* colmap;
00128 
00129   rowmap = new Epetra_Map (N[0], 0, Comm);
00130   colmap = new Epetra_Map (NZ[0], 0, Comm);
00131 
00132   errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, *rowmap, *colmap, Amat);
00133   Amat->OptimizeStorage();
00134 
00135   //begin epetra code
00136   Epetra_Vector x(Amat->DomainMap()); x.Random();
00137   Epetra_Vector w(Amat->DomainMap()); w.Random();
00138   Epetra_Vector z(Amat->RowMap()); z.Random();
00139   Epetra_Vector y(Amat->RowMap()); y.Random();
00140   Epetra_MultiVector X(Amat->DomainMap(), 5); X.Random();
00141   Epetra_MultiVector Y(Amat->RowMap(), 5); Y.Random();
00142 
00143   //begin example oski code
00144   
00145    
00146   Epetra_OskiUtils object; //Create an Oski utility object
00147   object.Init();  //Initialize Oski now we can make Oski matrices and vectors
00148 
00149   //Parameters to create a matrix based on.
00150   List.set("zerobased", true); //We are zero based in Epetra
00151   List.set("unique", true); //All entries are unique because optimize storage was called
00152 
00153   Epetra_OskiMatrix* OskiAmat = new Epetra_OskiMatrix(*Amat, List); //Create an OskiMatrix from an Epetra Matrix
00154   OskiAmat->Multiply(false, x, y);  //Perform y = A*x using OSKI multiply with Epetra Vectors
00155   OskiAmat->Multiply(true, y, x, 2, 1);  //Perform x = 2*A^T*y + x using OSKI multiply with Epetra Vectors
00156   
00157   //Create OskiVectors
00158   Epetra_OskiVector Oskix(x);
00159   Epetra_OskiVector Oskiy(y);
00160   Epetra_OskiVector Oskiw(w);
00161   Epetra_OskiVector Oskiz(z);
00162 
00163   OskiAmat->Multiply(false, Oskix, Oskiy);  //Perform y = A*x using OSKI multiply with Oski Vectors
00164   OskiAmat->Multiply(true, Oskiy, Oskix, 2, 1);  //Perform x = 2*A^T*y + x using OSKI multiply with Oski Vectors
00165   
00166   //Create OskiMultiVectors
00167   Epetra_OskiMultiVector OskiX(X);
00168   Epetra_OskiMultiVector OskiY(Y);
00169 
00170   OskiAmat->Multiply(false, OskiX, OskiY);  //Perform Y = A*X using OSKI multiply with Oski Vectors
00171   OskiAmat->Multiply(true, OskiY, OskiX, 2.0, 1.0);  //Perform X = 2*A^T*Y + X using OSKI multiply with Oski Vectors
00172 
00173   //Tune Multiply aggressively
00174   List2.set("singleblocksize", true); //Set machine specific hints here for blocks
00175   List2.set("row", 3); 
00176   List2.set("col", 3); 
00177   List2.set("alignedblocks", true); 
00178   OskiAmat->SetHintMultiply(false, 1.0, Oskix, 0.0, Oskiy, ALWAYS_TUNE_AGGRESSIVELY, List); //Pass routine specific hints
00179   OskiAmat->SetHint(List2); //Pass matrix specific hints
00180   OskiAmat->TuneMatrix();  //Tune matrix
00181   char* trans;  
00182   trans = OskiAmat->GetMatrixTransforms();  //Get and print out transforms performed
00183   std::cout << "Aggressive transforms performed are: " << trans << "\n";
00184   OskiAmat->Multiply(false, Oskix, Oskiy); //Perform the tuned multiply
00185 
00186   //Done for demonstration purposes
00187   delete OskiAmat;
00188   OskiAmat = new Epetra_OskiMatrix(*Amat, List); //Create an OskiMatrix from an Epetra Matrix
00189 
00190   //Tune MultiVec moderately
00191   //need tuning list params set here
00192   OskiAmat->SetHintMultiply(true, 2.0, OskiX, 1.0, OskiY, ALWAYS_TUNE, List); //Pass routine specific hints
00193   OskiAmat->SetHint(List2); //Pass matrix specific hints
00194   OskiAmat->TuneMatrix();  //Tune matrix
00195   trans = OskiAmat->GetMatrixTransforms();  //Get and print out transforms performed
00196   std::cout << "Moderate transforms performed are: " << trans << "\n";
00197   OskiAmat->Multiply(true, OskiX, OskiY, 2, 1); //Perform the tuned multiply
00198 
00199   //Done for demonstration purposes
00200   delete OskiAmat;
00201   OskiAmat = new Epetra_OskiMatrix(*Amat, List); //Create an OskiMatrix from an Epetra Matrix
00202 
00203   //Tune MultiVec based on calls
00204   OskiAmat->SetHintMultiply(true, 2.0, OskiX, 1.0, OskiY, 10, List); //Pass routine specific hints for 10 calls
00205   OskiAmat->SetHint(List2); //Pass matrix specific hints
00206   OskiAmat->TuneMatrix();  //Tune matrix
00207   trans = OskiAmat->GetMatrixTransforms();  //Get and print out transforms performed
00208   std::cout << "Moderate transforms performed are: " << trans << "\n";
00209   OskiAmat->Multiply(true, OskiX, OskiY, 2, 1); //Perform the tuned multiply
00210 
00211   //Composed multiplies
00212   //ATA
00213   OskiAmat->MatTransMatMultiply(true, Oskix, Oskiw, NULL);  //Perform the multi not saving the intermediate result.  This will not be tuned or composed by OSKI.
00214   List.set("invalidinter", true);  //Tell the tune function we're not storing the intermediate.
00215   OskiAmat->SetHintMatTransMatMultiply(true, 1.0, Oskix, 0.0, Oskiw, Oskiy, ALWAYS_TUNE, List); //Set our tuning hints.
00216   OskiAmat->TuneMatrix();  //Tune the matrix.
00217   OskiAmat->MatTransMatMultiply(true, Oskix, Oskiw, NULL); //Call the tuned matrix-vector multiply which will be composed.
00218 
00219   //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
00220   //     without communication in parallel.
00221   OskiAmat->MatTransMatMultiply(false, Oskiy, Oskiz, NULL); //Same as above
00222   OskiAmat->SetHintMatTransMatMultiply(false, 1.0, Oskiy, 0.0, Oskiz, Oskiy, ALWAYS_TUNE, List); //This time lets store the intermediate value.
00223   OskiAmat->TuneMatrix();  //Call the tune function
00224   OskiAmat->MatTransMatMultiply(false, Oskiy, Oskiz, &Oskix);  //Store the intermediate.  
00225 
00226   //2Mult
00227   OskiAmat->MultiplyAndMatTransMultiply(true, Oskix, Oskiy, Oskiz, Oskiw);  //Perform the two multiply routine.
00228   OskiAmat->SetHintMatTransMatMultiply(true, 1.0, Oskix, 0.0, Oskiy, Oskiz, ALWAYS_TUNE, List); //Tune the routine so we use the composed routine.
00229   OskiAmat->TuneMatrix();  //Tune the Matrix
00230   OskiAmat->MultiplyAndMatTransMultiply(true, Oskix, Oskiy, Oskiz, Oskiw); //Call the composed routine.
00231 
00232   OskiAmat->MultiplyAndMatTransMultiply(false, Oskix, Oskiy, Oskiw, Oskiz); //Don't transpose the second calculation.
00233 
00234   delete OskiAmat; 
00235   object.Close(); //close the OSKI object allows OSKI to do any garbage collection or freeing it needs.
00236   delete rowmap;
00237   delete [] trans;
00238   delete colmap;
00239   delete Amat;
00240 
00241 #ifdef HAVE_MPI
00242   MPI_Finalize();
00243 #endif
00244 
00245   return(EXIT_SUCCESS);
00246 } //main
00247 #endif
00248 #endif

Generated on Wed May 12 21:41:04 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7