|
Tpetra Matrix/Vector Services
Version of the Day
|
|
Tpetra Matrix/Vector Services
Version of the Day
|
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> ¶ms) 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
1.7.6.1
