Zoltan2
UserInputForTests.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //   Zoltan2: A package of combinatorial algorithms for scientific computing
00006 //                  Copyright 2012 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 Karen Devine      (kddevin@sandia.gov)
00039 //                    Erik Boman        (egboman@sandia.gov)
00040 //                    Siva Rajamanickam (srajama@sandia.gov)
00041 //
00042 // ***********************************************************************
00043 //
00044 // @HEADER
00045 
00050 #ifndef USERINPUTFORTESTS
00051 #define USERINPUTFORTESTS
00052 
00053 #include <Zoltan2_XpetraTraits.hpp>
00054 #include <Zoltan2_ChacoReader.hpp>
00055 
00056 #include <Tpetra_MultiVector.hpp>
00057 #include <Xpetra_Vector.hpp>
00058 #include <Xpetra_CrsMatrix.hpp>
00059 #include <Xpetra_CrsGraph.hpp>
00060 
00061 #include <MatrixMarket_Tpetra.hpp>
00062 #include <Galeri_XpetraProblemFactory.hpp>
00063 #include <Galeri_XpetraParameters.hpp>
00064 
00065 #include <Kokkos_DefaultNode.hpp>
00066 
00067 //#include <Xpetra_EpetraUtils.hpp>
00068 #ifdef HAVE_ZOLTAN2_MPI
00069 #include <Epetra_MpiComm.h>
00070 #else
00071 #include <Epetra_SerialComm.h>
00072 #endif
00073 
00074 using Teuchos::RCP;
00075 using Teuchos::ArrayRCP;
00076 using Teuchos::ArrayView;
00077 using Teuchos::Array;
00078 using Teuchos::Comm;
00079 using Teuchos::rcp;
00080 using Teuchos::arcp;
00081 using Teuchos::rcp_const_cast;
00082 using std::string;
00083 
00114 class UserInputForTests
00115 {
00116 public:
00117 
00118   typedef Tpetra::CrsMatrix<zscalar_t, zlno_t, zgno_t, znode_t> tcrsMatrix_t;
00119   typedef Tpetra::CrsGraph<zlno_t, zgno_t, znode_t> tcrsGraph_t;
00120   typedef Tpetra::Vector<zscalar_t, zlno_t, zgno_t, znode_t> tVector_t;
00121   typedef Tpetra::MultiVector<zscalar_t, zlno_t, zgno_t, znode_t> tMVector_t;
00122 
00123   typedef Xpetra::CrsMatrix<zscalar_t, zlno_t, zgno_t, znode_t> xcrsMatrix_t;
00124   typedef Xpetra::CrsGraph<zlno_t, zgno_t, znode_t> xcrsGraph_t;
00125   typedef Xpetra::Vector<zscalar_t, zlno_t, zgno_t, znode_t> xVector_t;
00126   typedef Xpetra::MultiVector<zscalar_t, zlno_t, zgno_t, znode_t> xMVector_t;
00127 
00128   typedef Tpetra::Map<zlno_t, zgno_t, znode_t> map_t;
00129   typedef Tpetra::Export<zlno_t, zgno_t, znode_t> export_t;
00130   typedef Tpetra::Import<zlno_t, zgno_t, znode_t> import_t;
00131   typedef KokkosClassic::DefaultNode::DefaultNodeType default_znode_t;
00132 
00152   UserInputForTests(string path, string testData, 
00153     const RCP<const Comm<int> > &c, bool debugInfo=false,
00154     bool distributeInput=true);
00155 
00172   UserInputForTests(int x, int y, int z, string matrixType,
00173     const RCP<const Comm<int> > &c, bool debugInfo=false,
00174     bool distributeInput=true);
00175 
00178   static void getUIRandomData(unsigned int seed, zlno_t length,
00179     zscalar_t min, zscalar_t max, ArrayView<ArrayRCP<zscalar_t > > data);
00180 
00181   RCP<tMVector_t> getUICoordinates();
00182 
00183   RCP<tMVector_t> getUIWeights();
00184 
00185   RCP<tMVector_t> getUIEdgeWeights();
00186 
00187   RCP<tcrsMatrix_t> getUITpetraCrsMatrix();
00188 
00189   RCP<tcrsGraph_t> getUITpetraCrsGraph();
00190 
00191   RCP<tVector_t> getUITpetraVector();
00192 
00193   RCP<tMVector_t> getUITpetraMultiVector(int nvec);
00194 
00195   RCP<xcrsMatrix_t> getUIXpetraCrsMatrix();
00196 
00197   RCP<xcrsGraph_t> getUIXpetraCrsGraph();
00198 
00199   RCP<xVector_t> getUIXpetraVector();
00200 
00201   RCP<xMVector_t> getUIXpetraMultiVector(int nvec);
00202 
00203 #ifdef HAVE_EPETRA_DATA_TYPES
00204   RCP<Epetra_CrsGraph> getUIEpetraCrsGraph();
00205 
00206   RCP<Epetra_CrsMatrix> getUIEpetraCrsMatrix();
00207 
00208   RCP<Epetra_Vector> getUIEpetraVector();
00209 
00210   RCP<Epetra_MultiVector> getUIEpetraMultiVector(int nvec);
00211 #endif
00212 
00213 private:
00214 
00215   bool verbose_;
00216 
00217   const RCP<const Comm<int> > tcomm_;
00218 
00219   RCP<tcrsMatrix_t> M_; 
00220   RCP<xcrsMatrix_t> xM_; 
00221 
00222   RCP<tMVector_t> xyz_;
00223   RCP<tMVector_t> vtxWeights_;
00224   RCP<tMVector_t> edgWeights_;
00225 
00226 #ifdef HAVE_EPETRA_DATA_TYPES
00227   RCP<const Epetra_Comm> ecomm_;
00228   RCP<Epetra_CrsMatrix> eM_; 
00229   RCP<Epetra_CrsGraph> eG_; 
00230 #endif
00231 
00232   // Read a Matrix Market file into M_
00233   // using Tpetra::MatrixMarket::Reader.
00234   // If there are "Tim Davis" style coordinates
00235   // that go with the file,  read those into xyz_.
00236 
00237   void readMatrixMarketFile(string path, string testData);
00238 
00239   // Build matrix M_ from a mesh and a problem type
00240   // with Galeri::Xpetra.
00241 
00242   void buildCrsMatrix(int xdim, int ydim, int zdim, string type,
00243                       bool distributeInput);
00244 
00245   // Read a Zoltan1 Chaco or Matrix Market file
00246   // into M_.  If it has geometric coordinates,
00247   // read them into xyz_.  If it has weights,
00248   // read those into vtxWeights_ and edgWeights_.
00249 
00250   void readZoltanTestData(string path, string testData,
00251                           bool distributeInput);
00252 
00253   // Read Zoltan data that is in a .graph file.
00254   void getUIChacoGraph(FILE *fptr, string name, bool distributeInput);
00255 
00256   // Read Zoltan data that is in a .coords file.
00257   void getUIChacoCoords(FILE *fptr, string name);
00258 };
00259 
00260 UserInputForTests::UserInputForTests(string path, string testData, 
00261   const RCP<const Comm<int> > &c, bool debugInfo, bool distributeInput):
00262     verbose_(debugInfo),
00263     tcomm_(c), M_(), xM_(), xyz_(), vtxWeights_(), edgWeights_()
00264 #ifdef HAVE_EPETRA_DATA_TYPES
00265     ,ecomm_(), eM_(), eG_()
00266 #endif
00267 {
00268   bool zoltan1 = false;
00269   string::size_type loc = path.find("/zoltan/test/");  // Zoltan1 data
00270   if (loc != string::npos)
00271     zoltan1 = true;
00272 
00273   if (zoltan1)
00274     readZoltanTestData(path, testData, distributeInput);
00275   else
00276     readMatrixMarketFile(path, testData);
00277   
00278 #ifdef HAVE_EPETRA_DATA_TYPES
00279   ecomm_ = Xpetra::toEpetra(c);
00280 #endif
00281 }
00282 
00283 UserInputForTests::UserInputForTests(int x, int y, int z, 
00284   string matrixType, const RCP<const Comm<int> > &c, bool debugInfo,
00285   bool distributeInput):
00286     verbose_(debugInfo),
00287     tcomm_(c), M_(), xM_(), xyz_(), vtxWeights_(), edgWeights_()
00288 #ifdef HAVE_EPETRA_DATA_TYPES
00289     ,ecomm_(), eM_(), eG_() 
00290 #endif
00291 {
00292   if (matrixType.size() == 0){
00293     int dim = 0;
00294     if (x > 0) dim++;
00295     if (y > 0) dim++;
00296     if (z > 0) dim++;
00297     if (dim == 1)
00298       matrixType = string("Laplace1D");
00299     else if (dim == 2)
00300       matrixType = string("Laplace2D");
00301     else if (dim == 3)
00302       matrixType = string("Laplace3D");
00303     else
00304       throw std::runtime_error("input");
00305 
00306     if (verbose_ && tcomm_->getRank() == 0)
00307       std::cout << "UserInputForTests, Matrix type : " << matrixType << std::endl;
00308   }
00309 
00310   buildCrsMatrix(x, y, z, matrixType, distributeInput);
00311   
00312 #ifdef HAVE_EPETRA_DATA_TYPES
00313   ecomm_ = Xpetra::toEpetra(c);
00314 #endif
00315 }
00316 
00317 RCP<UserInputForTests::tMVector_t> UserInputForTests::getUICoordinates()
00318 { 
00319   if (xyz_.is_null())
00320     throw std::runtime_error("could not read coord file");
00321   return xyz_;
00322 }
00323 
00324 RCP<UserInputForTests::tMVector_t> UserInputForTests::getUIWeights()
00325 { 
00326   return vtxWeights_;
00327 }
00328 
00329 RCP<UserInputForTests::tMVector_t> UserInputForTests::getUIEdgeWeights()
00330 { 
00331   return edgWeights_;
00332 }
00333 
00334 RCP<UserInputForTests::tcrsMatrix_t> UserInputForTests::getUITpetraCrsMatrix() 
00335 { 
00336   if (M_.is_null())
00337     throw std::runtime_error("could not read mtx file");
00338   return M_;
00339 }
00340 
00341 RCP<UserInputForTests::tcrsGraph_t> UserInputForTests::getUITpetraCrsGraph() 
00342 { 
00343   if (M_.is_null())
00344     throw std::runtime_error("could not read mtx file");
00345   return rcp_const_cast<tcrsGraph_t>(M_->getCrsGraph());
00346 }
00347 
00348 RCP<UserInputForTests::tVector_t> UserInputForTests::getUITpetraVector() 
00349 { 
00350   RCP<tVector_t> V = rcp(new tVector_t(M_->getRowMap(),  1));
00351   V->randomize();
00352   
00353   return V;
00354 }
00355 
00356 RCP<UserInputForTests::tMVector_t> UserInputForTests::getUITpetraMultiVector(int nvec) 
00357 { 
00358   RCP<tMVector_t> mV = rcp(new tMVector_t(M_->getRowMap(), nvec));
00359   mV->randomize();
00360   
00361   return mV;
00362 }
00363 
00364 RCP<UserInputForTests::xcrsMatrix_t> UserInputForTests::getUIXpetraCrsMatrix() 
00365 { 
00366   if (M_.is_null())
00367     throw std::runtime_error("could not read mtx file");
00368   return xM_;
00369 }
00370 
00371 RCP<UserInputForTests::xcrsGraph_t> UserInputForTests::getUIXpetraCrsGraph() 
00372 { 
00373   if (M_.is_null())
00374     throw std::runtime_error("could not read mtx file");
00375   return rcp_const_cast<xcrsGraph_t>(xM_->getCrsGraph());
00376 }
00377 
00378 RCP<UserInputForTests::xVector_t> UserInputForTests::getUIXpetraVector() 
00379 { 
00380   RCP<const tVector_t> tV = getUITpetraVector();
00381   RCP<const xVector_t> xV =
00382     Zoltan2::XpetraTraits<tVector_t>::convertToXpetra(tV);
00383   return rcp_const_cast<xVector_t>(xV);
00384 }
00385 
00386 RCP<UserInputForTests::xMVector_t> UserInputForTests::getUIXpetraMultiVector(int nvec) 
00387 { 
00388   RCP<const tMVector_t> tMV = getUITpetraMultiVector(nvec);
00389   RCP<const xMVector_t> xMV =
00390     Zoltan2::XpetraTraits<tMVector_t>::convertToXpetra(tMV);
00391   return rcp_const_cast<xMVector_t>(xMV);
00392 }
00393 
00394 #ifdef HAVE_EPETRA_DATA_TYPES
00395 RCP<Epetra_CrsGraph> UserInputForTests::getUIEpetraCrsGraph()
00396 {
00397   if (M_.is_null())
00398     throw std::runtime_error("could not read mtx file");
00399   RCP<const tcrsGraph_t> tgraph = M_->getCrsGraph();
00400   RCP<const Tpetra::Map<zlno_t, zgno_t> > trowMap = tgraph->getRowMap();
00401   RCP<const Tpetra::Map<zlno_t, zgno_t> > tcolMap = tgraph->getColMap();
00402 
00403   int nElts = static_cast<int>(trowMap->getGlobalNumElements());
00404   int nMyElts = static_cast<int>(trowMap->getNodeNumElements());
00405   int base = trowMap->getIndexBase();
00406   ArrayView<const int> gids = trowMap->getNodeElementList();
00407 
00408   Epetra_BlockMap erowMap(nElts, nMyElts,
00409     gids.getRawPtr(), 1, base, *ecomm_);
00410 
00411   Array<int> rowSize(nMyElts);
00412   for (int i=0; i < nMyElts; i++){
00413     rowSize[i] = static_cast<int>(M_->getNumEntriesInLocalRow(i+base));
00414   }
00415 
00416   size_t maxRow = M_->getNodeMaxNumRowEntries();
00417   Array<int> colGids(maxRow);
00418   ArrayView<const int> colLid;
00419 
00420   eG_ = rcp(new Epetra_CrsGraph(Copy, erowMap, 
00421     rowSize.getRawPtr(), true));
00422 
00423   for (int i=0; i < nMyElts; i++){
00424     tgraph->getLocalRowView(i+base, colLid);
00425     for (int j=0; j < colLid.size(); j++)
00426       colGids[j] = tcolMap->getGlobalElement(colLid[j]);
00427     eG_->InsertGlobalIndices(gids[i], rowSize[i], colGids.getRawPtr());
00428   }
00429   eG_->FillComplete();
00430   return eG_;
00431 }
00432 
00433 RCP<Epetra_CrsMatrix> UserInputForTests::getUIEpetraCrsMatrix()
00434 {
00435   if (M_.is_null())
00436     throw std::runtime_error("could not read mtx file");
00437   RCP<Epetra_CrsGraph> egraph = getUIEpetraCrsGraph();
00438   eM_ = rcp(new Epetra_CrsMatrix(Copy, *egraph));
00439 
00440   size_t maxRow = M_->getNodeMaxNumRowEntries();
00441   int nrows = egraph->NumMyRows();
00442   int base = egraph->IndexBase();
00443   const Epetra_BlockMap &rowMap = egraph->RowMap();
00444   const Epetra_BlockMap &colMap = egraph->ColMap();
00445   Array<int> colGid(maxRow);
00446 
00447   for (int i=0; i < nrows; i++){
00448     ArrayView<const int> colLid;
00449     ArrayView<const zscalar_t> nz;
00450     M_->getLocalRowView(i+base, colLid, nz);
00451     size_t rowSize = colLid.size();
00452     int rowGid = rowMap.GID(i+base);
00453     for (size_t j=0; j < rowSize; j++){
00454       colGid[j] = colMap.GID(colLid[j]);
00455     }
00456     eM_->InsertGlobalValues(
00457       rowGid, rowSize, nz.getRawPtr(), colGid.getRawPtr());
00458   }
00459   eM_->FillComplete();
00460   return eM_;
00461 }
00462 
00463 RCP<Epetra_Vector> UserInputForTests::getUIEpetraVector() 
00464 { 
00465   RCP<Epetra_CrsGraph> egraph = getUIEpetraCrsGraph();
00466   RCP<Epetra_Vector> V = rcp(new Epetra_Vector(egraph->RowMap()));
00467   V->Random();
00468   return V;
00469 }
00470 
00471 RCP<Epetra_MultiVector> UserInputForTests::getUIEpetraMultiVector(int nvec) 
00472 { 
00473   RCP<Epetra_CrsGraph> egraph = getUIEpetraCrsGraph();
00474   RCP<Epetra_MultiVector> mV = 
00475     rcp(new Epetra_MultiVector(egraph->RowMap(), nvec));
00476   mV->Random();
00477   return mV;
00478 }
00479 #endif
00480 
00481 void UserInputForTests::getUIRandomData(unsigned int seed, zlno_t length, 
00482     zscalar_t min, zscalar_t max,
00483     ArrayView<ArrayRCP<zscalar_t > > data)
00484 {
00485   if (length < 1)
00486     return;
00487 
00488   size_t dim = data.size();
00489   for (size_t i=0; i < dim; i++){
00490     zscalar_t *tmp = new zscalar_t [length];
00491     if (!tmp)
00492        throw (std::bad_alloc());
00493     data[i] = Teuchos::arcp(tmp, 0, length, true);
00494   }
00495 
00496   zscalar_t scalingFactor = (max-min) / RAND_MAX;
00497   srand(seed);
00498   for (size_t i=0; i < dim; i++){
00499     zscalar_t *x = data[i].getRawPtr();
00500     for (zlno_t j=0; j < length; j++)
00501       *x++ = min + (zscalar_t(rand()) * scalingFactor);
00502   }
00503 }
00504 
00505 void UserInputForTests::readMatrixMarketFile(string path, string testData)
00506 {
00507   std::ostringstream fname;
00508   fname << path << "/" << testData << ".mtx";
00509   RCP<KokkosClassic::DefaultNode::DefaultNodeType> dnode 
00510     = KokkosClassic::DefaultNode::getDefaultNode();
00511 
00512   if (verbose_ && tcomm_->getRank() == 0)
00513     std::cout << "UserInputForTests, Read: " << fname.str() << std::endl;
00514 
00515   // This reader has some problems.  "Pattern" matrices
00516   // cannot be read.  Until the
00517   // reader is fixed, we'll have to get inputs that are consistent with
00518   // the reader. (Tpetra bug 5611 and 5624)
00519  
00520   bool aok = true;
00521   try{
00522     M_ = Tpetra::MatrixMarket::Reader<tcrsMatrix_t>::readSparseFile(
00523       fname.str(), tcomm_, dnode, true, true, false);
00524   }
00525   catch (std::exception &e) {
00526     //TEST_FAIL_AND_THROW(*tcomm_, 1, e.what());
00527     aok = false;
00528   }
00529 
00530   if (aok){
00531     RCP<const xcrsMatrix_t> xm = 
00532       Zoltan2::XpetraTraits<tcrsMatrix_t>::convertToXpetra(M_);
00533     xM_ = rcp_const_cast<xcrsMatrix_t>(xm);
00534   }
00535   else{
00536     if (tcomm_->getRank() == 0)
00537       std::cout << "UserInputForTests unable to read matrix." << std::endl;
00538   }
00539 
00540   // Open the coordinate file.
00541 
00542   fname.str("");
00543   fname << path << "/" << testData << "_coord.mtx";
00544 
00545   size_t coordDim = 0, numGlobalCoords = 0;
00546   size_t msg[2]={0,0};
00547   ArrayRCP<ArrayRCP<zscalar_t> > xyz;
00548   std::ifstream coordFile;
00549 
00550   if (tcomm_->getRank() == 0){
00551 
00552     if (verbose_)
00553       std::cout << "UserInputForTests, Read: " << 
00554          fname.str() << std::endl;
00555   
00556     int fail = 0;
00557     try{
00558       coordFile.open(fname.str().c_str());
00559     }
00560     catch (std::exception &e){ // there is no coordinate file
00561       fail = 1;
00562     }
00563   
00564     if (!fail){
00565   
00566       // Read past banner to number and dimension of coordinates.
00567     
00568       char c[256];
00569       bool done=false;
00570     
00571       while (!done && !fail && coordFile.good()){
00572         coordFile.getline(c, 256);
00573         if (!c[0])
00574           fail = 1;
00575         else if (c[0] == '%')
00576           continue;
00577         else {
00578           done=true;
00579           std::istringstream s(c);
00580           s >> numGlobalCoords >> coordDim;
00581           if (!s.eof() || numGlobalCoords < 1 || coordDim < 1)
00582             fail=1;
00583         }
00584       }
00585 
00586       if (done){
00587     
00588         // Read in the coordinates.
00589       
00590         xyz = Teuchos::arcp(new ArrayRCP<zscalar_t> [coordDim], 0, coordDim);
00591       
00592         for (size_t dim=0; !fail && dim < coordDim; dim++){
00593           size_t idx;
00594           zscalar_t *tmp = new zscalar_t [numGlobalCoords];
00595           if (!tmp)
00596             fail = 1;
00597           else{
00598             xyz[dim] = Teuchos::arcp(tmp, 0, numGlobalCoords);
00599         
00600             for (idx=0; !coordFile.eof() && idx < numGlobalCoords; idx++){
00601               coordFile.getline(c, 256);
00602               std::istringstream s(c);
00603               s >> tmp[idx];
00604             }
00605       
00606             if (idx < numGlobalCoords)
00607               fail = 1;
00608           }
00609         }
00610 
00611         if (fail){
00612           ArrayRCP<zscalar_t> emptyArray;
00613           for (size_t dim=0; dim < coordDim; dim++)
00614             xyz[dim] = emptyArray;   // free the memory
00615 
00616           coordDim = 0;
00617         }
00618       }
00619       else{
00620         fail = 1;
00621       }
00622     
00623       coordFile.close();
00624     }
00625 
00626     msg[0] = coordDim;
00627     msg[1] = numGlobalCoords;
00628   }
00629 
00630   // Broadcast coordinate dimension
00631   Teuchos::broadcast<int, size_t>(*tcomm_, 0, 2, msg);
00632 
00633   coordDim = msg[0];
00634   numGlobalCoords= msg[1];
00635 
00636   if (coordDim == 0)
00637     return;
00638 
00639   zgno_t base;
00640   RCP<const map_t> toMap;
00641 
00642   if (!M_.is_null()){
00643     base = M_->getIndexBase();
00644     const RCP<const map_t> &mapM = M_->getRowMap();
00645     toMap = mapM;
00646   }
00647   else{
00648     base = 0;
00649     toMap = rcp(new map_t(numGlobalCoords, base, tcomm_));
00650   }
00651 
00652   // Export coordinates to their owners
00653 
00654   xyz_ = rcp(new tMVector_t(toMap, coordDim));
00655 
00656   ArrayRCP<ArrayView<const zscalar_t> > coordLists(coordDim);
00657 
00658   if (tcomm_->getRank() == 0){
00659 
00660     for (size_t dim=0; dim < coordDim; dim++)
00661       coordLists[dim] = xyz[dim].view(0, numGlobalCoords);
00662 
00663     zgno_t *tmp = new zgno_t [numGlobalCoords];
00664     if (!tmp)
00665       throw std::bad_alloc();
00666 
00667     ArrayRCP<const zgno_t> rowIds = Teuchos::arcp(tmp, 0, numGlobalCoords);
00668 
00669     zgno_t basePlusNumGlobalCoords = base+numGlobalCoords;
00670     for (zgno_t id=base; id < basePlusNumGlobalCoords; id++)
00671       *tmp++ = id;
00672 
00673     RCP<const map_t> fromMap = rcp(new map_t(numGlobalCoords, 
00674       rowIds.view(0, numGlobalCoords), base, tcomm_));
00675 
00676     tMVector_t allCoords(fromMap, coordLists.view(0, coordDim), coordDim);
00677 
00678     export_t exporter(fromMap, toMap);
00679 
00680     xyz_->doExport(allCoords, exporter, Tpetra::INSERT);
00681   }
00682   else{
00683 
00684     RCP<const map_t> fromMap = rcp(new map_t(numGlobalCoords, 
00685       ArrayView<zgno_t>(), base, tcomm_));
00686 
00687     tMVector_t allCoords(fromMap, coordLists.view(0, coordDim), coordDim);
00688 
00689     export_t exporter(fromMap, toMap);
00690 
00691     xyz_->doExport(allCoords, exporter, Tpetra::INSERT);
00692   }
00693 }
00694 
00695 void UserInputForTests::buildCrsMatrix(int xdim, int ydim, int zdim, 
00696   string problemType, bool distributeInput)
00697 {
00698   Teuchos::CommandLineProcessor tclp;
00699   Galeri::Xpetra::Parameters<zgno_t> params(tclp,
00700      xdim, ydim, zdim, problemType);
00701 
00702   RCP<const Tpetra::Map<zlno_t, zgno_t> > map;
00703   if (distributeInput)
00704     map = rcp(new Tpetra::Map<zlno_t, zgno_t>(params.GetNumGlobalElements(),
00705                                             0, tcomm_));
00706   else {
00707     // All data initially on rank 0
00708     size_t nGlobalElements = params.GetNumGlobalElements();
00709     size_t nLocalElements = ((tcomm_->getRank() == 0) ? nGlobalElements : 0);
00710     map = rcp(new Tpetra::Map<zlno_t, zgno_t>(nGlobalElements, nLocalElements, 0,
00711                                             tcomm_));
00712   }
00713 
00714   if (verbose_ && tcomm_->getRank() == 0){
00715     std::cout << "UserInputForTests, Create matrix with " << problemType;
00716     std::cout << " (and " << xdim;
00717     if (zdim > 0)
00718       std::cout << " x " << ydim << " x " << zdim;
00719     else if (ydim > 0)
00720       std::cout << " x"  << ydim << " x 1";
00721     else
00722       std::cout << "x 1 x 1";
00723 
00724     std::cout << " mesh)" << std::endl;
00725   }
00726 
00727   try{
00728     RCP<Galeri::Xpetra::Problem<Tpetra::Map<zlno_t, zgno_t>, Tpetra::CrsMatrix<zscalar_t, zlno_t, zgno_t>, Tpetra::MultiVector<zscalar_t, zlno_t, zgno_t> > > Pr =
00729         Galeri::Xpetra::BuildProblem<zscalar_t, zlno_t, zgno_t, Tpetra::Map<zlno_t, zgno_t>, Tpetra::CrsMatrix<zscalar_t, zlno_t, zgno_t>, Tpetra::MultiVector<zscalar_t, zlno_t, zgno_t> >
00730         (params.GetMatrixType(), map, params.GetParameterList());
00731     M_ = Pr->BuildMatrix();
00732   }
00733   catch (std::exception &e) {    // Probably not enough memory
00734     TEST_FAIL_AND_THROW(*tcomm_, 1, e.what());
00735   }
00736 
00737   RCP<const xcrsMatrix_t> xm = 
00738     Zoltan2::XpetraTraits<tcrsMatrix_t>::convertToXpetra(M_);
00739   xM_ = rcp_const_cast<xcrsMatrix_t>(xm);
00740 
00741   // Compute the coordinates for the matrix rows.
00742 
00743   if (verbose_ && tcomm_->getRank() == 0)
00744     std::cout << 
00745      "UserInputForTests, Implied matrix row coordinates computed" << 
00746        std::endl;
00747   
00748   ArrayView<const zgno_t> gids = map->getNodeElementList();
00749   zlno_t count = gids.size();
00750   int dim = 3;
00751   size_t pos = problemType.find("2D");
00752   if (pos != string::npos)
00753     dim = 2;
00754   else if (problemType == string("Laplace1D") || 
00755            problemType == string("Identity"))
00756     dim = 1;
00757 
00758   Array<ArrayRCP<zscalar_t> > coordinates(dim);
00759 
00760   if (count > 0){
00761     for (int i=0; i < dim; i++){
00762       zscalar_t *c = new zscalar_t [count];
00763       if (!c)
00764         throw(std::bad_alloc());
00765       coordinates[i] = Teuchos::arcp(c, 0, count, true);
00766     }
00767   
00768     if (dim==3){
00769       zscalar_t *x = coordinates[0].getRawPtr();
00770       zscalar_t *y = coordinates[1].getRawPtr();
00771       zscalar_t *z = coordinates[2].getRawPtr();
00772       zgno_t xySize = xdim * ydim;
00773       for (zlno_t i=0; i < count; i++){
00774         zgno_t iz = gids[i] / xySize;
00775         zgno_t xy = gids[i] - iz*xySize;
00776         z[i] = zscalar_t(iz);
00777         y[i] = zscalar_t(xy / xdim);
00778         x[i] = zscalar_t(xy % xdim);
00779       }
00780     }
00781     else if (dim==2){
00782       zscalar_t *x = coordinates[0].getRawPtr();
00783       zscalar_t *y = coordinates[1].getRawPtr();
00784       for (zlno_t i=0; i < count; i++){
00785         y[i] = zscalar_t(gids[i] / xdim);
00786         x[i] = zscalar_t(gids[i] % xdim);
00787       }
00788     }
00789     else{
00790       zscalar_t *x = coordinates[0].getRawPtr();
00791       for (zlno_t i=0; i < count; i++)
00792         x[i] = zscalar_t(gids[i]);
00793     }
00794   }
00795 
00796   Array<ArrayView<const zscalar_t> > coordView(dim);
00797   if (count > 0)
00798     for (int i=0; i < dim; i++)
00799       coordView[i] = coordinates[i].view(0,count);
00800 
00801   xyz_ = rcp(new tMVector_t(map, coordView.view(0, dim), dim));
00802 }
00803 
00804 void UserInputForTests::readZoltanTestData(string path, string testData,
00805                                            bool distributeInput)
00806 {
00807   int rank = tcomm_->getRank();
00808   FILE *graphFile = NULL;
00809   FILE *coordFile = NULL;
00810   int fileInfo[2];
00811 
00812   if (rank == 0){
00813     std::ostringstream chGraphFileName;
00814     chGraphFileName << path << "/ch_" << testData << "/" << testData << ".graph";
00815     std::ostringstream chCoordFileName;
00816     chCoordFileName << path << "/ch_" << testData << "/" << testData << ".coords";
00817     memset(fileInfo, 0, sizeof(int) * 2);
00818 
00819     graphFile = fopen(chGraphFileName.str().c_str(), "r");
00820     if (graphFile){
00821       fileInfo[0] = 1;
00822       if (verbose_ && tcomm_->getRank() == 0)
00823         std::cout << "UserInputForTests, open " << 
00824           chGraphFileName << std::endl;
00825       
00826       coordFile = fopen(chCoordFileName.str().c_str(), "r");
00827       if (coordFile){
00828         fileInfo[1] = 1;
00829         if (verbose_ && tcomm_->getRank() == 0)
00830           std::cout << "UserInputForTests, open " << 
00831             chCoordFileName << std::endl;
00832       }
00833     }
00834   }
00835 
00836   Teuchos::broadcast<int, int>(*tcomm_, 0, 2, fileInfo);
00837 
00838   bool haveGraph = (fileInfo[0] == 1);
00839   bool haveCoords = (fileInfo[1] == 1);
00840 
00841   if (haveGraph){
00842     // Writes M_, vtxWeights_, and edgWeights_ and closes file.
00843     try{
00844       getUIChacoGraph(graphFile, testData, distributeInput);
00845     }
00846     Z2_FORWARD_EXCEPTIONS
00847 
00848     if (haveCoords){
00849       // Writes xyz_ and closes the file.
00850       try{
00851         getUIChacoCoords(coordFile, testData);
00852       }
00853       Z2_FORWARD_EXCEPTIONS
00854     }
00855   }
00856 
00857   RCP<const xcrsMatrix_t> xm =
00858     Zoltan2::XpetraTraits<tcrsMatrix_t>::convertToXpetra(M_);
00859   xM_ = rcp_const_cast<xcrsMatrix_t>(xm);
00860 }
00861 
00862 void UserInputForTests::getUIChacoGraph(FILE *fptr, string fname,
00863                                       bool distributeInput)
00864 {
00865   int rank = tcomm_->getRank();
00866   int graphCounts[5];
00867   int nvtxs=0, nedges=0;
00868   int nVwgts=0, nEwgts=0;
00869   int *start = NULL, *adj = NULL;
00870   float *ewgts = NULL, *vwgts = NULL;
00871   size_t *nzPerRow = NULL;
00872   size_t maxRowLen = 0;
00873   zgno_t base = 0;
00874   ArrayRCP<const size_t> rowSizes;
00875   int fail = 0;
00876   bool haveEdges = true;
00877 
00878   if (rank == 0){
00879 
00880     memset(graphCounts, 0, 5*sizeof(int));
00881   
00882     // This function is in the Zoltan C-library.
00883   
00884     // Reads in the file and closes it when done.
00885     char *nonConstName = new char [fname.size() + 1];
00886     strcpy(nonConstName, fname.c_str());
00887 
00888     fail = Zoltan2::chaco_input_graph(fptr, nonConstName,
00889       &start, &adj, &nvtxs, &nVwgts, &vwgts, &nEwgts, &ewgts);
00890     delete [] nonConstName;
00891 
00892     // There are Zoltan2 test graphs that have no edges.
00893 
00894     if (start == NULL)
00895       haveEdges = false;
00896 
00897     if (verbose_){
00898       std::cout << "UserInputForTests, " << nvtxs << " vertices,";
00899       if (haveEdges)
00900         std::cout << start[nvtxs] << " edges,";
00901       else
00902         std::cout << "no edges,";
00903       std::cout << nVwgts << " vertex weights, ";
00904       std::cout << nEwgts << " edge weights" << std::endl;
00905     }
00906 
00907     if (nvtxs==0)
00908       fail = true;
00909 
00910     if (fail){
00911       Teuchos::broadcast<int, int>(*tcomm_, 0, 5, graphCounts);
00912       throw std::runtime_error("Unable to read chaco file");
00913     }
00914 
00915     if (haveEdges)
00916       nedges = start[nvtxs];
00917 
00918     nzPerRow = new size_t [nvtxs];
00919     if (!nzPerRow)
00920       throw std::bad_alloc();
00921     rowSizes = arcp(nzPerRow, 0, nvtxs, true);
00922 
00923     if (haveEdges){
00924       for (int i=0; i < nvtxs; i++){
00925         nzPerRow[i] = start[i+1] - start[i];
00926         if (nzPerRow[i] > maxRowLen)
00927           maxRowLen = nzPerRow[i];
00928       }
00929     }
00930     else{
00931       memset(nzPerRow, 0, sizeof(size_t) * nvtxs);
00932     }
00933 
00934     if (haveEdges){
00935       free(start);
00936       start = NULL;
00937     }
00938   
00939     // Make sure base gid is zero.
00940 
00941     if (nedges){
00942       int chbase = adj[0];
00943       for (int i=1; i < nedges; i++)
00944         if (adj[i] < chbase)
00945           chbase = adj[i];
00946   
00947       if (chbase > 0){
00948         for (int i=0; i < nedges; i++)
00949           adj[i] -= chbase;
00950       }
00951     }
00952 
00953     graphCounts[0] = nvtxs;
00954     graphCounts[1] = nedges;
00955     graphCounts[2] = nVwgts;
00956     graphCounts[3] = nEwgts;
00957     graphCounts[4] = maxRowLen; // size_t maxRowLen will fit; it is <= (int-int)
00958   }
00959   
00960   Teuchos::broadcast<int, int>(*tcomm_, 0, 5, graphCounts);
00961 
00962   if (graphCounts[0] == 0)
00963     throw std::runtime_error("Unable to read chaco file");
00964 
00965   haveEdges = (graphCounts[1] > 0);
00966 
00967   RCP<tcrsMatrix_t> fromMatrix;
00968   RCP<const map_t> fromMap;
00969 
00970   // Create a Tpetra::CrsMatrix where rank 0 has entire matrix.
00971 
00972   if (rank == 0){
00973     fromMap = rcp(new map_t(nvtxs, nvtxs, base, tcomm_));
00974 
00975     fromMatrix = 
00976       rcp(new tcrsMatrix_t(fromMap, rowSizes, Tpetra::StaticProfile));
00977 
00978     if (haveEdges){
00979 
00980       zgno_t *edgeIds = new zgno_t [nedges];
00981       if (nedges && !edgeIds)
00982         throw std::bad_alloc();
00983       for (int i=0; i < nedges; i++)
00984          edgeIds[i] = adj[i];
00985   
00986       free(adj);
00987       adj = NULL; 
00988 
00989       zgno_t *nextId = edgeIds;
00990       Array<zscalar_t> values(maxRowLen, 1.0);
00991   
00992       for (int i=0; i < nvtxs; i++){
00993         if (nzPerRow[i] > 0){
00994           ArrayView<const zgno_t> rowNz(nextId, nzPerRow[i]);
00995           fromMatrix->insertGlobalValues(i, rowNz, values.view(0,nzPerRow[i]));
00996           nextId += nzPerRow[i];
00997         }
00998       }
00999   
01000       delete [] edgeIds;
01001       edgeIds = NULL;
01002     }
01003 
01004     fromMatrix->fillComplete();
01005   }
01006   else{
01007     nvtxs = graphCounts[0];
01008     nedges = graphCounts[1];
01009     nVwgts = graphCounts[2];
01010     nEwgts  = graphCounts[3];
01011     maxRowLen = graphCounts[4];
01012 
01013     // Create a Tpetra::CrsMatrix where rank 0 has entire matrix.
01014 
01015     fromMap = rcp(new map_t(nvtxs, 0, base, tcomm_));
01016 
01017     fromMatrix = 
01018       rcp(new tcrsMatrix_t(fromMap, rowSizes, Tpetra::StaticProfile));
01019 
01020     fromMatrix->fillComplete();
01021   }
01022 
01023   
01024   RCP<const map_t> toMap;
01025   RCP<tcrsMatrix_t> toMatrix;
01026   RCP<import_t> importer;
01027 
01028   if (distributeInput) {
01029     // Create a Tpetra::CrsMatrix with default row distribution.
01030     toMap = rcp(new map_t(nvtxs, base, tcomm_));
01031     toMatrix = rcp(new tcrsMatrix_t(toMap, maxRowLen));
01032 
01033     // Import the data.
01034     importer = rcp(new import_t(fromMap, toMap));
01035     toMatrix->doImport(*fromMatrix, *importer, Tpetra::INSERT);
01036     toMatrix->fillComplete();
01037   }
01038   else {
01039     toMap = fromMap;
01040     toMatrix = fromMatrix;
01041   }
01042 
01043   M_ = toMatrix;
01044 
01045   // Vertex weights, if any
01046 
01047   typedef ArrayRCP<const ArrayView<const zscalar_t> > arrayArray_t;
01048 
01049   if (nVwgts > 0){
01050 
01051     ArrayRCP<zscalar_t> weightBuf;
01052     ArrayView<const zscalar_t> *wgts = new ArrayView<const zscalar_t> [nVwgts];
01053 
01054     if (rank == 0){
01055       size_t len = nVwgts * nvtxs;
01056       zscalar_t *buf = new zscalar_t [len];
01057       if (!buf) throw std::bad_alloc();
01058       weightBuf = arcp(buf, 0, len, true);
01059 
01060       for (int widx=0; widx < nVwgts; widx++){
01061         wgts[widx] = ArrayView<const zscalar_t>(buf, nvtxs);
01062         float *vw = vwgts + widx;
01063         for (int i=0; i < nvtxs; i++, vw += nVwgts)
01064           buf[i] = *vw;
01065         buf += nvtxs;
01066       }
01067 
01068       free(vwgts);
01069       vwgts = NULL;
01070     }
01071 
01072     arrayArray_t vweights = arcp(wgts, 0, nVwgts, true);
01073 
01074     RCP<tMVector_t> fromVertexWeights = 
01075       rcp(new tMVector_t(fromMap, vweights.view(0, nVwgts), nVwgts));
01076 
01077     RCP<tMVector_t> toVertexWeights;
01078     if (distributeInput) {
01079       toVertexWeights = rcp(new tMVector_t(toMap, nVwgts));
01080       toVertexWeights->doImport(*fromVertexWeights, *importer, Tpetra::INSERT);
01081     }
01082     else
01083       toVertexWeights = fromVertexWeights;
01084 
01085     vtxWeights_ = toVertexWeights;
01086   }
01087 
01088   // Edge weights, if any
01089 
01090   if (haveEdges && nEwgts > 0){
01091 
01092     ArrayRCP<zscalar_t> weightBuf;
01093     ArrayView<const zscalar_t> *wgts = new ArrayView<const zscalar_t> [nEwgts];
01094 
01095     toMap = rcp(new map_t(nedges, M_->getNodeNumEntries(), base, tcomm_));
01096 
01097     if (rank == 0){
01098       size_t len = nEwgts * nedges;
01099       zscalar_t *buf = new zscalar_t [len];
01100       if (!buf) throw std::bad_alloc();
01101       weightBuf = arcp(buf, 0, len, true);
01102 
01103       for (int widx=0; widx < nEwgts; widx++){
01104         wgts[widx] = ArrayView<const zscalar_t>(buf, nedges);
01105         float *ew = ewgts + widx;
01106         for (int i=0; i < nedges; i++, ew += nEwgts)
01107           buf[i] = *ew;
01108         buf += nedges;
01109       }
01110 
01111       free(ewgts);
01112       ewgts = NULL;
01113       fromMap = rcp(new map_t(nedges, nedges, base, tcomm_));
01114     }
01115     else{
01116       fromMap = rcp(new map_t(nedges, 0, base, tcomm_));
01117     }
01118 
01119     arrayArray_t eweights = arcp(wgts, 0, nEwgts, true);
01120 
01121     RCP<tMVector_t> fromEdgeWeights;
01122     RCP<tMVector_t> toEdgeWeights;
01123     RCP<import_t> edgeImporter;
01124     if (distributeInput) {
01125       fromEdgeWeights = 
01126         rcp(new tMVector_t(fromMap, eweights.view(0, nEwgts), nEwgts));
01127       toEdgeWeights = rcp(new tMVector_t(toMap, nEwgts));
01128       edgeImporter = rcp(new import_t(fromMap, toMap));
01129       toEdgeWeights->doImport(*fromEdgeWeights, *edgeImporter, Tpetra::INSERT);
01130     }
01131     else
01132       toEdgeWeights = fromEdgeWeights;
01133 
01134     edgWeights_ = toEdgeWeights;
01135   }
01136 }
01137 
01138 void UserInputForTests::getUIChacoCoords(FILE *fptr, string fname)
01139 {
01140   int rank = tcomm_->getRank();
01141   int ndim=0;
01142   float *x=NULL, *y=NULL, *z=NULL;
01143   int fail = 0;
01144 
01145   size_t localNumVtx = M_->getNodeNumRows();
01146   size_t globalNumVtx = M_->getGlobalNumRows();
01147 
01148   if (rank == 0){
01149   
01150     // This function is in the Zoltan C-library.
01151   
01152     // Reads in the file and closes it when done.
01153     char *nonConstName = new char [fname.size() + 1];
01154     strcpy(nonConstName, fname.c_str());
01155     fail = Zoltan2::chaco_input_geom(fptr, nonConstName, globalNumVtx,
01156       &ndim, &x, &y, &z);
01157     delete [] nonConstName;
01158 
01159     if (fail)
01160       ndim = 0;
01161 
01162     if (verbose_){
01163       std::cout << "UserInputForTests, read " << globalNumVtx;
01164       std::cout << " " << ndim << "-dimensional coordinates." << std::endl;
01165     }
01166   }
01167   
01168   Teuchos::broadcast<int, int>(*tcomm_, 0, 1, &ndim);
01169 
01170   if (ndim == 0)
01171     throw std::runtime_error("Can't read coordinate file");
01172 
01173   ArrayRCP<ArrayRCP<const zscalar_t> > coords(ndim);
01174   zlno_t len = 0;
01175 
01176   if (rank == 0){
01177 
01178     for (int dim=0; dim < ndim; dim++){
01179       zscalar_t *v = new zscalar_t [globalNumVtx];
01180       if (!v)
01181         throw std::bad_alloc();
01182       coords[dim] = arcp<const zscalar_t>(v, 0, globalNumVtx, true);
01183       float *val = (dim==0 ? x : (dim==1 ? y : z));
01184       for (size_t i=0; i < globalNumVtx; i++)
01185         v[i] = zscalar_t(val[i]);
01186 
01187       free(val);
01188     }
01189 
01190     len = globalNumVtx;;
01191   }
01192 
01193   RCP<const map_t> fromMap = rcp(new map_t(globalNumVtx, len, 0, tcomm_));
01194   RCP<const map_t> toMap = rcp(new map_t(globalNumVtx, localNumVtx, 0, tcomm_));
01195   RCP<import_t> importer = rcp(new import_t(fromMap, toMap));
01196 
01197   Array<ArrayView<const zscalar_t> > coordData;
01198   for (int dim=0; dim < ndim; dim++)
01199     coordData.push_back(coords[dim].view(0, len));
01200 
01201   RCP<tMVector_t> fromCoords = 
01202     rcp(new tMVector_t(fromMap, coordData.view(0, ndim), ndim));
01203 
01204   RCP<tMVector_t> toCoords = rcp(new tMVector_t(toMap, ndim));
01205 
01206   toCoords->doImport(*fromCoords, *importer, Tpetra::INSERT);
01207 
01208   xyz_ = toCoords;
01209 
01210 }
01211 
01212 #endif