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