Tpetra Matrix/Vector Services Version of the Day
Tpetra_MatrixIO_def.hpp
00001 #ifndef TPETRA_MATRIX_IO_DEF
00002 #define TPETRA_MATRIX_IO_DEF
00003 
00004 #include "Tpetra_CrsMatrix.hpp"
00005 #include "Tpetra_MatrixIO.hpp"
00006 #include <iostream>
00007 
00008 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00009 void
00010 Tpetra::Utils::generateMatrix(const Teuchos::RCP<Teuchos::ParameterList> &plist,
00011                               const Teuchos::RCP<const Teuchos::Comm<int> > &comm, 
00012                               const Teuchos::RCP<Node> &node,
00013                               Teuchos::RCP< Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > &A) 
00014 {
00015   typedef Teuchos::ScalarTraits<Scalar> ST;
00016   TEST_FOR_EXCEPTION( plist == Teuchos::null, std::runtime_error,
00017       "Tpetra::Utils::generateMatrix(): ParameterList is null.");
00018   TEST_FOR_EXCEPTION( Teuchos::isParameterType<std::string>(*plist,"mat_type") == false, std::runtime_error,
00019       "Tpetra::Utils::generateMatrix(): ParameterList did not contain string parameter ""mat_type"".");
00020   std::string mat_type = plist->get<std::string>("mat_type");
00021   if (mat_type == "Lap3D") {
00022     // 3D Laplacian, grid is a cube with dimension gridSize x gridSize x gridSize
00023     const int gridSize = plist->get<int>("gridSize",100);
00024     const GlobalOrdinal gS2 = (GlobalOrdinal)gridSize*(GlobalOrdinal)gridSize;
00025     const GlobalOrdinal numRows = gS2*(GlobalOrdinal)gridSize;
00026     Teuchos::RCP<Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap;
00027     rowMap = Teuchos::rcp(new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>((global_size_t)numRows,(GlobalOrdinal)0,comm,GloballyDistributed,node));
00028     // create with DynamicProfile, so that the fillComplete(DoOptimizeStorage) can do first-touch reallocation 
00029     A = rcp(new Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rowMap,7,Tpetra::DynamicProfile));
00030     // fill matrix, one row at a time
00031     Teuchos::Array<GlobalOrdinal> neighbors;
00032     Teuchos::Array<Scalar> values(7, -ST::one());
00033     values[0] = (Scalar)6;
00034     for (GlobalOrdinal r = rowMap->getMinGlobalIndex(); r <= rowMap->getMaxGlobalIndex(); ++r) {
00035       neighbors.clear();
00036       neighbors.push_back(r); // add diagonal
00037       GlobalOrdinal ixy, iz, ix, iy;  // (x,y,z) coords and index in xy plane
00038       ixy = r%gS2;
00039       iz = (r - ixy)/gS2;
00040       ix = ixy%gridSize;
00041       iy = (ixy - ix)/gridSize;
00042       //
00043       if ( ix != 0 )          neighbors.push_back( r-1 );
00044       if ( ix != gridSize-1 ) neighbors.push_back( r+1 );
00045       if ( iy != 0 )          neighbors.push_back( r-gridSize );
00046       if ( iy != gridSize-1 ) neighbors.push_back( r+gridSize );
00047       if ( iz != 0 )          neighbors.push_back( r-gS2 );
00048       if ( iz != gridSize-1 ) neighbors.push_back( r+gS2 );
00049       A->insertGlobalValues( r, neighbors(), values(0,neighbors.size()) );
00050     }
00051     A->fillComplete(DoOptimizeStorage);
00052   }
00053   else {
00054     TEST_FOR_EXCEPTION( true, std::runtime_error, 
00055         "Tpetra::Utils::generateMatrix(): ParameterList specified unsupported ""mat_type"".");
00056   }
00057 }
00058 
00059 
00060 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00061 void
00062 Tpetra::Utils::readHBMatrix(const std::string &filename, 
00063                              const Teuchos::RCP<const Teuchos::Comm<int> > &comm, 
00064                              const Teuchos::RCP<Node> &node,
00065                              Teuchos::RCP< Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > &A,
00066                              Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap)
00067 {
00068   const int myRank = comm->getRank();
00069   int numRows,numCols,numNZ;
00070   Teuchos::ArrayRCP<Scalar> svals;
00071   Teuchos::ArrayRCP<GlobalOrdinal> colinds;
00072   Teuchos::ArrayRCP<int>           rowptrs;
00073   Teuchos::ArrayRCP<size_t>        nnzPerRow;
00074   //Teuchos::ArrayRCP<char>          type;
00075   int fail = 0;
00076   if (myRank == 0) {
00077     bool isSymmetric=false;
00078     Teuchos::ArrayRCP<double> dvals;
00079     Teuchos::ArrayRCP<int> colptrs, rowinds;
00080     std::string type;
00081     Tpetra::Utils::readHBMatDouble(filename,numRows,numCols,numNZ,type,colptrs,rowinds,dvals);
00082     TEST_FOR_EXCEPT(type.size() != 3);
00083     if (type[0] != 'R' && type[0] != 'r') {
00084       // only real matrices right now
00085       fail = 1;
00086     }
00087     if (fail == 0 && numNZ > 0) {
00088       if (type[1] == 'S' || type[1] == 's') {
00089         isSymmetric = true;
00090       }
00091       else {
00092         isSymmetric = false;
00093       }
00094     }
00095     if (fail == 0 && numNZ > 0) {
00096       // find num non-zero per row
00097       nnzPerRow = Teuchos::arcp<size_t>(numRows);
00098       std::fill(nnzPerRow.begin(), nnzPerRow.end(), 0);
00099       for (Teuchos::ArrayRCP<int>::const_iterator ri=rowinds.begin(); ri != rowinds.end(); ++ri) {
00100         // count each row index towards its row
00101         ++nnzPerRow[*ri-1];
00102       }
00103       if (isSymmetric) {
00104         // count each column toward the corresponding row as well
00105         for (int c=0; c < numCols; ++c) {
00106           // the diagonal was already counted; neglect it, if it exists
00107           for (int i=colptrs[c]-1; i != colptrs[c+1]-1; ++i) {
00108             if (rowinds[i] != c+1) {
00109               ++nnzPerRow[c];
00110               ++numNZ;
00111             }
00112           }
00113         }
00114       }
00115       // allocate/set new matrix data
00116       svals = Teuchos::arcp<Scalar>(numNZ);
00117       colinds = Teuchos::arcp<GlobalOrdinal>(numNZ);
00118       rowptrs = Teuchos::arcp<int>(numRows+1);
00119       rowptrs[0] = 0;
00120 #ifdef HAVE_TPETRA_DEBUG
00121       Teuchos::ArrayRCP<size_t> nnzPerRow_debug(nnzPerRow.size());
00122       std::copy(nnzPerRow.begin(), nnzPerRow.end(), nnzPerRow_debug.begin());
00123 #endif
00124       for (int j=1; j <= numRows; ++j) {
00125         rowptrs[j] = rowptrs[j-1] + nnzPerRow[j-1];
00126         nnzPerRow[j-1] = 0;
00127       }
00128       // translate from column-oriented to row-oriented
00129       for (int col=0; col<numCols; ++col) {
00130         for (int i=colptrs[col]-1; i != colptrs[col+1]-1; ++i) {
00131           const int row = rowinds[i]-1;
00132           // add entry to (row,col), with value dvals[i]
00133           const size_t entry = rowptrs[row] + nnzPerRow[row];
00134           svals[entry] = Teuchos::as<Scalar>(dvals[i]);
00135           colinds[entry] = Teuchos::as<GlobalOrdinal>(col);
00136           ++nnzPerRow[row];
00137           if (isSymmetric && row != col) {
00138             // add entry to (col,row), with value dvals[i]
00139             const size_t symentry = rowptrs[col] + nnzPerRow[col];
00140             svals[symentry] = Teuchos::as<Scalar>(dvals[i]);
00141             colinds[symentry] = Teuchos::as<GlobalOrdinal>(row);
00142             ++nnzPerRow[col];
00143           }
00144         }
00145       }
00146 #ifdef HAVE_TPETRA_DEBUG
00147       {
00148         bool isequal = true;
00149         Teuchos::ArrayRCP<size_t>::const_iterator it1, it2;
00150         for (it1 = nnzPerRow.begin(), it2 = nnzPerRow_debug.begin(); it1 != nnzPerRow.end(); ++it1, ++it2) {
00151           if (*it1 != *it2) {
00152             isequal = false; 
00153             break;
00154           }
00155         }
00156         TEST_FOR_EXCEPTION(!isequal || nnzPerRow.size() != nnzPerRow_debug.size(), std::logic_error,
00157             "Tpetra::Utils::readHBMatrix(): Logic error.");
00158       }
00159 #endif
00160     }
00161     // std::cout << "Matrix " << filename << " of type " << type << ": " << numRows << " by " << numCols << ", " << numNZ << " nonzeros" << std::endl;
00162   }
00163   // check for read errors
00164   broadcast(*comm,0,&fail);
00165   TEST_FOR_EXCEPTION(fail == 1, std::runtime_error, "Tpetra::Utils::readHBMatrix() can only read Real matrices.");
00166   // distribute global matrix info
00167   broadcast(*comm,0,&numRows);
00168   broadcast(*comm,0,&numCols);
00169   // create map with uniform partitioning
00170   if (rowMap == Teuchos::null) {
00171     rowMap = Teuchos::rcp(new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>((global_size_t)numRows,(GlobalOrdinal)0,comm,GloballyDistributed,node));
00172   }
00173   else {
00174     TEST_FOR_EXCEPTION( rowMap->getGlobalNumElements() != (global_size_t)numRows, std::runtime_error,
00175         "Tpetra::Utils::readHBMatrix(): specified map has incorrect number of elements.");
00176     TEST_FOR_EXCEPTION( rowMap->isDistributed() == false && comm->getSize() > 1, std::runtime_error,
00177         "Tpetra::Utils::readHBMatrix(): specified map is not distributed.");
00178   }
00179   Teuchos::ArrayRCP<size_t> myNNZ;
00180   if (rowMap->getNodeNumElements()) {
00181     myNNZ = Teuchos::arcp<size_t>(rowMap->getNodeNumElements());
00182   }
00183   if (myRank == 0) {
00184     size_t numRowsAlreadyDistributed = rowMap->getNodeNumElements();
00185     std::copy(nnzPerRow.begin(), nnzPerRow.begin()+numRowsAlreadyDistributed,myNNZ);
00186     for (int p=1; p < Teuchos::size(*comm); ++p) {
00187       size_t numRowsForP;
00188       Teuchos::receive(*comm,p,&numRowsForP);
00189       if (numRowsForP) {
00190         Teuchos::send<int,size_t>(*comm,numRowsForP,nnzPerRow(numRowsAlreadyDistributed,numRowsForP).getRawPtr(),p);
00191         numRowsAlreadyDistributed += numRowsForP;
00192       }
00193     }
00194   }
00195   else {
00196     const size_t numMyRows = rowMap->getNodeNumElements();
00197     Teuchos::send(*comm,numMyRows,0);
00198     if (numMyRows) {
00199       Teuchos::receive<int,size_t>(*comm,0,numMyRows,myNNZ(0,numMyRows).getRawPtr());
00200     }
00201   }
00202   nnzPerRow = Teuchos::null;
00203   // create column map
00204   Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap;
00205   if (numRows == numCols) {
00206     domMap = rowMap;
00207   }
00208   else {
00209     domMap = createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(numCols,comm,node);
00210   }
00211   // create with DynamicProfile, so that the fillComplete(DoOptimizeStorage) can do first-touch reallocation 
00212   A = rcp(new Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rowMap,myNNZ,Tpetra::DynamicProfile));
00213   // free this locally, A will keep it allocated as long as it is needed by A (up until allocation of nonzeros)
00214   myNNZ = Teuchos::null;
00215   if (myRank == 0 && numNZ > 0) {
00216     for (int r=0; r < numRows; ++r) {
00217       const size_t nnz = rowptrs[r+1] - rowptrs[r];
00218       if (nnz > 0) {
00219         Teuchos::ArrayView<const GlobalOrdinal> inds = colinds(rowptrs[r],nnz);
00220         Teuchos::ArrayView<const        Scalar> vals = svals(  rowptrs[r],nnz);
00221         A->insertGlobalValues(r, inds, vals);
00222       }
00223     }
00224   }
00225   // don't need these anymore
00226   colinds = Teuchos::null;
00227   svals   = Teuchos::null;
00228   rowptrs = Teuchos::null;
00229   A->fillComplete(domMap,rowMap,Tpetra::DoOptimizeStorage);
00230 }
00231 
00232 
00233 
00234 //
00235 // Explicit instantiation macro
00236 //
00237 // Must be expanded from within the Tpetra::Utils namespace!
00238 //
00239 
00240 #define TPETRA_MATRIXIO_INSTANT(SCALAR,LO,GO,NODE)                                                                                                      \
00241   template                                                                                                                                              \
00242   void                                                                                                                                                  \
00243   readHBMatrix<SCALAR,LO,GO,NODE,Kokkos::DefaultKernels<SCALAR,LO,NODE>::SparseOps>(                                                                    \
00244           const std::string &, const Teuchos::RCP<const Teuchos::Comm<int> > &, const Teuchos::RCP<NODE > &,                                            \
00245           Teuchos::RCP< CrsMatrix<SCALAR,LO,GO,NODE,Kokkos::DefaultKernels<SCALAR,LO,NODE>::SparseOps > > &,                                            \
00246           Teuchos::RCP< const Tpetra::Map<LO,GO,NODE> >);                                                                                               \
00247                                                                                                                                                         \
00248   template                                                                                                                                              \
00249   void                                                                                                                                                  \
00250   generateMatrix<SCALAR,LO,GO,NODE,Kokkos::DefaultKernels<SCALAR,LO,NODE>::SparseOps>(                                                                  \
00251           const Teuchos::RCP<Teuchos::ParameterList> &plist, const Teuchos::RCP<const Teuchos::Comm<int> > &, const Teuchos::RCP<NODE > &,              \
00252           Teuchos::RCP< CrsMatrix<SCALAR,LO,GO,NODE,Kokkos::DefaultKernels<SCALAR,LO,NODE>::SparseOps > > &);
00253 
00254 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines