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