test/OSKI/cxx_main.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               ML: A Multilevel Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 int main(int argc, char **argv){
00030   return 0;
00031 }
00032 
00033 #ifdef THIS_SHOULD_NOT_BE_DEFINED
00034 
00035 #include "Epetra_Time.h"
00036 #include "Epetra_OskiMatrix.h"
00037 #include "Epetra_OskiVector.h"
00038 #include "Epetra_OskiUtils.h"
00039 #include "Epetra_OskiPermutation.h"
00040 
00041 #ifdef EPETRA_MPI
00042 #include "mpi.h"
00043 #include "Epetra_MpiComm.h"
00044 #else
00045 #include "Epetra_SerialComm.h"
00046 #endif
00047 #include "Epetra_Map.h"
00048 #include "Epetra_Vector.h"
00049 #include "Epetra_CrsMatrix.h"
00050 #include "Epetra_LinearProblem.h"
00051 #include "Teuchos_FileInputSource.hpp"
00052 #include "Teuchos_XMLObject.hpp"
00053 #include "Teuchos_XMLParameterListReader.hpp"
00054 #include "EpetraExt_BlockMapIn.h"
00055 #include "EpetraExt_CrsMatrixIn.h"
00056 #include "EpetraExt_VectorIn.h"
00057 #include "EpetraExt_MultiVectorIn.h"
00058 
00059 using namespace Teuchos;
00060 
00061 // Turn on timing
00062 //#define ML_SCALING
00063 
00064 // prototypes defined after main
00065 void ML_Exit(int mypid, const char *msg, int code);
00066 void ML_Print_Help();
00067 void ML_Read_Matrix_Dimensions(const char *filename, int *numGlobalRows, Epetra_Comm &Comm);
00068 
00069 int main(int argc, char *argv[])
00070 {
00071 #ifdef HAVE_MPI
00072   MPI_Init(&argc,&argv);
00073   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00074   int mypid = Comm.MyPID();
00075 #else
00076   Epetra_SerialComm Comm;
00077   int mypid = 0;
00078 #endif
00079 
00080   // Read XML input deck
00081   ParameterList masterList;
00082   if (argc > 1) {
00083     if (strncmp("-h",argv[1],2) == 0) {
00084       cout << "help" << endl;
00085       ML_Print_Help();
00086       ML_Exit(mypid,0,EXIT_SUCCESS);
00087     }
00088     else {
00089       int i=0,j;
00090       FILE* fid = fopen(argv[1],"r");
00091       if (fid) {
00092         i++;
00093         fclose(fid);
00094       }
00095       Comm.SumAll(&i, &j, 1);
00096       if (j!=Comm.NumProc()) {
00097         cout << "Could not open input file." << endl;
00098         ML_Print_Help();
00099         ML_Exit(mypid,0,EXIT_FAILURE);
00100       }
00101       FileInputSource fileSrc(argv[1]);
00102       XMLObject fileXML = fileSrc.getObject();
00103       XMLParameterListReader ListReader;
00104       masterList = ListReader.toParameterList(fileXML);
00105     }
00106   } else {
00107     cout << "No input file specified." << endl;
00108     ML_Print_Help();
00109     ML_Exit(mypid,0,EXIT_SUCCESS);
00110 }
00111 
00112   ParameterList *fileList, *AztecOOList;
00113   try {fileList = &(masterList.sublist("data files",true));}
00114   catch(...) {ML_Exit(mypid,"Missing \"data files\" sublist.",EXIT_FAILURE);}
00115   try {AztecOOList = &(masterList.sublist("AztecOO"));}
00116   catch(...) {ML_Exit(mypid,"Missing \"AztecOO\" sublist.",EXIT_FAILURE);}
00117 
00118 #ifdef ML_SCALING
00119    const int ntimers=4;
00120    enum {total, probBuild, precBuild, solve};
00121    ml_DblLoc timeVec[ntimers], maxTime[ntimers], minTime[ntimers];
00122 
00123   for (int i=0; i<ntimers; i++) timeVec[i].rank = Comm.MyPID();
00124   timeVec[total].value = MPI_Wtime();
00125 #endif
00126   string matrixfile = fileList->get("matrix input file","A.dat");
00127   const char *datafile = matrixfile.c_str();
00128   int numGlobalRows;
00129   ML_Read_Matrix_Dimensions(datafile, &numGlobalRows, Comm);
00130 
00131 #ifdef ML_SCALING
00132   timeVec[probBuild].value = MPI_Wtime();
00133 #endif
00134 
00135   // ===================================================== //
00136   // READ IN MATRICES FROM FILE                            //
00137   // ===================================================== //
00138 
00139   if (!mypid) printf("reading %s\n",datafile); fflush(stdout);
00140   Epetra_CrsMatrix *Amat=NULL;
00141   //Epetra_Map *RowMap=NULL;
00142   int errCode=0;
00143   //if (RowMap) errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, *RowMap, Amat);
00144   //else        errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, Comm, Amat);
00145   errCode=EpetraExt::MatrixMarketFileToCrsMatrix(datafile, Comm, Amat);
00146   if (errCode) ML_Exit(mypid,"error reading matrix", EXIT_FAILURE);
00147   Amat->OptimizeStorage();
00148     
00149   Epetra_Vector LHS(Amat->RowMap()); LHS.Random();
00150   Epetra_Vector RHS(Amat->RowMap()); RHS.PutScalar(0.0);
00151   Epetra_LinearProblem Problem(Amat, &LHS, &RHS);
00152 
00153 #ifdef ML_SCALING
00154   timeVec[probBuild].value = MPI_Wtime() - timeVec[probBuild].value;
00155 #endif
00156 
00157   // =========================== build preconditioner ===========================
00158   
00159 #ifdef ML_SCALING
00160   timeVec[precBuild].value = MPI_Wtime();
00161 #endif
00162 
00163   // no preconditioner right now
00164 
00165 #ifdef ML_SCALING
00166   timeVec[precBuild].value = MPI_Wtime() - timeVec[precBuild].value;
00167 #endif
00168 
00169   // =========================== outer solver =============================
00170 
00171 #ifdef ML_SCALING
00172   timeVec[solve].value = MPI_Wtime();
00173 #endif
00174   solver.SetParameters(*AztecOOList);
00175   int maxits = AztecOOList->get("Aztec iterations",250);
00176   double tol = AztecOOList->get("Aztec tolerance",1e-10);
00177 #ifdef ML_SCALING
00178   timeVec[solve].value = MPI_Wtime() - timeVec[solve].value;
00179 #endif
00180 
00181   // compute the real residual
00182   double residual;
00183   LHS.Norm2(&residual);
00184   
00185   if( Comm.MyPID()==0 ) {
00186     cout << "||b-Ax||_2 = " << residual << endl;
00187   }
00188 
00189   delete Amat;
00190   //delete RowMap;
00191 
00192 #ifdef ML_SCALING
00193   timeVec[total].value = MPI_Wtime() - timeVec[total].value;
00194 
00195   //avg
00196   double dupTime[ntimers],avgTime[ntimers];
00197   for (int i=0; i<ntimers; i++) dupTime[i] = timeVec[i].value;
00198   MPI_Reduce(dupTime,avgTime,ntimers,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
00199   for (int i=0; i<ntimers; i++) avgTime[i] = avgTime[i]/Comm.NumProc();
00200   //min
00201   MPI_Reduce(timeVec,minTime,ntimers,MPI_DOUBLE_INT,MPI_MINLOC,0,MPI_COMM_WORLD);
00202   //max
00203   MPI_Reduce(timeVec,maxTime,ntimers,MPI_DOUBLE_INT,MPI_MAXLOC,0,MPI_COMM_WORLD);
00204 
00205   if (Comm.MyPID() == 0) {
00206     printf("timing :  max (pid)  min (pid)  avg\n");
00207     printf("Problem build         :   %2.3e (%d)  %2.3e (%d)  %2.3e \n",
00208              maxTime[probBuild].value,maxTime[probBuild].rank,
00209              minTime[probBuild].value,minTime[probBuild].rank,
00210              avgTime[probBuild]);
00211     printf("Preconditioner build  :   %2.3e (%d)  %2.3e (%d)  %2.3e \n",
00212              maxTime[precBuild].value,maxTime[precBuild].rank,
00213              minTime[precBuild].value,minTime[precBuild].rank,
00214              avgTime[precBuild]);
00215     printf("Solve                 :   %2.3e (%d)  %2.3e (%d)  %2.3e \n",
00216              maxTime[solve].value,maxTime[solve].rank,
00217              minTime[solve].value,minTime[solve].rank,
00218              avgTime[solve]);
00219     printf("Total                 :   %2.3e (%d)  %2.3e (%d)  %2.3e \n",
00220              maxTime[total].value,maxTime[total].rank,
00221              minTime[total].value,minTime[total].rank,
00222              avgTime[total]);
00223   }
00224 #endif
00225 
00226 #ifdef HAVE_MPI
00227   MPI_Finalize();
00228 #endif
00229 
00230 
00231   return(EXIT_SUCCESS);
00232 } //main
00233 
00234 void ML_Exit(int mypid, const char *msg, int code)
00235 {
00236   if (!mypid && msg != 0)
00237     printf("%s\n",msg);
00238 #ifdef ML_MPI
00239   MPI_Finalize();
00240 #endif
00241   exit(code);
00242 }
00243 
00244 void ML_Print_Help()
00245 {
00246   printf("Usage: ml_scaling.exe [XML input file]\n");
00247   printf("  The XML input file must have three sublists:\n");
00248   printf("    1) \"data files\"               -- input file names\n");
00249   printf("    2) \"AztecOO\"                  -- outer Krylov options\n");
00250 }
00251 
00252 void ML_Read_Matrix_Dimensions(const char *filename, int *numGlobalRows, Epetra_Comm &Comm)
00253 {
00254     char line[35], token1[35], token2[35], token3[35], token4[35], token5[35];
00255     int lineLength = 1025;
00256     FILE *fid = fopen(filename,"r");
00257     int N, NZ;
00258     if(fgets(line, lineLength, fid)==0) {
00259       if (fid!=0) fclose(fid);
00260       ML_Exit(Comm.MyPID(),"error opening matrix file", EXIT_FAILURE);
00261     }
00262     if(sscanf(line, "%s %s %s %s %s", token1, token2, token3, token4, token5 )==0) {
00263       if (fid!=0) fclose(fid);
00264       ML_Exit(Comm.MyPID(),"error reading matrix file header", EXIT_FAILURE);
00265     }
00266     if (strcmp(token1, "%%MatrixMarket") || strcmp(token2, "matrix") ||
00267         strcmp(token3, "coordinate") || strcmp(token4, "real") ||
00268         strcmp(token5, "general"))
00269     {
00270       if (fid!=0) fclose(fid);
00271       ML_Exit(Comm.MyPID(),"error reading matrix file header", EXIT_FAILURE);
00272     }
00273     // Next, strip off header lines (which start with "%")
00274     do {
00275       if(fgets(line, lineLength, fid)==0) {
00276         if (fid!=0) fclose(fid);
00277         ML_Exit(Comm.MyPID(),"error reading matrix file comments", EXIT_FAILURE);
00278       }
00279     } while (line[0] == '%');
00280 
00281     // Next get problem dimensions: M, N, NZ
00282     if(sscanf(line, "%d %d %d", numGlobalRows, &N, &NZ)==0) {
00283       if (fid!=0) fclose(fid);
00284       ML_Exit(Comm.MyPID(),"error reading matrix file dimensions", EXIT_FAILURE);
00285     }
00286 } //ML_Read_Matrix_Dimensions()
00287 
00288 #endif

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