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