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