|
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 // create with DynamicProfile, so that the fillComplete(DoOptimizeStorage) can do first-touch reallocation 00072 A = rcp(new Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rowMap,7,Tpetra::DynamicProfile)); 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(DoOptimizeStorage); 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 { 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 //Teuchos::ArrayRCP<char> type; 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 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 size_t 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 // create with DynamicProfile, so that the fillComplete(DoOptimizeStorage) can do first-touch reallocation 00255 A = rcp(new Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rowMap,myNNZ,Tpetra::DynamicProfile)); 00256 // free this locally, A will keep it allocated as long as it is needed by A (up until allocation of nonzeros) 00257 myNNZ = Teuchos::null; 00258 if (myRank == 0 && numNZ > 0) { 00259 for (int r=0; r < numRows; ++r) { 00260 const size_t nnz = rowptrs[r+1] - rowptrs[r]; 00261 if (nnz > 0) { 00262 Teuchos::ArrayView<const GlobalOrdinal> inds = colinds(rowptrs[r],nnz); 00263 Teuchos::ArrayView<const Scalar> vals = svals( rowptrs[r],nnz); 00264 A->insertGlobalValues(r, inds, vals); 00265 } 00266 } 00267 } 00268 // don't need these anymore 00269 colinds = Teuchos::null; 00270 svals = Teuchos::null; 00271 rowptrs = Teuchos::null; 00272 A->fillComplete(domMap,rowMap,Tpetra::DoOptimizeStorage); 00273 } 00274 00275 00276 00277 // 00278 // Explicit instantiation macro 00279 // 00280 // Must be expanded from within the Tpetra::Utils namespace! 00281 // 00282 00283 #define TPETRA_MATRIXIO_INSTANT(SCALAR,LO,GO,NODE) \ 00284 template \ 00285 void \ 00286 readHBMatrix<SCALAR,LO,GO,NODE,Kokkos::DefaultKernels<SCALAR,LO,NODE>::SparseOps>( \ 00287 const std::string &, const Teuchos::RCP<const Teuchos::Comm<int> > &, const Teuchos::RCP<NODE > &, \ 00288 Teuchos::RCP< CrsMatrix<SCALAR,LO,GO,NODE,Kokkos::DefaultKernels<SCALAR,LO,NODE>::SparseOps > > &, \ 00289 Teuchos::RCP< const Tpetra::Map<LO,GO,NODE> >); \ 00290 \ 00291 template \ 00292 void \ 00293 generateMatrix<SCALAR,LO,GO,NODE,Kokkos::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,Kokkos::DefaultKernels<SCALAR,LO,NODE>::SparseOps > > &); 00296 00297 #endif
1.7.4
