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 void
00310 Matrix_Local::setCommSizes()
00311 {
00312 }
00313 
00314 int
00315 Matrix_Local::gatherFromOverlap(bool accumulate)
00316 {
00317   (void)accumulate;
00318   return(0);
00319 }
00320 
00321 int
00322 Matrix_Local::writeToFile(const char* filename,
00323                           bool matrixMarketFormat)
00324 {
00325   fei::SharedPtr<fei::VectorSpace> vspace = matrixGraph_->getRowSpace();
00326 
00327   MPI_Comm comm = vspace->getCommunicator();
00328 
00329   int localProc = fei::localProc(comm);
00330 
00331   FEI_OSTRINGSTREAM osstr;
00332   osstr << filename << "." << localProc << ".mtx";
00333   std::string fullname = osstr.str();
00334 
00335   FEI_OFSTREAM ofstr(fullname.c_str(), IOS_OUT);
00336 
00337   return( writeToStream(ofstr, matrixMarketFormat) );
00338 }
00339 
00340 int
00341 Matrix_Local::writeToStream(FEI_OSTREAM& ostrm,
00342                             bool matrixMarketFormat)
00343 {
00344   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
00345 
00346   fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
00347   fei::SharedPtr<fei::VectorSpace> cspace = matrixGraph_->getColSpace();
00348 
00349   int numRows = getLocalNumRows();
00350   int numCols = cspace->getEqnNumbers().size();
00351   int nnz = coefs_.size();
00352 
00353   if (matrixMarketFormat) {
00354     ostrm << mmbanner << FEI_ENDL;
00355     ostrm << numRows << " " << numCols << " " << nnz << FEI_ENDL;
00356   }
00357   else {
00358     ostrm << numRows << " " << numCols << " "<< FEI_ENDL;
00359   }
00360 
00361   std::vector<int>& rowNumbers = sparseRowGraph_->rowNumbers;
00362   std::vector<int>& rowOffsets = sparseRowGraph_->rowOffsets;
00363   std::vector<int>& colIndices = sparseRowGraph_->packedColumnIndices;
00364 
00365   ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00366   ostrm.precision(13);
00367 
00368   int offset = 0;
00369   for(unsigned i=0; i<rowNumbers.size(); ++i) {
00370     int rowlen = rowOffsets[i+1]-rowOffsets[i];
00371 
00372     for(int j=0; j<rowlen; ++j) {
00373       if (matrixMarketFormat) {
00374         ostrm << rowNumbers[i]+1 << " " << colIndices[offset]+1
00375            << " " << coefs_[offset] << FEI_ENDL;
00376       }
00377       else {
00378         ostrm << rowNumbers[i] << " " << colIndices[offset]
00379            << " " << coefs_[offset] << FEI_ENDL;
00380       }
00381       ++offset;
00382     }
00383   }
00384 
00385   return(0);
00386 }
00387 
00388 bool
00389 Matrix_Local::usingBlockEntryStorage()
00390 { return( false ); }
00391 
00392 void
00393 Matrix_Local::markState()
00394 {
00395   stateChanged_ = false;
00396 }
00397 
00398 bool
00399 Matrix_Local::changedSinceMark()
00400 { return(stateChanged_); }
00401 
00402 const std::vector<int>&
00403 Matrix_Local::getRowNumbers() const
00404 { return( sparseRowGraph_->rowNumbers ); }
00405 
00406 const std::vector<int>&
00407 Matrix_Local::getRowOffsets() const
00408 { return( sparseRowGraph_->rowOffsets ); }
00409 
00410 const std::vector<int>&
00411 Matrix_Local::getColumnIndices() const
00412 { return( sparseRowGraph_->packedColumnIndices ); }
00413 
00414 const std::vector<double>&
00415 Matrix_Local::getCoefs() const
00416 { return( coefs_ ); }
00417 
00418 }//namespace fei
00419 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends