Epetra Package Browser (Single Doxygen Collection) Development
test/OSKI/cxx_main.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright 2011 Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 
00043 int main(int argc, char **argv){
00044   return 0;
00045 }
00046 
00047 #ifdef THIS_SHOULD_NOT_BE_DEFINED
00048 
00049 #include "Epetra_Time.h"
00050 #include "Epetra_OskiMatrix.h"
00051 #include "Epetra_OskiVector.h"
00052 #include "Epetra_OskiUtils.h"
00053 #include "Epetra_OskiPermutation.h"
00054 
00055 #ifdef EPETRA_MPI
00056 #include "mpi.h"
00057 #include "Epetra_MpiComm.h"
00058 #else
00059 #include "Epetra_SerialComm.h"
00060 #endif
00061 #include "Epetra_Map.h"
00062 #include "Epetra_Vector.h"
00063 #include "Epetra_CrsMatrix.h"
00064 #include "Epetra_LinearProblem.h"
00065 #include "Teuchos_FileInputSource.hpp"
00066 #include "Teuchos_XMLObject.hpp"
00067 #include "Teuchos_XMLParameterListReader.hpp"
00068 #include "EpetraExt_BlockMapIn.h"
00069 #include "EpetraExt_CrsMatrixIn.h"
00070 #include "EpetraExt_VectorIn.h"
00071 #include "EpetraExt_MultiVectorIn.h"
00072 
00073 using namespace Teuchos;
00074 
00075 // Turn on timing
00076 //#define ML_SCALING
00077 
00078 // prototypes defined after main
00079 void ML_Exit(int mypid, const char *msg, int code);
00080 void ML_Print_Help();
00081 void ML_Read_Matrix_Dimensions(const char *filename, int *numGlobalRows, Epetra_Comm &Comm);
00082 
00083 int main(int argc, char *argv[])
00084 {
00085 #ifdef HAVE_MPI
00086   MPI_Init(&argc,&argv);
00087   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00088   int mypid = Comm.MyPID();
00089 #else
00090   Epetra_SerialComm Comm;
00091   int mypid = 0;
00092 #endif
00093 
00094   // Read XML input deck
00095   ParameterList masterList;
00096   if (argc > 1) {
00097     if (strncmp("-h",argv[1],2) == 0) {
00098       cout << "help" << endl;
00099       ML_Print_Help();
00100       ML_Exit(mypid,0,EXIT_SUCCESS);
00101     }
00102     else {
00103       int i=0,j;
00104       FILE* fid = fopen(argv[1],"r");
00105       if (fid) {
00106         i++;
00107         fclose(fid);
00108       }
00109       Comm.SumAll(&i, &j, 1);
00110       if (j!=Comm.NumProc()) {
00111         cout << "Could not open input file." << endl;
00112         ML_Print_Help();
00113         ML_Exit(mypid,0,EXIT_FAILURE);
00114       }
00115       FileInputSource fileSrc(argv[1]);
00116       XMLObject fileXML = fileSrc.getObject();
00117       XMLParameterListReader ListReader;
00118       masterList = ListReader.toParameterList(fileXML);
00119     }
00120   } else {
00121     cout << "No input file specified." << endl;
00122     ML_Print_Help();
00123     ML_Exit(mypid,0,EXIT_SUCCESS);
00124 }
00125 
00126   ParameterList *fileList, *AztecOOList;
00127   try {fileList = &(masterList.sublist("data files",true));}
00128   catch(...) {ML_Exit(mypid,"Missing \"data files\" sublist.",EXIT_FAILURE);}
00129   try {AztecOOList = &(masterList.sublist("AztecOO"));}
00130   catch(...) {ML_Exit(mypid,"Missing \"AztecOO\" sublist.",EXIT_FAILURE);}
00131 
00132 #ifdef ML_SCALING
00133    const int ntimers=4;
00134    enum {total, probBuild, precBuild, solve};
00135    ml_DblLoc timeVec[ntimers], maxTime[ntimers], minTime[ntimers];
00136 
00137   for (int i=0; i<ntimers; i++) timeVec[i].rank = Comm.MyPID();
00138   timeVec[total].value = MPI_Wtime();
00139 #endif
00140   string matrixfile = fileList->get("matrix input file","A.dat");
00141   const char *datafile = matrixfile.c_str();
00142   int numGlobalRows;
00143   ML_Read_Matrix_Dimensions(datafile, &numGlobalRows, Comm);
00144 
00145 #ifdef ML_SCALING
00146   timeVec[probBuild].value = MPI_Wtime();
00147 #endif
00148 
00149   // ===================================================== //
00150   // READ IN MATRICES FROM FILE                            //
00151   // ===================================================== //
00152 
00153   if (!mypid) printf("reading %s\n",datafile); fflush(stdout);
00154   Epetra_CrsMatrix *Amat=NULL;
00155   //Epetra_Map *RowMap=NULL;
00156   int errCode=0;
00157   //if (RowMap) errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, *RowMap, Amat);
00158   //else        errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, Comm, Amat);
00159   errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, Comm, Amat);
00160   if (errCode) ML_Exit(mypid,"error reading matrix", EXIT_FAILURE);
00161   Amat->OptimizeStorage();
00162     
00163   Epetra_Vector LHS(Amat->RowMap()); LHS.Random();
00164   Epetra_Vector RHS(Amat->RowMap()); RHS.PutScalar(0.0);
00165   Epetra_LinearProblem Problem(Amat, &LHS, &RHS);
00166 
00167 #ifdef ML_SCALING
00168   timeVec[probBuild].value = MPI_Wtime() - timeVec[probBuild].value;
00169 #endif
00170 
00171   // =========================== build preconditioner ===========================
00172   
00173 #ifdef ML_SCALING
00174   timeVec[precBuild].value = MPI_Wtime();
00175 #endif
00176 
00177   // no preconditioner right now
00178 
00179 #ifdef ML_SCALING
00180   timeVec[precBuild].value = MPI_Wtime() - timeVec[precBuild].value;
00181 #endif
00182 
00183   // =========================== outer solver =============================
00184 
00185 #ifdef ML_SCALING
00186   timeVec[solve].value = MPI_Wtime();
00187 #endif
00188   solver.SetParameters(*AztecOOList);
00189   int maxits = AztecOOList->get("Aztec iterations",250);
00190   double tol = AztecOOList->get("Aztec tolerance",1e-10);
00191 #ifdef ML_SCALING
00192   timeVec[solve].value = MPI_Wtime() - timeVec[solve].value;
00193 #endif
00194 
00195   // compute the real residual
00196   double residual;
00197   LHS.Norm2(&residual);
00198   
00199   if( Comm.MyPID()==0 ) {
00200     cout << "||b-Ax||_2 = " << residual << endl;
00201   }
00202 
00203   delete Amat;
00204   //delete RowMap;
00205 
00206 #ifdef ML_SCALING
00207   timeVec[total].value = MPI_Wtime() - timeVec[total].value;
00208 
00209   //avg
00210   double dupTime[ntimers],avgTime[ntimers];
00211   for (int i=0; i<ntimers; i++) dupTime[i] = timeVec[i].value;
00212   MPI_Reduce(dupTime,avgTime,ntimers,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
00213   for (int i=0; i<ntimers; i++) avgTime[i] = avgTime[i]/Comm.NumProc();
00214   //min
00215   MPI_Reduce(timeVec,minTime,ntimers,MPI_DOUBLE_INT,MPI_MINLOC,0,MPI_COMM_WORLD);
00216   //max
00217   MPI_Reduce(timeVec,maxTime,ntimers,MPI_DOUBLE_INT,MPI_MAXLOC,0,MPI_COMM_WORLD);
00218 
00219   if (Comm.MyPID() == 0) {
00220     printf("timing :  max (pid)  min (pid)  avg\n");
00221     printf("Problem build         :   %2.3e (%d)  %2.3e (%d)  %2.3e \n",
00222              maxTime[probBuild].value,maxTime[probBuild].rank,
00223              minTime[probBuild].value,minTime[probBuild].rank,
00224              avgTime[probBuild]);
00225     printf("Preconditioner build  :   %2.3e (%d)  %2.3e (%d)  %2.3e \n",
00226              maxTime[precBuild].value,maxTime[precBuild].rank,
00227              minTime[precBuild].value,minTime[precBuild].rank,
00228              avgTime[precBuild]);
00229     printf("Solve                 :   %2.3e (%d)  %2.3e (%d)  %2.3e \n",
00230              maxTime[solve].value,maxTime[solve].rank,
00231              minTime[solve].value,minTime[solve].rank,
00232              avgTime[solve]);
00233     printf("Total                 :   %2.3e (%d)  %2.3e (%d)  %2.3e \n",
00234              maxTime[total].value,maxTime[total].rank,
00235              minTime[total].value,minTime[total].rank,
00236              avgTime[total]);
00237   }
00238 #endif
00239 
00240 #ifdef HAVE_MPI
00241   MPI_Finalize();
00242 #endif
00243 
00244 
00245   return(EXIT_SUCCESS);
00246 } //main
00247 
00248 void ML_Exit(int mypid, const char *msg, int code)
00249 {
00250   if (!mypid && msg != 0)
00251     printf("%s\n",msg);
00252 #ifdef ML_MPI
00253   MPI_Finalize();
00254 #endif
00255   exit(code);
00256 }
00257 
00258 void ML_Print_Help()
00259 {
00260   printf("Usage: ml_scaling.exe [XML input file]\n");
00261   printf("  The XML input file must have three sublists:\n");
00262   printf("    1) \"data files\"               -- input file names\n");
00263   printf("    2) \"AztecOO\"                  -- outer Krylov options\n");
00264 }
00265 
00266 void ML_Read_Matrix_Dimensions(const char *filename, int *numGlobalRows, Epetra_Comm &Comm)
00267 {
00268     char line[35], token1[35], token2[35], token3[35], token4[35], token5[35];
00269     int lineLength = 1025;
00270     FILE *fid = fopen(filename,"r");
00271     int N, NZ;
00272     if(fgets(line, lineLength, fid)==0) {
00273       if (fid!=0) fclose(fid);
00274       ML_Exit(Comm.MyPID(),"error opening matrix file", EXIT_FAILURE);
00275     }
00276     if(sscanf(line, "%s %s %s %s %s", token1, token2, token3, token4, token5 )==0) {
00277       if (fid!=0) fclose(fid);
00278       ML_Exit(Comm.MyPID(),"error reading matrix file header", EXIT_FAILURE);
00279     }
00280     if (strcmp(token1, "%%MatrixMarket") || strcmp(token2, "matrix") ||
00281         strcmp(token3, "coordinate") || strcmp(token4, "real") ||
00282         strcmp(token5, "general"))
00283     {
00284       if (fid!=0) fclose(fid);
00285       ML_Exit(Comm.MyPID(),"error reading matrix file header", EXIT_FAILURE);
00286     }
00287     // Next, strip off header lines (which start with "%")
00288     do {
00289       if(fgets(line, lineLength, fid)==0) {
00290         if (fid!=0) fclose(fid);
00291         ML_Exit(Comm.MyPID(),"error reading matrix file comments", EXIT_FAILURE);
00292       }
00293     } while (line[0] == '%');
00294 
00295     // Next get problem dimensions: M, N, NZ
00296     if(sscanf(line, "%d %d %d", numGlobalRows, &N, &NZ)==0) {
00297       if (fid!=0) fclose(fid);
00298       ML_Exit(Comm.MyPID(),"error reading matrix file dimensions", EXIT_FAILURE);
00299     }
00300 } //ML_Read_Matrix_Dimensions()
00301 
00302 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines