EpetraExt Package Browser (Single Doxygen Collection) Development
test/MatrixMatrix/cxx_main.cpp
Go to the documentation of this file.
00001 
00002 #include <cstring>
00003 #include <cstdio>
00004 #include <iostream>
00005 #include <fstream>
00006 
00007 #include <Epetra_ConfigDefs.h>
00008 
00009 #ifdef EPETRA_MPI
00010 #include <mpi.h>
00011 #include <Epetra_MpiComm.h>
00012 #endif
00013 
00014 #include <Epetra_SerialComm.h>
00015 #include <Epetra_Time.h>
00016 #include <Epetra_Import.h>
00017 #include <Epetra_Map.h>
00018 #include <Epetra_LocalMap.h>
00019 #include <Epetra_CrsGraph.h>
00020 #include <Epetra_CrsMatrix.h>
00021 #include <Epetra_Vector.h>
00022 #include <EpetraExt_MatrixMatrix.h>
00023 
00024 #include <EpetraExt_BlockMapIn.h>
00025 #include <EpetraExt_CrsMatrixIn.h>
00026 #include <EpetraExt_RowMatrixOut.h>
00027 
00028 namespace EpetraExt {
00029 extern
00030 Epetra_Map* find_rows_containing_cols(const Epetra_CrsMatrix& M,
00031                                       const Epetra_Map* colmap);
00032 }
00033 
00034 int read_input_file(Epetra_Comm& Comm,
00035                     const char* input_file_name,
00036                     const char**& filenames,
00037                     int& numfiles,
00038                     int& numfilenames_allocated);
00039 
00040 int read_matrix_file_names(Epetra_Comm& Comm,
00041                            const char* input_file_name,
00042                            char*& A_file,
00043                            bool& transA,
00044                            char*& B_file,
00045                            bool& transB,
00046                            char*& C_file);
00047 
00048 int broadcast_name(Epetra_Comm& Comm, const char*& name);
00049 
00050 int create_maps(Epetra_Comm& Comm,
00051                 const char* input_file_name,
00052                 Epetra_Map*& row_map,
00053                 Epetra_Map*& col_map,
00054                 Epetra_Map*& range_map,
00055                 Epetra_Map*& domain_map);
00056 
00057 int read_matrix(const char* filename,
00058                 Epetra_Comm& Comm,
00059                 const Epetra_Map* rowmap,
00060                 Epetra_Map* colmap,
00061                 const Epetra_Map* rangemap,
00062                 const Epetra_Map* domainmap,
00063                 Epetra_CrsMatrix*& mat);
00064 
00065 int run_test(Epetra_Comm& Comm,
00066              const char* filename,
00067              bool result_mtx_to_file=false,
00068              bool verbose=false);
00069 
00070 int two_proc_test(Epetra_Comm& Comm,
00071                   bool verbose=false);
00072 
00073 int test_find_rows(Epetra_Comm& Comm);
00074 
00075 Epetra_CrsMatrix* create_epetra_crsmatrix(int numProcs,
00076                                           int localProc,
00077                                           int local_n,
00078                                           bool callFillComplete = true,
00079                                           bool symmetric = true);
00080 
00081 int time_matrix_matrix_multiply(Epetra_Comm& Comm,
00082                                 bool verbose);
00083 
00085 //Global variable!!!!
00086 char* path;
00088 
00089 int main(int argc, char** argv) {
00090 
00091 #ifdef EPETRA_MPI
00092   MPI_Init(&argc,&argv);
00093   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00094 #else
00095   Epetra_SerialComm Comm;
00096 #endif
00097 
00098   bool write_result_mtx = false;
00099   bool verbose = false;
00100   int write = 0;
00101   bool path_specified = false;
00102   char* input_file = NULL;
00103   bool input_file_specified = false;
00104 
00105   if (Comm.MyPID()==0) {
00106     for(int ii=0; ii<argc; ++ii) {
00107       if (!strcmp("-write_result", argv[ii])) write_result_mtx = true;
00108       if (!strcmp("-v", argv[ii])) verbose = true;
00109       if (!strcmp("-i", argv[ii])) {
00110         input_file = argv[ii+1];
00111         input_file_specified = true;
00112       }
00113       if (!strcmp("-d",argv[ii])) {
00114         path = argv[ii+1];
00115         path_specified = true;
00116       }
00117     }
00118     write = write_result_mtx ? 1 : 0;
00119   }
00120 #ifdef EPETRA_MPI
00121   MPI_Bcast(&write, 1, MPI_INT, 0, MPI_COMM_WORLD);
00122   if (write) write_result_mtx = true;
00123 #endif
00124 
00125   if (!path_specified) {
00126     path = new char[32];
00127     sprintf(path, "%s", ".");
00128   }
00129 
00130   int err = two_proc_test(Comm, verbose);
00131   if (err != 0) {
00132     std::cout << "two_proc_test returned err=="<<err<<std::endl;
00133     return(err);
00134   }
00135 
00136   if (!input_file_specified) {
00137     input_file = new char[64];
00138     sprintf(input_file, "%s", "infiles");
00139   }
00140 
00141   const char** filenames = NULL;
00142   int numfiles = 0;
00143   int numfilenames_allocated = 0;
00144 
00145   err = read_input_file(Comm, input_file,
00146                         filenames, numfiles, numfilenames_allocated);
00147   if (err != 0) {
00148     if (path_specified) path_specified = false;
00149     sprintf(path, "%s", "./MatrixMatrix");
00150     read_input_file(Comm, input_file,
00151                     filenames, numfiles, numfilenames_allocated);
00152   }
00153 
00154   err = test_find_rows(Comm);
00155   if (err != 0) {
00156     std::cout << "test_find_rows returned err=="<<err<<std::endl;
00157     return(err);
00158   }
00159 
00160   for(int i=0; i<numfiles; ++i) {
00161     err = run_test(Comm, filenames[i], write_result_mtx, verbose);
00162     delete [] filenames[i];
00163     if (err != 0) break;
00164   }
00165 
00166   for(int j=numfiles; j<numfilenames_allocated; ++j) {
00167     delete [] filenames[j];
00168   }
00169 
00170   delete [] filenames;
00171 
00172   if (!input_file_specified) delete [] input_file;
00173   if (!path_specified) delete [] path;
00174 
00175   Epetra_CrsMatrix* D = create_epetra_crsmatrix(Comm.NumProc(),
00176                                                 Comm.MyPID(), 10,
00177                                                 true, false);
00178 
00179   std::cout << "D: \n"  << *D << std::endl;
00180 
00181   EpetraExt::MatrixMatrix::Add(*D, true, 0.5, *D, 0.5);
00182 
00183 //  std::cout << "symm D: \n"  << *D << std::endl;
00184 
00185   delete D;
00186 
00187   if (err == 0) {
00188     err = time_matrix_matrix_multiply(Comm, verbose);
00189   }
00190 
00191   int global_err = err;
00192 
00193 #ifdef EPETRA_MPI
00194   MPI_Allreduce(&err, &global_err, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
00195   MPI_Finalize();
00196 #endif
00197 
00198   return(global_err);
00199 
00200 }
00201 
00202 int test_find_rows(Epetra_Comm& Comm)
00203 {
00204   int numprocs = Comm.NumProc();
00205   int localproc = Comm.MyPID();
00206   int numlocalrows = 2;
00207   int numglobalrows = numprocs*numlocalrows;
00208   Epetra_Map rowmap(numlocalrows*numprocs, 0, Comm);
00209   Epetra_CrsMatrix matrix(Copy, rowmap, numglobalrows);
00210 
00211   int err = 0;
00212   int* cols = new int[numglobalrows];
00213   double*vals = new double[numglobalrows];
00214 
00215   for(int j=0; j<numglobalrows; ++j) {
00216     cols[j] = j;
00217     vals[j] = 1.0;
00218   }
00219 
00220   Epetra_Map colmap(-1, numglobalrows, cols, 0, Comm);
00221 
00222   for(int i=0; i<numlocalrows; ++i) {
00223     int row = localproc*numlocalrows+i;
00224     err = matrix.InsertGlobalValues(row, 1, &(vals[i]), &row);
00225     if (err != 0) {
00226       return(err);
00227     }
00228   }
00229 
00230   err = matrix.FillComplete();
00231   if (err != 0) {
00232     return(err);
00233   }
00234 
00235   Epetra_Map* map_rows = EpetraExt::find_rows_containing_cols(matrix,
00236                     &colmap);
00237 
00238   if (map_rows->NumMyElements() != numglobalrows) {
00239     return(-1);
00240   }
00241 
00242   delete map_rows;
00243   delete [] cols;
00244   delete [] vals;
00245 
00246   return(0);
00247 }
00248 
00249 int expand_name_list(const char* newname,
00250                      const char**& names,
00251                      int& alloc_len,
00252                      int& num_names)
00253 {
00254   int offset = num_names;
00255   if (offset >= alloc_len) {
00256     int alloc_increment = 8;
00257     const char** newlist = new const char*[alloc_len+alloc_increment];
00258     for(int i=0; i<offset; ++i) {
00259       newlist[i] = names[i];
00260     }
00261     delete [] names;
00262     names = newlist;
00263     alloc_len += alloc_increment;
00264     for(int i=offset; i<alloc_len; ++i) {
00265       names[i] = NULL;
00266     }
00267   }
00268 
00269   names[offset] = newname;
00270   ++num_names;
00271   return(0);
00272 }
00273 
00274 int broadcast_name(Epetra_Comm& Comm, const char*& name)
00275 {
00276   if (Comm.NumProc() < 2) return(0);
00277 
00278 #ifdef EPETRA_MPI
00279   int len;
00280   int localProc = Comm.MyPID();
00281   if (localProc == 0) {
00282     len = (int)strlen(name)+1;
00283     
00284     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00285     MPI_Bcast((void*)name, len, MPI_CHAR, 0, MPI_COMM_WORLD);
00286   }
00287   else {
00288     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00289     name = new char[len];
00290     MPI_Bcast((void*)name, len, MPI_CHAR, 0, MPI_COMM_WORLD);
00291   }
00292 
00293 #endif
00294   return(0);
00295 }
00296 
00297 int read_input_file(Epetra_Comm& Comm,
00298                     const char* input_file_name,
00299                     const char**& filenames,
00300                     int& numfiles,
00301                     int& numfilenames_allocated)
00302 {
00303   int local_err = 0, global_err = 0;
00304   std::ifstream* infile = NULL;
00305   int pathlen = path != 0 ? (int)strlen(path): 0;
00306 
00307   if (Comm.MyPID() == 0) {
00308     char* full_name = NULL;
00309     int filenamelen = input_file_name != 0 ? (int)strlen(input_file_name) : 0;
00310 
00311     full_name = new char[pathlen+filenamelen+2];
00312     if (path != 0) {
00313       sprintf(full_name, "%s/%s",path,input_file_name);
00314     }
00315     else {
00316       sprintf(full_name, "%s", input_file_name);
00317     }
00318 
00319     infile = new std::ifstream(full_name);
00320     if (!(*infile)) {
00321       local_err = -1;
00322       delete infile;
00323     }
00324     delete [] full_name;
00325   }
00326 
00327   Comm.SumAll(&local_err, &global_err, 1);
00328 
00329   if (global_err != 0) {
00330     return(global_err);
00331   }
00332 
00333 
00334   if (Comm.MyPID() == 0) {
00335     int linelen = 512;
00336     char* line = NULL;
00337 
00338     std::ifstream& ifile = *infile;
00339     while(!ifile.eof()) {
00340       line = new char[pathlen+1+linelen];
00341       if (pathlen>0) {
00342         sprintf(line,"%s/",path);
00343         ifile.getline(&(line[pathlen+1]), linelen);
00344       }
00345       else {
00346         ifile.getline(line, linelen);
00347       }
00348 
00349       if (ifile.fail()) {
00350   delete [] line;
00351         break;
00352       }
00353       if (strchr(line, '#') == NULL) {
00354         expand_name_list(line, filenames, numfilenames_allocated, numfiles);
00355       }
00356       else {
00357         delete [] line;
00358       }
00359     }
00360 
00361 #ifdef EPETRA_MPI
00362     MPI_Bcast(&numfiles, 1, MPI_INT, 0, MPI_COMM_WORLD);
00363 #endif
00364     for(int i=0; i<numfiles; ++i) {
00365       broadcast_name(Comm, filenames[i]);
00366     }
00367 
00368     delete infile;
00369   }
00370   else {
00371 #ifdef EPETRA_MPI
00372     MPI_Bcast(&numfiles, 1, MPI_INT, 0, MPI_COMM_WORLD);
00373 #endif
00374     filenames = new const char*[numfiles];
00375     numfilenames_allocated = numfiles;
00376     for(int i=0; i<numfiles; ++i) {
00377       broadcast_name(Comm, filenames[i]);
00378     }
00379   }
00380   
00381   return(0);
00382 }
00383 
00384 int run_test(Epetra_Comm& Comm,
00385              const char* filename,
00386              bool result_mtx_to_file,
00387              bool verbose)
00388 {
00389   char* A_file = NULL;
00390   char AT[3]; AT[0] = '^'; AT[1] = 'T'; AT[2] = '\0';
00391   char* B_file = NULL;
00392   char BT[3]; BT[0] = '^'; BT[1] = 'T'; BT[2] = '\0';
00393   char* C_file = NULL;
00394   bool transA, transB;
00395 
00396   int err = read_matrix_file_names(Comm, filename, A_file, transA,
00397                                    B_file, transB, C_file);
00398   if (err != 0) {
00399     std::cout << "Error, read_matrix_file_names returned " << err << std::endl;
00400     return(err);
00401   }
00402 
00403   if (!transA) AT[0] = '\0';
00404   if (!transB) BT[0] = '\0';
00405 
00406   int localProc = Comm.MyPID();
00407 
00408   if (localProc == 0 && verbose) {
00409     std::cout << "Testing C=A"<<AT<<"*B"<<BT<< "; A:" << A_file
00410               << ", B:" << B_file << ", C:" << C_file << std::endl;
00411   }
00412 
00413   Epetra_CrsMatrix* A = NULL;
00414   Epetra_CrsMatrix* B = NULL;
00415   Epetra_CrsMatrix* C = NULL;
00416   Epetra_CrsMatrix* C_check = NULL;
00417 
00418   Epetra_Map* A_row_map = NULL;
00419   Epetra_Map* A_col_map = NULL;
00420   Epetra_Map* A_range_map = NULL;
00421   Epetra_Map* A_domain_map = NULL;
00422   err = create_maps(Comm, A_file, A_row_map, A_col_map, A_range_map, A_domain_map);
00423   if (err != 0) {
00424     std::cout << "create_maps A returned err=="<<err<<std::endl;
00425     return(err);
00426   }
00427 
00428   Epetra_Map* B_row_map = NULL;
00429   Epetra_Map* B_col_map = NULL;
00430   Epetra_Map* B_range_map = NULL;
00431   Epetra_Map* B_domain_map = NULL;
00432   err = create_maps(Comm, B_file, B_row_map, B_col_map, B_range_map, B_domain_map);
00433   if (err != 0) {
00434     std::cout << "create_maps A returned err=="<<err<<std::endl;
00435     return(err);
00436   }
00437 
00438   err = read_matrix(A_file, Comm, A_row_map, A_col_map,
00439                     A_range_map, A_domain_map, A);
00440   delete [] A_file;
00441   if (err != 0) {
00442     std::cout << "read_matrix A returned err=="<<err<<std::endl;
00443     return(err);
00444   }
00445 
00446   err = read_matrix(B_file, Comm, B_row_map, B_col_map,
00447                     B_range_map, B_domain_map, B);
00448   delete [] B_file;
00449   if (err != 0) {
00450     std::cout << "read_matrix B returned err=="<<err<<std::endl;
00451     return(-1);
00452   }
00453 
00454   const Epetra_Map* rowmap = transA ? &(A->DomainMap()) : &(A->RowMap());
00455 
00456   C = new Epetra_CrsMatrix(Copy, *rowmap, 1);
00457 
00458   err = EpetraExt::MatrixMatrix::Multiply(*A, transA, *B, transB, *C);
00459   if (err != 0) {
00460     std::cout << "err "<<err<<" from MatrixMatrix::Multiply"<<std::endl;
00461     return(err);
00462   }
00463 
00464 //  std::cout << "A: " << *A << std::endl << "B: "<<*B<<std::endl<<"C: "<<*C<<std::endl;
00465   if (result_mtx_to_file) {
00466     EpetraExt::RowMatrixToMatrixMarketFile("result.mtx", *C);
00467   }
00468 
00469   Epetra_Map* Cck_row_map = NULL;
00470   Epetra_Map* Cck_col_map = NULL;
00471   Epetra_Map* Cck_range_map = NULL;
00472   Epetra_Map* Cck_domain_map = NULL;
00473   err = create_maps(Comm, C_file, Cck_row_map, Cck_col_map,
00474                     Cck_range_map, Cck_domain_map);
00475   if (err != 0) {
00476     std::cout << "create_maps C returned err=="<<err<<std::endl;
00477     return(err);
00478   }
00479 
00480   err = read_matrix(C_file, Comm, Cck_row_map, Cck_col_map,
00481                      Cck_range_map, Cck_domain_map, C_check);
00482   delete [] C_file;
00483   if (err != 0) {
00484     std::cout << "read_matrix C returned err=="<<err<<std::endl;
00485     return(-1);
00486   }
00487 
00488   err = EpetraExt::MatrixMatrix::Multiply(*A, transA, *B, transB, *C);
00489   if (err != 0) {
00490     std::cout << "err "<<err<<" from MatrixMatrix::Multiply"<<std::endl;
00491     return(err);
00492   }
00493 
00494   EpetraExt::MatrixMatrix::Add(*C, false, -1.0, *C_check, 1.0);
00495 
00496   double inf_norm = C_check->NormInf();
00497 
00498   int return_code = 0;
00499 
00500   if (inf_norm < 1.e-13) {
00501     if (localProc == 0 && verbose) {
00502       std::cout << "Test Passed" << std::endl;
00503     }
00504   }
00505   else {
00506     return_code = -1;
00507     if (localProc == 0) {
00508       std::cout << "Test Failed, inf_norm = " << inf_norm << std::endl;
00509     }
00510   }
00511 
00512   delete A;
00513   delete B;
00514   delete C;
00515   delete C_check;
00516 
00517   delete A_row_map;
00518   delete A_col_map;
00519   delete A_range_map;
00520   delete A_domain_map;
00521 
00522   delete B_row_map;
00523   delete B_col_map;
00524   delete B_range_map;
00525   delete B_domain_map;
00526 
00527   delete Cck_row_map;
00528   delete Cck_col_map;
00529   delete Cck_range_map;
00530   delete Cck_domain_map;
00531 
00532   return(return_code);
00533 }
00534 
00535 int read_matrix_file_names(Epetra_Comm& Comm,
00536                            const char* input_file_name,
00537                            char*& A_file,
00538                            bool& transA,
00539                            char*& B_file,
00540                            bool& transB,
00541                            char*& C_file)
00542 {
00543   int pathlen = path!=0 ? (int)strlen(path) : 0;
00544 
00545   if (Comm.MyPID()==0) {
00546     std::ifstream infile(input_file_name);
00547     if (!infile) {
00548       std::cout << "error opening input file " << input_file_name << std::endl;
00549       return(-1);
00550     }
00551 
00552     int linelen = 512;
00553     char line[512];
00554 
00555     infile.getline(line, linelen);
00556     if (!infile.eof()) {
00557       if (strchr(line, '#') == NULL) {
00558         A_file = new char[pathlen+strlen(line)+2];
00559         sprintf(A_file, "%s/%s",path,line);
00560       }
00561     }
00562 
00563     infile.getline(line, linelen);
00564     if (!infile.eof()) {
00565       if (!strcmp(line, "TRANSPOSE")) {
00566         transA = true;
00567       }
00568       else transA = false;
00569     }
00570 
00571     infile.getline(line, linelen);
00572     if (!infile.eof()) {
00573       if (strchr(line, '#') == NULL) {
00574         B_file = new char[pathlen+strlen(line)+2];
00575         sprintf(B_file, "%s/%s",path,line);
00576       }
00577     }
00578 
00579     infile.getline(line, linelen);
00580     if (!infile.eof()) {
00581       if (!strcmp(line, "TRANSPOSE")) {
00582         transB = true;
00583       }
00584       else transB = false;
00585     }
00586 
00587     infile.getline(line, linelen);
00588     if (!infile.eof()) {
00589       if (strchr(line, '#') == NULL) {
00590         C_file = new char[pathlen+strlen(line)+2];
00591         sprintf(C_file, "%s/%s", path, line);
00592       }
00593     }
00594 
00595     broadcast_name(Comm, (const char*&)A_file);
00596     broadcast_name(Comm, (const char*&)B_file);
00597     broadcast_name(Comm, (const char*&)C_file);
00598 #ifdef EPETRA_MPI
00599     int len = transA ? 1 : 0;
00600     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00601     len = transB ? 1 : 0;
00602     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00603 #endif
00604   }
00605   else {
00606     broadcast_name(Comm, (const char*&)A_file);
00607     broadcast_name(Comm, (const char*&)B_file);
00608     broadcast_name(Comm, (const char*&)C_file);
00609 #ifdef EPETRA_MPI
00610     int len = 0;
00611     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00612     transA = len==1 ? true : false;
00613     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00614     transB = len==1 ? true : false;
00615 #endif
00616   }
00617 
00618   return(0);
00619 }
00620 
00621 int create_maps(Epetra_Comm& Comm,
00622                 const char* input_file_name,
00623                 Epetra_Map*& row_map,
00624                 Epetra_Map*& col_map,
00625                 Epetra_Map*& range_map,
00626                 Epetra_Map*& domain_map)
00627 {
00628   return( EpetraExt::MatrixMarketFileToBlockMaps(input_file_name,
00629                                          Comm,
00630                                          (Epetra_BlockMap*&)row_map,
00631                                          (Epetra_BlockMap*&)col_map,
00632                                          (Epetra_BlockMap*&)range_map,
00633                                          (Epetra_BlockMap*&)domain_map)
00634   );
00635 }
00636 
00637 int read_matrix(const char* filename,
00638                 Epetra_Comm& Comm,
00639                 const Epetra_Map* rowmap,
00640                 Epetra_Map* colmap,
00641                 const Epetra_Map* rangemap,
00642                 const Epetra_Map* domainmap,
00643                 Epetra_CrsMatrix*& mat)
00644 {
00645   (void)Comm;
00646   int err = EpetraExt::MatrixMarketFileToCrsMatrix(filename, *rowmap, *colmap,
00647                                                    *rangemap, *domainmap, mat);
00648 
00649   return(err);
00650 }
00651 
00652 int two_proc_test(Epetra_Comm& Comm,
00653                   bool verbose)
00654 {
00655   (void)verbose;
00656   int thisproc = Comm.MyPID();
00657   int numprocs = Comm.NumProc();
00658 
00659   //only run this test on 2 procs
00660   if (numprocs != 2) return(0);
00661 
00662   //set up a row-std::map with 2 global elements,
00663   //1 on each proc.
00664   int numGlobalRows = 2;
00665   int numMyRows = 1;
00666   int myrow = 3;
00667   if (thisproc == 1) myrow = 7;
00668   Epetra_Map rowmap(numGlobalRows, numMyRows, &myrow, 0, Comm);
00669 
00670   //set up a domain-std::map with columns 0 - 4 on proc 0,
00671   //and columns 5 - 9 on proc 1.
00672   int numGlobalCols = 10;
00673   int numMyCols = 5;
00674   int* mycols = new int[numGlobalCols];
00675   int i;
00676   for(i=0; i<numGlobalCols; ++i) {
00677     mycols[i] = i;
00678   }
00679 
00680   Epetra_Map domainmap(numGlobalCols, numMyCols, &(mycols[thisproc*numMyCols]),
00681                        0, Comm);
00682 
00683   //now create matrices A, B and C with rowmap.
00684   Epetra_CrsMatrix A(Copy, rowmap, 10);
00685   Epetra_CrsMatrix B(Copy, rowmap, 10);
00686   Epetra_CrsMatrix C(Copy, rowmap, 10);
00687 
00688   double* coefs = new double[numGlobalCols];
00689   for(i=0; i<numGlobalCols; ++i) {
00690     coefs[i] = 1.0*i;
00691   }
00692 
00693   int err = A.InsertGlobalValues(myrow, numGlobalCols, coefs, mycols);
00694 
00695   err += B.InsertGlobalValues(myrow, numMyCols, &(coefs[thisproc*numMyCols]),
00696                        &(mycols[thisproc*numMyCols]));
00697 
00698   err += A.FillComplete(domainmap, rowmap);
00699   err += B.FillComplete(domainmap, rowmap);
00700 
00701   err += EpetraExt::MatrixMatrix::Multiply(A, false, B, true, C);
00702 
00703   //std::cout << "two_proc_test, A: "<<std::endl;
00704   //std::cout << A << std::endl;
00705 
00706   //std::cout << "two_proc_test, B: "<<std::endl;
00707   //std::cout << B << std::endl;
00708 
00709   //std::cout << "two_proc_test, C: "<<std::endl;
00710   //std::cout << C << std::endl;
00711 
00712   if (C.NumGlobalNonzeros() != 4) {
00713     err += 1;
00714   }
00715 
00716   delete [] coefs;
00717   delete [] mycols;
00718 
00719   return(err);
00720 }
00721 
00722 int time_matrix_matrix_multiply(Epetra_Comm& Comm, bool verbose)
00723 {
00724 
00725   const int magic_num = 3000;
00726   // 2009/02/23: rabartl: If you are going to do a timing test you need to
00727   // make this number adjustable form the command-line and you need to put in
00728   // a real test that compares against hard numbers for pass/fail.
00729 
00730   int localn = magic_num/Comm.NumProc();
00731 
00732   Epetra_CrsMatrix* A = create_epetra_crsmatrix(Comm.NumProc(),
00733                                                 Comm.MyPID(),
00734                                                 localn);
00735 
00736   Epetra_CrsMatrix* C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
00737 
00738   Epetra_Time timer(Comm);
00739   double start_time = timer.WallTime();
00740 
00741   int err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, false, *C);
00742 
00743   int globaln = localn*Comm.NumProc();
00744   if (verbose) {
00745     std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A, time: "
00746        << timer.WallTime()-start_time << std::endl;
00747   }
00748 
00749   C->FillComplete();
00750 
00751   start_time = timer.WallTime();
00752 
00753   err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, false, *C);
00754 
00755   if (verbose) {
00756     std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A, time: "
00757        << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
00758   }
00759 
00760   delete C;
00761 
00762   C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
00763 
00764   start_time = timer.WallTime();
00765 
00766   err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, true, *C);
00767 
00768   if (verbose) {
00769     std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A^T, time: "
00770        << timer.WallTime()-start_time << std::endl;
00771   }
00772 
00773   C->FillComplete();
00774 
00775   start_time = timer.WallTime();
00776 
00777   err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, true, *C);
00778 
00779   if (verbose) {
00780     std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A^T, time: "
00781        << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
00782   }
00783 
00784   delete C;
00785 
00786   C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
00787 
00788   start_time = timer.WallTime();
00789 
00790   err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, false, *C);
00791 
00792   if (verbose) {
00793     std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A, time: "
00794        << timer.WallTime()-start_time << std::endl;
00795   }
00796 
00797   C->FillComplete();
00798 
00799   start_time = timer.WallTime();
00800 
00801   err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, false, *C);
00802 
00803   if (verbose) {
00804     std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A, time: "
00805        << timer.WallTime()-start_time << " (C already Filled)\n"<<std::endl;
00806   }
00807 
00808   delete C;
00809 
00810   C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
00811 
00812   start_time = timer.WallTime();
00813 
00814   err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, true, *C);
00815 
00816   if (verbose) {
00817     std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A^T, time: "
00818        << timer.WallTime()-start_time << std::endl;
00819   }
00820 
00821   C->FillComplete();
00822 
00823   start_time = timer.WallTime();
00824 
00825   err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, true, *C);
00826 
00827   if (verbose) {
00828     std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A^T, time: "
00829        << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
00830   }
00831 
00832   delete C;
00833 
00834   delete A;
00835 
00836   return(err);
00837 }
00838 
00839 Epetra_CrsMatrix* create_epetra_crsmatrix(int numProcs,
00840                                           int localProc,
00841                                           int local_n,
00842                                           bool callFillComplete,
00843                                           bool symmetric)
00844 {
00845   (void)localProc;
00846 #ifdef EPETRA_MPI
00847   Epetra_MpiComm comm(MPI_COMM_WORLD);
00848 #else
00849   Epetra_SerialComm comm;
00850 #endif
00851   int global_num_rows = numProcs*local_n;
00852 
00853   Epetra_Map rowmap(global_num_rows, local_n, 0, comm);
00854 
00855   int nnz_per_row = 9;
00856   Epetra_CrsMatrix* matrix =
00857     new Epetra_CrsMatrix(Copy, rowmap, nnz_per_row);
00858 
00859   // Add  rows one-at-a-time
00860   double negOne = -1.0;
00861   double posTwo = 2.0;
00862   double val_L = symmetric ? negOne : -0.5;
00863 
00864   for (int i=0; i<local_n; i++) {
00865     int GlobalRow = matrix->GRID(i);
00866     int RowLess1 = GlobalRow - 1;
00867     int RowPlus1 = GlobalRow + 1;
00868     int RowLess5 = GlobalRow - 5;
00869     int RowPlus5 = GlobalRow + 5;
00870     int RowLess9 = GlobalRow - 9;
00871     int RowPlus9 = GlobalRow + 9;
00872     int RowLess24 = GlobalRow - 24;
00873     int RowPlus24 = GlobalRow + 24;
00874     int RowLess48 = GlobalRow - 48;
00875     int RowPlus48 = GlobalRow + 48;
00876 
00877 //    if (!symmetric) RowLess5 -= 2;
00878 
00879     if (RowLess48>=0) {
00880       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess48);
00881     }
00882     if (RowLess24>=0) {
00883       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess24);
00884     }
00885     if (RowLess9>=0) {
00886       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess9);
00887     }
00888     if (RowLess5>=0) {
00889       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess5);
00890     }
00891     if (RowLess1>=0) {
00892       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess1);
00893     }
00894     if (RowPlus1<global_num_rows) {
00895       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
00896     }
00897     if (RowPlus5<global_num_rows) {
00898       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus5);
00899     }
00900     if (RowPlus9<global_num_rows) {
00901       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus9);
00902     }
00903     if (RowPlus24<global_num_rows) {
00904       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus24);
00905     }
00906     if (RowPlus48<global_num_rows) {
00907       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus48);
00908     }
00909 
00910     matrix->InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
00911   }
00912 
00913   if (callFillComplete) {
00914     int err = matrix->FillComplete();
00915     if (err != 0) {
00916       std::cout << "create_epetra_matrix: error in matrix.FillComplete()"
00917          <<std::endl;
00918     }
00919   }
00920 
00921   return(matrix);
00922 }
00923 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines