FEI Version of the Day
fei_Matrix_Local.cpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 #include <fei_ParameterSet.hpp>
00045 #include "fei_Matrix_Local.hpp"
00046 #include "fei_Matrix_core.hpp"
00047 #include "fei_sstream.hpp"
00048 #include "fei_fstream.hpp"
00049 
00050 namespace fei {
00051 
00052 Matrix_Local::Matrix_Local(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00053                fei::SharedPtr<fei::SparseRowGraph> sparseRowGraph)
00054  : matrixGraph_(matrixGraph),
00055    sparseRowGraph_(sparseRowGraph),
00056    coefs_(sparseRowGraph->packedColumnIndices.size(), 0.0),
00057    stateChanged_(false),
00058    work_data1D_(),
00059    work_data2D_()
00060 {
00061 }
00062 
00063 Matrix_Local::~Matrix_Local()
00064 {
00065 }
00066 
00067 //----------------------------------------------------------------------------
00068 fei::SharedPtr<fei::Matrix>
00069 Matrix_Local::create_Matrix_Local(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00070                                 bool blockEntry)
00071 {
00072   fei::SharedPtr<fei::SparseRowGraph> srg =
00073     matrixGraph->createGraph(blockEntry, true);
00074   fei::SharedPtr<fei::Matrix> mat(new fei::Matrix_Local(matrixGraph, srg));
00075   return(mat);
00076 }
00077 
00078 //----------------------------------------------------------------------------
00079 const char*
00080 Matrix_Local::typeName()
00081 { return( "fei::Matrix_Local" ); }
00082 
00083 //----------------------------------------------------------------------------
00084 int
00085 Matrix_Local::parameters(const fei::ParameterSet& /*paramset*/)
00086 {
00087   return(0);
00088 }
00089 
00090 //----------------------------------------------------------------------------
00091 int
00092 Matrix_Local::parameters(int /*numParams*/, const char* const* /*paramStrings*/)
00093 {
00094   return(0);
00095 }
00096 
00097 fei::SharedPtr<fei::MatrixGraph>
00098 Matrix_Local::getMatrixGraph() const
00099 { return( matrixGraph_ ); }
00100 
00101 void
00102 Matrix_Local::setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00103 { matrixGraph_ = matrixGraph; }
00104 
00105 int
00106 Matrix_Local::getGlobalNumRows() const
00107 { return( sparseRowGraph_->rowNumbers.size() ); }
00108 
00109 int
00110 Matrix_Local::getLocalNumRows() const
00111 { return( getGlobalNumRows() ); }
00112 
00113 int
00114 Matrix_Local::getRowIndex(int rowNumber) const
00115 {
00116   int* rows = &(sparseRowGraph_->rowNumbers[0]);
00117   int numRows = getLocalNumRows();
00118   return( fei::binarySearch(rowNumber, rows, numRows) );
00119 }
00120 
00121 int
00122 Matrix_Local::getRowLength(int row, int& length) const
00123 {
00124   int idx = getRowIndex(row);
00125   if (idx < 0) return(idx);
00126 
00127   length = sparseRowGraph_->rowOffsets[idx+1] -
00128            sparseRowGraph_->rowOffsets[idx];
00129   return(0);
00130 }
00131 
00132 int
00133 Matrix_Local::putScalar(double scalar)
00134 {
00135   for(unsigned i=0; i<coefs_.size(); ++i) coefs_[i] = scalar;
00136   stateChanged_ = true;
00137   return(0);
00138 }
00139 
00140 int
00141 Matrix_Local::copyOutRow(int row, int len, double* coefs, int* indices) const
00142 {
00143   int idx = getRowIndex(row);
00144   if (idx < 0) return(idx);
00145 
00146   int offset = sparseRowGraph_->rowOffsets[idx];
00147   int length = sparseRowGraph_->rowOffsets[idx+1]-offset;
00148   if (length > len) length = len;
00149 
00150   for(int i=0; i<length; ++i) {
00151     indices[i] = sparseRowGraph_->packedColumnIndices[offset+i];
00152     coefs[i] = coefs_[offset+i];
00153   }
00154 
00155   return(0);
00156 }
00157 
00158 int
00159 Matrix_Local::giveToMatrix(int numRows, const int* rows,
00160                       int numCols, const int* cols,
00161                       const double* const* values,
00162                       bool sumInto,
00163                       int format)
00164 {
00165   if (numRows == 0 || numCols == 0) {
00166     return(0);
00167   }
00168 
00169   if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
00170     return(-1);
00171   }
00172 
00173   const double** myvalues = const_cast<const double**>(values);
00174   if (format != FEI_DENSE_ROW) {
00175     fei::Matrix_core::copyTransposeToWorkArrays(numRows, numCols, values,
00176                               work_data1D_, work_data2D_);
00177     myvalues = &work_data2D_[0];
00178   }
00179 
00180   for(int i=0; i<numRows; ++i) {
00181     int idx = getRowIndex(rows[i]);
00182     if (idx < 0) {
00183       throw std::runtime_error("fei::Matrix_Local::sumIn ERROR, row not found.");
00184     }
00185 
00186     int offset = sparseRowGraph_->rowOffsets[idx];
00187     int len = sparseRowGraph_->rowOffsets[idx+1] - offset;
00188 
00189     int* colInds = &(sparseRowGraph_->packedColumnIndices[offset]);
00190     double* coefs   = &(coefs_[offset]);
00191 
00192     for(int j=0; j<numCols; ++j) {
00193       int idx2 = fei::binarySearch(cols[j], colInds, len);
00194       if (idx2 < 0) {
00195         throw std::runtime_error("fei::Matrix_Local::sumIn ERROR, col not found.");
00196       }
00197 
00198       if (sumInto) {
00199         coefs[idx2] += myvalues[i][j];
00200       }
00201       else {
00202         coefs[idx2] = myvalues[i][j];
00203       }
00204     }
00205   }
00206 
00207   stateChanged_ = true;
00208   return(0);
00209 }
00210 
00211 int
00212 Matrix_Local::sumIn(int numRows, const int* rows,
00213                     int numCols, const int* cols,
00214                     const double* const* values,
00215                     int format)
00216 {
00217   return( giveToMatrix(numRows, rows, numCols, cols, values,
00218                        true, format) );
00219 }
00220 
00221 int
00222 Matrix_Local::copyIn(int numRows, const int* rows,
00223                        int numCols, const int* cols,
00224                        const double* const* values,
00225                       int format)
00226 {
00227   return( giveToMatrix(numRows, rows, numCols, cols, values,
00228                        false, format) );
00229 }
00230 
00231 int
00232 Matrix_Local::sumInFieldData(int fieldID,
00233                                int idType,
00234                                int rowID,
00235                                int colID,
00236                                const double* const* data,
00237                                int format)
00238 {
00239   fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
00240   fei::SharedPtr<fei::VectorSpace> cspace = matrixGraph_->getColSpace();
00241 
00242   int fieldSize = (int)rspace->getFieldSize(fieldID);
00243   std::vector<int> indices(2*fieldSize);
00244 
00245   rspace->getGlobalIndex(idType, rowID, fieldID, indices[0]);
00246   for(int i=1; i<fieldSize; ++i) {
00247     indices[i] = indices[0]+i;
00248   }
00249 
00250   cspace->getGlobalIndex(idType, colID, fieldID, indices[fieldSize]);
00251   for(int i=1; i<fieldSize; ++i) {
00252     indices[fieldSize+i] = indices[fieldSize]+i;
00253   }
00254 
00255   return( giveToMatrix(fieldSize, &indices[0], fieldSize, &indices[fieldSize],
00256                        data, true, format) );
00257 }
00258 
00259 int
00260 Matrix_Local::sumInFieldData(int fieldID,
00261                                int idType,
00262                                int rowID,
00263                                int colID,
00264                                const double* data,
00265                                int format)
00266 {
00267   fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
00268 
00269   int fieldSize = (int)rspace->getFieldSize(fieldID);
00270   std::vector<const double*> data2D(fieldSize);
00271 
00272   int offset = 0;
00273   for(int i=0; i<fieldSize; ++i) {
00274     data2D[i] = &data[offset];
00275     offset += fieldSize;
00276   }
00277 
00278   return( sumInFieldData(fieldID, idType, rowID, colID,
00279                          &data2D[0], format) );
00280 }
00281 
00282 int
00283 Matrix_Local::sumIn(int blockID, int connectivityID,
00284                     const double* const* values,
00285                     int format)
00286 {
00287   int numIndices = matrixGraph_->getConnectivityNumIndices(blockID);
00288   std::vector<int> indices(numIndices);
00289 
00290   matrixGraph_->getConnectivityIndices(blockID, connectivityID,
00291                                        numIndices, &indices[0], numIndices);
00292 
00293   return( giveToMatrix(numIndices, &indices[0], numIndices, &indices[0],
00294                        values, true, format) );
00295 }
00296 
00297 int
00298 Matrix_Local::globalAssemble()
00299 { return(0); }
00300 
00301 int
00302 Matrix_Local::multiply(fei::Vector* x,
00303                        fei::Vector* y)
00304 {
00305   FEI_COUT << "fei::Matrix_Local::multiply NOT IMPLEMENTED."<<FEI_ENDL;
00306   return(-1);
00307 }
00308 
00309 int
00310 Matrix_Local::gatherFromOverlap(bool accumulate)
00311 {
00312   (void)accumulate;
00313   return(0);
00314 }
00315 
00316 int
00317 Matrix_Local::writeToFile(const char* filename,
00318                           bool matrixMarketFormat)
00319 {
00320   fei::SharedPtr<fei::VectorSpace> vspace = matrixGraph_->getRowSpace();
00321 
00322   MPI_Comm comm = vspace->getCommunicator();
00323 
00324   int localProc = fei::localProc(comm);
00325 
00326   FEI_OSTRINGSTREAM osstr;
00327   osstr << filename << "." << localProc << ".mtx";
00328   std::string fullname = osstr.str();
00329 
00330   FEI_OFSTREAM ofstr(fullname.c_str(), IOS_OUT);
00331 
00332   return( writeToStream(ofstr, matrixMarketFormat) );
00333 }
00334 
00335 int
00336 Matrix_Local::writeToStream(FEI_OSTREAM& ostrm,
00337                             bool matrixMarketFormat)
00338 {
00339   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
00340 
00341   fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
00342   fei::SharedPtr<fei::VectorSpace> cspace = matrixGraph_->getColSpace();
00343 
00344   int numRows = getLocalNumRows();
00345   int numCols = cspace->getEqnNumbers().size();
00346   int nnz = coefs_.size();
00347 
00348   if (matrixMarketFormat) {
00349     ostrm << mmbanner << FEI_ENDL;
00350     ostrm << numRows << " " << numCols << " " << nnz << FEI_ENDL;
00351   }
00352   else {
00353     ostrm << numRows << " " << numCols << " "<< FEI_ENDL;
00354   }
00355 
00356   std::vector<int>& rowNumbers = sparseRowGraph_->rowNumbers;
00357   std::vector<int>& rowOffsets = sparseRowGraph_->rowOffsets;
00358   std::vector<int>& colIndices = sparseRowGraph_->packedColumnIndices;
00359 
00360   ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00361   ostrm.precision(13);
00362 
00363   int offset = 0;
00364   for(unsigned i=0; i<rowNumbers.size(); ++i) {
00365     int rowlen = rowOffsets[i+1]-rowOffsets[i];
00366 
00367     for(int j=0; j<rowlen; ++j) {
00368       if (matrixMarketFormat) {
00369         ostrm << rowNumbers[i]+1 << " " << colIndices[offset]+1
00370            << " " << coefs_[offset] << FEI_ENDL;
00371       }
00372       else {
00373         ostrm << rowNumbers[i] << " " << colIndices[offset]
00374            << " " << coefs_[offset] << FEI_ENDL;
00375       }
00376       ++offset;
00377     }
00378   }
00379 
00380   return(0);
00381 }
00382 
00383 bool
00384 Matrix_Local::usingBlockEntryStorage()
00385 { return( false ); }
00386 
00387 void
00388 Matrix_Local::markState()
00389 {
00390   stateChanged_ = false;
00391 }
00392 
00393 bool
00394 Matrix_Local::changedSinceMark()
00395 { return(stateChanged_); }
00396 
00397 const std::vector<int>&
00398 Matrix_Local::getRowNumbers() const
00399 { return( sparseRowGraph_->rowNumbers ); }
00400 
00401 const std::vector<int>&
00402 Matrix_Local::getRowOffsets() const
00403 { return( sparseRowGraph_->rowOffsets ); }
00404 
00405 const std::vector<int>&
00406 Matrix_Local::getColumnIndices() const
00407 { return( sparseRowGraph_->packedColumnIndices ); }
00408 
00409 const std::vector<double>&
00410 Matrix_Local::getCoefs() const
00411 { return( coefs_ ); }
00412 
00413 }//namespace fei
00414 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends