EpetraExt Package Browser (Single Doxygen Collection) Development
test/MatrixMatrix/cxx_main.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00006 //                 Copyright (2011) Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //@HEADER
00042 */
00043 
00044 
00045 #include <string>
00046 #include <vector>
00047 #include <iostream>
00048 #include <fstream>
00049 
00050 #include <Epetra_ConfigDefs.h>
00051 
00052 #ifdef EPETRA_MPI
00053 #include <mpi.h>
00054 #include <Epetra_MpiComm.h>
00055 #endif
00056 
00057 #include <Epetra_SerialComm.h>
00058 #include <Epetra_Time.h>
00059 #include <Epetra_Import.h>
00060 #include <Epetra_Export.h>
00061 #include <Epetra_Map.h>
00062 #include <Epetra_LocalMap.h>
00063 #include <Epetra_CrsGraph.h>
00064 #include <Epetra_CrsMatrix.h>
00065 #include <Epetra_Vector.h>
00066 #include <EpetraExt_MMHelpers.h>
00067 #include <EpetraExt_MatrixMatrix.h>
00068 
00069 #include <EpetraExt_BlockMapIn.h>
00070 #include <EpetraExt_CrsMatrixIn.h>
00071 #include <EpetraExt_RowMatrixOut.h>
00072 
00073 namespace EpetraExt {
00074 extern
00075 Epetra_Map* find_rows_containing_cols(const Epetra_CrsMatrix& M,
00076                                       const Epetra_Map& column_map);
00077 }
00078 
00079 int read_input_file(Epetra_Comm& Comm,
00080                     const std::string& input_file_name,
00081                     std::vector<std::string>& filenames);
00082 
00083 int read_matrix_file_names(Epetra_Comm& Comm,
00084                            const std::string& input_file_name,
00085                            std::string& A_file,
00086                            bool& transA,
00087                            std::string& B_file,
00088                            bool& transB,
00089                            std::string& C_file);
00090 
00091 int broadcast_name(Epetra_Comm& Comm, std::string& name);
00092 
00093 int create_maps(Epetra_Comm& Comm,
00094                 const std::string& input_file_name,
00095                 Epetra_Map*& row_map,
00096                 Epetra_Map*& col_map,
00097                 Epetra_Map*& range_map,
00098                 Epetra_Map*& domain_map);
00099 
00100 int read_matrix(const std::string& filename,
00101                 Epetra_Comm& Comm,
00102                 const Epetra_Map* rowmap,
00103                 Epetra_Map* colmap,
00104                 const Epetra_Map* rangemap,
00105                 const Epetra_Map* domainmap,
00106                 Epetra_CrsMatrix*& mat);
00107 
00108 int run_test(Epetra_Comm& Comm,
00109              const std::string& filename,
00110        bool do_FillComplete,
00111              bool result_mtx_to_file=false,
00112              bool verbose=false);
00113 
00114 int two_proc_test(Epetra_Comm& Comm,
00115                   bool verbose=false);
00116 
00117 int test_find_rows(Epetra_Comm& Comm);
00118 
00119 Epetra_CrsMatrix* create_epetra_crsmatrix(int numProcs,
00120                                           int localProc,
00121                                           int local_n,
00122                                           bool callFillComplete = true,
00123                                           bool symmetric = true);
00124 
00125 int time_matrix_matrix_multiply(Epetra_Comm& Comm,
00126                                 bool verbose);
00127 
00128 int test_drumm1(Epetra_Comm& Comm);
00129 
00131 //Global variable!!!!
00132 std::string path;
00134 
00135 int main(int argc, char** argv) {
00136 
00137 #ifdef EPETRA_MPI
00138   MPI_Init(&argc,&argv);
00139   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00140 #else
00141   Epetra_SerialComm Comm;
00142 #endif
00143 
00144   bool write_result_mtx = false;
00145   bool verbose = false;
00146   int write = 0, inp_specified = 0;
00147   bool path_specified = false;
00148   std::string input_file;
00149   bool input_file_specified = false;
00150 
00151   if (Comm.MyPID()==0) {
00152     for(int ii=0; ii<argc; ++ii) {
00153       if (!strcmp("-write_result", argv[ii])) write_result_mtx = true;
00154       if (!strcmp("-v", argv[ii])) verbose = true;
00155       if (!strcmp("-i", argv[ii])) {
00156         input_file = argv[ii+1];
00157         input_file_specified = true;
00158       }
00159       if (!strcmp("-d",argv[ii])) {
00160         path = argv[ii+1];
00161         path_specified = true;
00162       }
00163     }
00164     write = write_result_mtx ? 1 : 0;
00165     inp_specified = input_file_specified ? 1 : 0;
00166   }
00167 #ifdef EPETRA_MPI
00168   MPI_Bcast(&write, 1, MPI_INT, 0, MPI_COMM_WORLD);
00169   if (write) write_result_mtx = true;
00170   MPI_Bcast(&inp_specified, 1, MPI_INT, 0, MPI_COMM_WORLD);
00171   if (inp_specified) input_file_specified = true;
00172 #endif
00173 
00174   if (!path_specified) {
00175     path = ".";
00176   }
00177 
00178   int err = two_proc_test(Comm, verbose);
00179   if (err != 0) {
00180     std::cout << "two_proc_test returned err=="<<err<<std::endl;
00181     return(err);
00182   }
00183 
00184   std::vector<std::string> filenames;
00185 
00186   if (input_file_specified) {
00187     filenames.push_back(std::string(input_file));
00188   }
00189   else {
00190     input_file = "infiles";
00191     std::cout << "specific input-file not specified, getting list of input-files from file 'infiles'." << std::endl;
00192 
00193     err = read_input_file(Comm, input_file, filenames);
00194     if (err != 0) {
00195       if (path_specified) path_specified = false;
00196       path = "./MatrixMatrix";
00197       read_input_file(Comm, input_file, filenames);
00198     }
00199   }
00200 
00201   err = test_find_rows(Comm);
00202   if (err != 0) {
00203     std::cout << "test_find_rows returned err=="<<err<<std::endl;
00204     return(err);
00205   }
00206 
00207   err = test_drumm1(Comm);
00208   if (err != 0) {
00209     std::cout << "test_drumm1 returned err=="<<err<<std::endl;
00210     return(err);
00211   }
00212 
00213   for(size_t i=0; i<filenames.size(); ++i) {
00214     err = run_test(Comm, filenames[i], true, write_result_mtx, verbose);
00215     if (err != 0) break;
00216     err = run_test(Comm, filenames[i], false, write_result_mtx, verbose);
00217     if (err != 0) break;
00218   }
00219 
00220   Epetra_CrsMatrix* D = create_epetra_crsmatrix(Comm.NumProc(),
00221                                                 Comm.MyPID(), 10,
00222                                                 true, false);
00223 
00224 //  std::cout << "D: \n"  << *D << std::endl;
00225 
00226   EpetraExt::MatrixMatrix::Add(*D, true, 0.5, *D, 0.5);
00227 
00228 //  std::cout << "symm D: \n"  << *D << std::endl;
00229 
00230   delete D;
00231 
00232   if (err == 0) {
00233     err = time_matrix_matrix_multiply(Comm, verbose);
00234   }
00235 
00236   int global_err = err;
00237 
00238 #ifdef EPETRA_MPI
00239   MPI_Allreduce(&err, &global_err, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
00240   MPI_Finalize();
00241 #endif
00242 
00243   return(global_err);
00244 }
00245 
00246 int test_find_rows(Epetra_Comm& Comm)
00247 {
00248   int numprocs = Comm.NumProc();
00249   int localproc = Comm.MyPID();
00250   int numlocalrows = 2;
00251   int numglobalrows = numprocs*numlocalrows;
00252   Epetra_Map rowmap(numlocalrows*numprocs, 0, Comm);
00253   Epetra_CrsMatrix matrix(Copy, rowmap, numglobalrows);
00254 
00255   int err = 0;
00256   int* cols = new int[numglobalrows];
00257   double*vals = new double[numglobalrows];
00258 
00259   for(int j=0; j<numglobalrows; ++j) {
00260     cols[j] = j;
00261     vals[j] = 1.0;
00262   }
00263 
00264   Epetra_Map colmap(-1, numglobalrows, cols, 0, Comm);
00265 
00266   for(int i=0; i<numlocalrows; ++i) {
00267     int row = localproc*numlocalrows+i;
00268     err = matrix.InsertGlobalValues(row, 1, &(vals[i]), &row);
00269     if (err != 0) {
00270       return(err);
00271     }
00272   }
00273 
00274   err = matrix.FillComplete();
00275   if (err != 0) {
00276     return(err);
00277   }
00278 
00279   Epetra_Map* map_rows = EpetraExt::find_rows_containing_cols(matrix, colmap);
00280 
00281   if (map_rows->NumMyElements() != numglobalrows) {
00282     return(-1);
00283   }
00284 
00285   delete map_rows;
00286   delete [] cols;
00287   delete [] vals;
00288 
00289   return(0);
00290 }
00291 
00292 int expand_name_list(const char* newname,
00293                      const char**& names,
00294                      int& alloc_len,
00295                      int& num_names)
00296 {
00297   int offset = num_names;
00298   if (offset >= alloc_len) {
00299     int alloc_increment = 8;
00300     const char** newlist = new const char*[alloc_len+alloc_increment];
00301     for(int i=0; i<offset; ++i) {
00302       newlist[i] = names[i];
00303     }
00304     delete [] names;
00305     names = newlist;
00306     alloc_len += alloc_increment;
00307     for(int i=offset; i<alloc_len; ++i) {
00308       names[i] = NULL;
00309     }
00310   }
00311 
00312   names[offset] = newname;
00313   ++num_names;
00314   return(0);
00315 }
00316 
00317 int broadcast_name(Epetra_Comm& Comm, std::string& name)
00318 {
00319   if (Comm.NumProc() < 2) return(0);
00320 
00321 #ifdef EPETRA_MPI
00322   int len = name.size();
00323   int localProc = Comm.MyPID();
00324   if (localProc == 0) {
00325     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00326     MPI_Bcast((void*)name.c_str(), len+1, MPI_CHAR, 0, MPI_COMM_WORLD);
00327   }
00328   else {
00329     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00330     name.resize(len+1, ' ');
00331     MPI_Bcast((void*)name.c_str(), len+1, MPI_CHAR, 0, MPI_COMM_WORLD);
00332   }
00333 
00334 #endif
00335   return(0);
00336 }
00337 
00338 int read_input_file(Epetra_Comm& Comm,
00339                     const std::string& input_file_name,
00340                     std::vector<std::string>& filenames)
00341 {
00342   int local_err = 0, global_err = 0;
00343   std::ifstream* infile = NULL;
00344   int pathlen = path.size();
00345 
00346   if (Comm.MyPID() == 0) {
00347     std::string full_name = path.empty() ? input_file_name : path+"/"+input_file_name;
00348 
00349     infile = new std::ifstream(full_name.c_str());
00350     if (!(*infile)) {
00351       local_err = -1;
00352       delete infile;
00353     }
00354   }
00355 
00356   Comm.SumAll(&local_err, &global_err, 1);
00357 
00358   if (global_err != 0) {
00359     return(global_err);
00360   }
00361 
00362 
00363   if (Comm.MyPID() == 0) {
00364     const int linelen = 512;
00365     char line[linelen];
00366 
00367     std::ifstream& ifile = *infile;
00368     while(!ifile.eof()) {
00369       if (pathlen>0) {
00370         sprintf(line,"%s/",path.c_str());
00371         ifile.getline(&(line[pathlen+1]), linelen);
00372       }
00373       else {
00374         ifile.getline(line, linelen);
00375       }
00376 
00377       if (ifile.fail()) {
00378         break;
00379       }
00380       if (strchr(line, '#') == NULL) {
00381         filenames.push_back(std::string(line));
00382       }
00383     }
00384 
00385     int numfiles = filenames.size();
00386 #ifdef EPETRA_MPI
00387     MPI_Bcast(&numfiles, 1, MPI_INT, 0, MPI_COMM_WORLD);
00388 #endif
00389     for(int i=0; i<numfiles; ++i) {
00390       broadcast_name(Comm, filenames[i]);
00391     }
00392 
00393     delete infile;
00394   }
00395   else {
00396     int numfiles = 0;
00397 #ifdef EPETRA_MPI
00398     MPI_Bcast(&numfiles, 1, MPI_INT, 0, MPI_COMM_WORLD);
00399 #endif
00400     filenames.resize(numfiles);
00401     for(int i=0; i<numfiles; ++i) {
00402       broadcast_name(Comm, filenames[i]);
00403     }
00404   }
00405   
00406   return(0);
00407 }
00408 
00409 int run_test(Epetra_Comm& Comm,
00410              const std::string& filename,
00411        bool do_FillComplete,
00412              bool result_mtx_to_file,
00413              bool verbose)
00414 {
00415   std::string A_file;
00416   char AT[3]; AT[0] = '^'; AT[1] = 'T'; AT[2] = '\0';
00417   std::string B_file;
00418   char BT[3]; BT[0] = '^'; BT[1] = 'T'; BT[2] = '\0';
00419   std::string C_file;
00420   bool transA, transB;
00421 
00422   if(!Comm.MyPID()) std::cout<<"Testing: "<<filename<<std::endl;
00423 
00424   int err = read_matrix_file_names(Comm, filename, A_file, transA,
00425                                    B_file, transB, C_file);
00426   if (err != 0) {
00427     std::cout << "Error, read_matrix_file_names returned " << err << std::endl;
00428     return(err);
00429   }
00430 
00431   if (!transA) AT[0] = '\0';
00432   if (!transB) BT[0] = '\0';
00433 
00434   int localProc = Comm.MyPID();
00435 
00436   if (localProc == 0 && verbose) {
00437     std::cout << "Testing C=A"<<AT<<"*B"<<BT<< "; A:" << A_file
00438               << ", B:" << B_file << ", C:" << C_file << std::endl;
00439   }
00440 
00441   Epetra_CrsMatrix* A = NULL;
00442   Epetra_CrsMatrix* B = NULL;
00443   Epetra_CrsMatrix* C = NULL;
00444   Epetra_CrsMatrix* C_check = NULL;
00445 
00446   Epetra_Map* A_row_map = NULL;
00447   Epetra_Map* A_col_map = NULL;
00448   Epetra_Map* A_range_map = NULL;
00449   Epetra_Map* A_domain_map = NULL;
00450   err = create_maps(Comm, A_file, A_row_map, A_col_map, A_range_map, A_domain_map);
00451   if (err != 0) {
00452     std::cout << "create_maps A returned err=="<<err<<std::endl;
00453     return(err);
00454   }
00455 
00456   Epetra_Map* B_row_map = NULL;
00457   Epetra_Map* B_col_map = NULL;
00458   Epetra_Map* B_range_map = NULL;
00459   Epetra_Map* B_domain_map = NULL;
00460   err = create_maps(Comm, B_file, B_row_map, B_col_map, B_range_map, B_domain_map);
00461   if (err != 0) {
00462     std::cout << "create_maps B returned err=="<<err<<std::endl;
00463     return(err);
00464   }
00465 
00466   err = read_matrix(A_file, Comm, A_row_map, A_col_map,
00467                     A_range_map, A_domain_map, A);
00468   if (err != 0) {
00469     std::cout << "read_matrix A returned err=="<<err<<std::endl;
00470     return(err);
00471   }
00472 
00473   err = read_matrix(B_file, Comm, B_row_map, B_col_map,
00474                     B_range_map, B_domain_map, B);
00475   if (err != 0) {
00476     std::cout << "read_matrix B returned err=="<<err<<std::endl;
00477     return(-1);
00478   }
00479 
00480   const Epetra_Map* rowmap    = transA ? &(A->DomainMap()) : &(A->RowMap());
00481   const Epetra_Map* domainMap = transB ? &(B->RangeMap()) : &(B->DomainMap());
00482   const Epetra_Map* rangeMap  = transA ? &(A->DomainMap()) : &(A->RangeMap());
00483 
00484 
00485   C = new Epetra_CrsMatrix(Copy, *rowmap, 1);
00486 
00487   if(C->Comm().MyPID()) printf("transA = %d transB = %d\n",(int)transA,(int)transB);
00488 
00489   err = EpetraExt::MatrixMatrix::Multiply(*A, transA, *B, transB, *C,do_FillComplete);
00490   if (err != 0) {
00491     std::cout << "err "<<err<<" from MatrixMatrix::Multiply"<<std::endl;
00492     return(err);
00493   }
00494 
00495   if(!do_FillComplete) C->FillComplete(*domainMap,*rangeMap);
00496 
00497 
00498 //  std::cout << "A: " << *A << std::endl << "B: "<<*B<<std::endl<<"C: "<<*C<<std::endl;
00499   if (result_mtx_to_file) {
00500     EpetraExt::RowMatrixToMatrixMarketFile("result.mtx", *C);
00501   }
00502 
00503   Epetra_Map* Cck_row_map = NULL;
00504   Epetra_Map* Cck_col_map = NULL;
00505   Epetra_Map* Cck_range_map = NULL;
00506   Epetra_Map* Cck_domain_map = NULL;
00507   err = create_maps(Comm, C_file, Cck_row_map, Cck_col_map,
00508                     Cck_range_map, Cck_domain_map);
00509   if (err != 0) {
00510     std::cout << "create_maps C returned err=="<<err<<std::endl;
00511     return(err);
00512   }
00513 
00514   err = read_matrix(C_file, Comm, Cck_row_map, Cck_col_map,
00515                      Cck_range_map, Cck_domain_map, C_check);
00516   if (err != 0) {
00517     std::cout << "read_matrix C returned err=="<<err<<std::endl;
00518     return(-1);
00519   }
00520 
00521   //now subtract our check-matrix from our result matrix (C) and
00522   //the infinity-norm of that should be zero.
00523   EpetraExt::MatrixMatrix::Add(*C, false, -1.0, *C_check, 1.0);
00524 
00525   double inf_norm = C_check->NormInf();
00526 
00527   int return_code = 0;
00528 
00529   if (inf_norm < 1.e-13) {
00530     if (localProc == 0 && verbose) {
00531       std::cout << "Test Passed" << std::endl;
00532     }
00533   }
00534   else {
00535     return_code = -1;
00536     if (localProc == 0) {
00537       std::cout << "Test Failed ("<<filename<<"), inf_norm = " << inf_norm << std::endl;
00538     }
00539 Comm.Barrier();
00540 std::cout << "C"<<std::endl;
00541 std::cout << *C<<std::endl;
00542   }
00543 
00544   // Test Matrix Jacobi for non-transpose matrices
00545   if(!transA && !transB && A->RowMap().SameAs(B->RowMap())) {
00546     if(!Comm.MyPID()) std::cout<<"--Testing: Jacobi for same matrices"<<std::endl;
00547     delete C; delete C_check;
00548 
00549     double omega=1.0;
00550     Epetra_Vector Dinv(B->RowMap());
00551     Dinv.PutScalar(1.0);
00552 
00553     // Jacobi version
00554     C = new Epetra_CrsMatrix(Copy,B->RowMap(),0);
00555     EpetraExt::MatrixMatrix::Jacobi(omega,Dinv,*A,*B,*C);
00556     
00557     // Multiply + Add version
00558     Dinv.PutScalar(omega);
00559     Epetra_CrsMatrix * AB = new Epetra_CrsMatrix(Copy,B->RowMap(),0);
00560     C_check = new Epetra_CrsMatrix(Copy,B->RowMap(),0);
00561     EpetraExt::MatrixMatrix::Multiply(*A,false,*B,false,*AB);
00562     AB->LeftScale(Dinv);
00563     EpetraExt::MatrixMatrix::Add(*AB,false,-1.0,*B,false,1.0,C_check);
00564 
00565     // Check the difference
00566     EpetraExt::MatrixMatrix::Add(*C, false, -1.0, *C_check, 1.0);    
00567     C_check->FillComplete(B->DomainMap(),B->RangeMap());
00568     
00569     // Error check
00570     inf_norm = C_check->NormInf();
00571     return_code = 0;    
00572     if (inf_norm < 1.e-13) {
00573       if (localProc == 0 && verbose) {
00574   std::cout << "Jacobi Test Passed" << std::endl;
00575       }
00576     }
00577     else {
00578       return_code = -1;
00579       if (localProc == 0) {
00580   std::cout << "Jacobi Test Failed ("<<filename<<"), inf_norm = " << inf_norm << std::endl;
00581     }
00582       Comm.Barrier();
00583       std::cout << "C"<<std::endl;
00584       std::cout << *C<<std::endl;
00585     }
00586        
00587     delete AB;
00588   }
00589 
00590 
00591   delete A;
00592   delete B;
00593   delete C;
00594   delete C_check;
00595 
00596   delete A_row_map;
00597   delete A_col_map;
00598   delete A_range_map;
00599   delete A_domain_map;
00600 
00601   delete B_row_map;
00602   delete B_col_map;
00603   delete B_range_map;
00604   delete B_domain_map;
00605 
00606   delete Cck_row_map;
00607   delete Cck_col_map;
00608   delete Cck_range_map;
00609   delete Cck_domain_map;
00610 
00611   return(return_code);
00612 }
00613 
00614 int read_matrix_file_names(Epetra_Comm& Comm,
00615                            const std::string& input_file_name,
00616                            std::string& A_file,
00617                            bool& transA,
00618                            std::string& B_file,
00619                            bool& transB,
00620                            std::string& C_file)
00621 {
00622   if (Comm.MyPID()==0) {
00623     std::ifstream infile(input_file_name.c_str());
00624     if (!infile) {
00625       std::cout << "error opening input file " << input_file_name << std::endl;
00626       return(-1);
00627     }
00628 
00629     const int linelen = 512;
00630     char line[linelen];
00631 
00632     infile.getline(line, linelen);
00633     if (!infile.eof()) {
00634       if (strchr(line, '#') == NULL) {
00635         A_file = path+"/"+line;
00636       }
00637     }
00638 
00639     infile.getline(line, linelen);
00640     if (!infile.eof()) {
00641       if (!strcmp(line, "TRANSPOSE")) {
00642         transA = true;
00643       }
00644       else transA = false;
00645     }
00646 
00647     infile.getline(line, linelen);
00648     if (!infile.eof()) {
00649       if (strchr(line, '#') == NULL) {
00650         B_file = path+"/"+line;
00651       }
00652     }
00653 
00654     infile.getline(line, linelen);
00655     if (!infile.eof()) {
00656       if (!strcmp(line, "TRANSPOSE")) {
00657         transB = true;
00658       }
00659       else transB = false;
00660     }
00661 
00662     infile.getline(line, linelen);
00663     if (!infile.eof()) {
00664       if (strchr(line, '#') == NULL) {
00665         C_file = path+"/"+line;
00666       }
00667     }
00668 
00669     broadcast_name(Comm, A_file);
00670     broadcast_name(Comm, B_file);
00671     broadcast_name(Comm, C_file);
00672 #ifdef EPETRA_MPI
00673     int len = transA ? 1 : 0;
00674     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00675     len = transB ? 1 : 0;
00676     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00677 #endif
00678   }
00679   else {
00680     broadcast_name(Comm, A_file);
00681     broadcast_name(Comm, B_file);
00682     broadcast_name(Comm, C_file);
00683 #ifdef EPETRA_MPI
00684     int len = 0;
00685     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00686     transA = len==1 ? true : false;
00687     MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
00688     transB = len==1 ? true : false;
00689 #endif
00690   }
00691 
00692   return(0);
00693 }
00694 
00695 int create_maps(Epetra_Comm& Comm,
00696                 const std::string& input_file_name,
00697                 Epetra_Map*& row_map,
00698                 Epetra_Map*& col_map,
00699                 Epetra_Map*& range_map,
00700                 Epetra_Map*& domain_map)
00701 {
00702   return( EpetraExt::MatrixMarketFileToBlockMaps(input_file_name.c_str(),
00703                                          Comm,
00704                                          (Epetra_BlockMap*&)row_map,
00705                                          (Epetra_BlockMap*&)col_map,
00706                                          (Epetra_BlockMap*&)range_map,
00707                                          (Epetra_BlockMap*&)domain_map)
00708   );
00709 }
00710 
00711 int read_matrix(const std::string& filename,
00712                 Epetra_Comm& Comm,
00713                 const Epetra_Map* rowmap,
00714                 Epetra_Map* colmap,
00715                 const Epetra_Map* rangemap,
00716                 const Epetra_Map* domainmap,
00717                 Epetra_CrsMatrix*& mat)
00718 {
00719   (void)Comm;
00720   int err = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), *rowmap,
00721                                                    *rangemap, *domainmap, mat);
00722 
00723   return(err);
00724 }
00725 
00726 int two_proc_test(Epetra_Comm& Comm,
00727                   bool verbose)
00728 {
00729   (void)verbose;
00730   int thisproc = Comm.MyPID();
00731   int numprocs = Comm.NumProc();
00732 
00733   //only run this test on 2 procs
00734   if (numprocs != 2) return(0);
00735 
00736   //set up a row-std::map with 2 global elements,
00737   //1 on each proc.
00738   int numGlobalRows = 2;
00739   int numMyRows = 1;
00740   int myrow = 3;
00741   if (thisproc == 1) myrow = 7;
00742   Epetra_Map rowmap(numGlobalRows, numMyRows, &myrow, 0, Comm);
00743 
00744   //set up a domain-map with columns 0 - 4 on proc 0,
00745   //and columns 5 - 9 on proc 1.
00746   int numGlobalCols = 10;
00747   int numMyCols = 5;
00748   int* mycols = new int[numGlobalCols];
00749   int i;
00750   for(i=0; i<numGlobalCols; ++i) {
00751     mycols[i] = i;
00752   }
00753 
00754   Epetra_Map domainmap(numGlobalCols, numMyCols, &(mycols[thisproc*numMyCols]),
00755                        0, Comm);
00756 
00757   //now create matrices A, B and C with rowmap.
00758   Epetra_CrsMatrix A(Copy, rowmap, 10);
00759   Epetra_CrsMatrix B(Copy, rowmap, 10);
00760   Epetra_CrsMatrix C(Copy, rowmap, 10);
00761 
00762   double* coefs = new double[numGlobalCols];
00763   for(i=0; i<numGlobalCols; ++i) {
00764     coefs[i] = 1.0*i;
00765   }
00766 
00767   int numCols = numGlobalCols - 2;
00768   int offset = 0;
00769   if (thisproc == 1) offset = 2;
00770 
00771   int err = A.InsertGlobalValues(myrow, numCols, &coefs[offset], &mycols[offset]);
00772 
00773   err += B.InsertGlobalValues(myrow, numMyCols, &(coefs[thisproc*numMyCols]),
00774                        &(mycols[thisproc*numMyCols]));
00775 
00776   err += A.FillComplete(domainmap, rowmap);
00777   err += B.FillComplete(domainmap, rowmap);
00778 
00779   err += EpetraExt::MatrixMatrix::Multiply(A, false, B, true, C);
00780 
00781   //std::cout << "two_proc_test, A: "<<std::endl;
00782   //std::cout << A << std::endl;
00783 
00784   //std::cout << "two_proc_test, B: "<<std::endl;
00785   //std::cout << B << std::endl;
00786 
00787   //std::cout << "two_proc_test, C: "<<std::endl;
00788   //std::cout << C << std::endl;
00789 
00790   if (C.NumGlobalNonzeros() != 4) {
00791     err += 1;
00792   }
00793 
00794   delete [] coefs;
00795   delete [] mycols;
00796 
00797   return(err);
00798 }
00799 
00800 int time_matrix_matrix_multiply(Epetra_Comm& Comm, bool verbose)
00801 {
00802 
00803   const int magic_num = 3000;
00804   // 2009/02/23: rabartl: If you are going to do a timing test you need to
00805   // make this number adjustable form the command-line and you need to put in
00806   // a real test that compares against hard numbers for pass/fail.
00807 
00808   int localn = magic_num/Comm.NumProc();
00809 
00810   Epetra_CrsMatrix* A = create_epetra_crsmatrix(Comm.NumProc(),
00811                                                 Comm.MyPID(),
00812                                                 localn);
00813 
00814   Epetra_CrsMatrix* C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
00815 
00816   Epetra_Time timer(Comm);
00817   double start_time = timer.WallTime();
00818 
00819   int err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, false, *C);
00820 
00821   int globaln = localn*Comm.NumProc();
00822   if (verbose) {
00823     std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A, time: "
00824        << timer.WallTime()-start_time << std::endl;
00825   }
00826 
00827   C->FillComplete();
00828 
00829   start_time = timer.WallTime();
00830 
00831   err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, false, *C);
00832 
00833   if (verbose) {
00834     std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A, time: "
00835        << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
00836   }
00837 
00838   delete C;
00839 
00840   C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
00841 
00842   start_time = timer.WallTime();
00843 
00844   err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, true, *C);
00845 
00846   if (verbose) {
00847     std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A^T, time: "
00848        << timer.WallTime()-start_time << std::endl;
00849   }
00850 
00851   C->FillComplete();
00852 
00853   start_time = timer.WallTime();
00854 
00855   err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, true, *C);
00856 
00857   if (verbose) {
00858     std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A^T, time: "
00859        << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
00860   }
00861 
00862   delete C;
00863 
00864   C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
00865 
00866   start_time = timer.WallTime();
00867 
00868   err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, false, *C);
00869 
00870   if (verbose) {
00871     std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A, time: "
00872        << timer.WallTime()-start_time << std::endl;
00873   }
00874 
00875   C->FillComplete();
00876 
00877   start_time = timer.WallTime();
00878 
00879   err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, false, *C);
00880 
00881   if (verbose) {
00882     std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A, time: "
00883        << timer.WallTime()-start_time << " (C already Filled)\n"<<std::endl;
00884   }
00885 
00886   delete C;
00887 
00888   C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
00889 
00890   start_time = timer.WallTime();
00891 
00892   err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, true, *C);
00893 
00894   if (verbose) {
00895     std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A^T, time: "
00896        << timer.WallTime()-start_time << std::endl;
00897   }
00898 
00899   C->FillComplete();
00900 
00901   start_time = timer.WallTime();
00902 
00903   err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, true, *C);
00904 
00905   if (verbose) {
00906     std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A^T, time: "
00907        << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
00908   }
00909 
00910   delete C;
00911 
00912   delete A;
00913 
00914   return(err);
00915 }
00916 
00917 Epetra_CrsMatrix* create_epetra_crsmatrix(int numProcs,
00918                                           int localProc,
00919                                           int local_n,
00920                                           bool callFillComplete,
00921                                           bool symmetric)
00922 {
00923   (void)localProc;
00924 #ifdef EPETRA_MPI
00925   Epetra_MpiComm comm(MPI_COMM_WORLD);
00926 #else
00927   Epetra_SerialComm comm;
00928 #endif
00929   int global_num_rows = numProcs*local_n;
00930 
00931   Epetra_Map rowmap(global_num_rows, local_n, 0, comm);
00932 
00933   int nnz_per_row = 9;
00934   Epetra_CrsMatrix* matrix =
00935     new Epetra_CrsMatrix(Copy, rowmap, nnz_per_row);
00936 
00937   // Add  rows one-at-a-time
00938   double negOne = -1.0;
00939   double posTwo = 2.0;
00940   double val_L = symmetric ? negOne : -0.5;
00941 
00942   for (int i=0; i<local_n; i++) {
00943     int GlobalRow = matrix->GRID(i);
00944     int RowLess1 = GlobalRow - 1;
00945     int RowPlus1 = GlobalRow + 1;
00946     int RowLess5 = GlobalRow - 5;
00947     int RowPlus5 = GlobalRow + 5;
00948     int RowLess9 = GlobalRow - 9;
00949     int RowPlus9 = GlobalRow + 9;
00950     int RowLess24 = GlobalRow - 24;
00951     int RowPlus24 = GlobalRow + 24;
00952     int RowLess48 = GlobalRow - 48;
00953     int RowPlus48 = GlobalRow + 48;
00954 
00955 //    if (!symmetric) RowLess5 -= 2;
00956 
00957     if (RowLess48>=0) {
00958       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess48);
00959     }
00960     if (RowLess24>=0) {
00961       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess24);
00962     }
00963     if (RowLess9>=0) {
00964       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess9);
00965     }
00966     if (RowLess5>=0) {
00967       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess5);
00968     }
00969     if (RowLess1>=0) {
00970       matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess1);
00971     }
00972     if (RowPlus1<global_num_rows) {
00973       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
00974     }
00975     if (RowPlus5<global_num_rows) {
00976       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus5);
00977     }
00978     if (RowPlus9<global_num_rows) {
00979       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus9);
00980     }
00981     if (RowPlus24<global_num_rows) {
00982       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus24);
00983     }
00984     if (RowPlus48<global_num_rows) {
00985       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus48);
00986     }
00987 
00988     matrix->InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
00989   }
00990 
00991   if (callFillComplete) {
00992     int err = matrix->FillComplete();
00993     if (err != 0) {
00994       std::cout << "create_epetra_matrix: error in matrix.FillComplete()"
00995          <<std::endl;
00996     }
00997   }
00998 
00999   return(matrix);
01000 }
01001 
01002 int test_drumm1(Epetra_Comm& Comm)
01003 {
01004   int size = Comm.NumProc();
01005   if (size != 2) return 0;
01006 
01007   int rank = Comm.MyPID();
01008 
01009   int indexBase = 0;
01010   int numGlobalElements = 2;
01011   
01012   Epetra_Map emap(numGlobalElements, indexBase, Comm);
01013 
01014   Epetra_CrsMatrix A(Copy, emap, 0);
01015 
01016   // 2x2 matrix:
01017   //   3 4
01018   //   1 2
01019   std::vector<std::vector<double> > vals(numGlobalElements);
01020   vals[0].push_back(3); vals[0].push_back(4);
01021   vals[1].push_back(1); vals[1].push_back(2);
01022 
01023   std::vector<int> indices;
01024   indices.push_back(0); indices.push_back(1);
01025 
01026   for (int row=0; row<numGlobalElements; ++row) {
01027     if ( A.MyGRID(row) )
01028       A.InsertGlobalValues(row, numGlobalElements, &(vals[row][0]), &indices[0]);
01029   }
01030 
01031   A.FillComplete();
01032 
01033   Epetra_CrsMatrix B(Copy, emap, 0);
01034   EpetraExt::MatrixMatrix::Multiply(A, true, A, false, B);
01035 
01036   // B = Transpose(A) x A should be
01037   //  10 14
01038   //  14 20
01039   int idx[2];
01040   int tmp;
01041   double val[2];
01042 
01043   //for this little test, global_row == rank:
01044   B.ExtractGlobalRowCopy(rank, 2, tmp, val, idx);
01045 
01046   int test_result = 0;
01047 
01048   if (rank == 0) {
01049     if (idx[0] == 0 && val[0] != 10.0) test_result = 1;
01050     if (idx[1] == 0 && val[1] != 10.0) test_result = 1;
01051     if (idx[0] == 1 && val[0] != 14.0) test_result = 1;
01052     if (idx[1] == 1 && val[1] != 14.0) test_result = 1;
01053   }
01054   else {
01055     if (idx[0] == 0 && val[0] != 14.0) test_result = 1;
01056     if (idx[1] == 0 && val[1] != 14.0) test_result = 1;
01057     if (idx[0] == 1 && val[0] != 20.0) test_result = 1;
01058     if (idx[1] == 1 && val[1] != 20.0) test_result = 1;
01059   }
01060 
01061   int global_test_result = 0;
01062   Comm.SumAll(&test_result, &global_test_result, 1);
01063 
01064   return global_test_result;
01065 }
01066 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines