FEI Version of the Day
fei_Matrix_Impl.hpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #ifndef _fei_Matrix_Impl_hpp_
00010 #define _fei_Matrix_Impl_hpp_
00011 
00012 #include <fei_fwd.hpp>
00013 #include <fei_mpi.h>
00014 #include <fei_defs.h>
00015 #include "fei_iostream.hpp"
00016 #include "fei_fstream.hpp"
00017 #include "fei_sstream.hpp"
00018 #include <fei_ArrayUtils.hpp>
00019 #include <fei_MatrixTraits.hpp>
00020 #include <fei_MatrixTraits_LinProbMgr.hpp>
00021 #include <fei_MatrixTraits_LinSysCore.hpp>
00022 #include <fei_MatrixTraits_FEData.hpp>
00023 #include <fei_MatrixTraits_FillableMat.hpp>
00024 #include <fei_FillableVec.hpp>
00025 #include <fei_FillableMat.hpp>
00026 
00027 #include <snl_fei_FEMatrixTraits.hpp>
00028 #include <snl_fei_FEMatrixTraits_FED.hpp>
00029 #include <snl_fei_BlockMatrixTraits.hpp>
00030 #include <fei_Matrix.hpp>
00031 #include <fei_MatrixGraph.hpp>
00032 #include <fei_Matrix_core.hpp>
00033 #include <snl_fei_Utils.hpp>
00034 
00035 #undef fei_file
00036 #define fei_file "fei_Matrix_Impl.hpp"
00037 #include <fei_ErrMacros.hpp>
00038 
00039 namespace fei {
00040 
00053   template<typename T>
00054   class Matrix_Impl : public fei::Matrix, public fei::Matrix_core {
00055   public:
00057     Matrix_Impl(fei::SharedPtr<T> matrix,
00058                fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00059                 int numLocalEqns);
00060 
00062     virtual ~Matrix_Impl();
00063 
00067     const char* typeName()
00068       {
00069         if (haveBlockMatrix()) {
00070           return(snl_fei::BlockMatrixTraits<T>::typeName());
00071         }
00072         else if (haveFEMatrix()) {
00073           return(snl_fei::FEMatrixTraits<T>::typeName());
00074         }
00075         else {
00076           return(fei::MatrixTraits<T>::typeName());
00077         }
00078       }
00079 
00082     int parameters(const fei::ParameterSet& paramset);
00083 
00087     fei::SharedPtr<T> getMatrix() { return( matrix_ ); }
00088     const fei::SharedPtr<T> getMatrix() const { return( matrix_ ); }
00089 
00090     fei::SharedPtr<fei::MatrixGraph> getMatrixGraph() const
00091       {return( Matrix_core::getMatrixGraph() ); }
00092 
00094     void setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph);
00095 
00098     int getGlobalNumRows() const
00099       {
00100         if ((int)(globalOffsets().size()) < numProcs()+1) return(-1);
00101         return(globalOffsets()[numProcs()]);
00102       }
00103 
00106     int getLocalNumRows() const
00107       {
00108         return(lastLocalOffset() - firstLocalOffset() + 1);
00109       }
00110 
00112     int putScalar(double scalar);
00113 
00119     int getRowLength(int row, int& length) const;
00120 
00130     int copyOutRow(int row, int len, double* coefs, int* indices) const;
00131 
00143     int sumIn(int numRows, const int* rows,
00144               int numCols, const int* cols,
00145               const double* const* values,
00146               int format=0);
00147 
00159     int copyIn(int numRows, const int* rows,
00160                int numCols, const int* cols,
00161                const double* const* values,
00162                int format=0);
00163 
00179     int sumInFieldData(int fieldID,
00180                        int idType,
00181                        int rowID,
00182                        int colID,
00183                        const double* const* data,
00184                        int format=0);
00185 
00203     int sumInFieldData(int fieldID,
00204                        int idType,
00205                        int rowID,
00206                        int colID,
00207                        const double* data,
00208                        int format=0);
00209 
00219     int sumIn(int blockID, int connectivityID,
00220               const double* const* values,
00221               int format=0);
00222 
00227     int globalAssemble();
00228 
00231     int multiply(fei::Vector* x,
00232                  fei::Vector* y);
00233 
00239     int gatherFromOverlap(bool accumulate = true);
00240 
00242     int writeToFile(const char* filename,
00243                     bool matrixMarketFormat=true);
00244 
00247     int writeToStream(FEI_OSTREAM& ostrm,
00248                       bool matrixMarketFormat=true);
00249 
00250     bool usingBlockEntryStorage()
00251       {
00252         return(haveBlockMatrix());
00253       }
00254 
00256     int giveToUnderlyingMatrix(int numRows, const int* rows,
00257                                int numCols, const int* cols,
00258                                const double* const* values,
00259                                bool sumInto,
00260                                int format);
00261 
00263     int giveToUnderlyingBlockMatrix(int row,
00264                                     int rowDim,
00265                                     int numCols,
00266                                     const int* cols,
00267                                     const int* LDAs,
00268                                     const int* colDims,
00269                                     const double* const* values,
00270                                     bool sumInto);
00271 
00272     void markState();
00273 
00274     bool changedSinceMark();
00275 
00276   private:
00277     int giveToMatrix(int numRows, const int* rows,
00278                      int numCols, const int* cols,
00279                      const double* const* values,
00280                      bool sumInto,
00281                      int format);
00282  
00283     int giveToBlockMatrix(int numRows, const int* rows,
00284                           int numCols, const int* cols,
00285                           const double* const* values,
00286                           bool sumInto);
00287 
00288     fei::SharedPtr<T> matrix_;
00289     bool globalAssembleCalled_;
00290     bool changedSinceMark_;
00291     std::string dbgprefix_;
00292   };//class Matrix_Impl
00293 }//namespace fei
00294 
00295 //----------------------------------------------------------------------------
00296 template<typename T>
00297 inline int fei::Matrix_Impl<T>::globalAssemble()
00298 {
00299   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00300     FEI_OSTREAM& os = *output_stream_;
00301     os << dbgprefix_<<"globalAssemble"<<FEI_ENDL;
00302   }
00303 
00304   globalAssembleCalled_ = true;
00305 
00306   if (haveBlockMatrix()) {
00307     return( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
00308   }
00309   return( fei::MatrixTraits<T>::globalAssemble(matrix_.get()) );
00310 }
00311 
00312 //----------------------------------------------------------------------------
00313 template<typename T>
00314 inline void fei::Matrix_Impl<T>::setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00315 {
00316   Matrix_core::setMatrixGraph(matrixGraph);
00317 }
00318 
00319 //----------------------------------------------------------------------------
00320 template<typename T>
00321 inline void fei::Matrix_Impl<T>::markState()
00322 {
00323   changedSinceMark_ = false;
00324 }
00325 
00326 //----------------------------------------------------------------------------
00327 template<typename T>
00328 inline bool fei::Matrix_Impl<T>::changedSinceMark()
00329 {
00330   return(changedSinceMark_);
00331 }
00332 
00333 //----------------------------------------------------------------------------
00334 template<typename T>
00335 int fei::Matrix_Impl<T>::giveToUnderlyingMatrix(int numRows, const int* rows,
00336                                                int numCols, const int* cols,
00337                                                const double* const* values,
00338                                                bool sumInto,
00339                                                int format)
00340 {
00341   if (format != FEI_DENSE_ROW) {
00342     ERReturn(-1);
00343   }
00344 
00345   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != 0) {
00346     FEI_OSTREAM& os = *output_stream_;
00347     for(int i=0; i<numRows; ++i) {
00348       os << dbgprefix_<<"giveToUnderlyingMatrix ";
00349       for(int j=0; j<numCols; ++j) {
00350         os << "("<<rows[i]<<","<<cols[j]<<","<<values[i][j]<<") ";
00351       }
00352       os << FEI_ENDL;
00353     }
00354   }
00355 
00356   int err = fei::MatrixTraits<T>::putValuesIn(matrix_.get(), numRows, rows,
00357                                         numCols, cols, values, sumInto);
00358   if (err != 0) {
00359     return(err);
00360   }
00361 
00362   changedSinceMark_ = true;
00363   return(0);
00364 }
00365 
00366 //----------------------------------------------------------------------------
00367 template<typename T>
00368 int fei::Matrix_Impl<T>::giveToUnderlyingBlockMatrix(int row,
00369                                                     int rowDim,
00370                                                     int numCols,
00371                                                     const int* cols,
00372                                                     const int* LDAs,
00373                                                     const int* colDims,
00374                                                     const double* const* values,
00375                                                     bool sumInto)
00376 {
00377   if (sumInto) {
00378     if ( snl_fei::BlockMatrixTraits<T>::sumIn(matrix_.get(),
00379                                          row, rowDim, numCols, cols,
00380                                          LDAs, colDims, values) != 0) {
00381       ERReturn(-1);
00382     }
00383   }
00384   else {
00385     if ( snl_fei::BlockMatrixTraits<T>::copyIn(matrix_.get(),
00386                                           row, rowDim, numCols, cols,
00387                                           LDAs, colDims, values) != 0) {
00388       ERReturn(-1);
00389     }
00390   }
00391 
00392   changedSinceMark_ = true;
00393 
00394   return(0);
00395 }
00396 
00397 #include <fei_macros.hpp>
00398 #include <fei_chk_mpi.hpp>
00399 
00400 #include <fei_TemplateUtils.hpp>
00401 
00402 #include <fei_Pattern.hpp>
00403 #include <fei_VectorSpace.hpp>
00404 
00405 #include <fei_ConnectivityBlock.hpp>
00406 
00407 #include <snl_fei_PointBlockMap.hpp>
00408 
00409 #include <fei_Record.hpp>
00410 
00411 #include <snl_fei_Utils.hpp>
00412 
00413 //----------------------------------------------------------------------------
00414 template<typename T>
00415 fei::Matrix_Impl<T>::Matrix_Impl(fei::SharedPtr<T> matrix,
00416                            fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00417                                 int numLocalEqns)
00418   : Matrix_core(matrixGraph, numLocalEqns),
00419     matrix_(matrix),
00420     globalAssembleCalled_(false),
00421     changedSinceMark_(true),
00422     dbgprefix_("MatImpl: ")
00423 {
00424   if (strcmp(snl_fei::FEMatrixTraits<T>::typeName(), "unsupported")) {
00425     setFEMatrix(true);
00426   }
00427   else {
00428     setFEMatrix(false);
00429   }
00430 
00431   if (strcmp(snl_fei::BlockMatrixTraits<T>::typeName(), "unsupported")) {
00432     setBlockMatrix(true);
00433   }
00434   else {
00435     setBlockMatrix(false);
00436   }
00437 }
00438 
00439 //----------------------------------------------------------------------------
00440 template<typename T>
00441 fei::Matrix_Impl<T>::~Matrix_Impl()
00442 {
00443 }
00444 
00445 //----------------------------------------------------------------------------
00446 template<typename T>
00447 int fei::Matrix_Impl<T>::parameters(const fei::ParameterSet& paramset)
00448 {
00449   Matrix_core::parameters(paramset);
00450   return 0;
00451 }
00452 
00453 //----------------------------------------------------------------------------
00454 template<typename T>
00455 int fei::Matrix_Impl<T>::getRowLength(int row, int& length) const
00456 {
00457   if (haveBlockMatrix()) {
00458     return( snl_fei::BlockMatrixTraits<T>::getPointRowLength(matrix_.get(), row, length) );
00459   }
00460   else {
00461     int code = fei::MatrixTraits<T>::getRowLength(matrix_.get(), row, length);
00462     if (code < 0) {
00463       //matrix row probably not local. maybe it's shared.
00464       int proc = getOwnerProc(row);
00465       const FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
00466       if (remote_mat->hasRow(row)) {
00467         const FillableVec* row_entries = remote_mat->getRow(row);
00468         length = row_entries->size();
00469       }
00470       else {
00471         length = 0;
00472       }
00473     }
00474   }
00475   return 0;
00476 }
00477 
00478 //----------------------------------------------------------------------------
00479 template<typename T>
00480 int fei::Matrix_Impl<T>::putScalar(double scalar)
00481 {
00482   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != 0) {
00483     FEI_OSTREAM& os = *output_stream_;
00484     os << dbgprefix_<<"putScalar("<<scalar<<")"<<FEI_ENDL;
00485   }
00486 
00487   if (haveFEMatrix()) {
00488     if (scalar != 0.0) return(-1);
00489     CHK_ERR( snl_fei::FEMatrixTraits<T>::reset(matrix_.get()) );
00490   }
00491   else if (haveBlockMatrix()) {
00492     if (globalAssembleCalled_ == true) {
00493       CHK_ERR( snl_fei::BlockMatrixTraits<T>::putScalar(matrix_.get(), scalar) );
00494     }
00495   }
00496   else {
00497     CHK_ERR( fei::MatrixTraits<T>::setValues(matrix_.get(), scalar) );
00498   }
00499 
00500   putScalar_remotelyOwned(scalar);
00501 
00502   changedSinceMark_ = true;
00503 
00504   return(0);
00505 }
00506 
00507 //----------------------------------------------------------------------------
00508 template<typename T>
00509 int fei::Matrix_Impl<T>::copyOutRow(int row, int len,
00510                                    double* coefs, int* indices) const
00511 {
00512   if (len < 1) {
00513     return 0;
00514   }
00515 
00516   if (haveBlockMatrix()) {
00517     int dummy;
00518     return( snl_fei::BlockMatrixTraits<T>::copyOutPointRow(matrix_.get(), firstLocalOffset(),
00519                                                   row, len,
00520                                                   coefs, indices, dummy));
00521   }
00522   else {
00523     int code = fei::MatrixTraits<T>::copyOutRow(matrix_.get(), row, len,
00524                                         coefs, indices);
00525     if (code < 0) {
00526       int proc = getOwnerProc(row);
00527       const FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
00528       if (remote_mat->hasRow(row)) {
00529         const FillableVec* row_entries = remote_mat->getRow(row);
00530         FillableVec::const_iterator it=row_entries->begin(); 
00531         for(int i=0; it!=row_entries->end(); ++it, ++i) {
00532           indices[i] = it->first;
00533           coefs[i] = it->second;
00534         }
00535       }
00536     }
00537   }
00538   return 0;
00539 }
00540 
00541 //----------------------------------------------------------------------------
00542 template<typename T>
00543 int fei::Matrix_Impl<T>::sumIn(int numRows, const int* rows,
00544                               int numCols, const int* cols,
00545                               const double* const* values,
00546                               int format)
00547 {
00548   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00549     FEI_OSTREAM& os = *output_stream_;
00550     os << dbgprefix_<<"sumIn"<<FEI_ENDL;
00551     if (output_level_ >= fei::FULL_LOGS) {
00552       for(int i=0; i<numRows; ++i) {
00553         for(int j=0; j<numCols; ++j) {
00554           os << "("<<rows[i]<<","<<cols[j]<<","<<values[i][j]<<") ";
00555         }
00556         os << FEI_ENDL;
00557       }
00558     }
00559   }
00560 
00561   return( giveToMatrix( numRows, rows, numCols, cols, values, true, format) );
00562 }
00563 
00564 //----------------------------------------------------------------------------
00565 template<typename T>
00566 int fei::Matrix_Impl<T>::copyIn(int numRows, const int* rows,
00567                                int numCols, const int* cols,
00568                                const double* const* values,
00569                                int format)
00570 {
00571   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
00572     FEI_OSTREAM& os = *output_stream_;
00573     os << "copyIn"<<FEI_ENDL;
00574     if (output_level_ >= fei::FULL_LOGS) {
00575       for(int i=0; i<numRows; ++i) {
00576         for(int j=0; j<numCols; ++j) {
00577           os << "("<<rows[i]<<","<<cols[j]<<","<<values[i][j]<<") ";
00578         }
00579         os << FEI_ENDL;
00580       }
00581     }
00582   }
00583   return( giveToMatrix( numRows, rows, numCols, cols, values, false, format) );
00584 }
00585 
00586 //----------------------------------------------------------------------------
00587 template<typename T>
00588 int fei::Matrix_Impl<T>::sumInFieldData(int fieldID,
00589                                        int idType,
00590                                        int rowID,
00591                                        int colID,
00592                                        const double* const* data,
00593                                        int format)
00594 {
00595   fei::SharedPtr<fei::VectorSpace> vspace = vecSpace();
00596 
00597   int fieldSize = vspace->getFieldSize(fieldID);
00598 
00599   if (fieldSize <= 0) ERReturn(-1);
00600 
00601   work_indices_.resize(fieldSize*2);
00602   int* indicesPtr = &work_indices_[0];
00603   int i;
00604 
00605   CHK_ERR( vspace->getGlobalIndices(1, &rowID, idType, fieldID, indicesPtr));
00606   for(i=1; i<fieldSize; ++i) {
00607     indicesPtr[i] = indicesPtr[i-1]+1;
00608   }
00609 
00610   CHK_ERR( vspace->getGlobalIndices(1, &colID, idType, fieldID,
00611                                         &(indicesPtr[fieldSize])) );
00612   for(i=fieldSize+1; i<2*fieldSize; ++i) {
00613     indicesPtr[i] = indicesPtr[i-1]+1;
00614   }
00615 
00616   CHK_ERR( sumIn(fieldSize, indicesPtr, fieldSize, &(indicesPtr[fieldSize]),
00617                  data, format) );
00618 
00619   return(0);
00620 }
00621 
00622 //----------------------------------------------------------------------------
00623 template<typename T>
00624 int fei::Matrix_Impl<T>::sumInFieldData(int fieldID,
00625                                        int idType,
00626                                        int rowID,
00627                                        int colID,
00628                                        const double* data,
00629                                        int format)
00630 {
00631   fei::SharedPtr<fei::VectorSpace> vspace = vecSpace();
00632 
00633   int fieldSize = vspace->getFieldSize(fieldID);
00634 
00635   if (fieldSize <= 0) ERReturn(-1);
00636 
00637   work_data2D_.resize(fieldSize);
00638 
00639   const double** data2DPtr = &work_data2D_[0];
00640   for(int i=0; i<fieldSize; ++i) {
00641     data2DPtr[i] = &(data[i*fieldSize]);
00642   }
00643 
00644   CHK_ERR( sumInFieldData(fieldID, idType, rowID, colID, data2DPtr, format) );
00645 
00646   return(0);
00647 }
00648 
00649 //----------------------------------------------------------------------------
00650 template<typename T>
00651 int fei::Matrix_Impl<T>::sumIn(int blockID, int connectivityID,
00652                               const double* const* values, int format)
00653 {
00654   if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
00655     FEI_OSTREAM& os = *output_stream_;
00656     os << dbgprefix_<<"sumIn blkID=" << blockID
00657        << ", connID=" << connectivityID << FEI_ENDL;
00658   }
00659 
00660   fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
00661 
00662   if (mgraph.get() == NULL) ERReturn(-1);
00663 
00664   const fei::ConnectivityBlock* cblock = mgraph->getConnectivityBlock(blockID);
00665   if (cblock==NULL) {
00666     FEI_OSTRINGSTREAM osstr;
00667     osstr << "fei::Matrix_Impl::sumIn ERROR, unable to "
00668           << "look up connectivity-block with ID "<<blockID;
00669     throw std::runtime_error(osstr.str());
00670   }
00671 
00672   bool symmetric = cblock->isSymmetric();
00673   const fei::Pattern* pattern = cblock->getRowPattern();
00674   const fei::Pattern* colpattern = symmetric ? NULL : cblock->getColPattern();
00675   const int* rowConn = cblock->getRowConnectivity(connectivityID);
00676 
00677   fei::SharedPtr<fei::VectorSpace> rspace = mgraph->getRowSpace();
00678   rspace->getGlobalIndicesL(pattern, rowConn, work_indices2_);
00679 
00680   int numRowIndices = work_indices2_.size();
00681   int* rowIndices = &work_indices2_[0];
00682 
00683   if (haveFEMatrix() || haveBlockMatrix()) {
00684     FieldDofMap<int>& fdofmap = rspace->getFieldDofMap();
00685     const std::map<int,int>& connIDs = cblock->getConnectivityIDs();
00686     std::map<int,int>::const_iterator
00687       iter = connIDs.find(connectivityID);
00688     if (iter == connIDs.end()) ERReturn(-1);
00689     int connOffset = iter->second;
00690     int numIDs = pattern->getNumIDs();
00691     const int* indicesPerID = pattern->getNumIndicesPerID();
00692     const int* fieldsPerID = pattern->getNumFieldsPerID();
00693     const int* fieldIDs = pattern->getFieldIDs();
00694 
00695     int numDofs = 0;
00696     for(int ii=0; ii<numIDs; ++ii) {
00697       numDofs += indicesPerID[ii];
00698     }
00699 
00700     const int* numIndicesPerID = pattern->getNumIndicesPerID();
00701     work_indices_.resize(numIDs+numDofs);
00702     int i, *nodeNumbers = &work_indices_[0];
00703     int* dof_ids = nodeNumbers+numIDs;
00704 
00705     int nodeType = 0;
00706     snl_fei::RecordCollection* records = NULL;
00707     rspace->getRecordCollection(nodeType, records);
00708     int foffset = 0;
00709     int doffset = 0;
00710     for(i=0; i<numIDs; ++i) {
00711       fei::Record<int>* rec = records->getRecordWithLocalID(rowConn[i]);
00712       nodeNumbers[i] = rec->getNumber();
00713       for(int ii=0; ii<fieldsPerID[i]; ++ii) {
00714         int fieldSize = rspace->getFieldSize(fieldIDs[foffset]);
00715         int dof_id = fdofmap.get_dof_id(fieldIDs[foffset++], 0);
00716         for(int j=0; j<fieldSize; ++j) {
00717           dof_ids[doffset++] = dof_id + j;
00718         }
00719       }
00720     }
00721 
00722     if (haveFEMatrix()) {
00723       CHK_ERR( snl_fei::FEMatrixTraits<T>::sumInElemMatrix(matrix_.get(), blockID, connOffset,
00724                                                   numIDs, nodeNumbers,
00725                                                   numIndicesPerID, dof_ids, values) );
00726       changedSinceMark_ = true;
00727     }
00728 
00729     if (haveBlockMatrix()) {
00730       if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL &&
00731           format != FEI_BLOCK_DIAGONAL_ROW) {
00732         fei::console_out() << "fei::Matrix_Impl::sumIn ERROR, for block-matrix, format must"
00733                  << " be FEI_DENSE_ROW or FEI_DENSE_COL."<<FEI_ENDL;
00734         ERReturn(-1);
00735       }
00736 
00737       int numPtIndices = pattern->getNumIndices();
00738       int* ptIndices = &work_indices2_[0];
00739 
00740       int numPtColIndices = symmetric ? numPtIndices : colpattern->getNumIndices();
00741 
00742       int len = numPtIndices*numPtColIndices;
00743 
00744       if (format == FEI_BLOCK_DIAGONAL_ROW) {
00745         len = 0;
00746         for(i=0; i<numIDs; ++i) {
00747           len += numIndicesPerID[i]*numIndicesPerID[i];
00748         }
00749       }
00750 
00751       double* ccvalues = new double[len];
00752       //Ouch! Data copy! Remember to optimize this later...
00753       if (format == FEI_BLOCK_DIAGONAL_ROW) {
00754         snl_fei::copy2DBlockDiagToColumnContig(numIDs, numIndicesPerID,
00755                                                values, format, ccvalues);
00756       }
00757       else {
00758         snl_fei::copy2DToColumnContig(numPtIndices, numPtColIndices, values, format,
00759                                       ccvalues);
00760       }
00761 
00762       int ptRowOffset = 0;
00763       for(i=0; i<numIDs; ++i) {
00764 
00765         fei::Record<int>* record = records->getRecordWithLocalID(rowConn[i]);
00766         if (record->getOwnerProc() == localProc()) {
00767 
00768           int numColIDs = numIDs;
00769           int* colNodeNums = nodeNumbers;
00770           const int* colDims = numIndicesPerID;
00771           int LDA = numPtColIndices;
00772           if (format == FEI_BLOCK_DIAGONAL_ROW) {
00773             numColIDs = 1;
00774             colNodeNums = &(nodeNumbers[i]);
00775             colDims = &(numIndicesPerID[i]);
00776             LDA = numIndicesPerID[i];
00777           }
00778 
00779           CHK_ERR( snl_fei::BlockMatrixTraits<T>::sumIn(matrix_.get(),
00780                                                nodeNumbers[i],
00781                                                numIndicesPerID[i],
00782                                                numColIDs, colNodeNums,
00783                                                colDims, LDA,
00784                                                &(ccvalues[ptRowOffset])) );
00785           changedSinceMark_ = true;
00786         }
00787         else {
00788           if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
00789             FEI_OSTREAM& os = *output_stream_;
00790             for(int ii=0; ii<numIndicesPerID[i]; ++ii) {
00791               os << "#   remote pt-row " << ptIndices[ptRowOffset]+ii <<" ";
00792               for(int jj=0; jj<numPtIndices; ++jj) {
00793                 os << "("<<ptIndices[jj]<<","<<values[ptRowOffset+ii][jj]<<") ";
00794               }
00795               os << FEI_ENDL;
00796             }
00797           }
00798 
00799           for(int ii=0; ii<numIndicesPerID[i]; ++ii) {
00800             int p=eqnComm_->getOwnerProc(ptIndices[ptRowOffset]+ii);
00801             FillableMat* remote_mat = getRemotelyOwnedMatrix(p);
00802             remote_mat->sumInRow(ptIndices[ptRowOffset]+ii,
00803                                       ptIndices,
00804                                       values[ptRowOffset+ii],
00805                                       numPtIndices);
00806             changedSinceMark_ = true;
00807           }
00808         }
00809 
00810         ptRowOffset += numIndicesPerID[i];
00811       }
00812 
00813       delete [] ccvalues;
00814     }
00815 
00816     return(0);
00817   }
00818 
00819   int numColIndices = symmetric ? numRowIndices : colpattern->getNumIndices();
00820   int* colIndices = rowIndices;
00821   const int* colConn = NULL;
00822 
00823   if (!symmetric) {
00824     colConn = cblock->getColConnectivity(connectivityID);
00825     mgraph->getColSpace()->getGlobalIndicesL(colpattern,
00826                                             colConn, work_indices_);
00827     colIndices = &work_indices_[0];
00828   }
00829 
00830   if (symmetric) {
00831     if (format == FEI_DENSE_ROW || format == FEI_DENSE_COL) {
00832       CHK_ERR( sumIn(numRowIndices, rowIndices, numRowIndices, rowIndices,
00833                      values, format) );
00834     }
00835     else if (format == FEI_BLOCK_DIAGONAL_ROW) {
00836       int totalnumfields = pattern->getTotalNumFields();
00837       const int* fieldIDs = pattern->getFieldIDs();
00838       fei::SharedPtr<fei::VectorSpace> rowspace = mgraph->getRowSpace();
00839       fei::VectorSpace* rowspaceptr = rowspace.get();
00840       int ioffset = 0;
00841       for(int i=0; i<totalnumfields; ++i) {
00842         int fieldsize = rowspaceptr->getFieldSize(fieldIDs[i]);
00843         if (ioffset+fieldsize > numRowIndices) {
00844           FEI_OSTRINGSTREAM osstr;
00845           osstr<<"snl_fei::sumIn, format=FEI_BLOCK_DIAGONAL_ROW, block-sizes"
00846                << " not consistent with total num-indices."<<FEI_ENDL;
00847           throw std::runtime_error(osstr.str());
00848         }
00849 
00850         CHK_ERR( sumIn(fieldsize, &(rowIndices[ioffset]),
00851                        fieldsize, &(rowIndices[ioffset]),
00852                        &(values[ioffset]), FEI_DENSE_ROW) );
00853         ioffset += fieldsize;
00854       }
00855     }
00856     else {
00857       FEI_OSTRINGSTREAM osstr;
00858       osstr << "fei::Matrix_Impl::sumIn, format="<<format<<" not supported."
00859             << FEI_ENDL;
00860       throw std::runtime_error(osstr.str());
00861     }
00862   }
00863   else {
00864     if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
00865       FEI_OSTRINGSTREAM osstr;
00866       osstr << "fei::Matrix_Impl::sumIn, format="<<format<<" not valid with"
00867             << " un-symmetric matrix contributions."<<FEI_ENDL;
00868       throw std::runtime_error(osstr.str());
00869     }
00870 
00871     CHK_ERR( sumIn(numRowIndices, rowIndices, numColIndices, colIndices,
00872                    values, format) );
00873   }
00874 
00875   return(0);
00876 }
00877 
00878 //----------------------------------------------------------------------------
00879 template<typename T>
00880 int fei::Matrix_Impl<T>::multiply(fei::Vector* x,
00881                                  fei::Vector* y)
00882 {
00883   return( fei::MatrixTraits<T>::matvec(matrix_.get(), x, y) );
00884 }
00885 
00886 //----------------------------------------------------------------------------
00887 template<typename T>
00888 int fei::Matrix_Impl<T>::gatherFromOverlap(bool accumulate)
00889 {
00890   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00891     (*output_stream_) << dbgprefix_<<"gatherFromOverlap"<<FEI_ENDL;
00892   }
00893 
00894   CHK_ERR( Matrix_core::gatherFromOverlap(accumulate) );
00895 
00896   return(0);
00897 }
00898 
00899 //----------------------------------------------------------------------------
00900 template<typename T>
00901 int fei::Matrix_Impl<T>::giveToMatrix(int numRows, const int* rows,
00902                                      int numCols, const int* cols,
00903                                      const double* const* values,
00904                                      bool sumInto,
00905                                      int format)
00906 {
00907   if (numRows == 0 || numCols == 0) {
00908     return(0);
00909   }
00910 
00911   if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
00912     return(-1);
00913   }
00914 
00915   const double** myvalues = const_cast<const double**>(values);
00916   if (format != FEI_DENSE_ROW) {
00917     copyTransposeToWorkArrays(numRows, numCols, values,
00918                               work_data1D_, work_data2D_);
00919     myvalues = &work_data2D_[0];
00920   }
00921 
00922   if ((int)work_ints_.size() < numRows) {
00923     work_ints_.resize(numRows);
00924   }
00925 
00926   if (haveBlockMatrix()) {
00927     return( giveToBlockMatrix(numRows, rows, numCols, cols,
00928                               myvalues, sumInto) );
00929   }
00930 
00931   int i; 
00932   int numRemote = 0;
00933   int* workIntPtr = &work_ints_[0];
00934   for(i=0; i<numRows; ++i) {
00935     int row = rows[i];
00936     if (row < firstLocalOffset() || row > lastLocalOffset()) {
00937       ++numRemote;
00938       workIntPtr[i] = 1;
00939     }
00940     else {
00941       workIntPtr[i] = 0;
00942     }
00943   }
00944 
00945   if (numRemote < 1) {
00946     int err = giveToUnderlyingMatrix(numRows, rows, numCols, cols, myvalues,
00947                                     sumInto, FEI_DENSE_ROW);
00948     if (err != 0) {
00949       FEI_OSTRINGSTREAM osstr;
00950       osstr << "fei::Matrix_Impl::giveToMatrix ERROR: err="<<err
00951         << " returned from giveToUnderlyingMatrix.";
00952       throw std::runtime_error(osstr.str());
00953     }
00954     return(0);
00955   }
00956 
00957   for(i=0; i<numRows; ++i) {
00958     int row = rows[i];
00959     const double*const rowvalues = myvalues[i];
00960 
00961     if (workIntPtr[i] > 0) {
00962       int proc = eqnComm_->getOwnerProc(row);
00963       FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
00964 
00965       if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
00966         FEI_OSTREAM& os = *output_stream_;
00967         os << dbgprefix_<<" remote_mat["<<proc<<"]: ";
00968         for(int jj=0; jj<numCols; ++jj) {
00969           os << "("<<row<<","<<cols[jj]<<","<<rowvalues[jj]<<") ";
00970         }
00971         os << FEI_ENDL;
00972       }
00973 
00974       if (sumInto) {
00975         remote_mat->sumInRow(row, cols, rowvalues, numCols);
00976         changedSinceMark_ = true;
00977       }
00978       else {
00979         remote_mat->putRow(row, cols, rowvalues, numCols);
00980         changedSinceMark_ = true;
00981       }
00982 
00983     }
00984     else {
00985       CHK_ERR( giveToUnderlyingMatrix(1, &row, numCols, cols, &rowvalues,
00986                                       sumInto, 0) );
00987     }
00988   }
00989 
00990   return(0);
00991 }
00992 
00993 //----------------------------------------------------------------------------
00994 template<typename T>
00995 int fei::Matrix_Impl<T>::giveToBlockMatrix(int numRows, const int* rows,
00996                                           int numCols, const int* cols,
00997                                           const double* const* values,
00998                                           bool sumInto)
00999 {
01000   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01001     FEI_OSTREAM& os = *output_stream_;
01002     os << "#  giveToBlockMatrix numRows: " << numRows
01003        << ", numCols: " << numCols << FEI_ENDL;
01004   }
01005 
01006   if (numRows == 0 || numCols == 0) {
01007     return(0);
01008   }
01009 
01010   snl_fei::PointBlockMap* pointBlockMap = vecSpace()->getPointBlockMap();
01011 
01012   if (sumInto && numProcs() == 1) {
01013     std::vector<int> temp((numRows+numCols)*2);
01014     int* tempPtr = &temp[0];
01015     int* blkRows = tempPtr;
01016     int* blkRowOffsets = tempPtr+numRows;
01017     int* blkCols = blkRowOffsets+numRows;
01018     int* blkColOffsets = blkCols+numCols;
01019 
01020     CHK_ERR( convertPtToBlk(numRows, rows, numCols, cols,
01021                             blkRows, blkRowOffsets,
01022                             blkCols, blkColOffsets) );
01023 
01024     std::vector<int> blockRows, blockRowSizes;
01025     std::vector<int> blockCols, blockColSizes;
01026     for(int i=0; i<numRows; ++i) fei::sortedListInsert(blkRows[i], blockRows);
01027     for(int i=0; i<numCols; ++i) fei::sortedListInsert(blkCols[i], blockCols);
01028 
01029     int rowSizeTotal = 0, colSizeTotal = 0;
01030 
01031     for(size_t i=0; i<blockRows.size(); ++i) {
01032       int size = pointBlockMap->getBlkEqnSize(blockRows[i]);
01033       blockRowSizes.push_back(size);
01034       rowSizeTotal += size;
01035     }
01036     for(size_t i=0; i<blockCols.size(); ++i) {
01037       int size = pointBlockMap->getBlkEqnSize(blockCols[i]);
01038       blockColSizes.push_back(size);
01039       colSizeTotal += size;
01040     }
01041     std::vector<double> coefs_1d(rowSizeTotal*colSizeTotal, 0.0);
01042     double* coefs1dPtr = &coefs_1d[0];
01043     std::vector<double*> coefs_2d(blockRows.size()*blockCols.size());
01044     double** coefs2dPtr = &coefs_2d[0];
01045 
01046     int blkCounter = 0;
01047     int offset = 0;
01048     for(size_t i=0; i<blockRows.size(); ++i) {
01049       for(size_t j=0; j<blockCols.size(); ++j) {
01050         coefs2dPtr[blkCounter++] = &(coefs1dPtr[offset]);
01051         offset += blockRowSizes[i]*blockColSizes[j];
01052       }
01053     }
01054 
01055     for(int i=0; i<numRows; ++i) {
01056       int rowind = fei::binarySearch(blkRows[i], blockRows);
01057       int rowsize = blockRowSizes[rowind];
01058 
01059       for(int j=0; j<numCols; ++j) {
01060         int colind = fei::binarySearch(blkCols[j], blockCols);
01061         int pos = blkColOffsets[j]*rowsize + blkRowOffsets[i];
01062 
01063         coefs2dPtr[rowind*blockCols.size()+colind][pos] += values[i][j];
01064       }
01065     }
01066 
01067     for(size_t i=0; i<blockRows.size(); ++i) {
01068       CHK_ERR( giveToUnderlyingBlockMatrix(blockRows[i],
01069                                            blockRowSizes[i],
01070                                            blockCols.size(),
01071                                            &blockCols[0],
01072                                            &blockColSizes[0],
01073                                            &blockColSizes[0],
01074                                            &coefs2dPtr[i*blockCols.size()],
01075                                            true) );
01076     }
01077 
01078     return(0);
01079   }
01080 
01081   int maxBlkEqnSize = pointBlockMap->getMaxBlkEqnSize();
01082   int coefBlkLen = maxBlkEqnSize*maxBlkEqnSize*2;
01083 
01084   for(int i=0; i<numRows; ++i) {
01085     int row = rows[i];
01086 
01087     if (row < firstLocalOffset() || row > lastLocalOffset()) {
01088       int proc = eqnComm_->getOwnerProc(row);
01089       FillableMat* remote_mat = getRemotelyOwnedMatrix(proc);
01090       if (sumInto) {
01091         remote_mat->sumInRow(row, cols, values[i], numCols);
01092         changedSinceMark_ = true;
01093       }
01094       else {
01095         remote_mat->putRow(row, cols, values[i], numCols);
01096         changedSinceMark_ = true;
01097       }
01098       continue;
01099     }
01100 
01101     int blockRow = pointBlockMap->eqnToBlkEqn(row);
01102     int blockRowSize = pointBlockMap->getBlkEqnSize(blockRow);
01103     int blockRowOffset = pointBlockMap->getBlkEqnOffset(blockRow, row);
01104 
01105     int blockRowLength = 0;
01106     CHK_ERR( snl_fei::BlockMatrixTraits<T>::getRowLength(matrix_.get(), blockRow,
01107                                                 blockRowLength) );
01108 
01109     std::vector<int> blkCols(blockRowLength);
01110     int* blkCols_ptr = &blkCols[0];
01111     std::vector<int> blkColDims(blockRowLength);
01112     int* blkColDims_ptr = &blkColDims[0];
01113     std::vector<double> coefs_1D(blockRowLength*coefBlkLen);
01114     double* coefs_1D_ptr = &coefs_1D[0];
01115     int coefsLen = coefs_1D.size();
01116     std::vector<double*> coefs_2D(blockRowLength);
01117     double** coefs_2D_ptr = &coefs_2D[0];
01118 
01119     std::vector<int> LDAs(blockRowLength);
01120     std::fill(LDAs.begin(), LDAs.end(), blockRowSize);
01121     std::fill(coefs_1D.begin(), coefs_1D.end(), 0.0);
01122 
01123     int checkRowLen = 0;
01124     CHK_ERR( snl_fei::BlockMatrixTraits<T>::copyOutRow(matrix_.get(),
01125                                               blockRow, blockRowLength,
01126                                               blockRowSize,
01127                                               blkCols_ptr,
01128                                               blkColDims_ptr,
01129                                               coefs_1D_ptr,
01130                                               coefsLen,
01131                                               checkRowLen) );
01132     int coefs_1D_offset = 0;
01133     for(int j=0; j<checkRowLen; ++j) {
01134       coefs_2D_ptr[j] = &(coefs_1D_ptr[coefs_1D_offset]);
01135       coefs_1D_offset += blockRowSize*blkColDims_ptr[j];
01136     }
01137 
01138     for(int j=0; j<numCols; ++j) {
01139       int blockCol, blkOffset;
01140       CHK_ERR( pointBlockMap->getPtEqnInfo(cols[j], blockCol, blkOffset) );
01141 
01142       for(int jj=0; jj<blockRowLength; ++jj) {
01143 
01144         if (blockCol == blkCols_ptr[jj]) {
01145           if (sumInto) {
01146             coefs_2D_ptr[jj][blkOffset*blockRowSize+blockRowOffset] += values[i][j];
01147           }
01148           else {
01149             coefs_2D_ptr[jj][blkOffset*blockRowSize+blockRowOffset] = values[i][j];
01150           }
01151 
01152           break;
01153         }
01154       }
01155     }
01156 
01157     //Now put the block-row back into the matrix
01158     CHK_ERR( giveToUnderlyingBlockMatrix(blockRow, blockRowSize,
01159                                          blockRowLength, blkCols_ptr,
01160                                          &LDAs[0],
01161                                          blkColDims_ptr,
01162                                          coefs_2D_ptr,
01163                                          false) );
01164   }
01165 
01166   return(0);
01167 }
01168 
01169 //----------------------------------------------------------------------------
01170 template<typename T>
01171 int fei::Matrix_Impl<T>::writeToFile(const char* filename,
01172                                     bool matrixMarketFormat)
01173 {
01174   fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
01175   if (mgraph.get() == NULL) {
01176     return(-1);
01177   }
01178 
01179   if (haveFEMatrix()) {
01180     return(0);//temporary no-op...
01181   }
01182 
01183   if (haveBlockMatrix()) {
01184     CHK_ERR( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
01185   }
01186 
01187   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
01188 
01189   fei::SharedPtr<fei::VectorSpace> vspace = mgraph->getRowSpace();
01190 
01191   int globalNNZ = 0;
01192   int globalNumRows = vspace->getGlobalNumIndices();
01193   int localNumRows = vspace->getNumIndices_Owned();
01194 
01195   fei::SharedPtr<fei::VectorSpace> cspace = mgraph->getColSpace();
01196   int globalNumCols = globalNumRows;
01197   if (cspace.get() != NULL) {
01198     globalNumCols = cspace->getGlobalNumIndices();
01199   }
01200 
01201   std::vector<int> indices_owned;
01202   int localNNZ = 0;
01203   CHK_ERR( vspace->getIndices_Owned(indices_owned) );
01204   int* rowsPtr = &indices_owned[0];
01205   for(int i=0; i<localNumRows; ++i) {
01206     int len;
01207     CHK_ERR( getRowLength(rowsPtr[i], len) );
01208     localNNZ += len;
01209   }
01210 
01211   CHK_MPI( GlobalSum(getCommunicator(), localNNZ, globalNNZ) );
01212 
01213   for(int p=0; p<numProcs(); ++p) {
01214     fei::Barrier(getCommunicator());
01215     if (p != localProc()) continue;
01216 
01217     FEI_OFSTREAM* outFile = NULL;
01218     if (p==0) {
01219       outFile = new FEI_OFSTREAM(filename, IOS_OUT);
01220       FEI_OFSTREAM& ofs = *outFile;
01221       if (matrixMarketFormat) {
01222         ofs << mmbanner << FEI_ENDL;
01223         ofs <<globalNumRows << " " <<globalNumCols << " " <<globalNNZ <<FEI_ENDL;
01224       }
01225       else {
01226         ofs <<globalNumRows << " " <<globalNumCols <<FEI_ENDL;
01227       }
01228     }
01229     else outFile = new FEI_OFSTREAM(filename, IOS_APP);
01230 
01231     outFile->setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
01232     outFile->precision(13);
01233     FEI_OFSTREAM& ofs = *outFile;
01234 
01235     int rowLength;
01236 
01237     for(int i=firstLocalOffset(); i<=lastLocalOffset(); ++i) {
01238       CHK_ERR( getRowLength(i, rowLength) );
01239 
01240       work_indices_.resize(rowLength);
01241       work_data1D_.resize(rowLength);
01242 
01243       int* indPtr = &work_indices_[0];
01244       double* coefPtr = &work_data1D_[0];
01245 
01246       CHK_ERR( copyOutRow(i, rowLength, coefPtr, indPtr) );
01247 
01248       fei::insertion_sort_with_companions<double>(rowLength,
01249                                             indPtr, coefPtr);
01250 
01251       for(int j=0; j<rowLength; ++j) {
01252         if (matrixMarketFormat) {
01253           ofs << i+1 <<" "<<indPtr[j]+1<<" "<<coefPtr[j]<<FEI_ENDL;
01254         }
01255         else {
01256           ofs << i <<" "<<indPtr[j]<<" "<<coefPtr[j]<<FEI_ENDL;
01257         }
01258       }
01259     }
01260 
01261     delete outFile;
01262   }
01263 
01264   return(0);
01265 }
01266 
01267 //----------------------------------------------------------------------------
01268 template<typename T>
01269 int fei::Matrix_Impl<T>::writeToStream(FEI_OSTREAM& ostrm,
01270                                       bool matrixMarketFormat)
01271 {
01272   fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
01273   if (mgraph.get() == NULL) {
01274     return(-1);
01275   }
01276 
01277   if (haveFEMatrix()) {
01278     return(0);//temporary no-op...
01279   }
01280 
01281   if (haveBlockMatrix()) {
01282     CHK_ERR( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
01283   }
01284 
01285   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
01286 
01287   fei::SharedPtr<fei::VectorSpace> vspace = mgraph->getRowSpace();
01288 
01289   int globalNNZ = 0;
01290   int localNumRows = vspace->getNumIndices_Owned();
01291   std::vector<int> indices_owned;
01292   int localNNZ = 0;
01293   CHK_ERR( vspace->getIndices_Owned(indices_owned));
01294   int* rowsPtr = &indices_owned[0];
01295   for(int i=0; i<localNumRows; ++i) {
01296     int len;
01297     CHK_ERR( getRowLength(rowsPtr[i], len) );
01298     localNNZ += len;
01299   }
01300 
01301   CHK_MPI( fei::GlobalSum(getCommunicator(), localNNZ, globalNNZ) );
01302 
01303   IOS_FMTFLAGS oldf = ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
01304 
01305   for(int p=0; p<numProcs(); ++p) {
01306     fei::Barrier(getCommunicator());
01307     if (p != localProc()) continue;
01308 
01309     if (p==0) {
01310       int globalSize = globalOffsets()[numProcs()]-1;
01311       if (matrixMarketFormat) {
01312         ostrm << mmbanner << FEI_ENDL;
01313         ostrm << globalSize << " " << globalSize << " " << globalNNZ << FEI_ENDL;
01314       }
01315       else {
01316         ostrm << globalSize << " " << globalSize << FEI_ENDL;
01317       }
01318     }
01319 
01320     int rowLength;
01321 
01322     for(int i=firstLocalOffset(); i<=lastLocalOffset(); ++i) {
01323       CHK_ERR( getRowLength(i, rowLength) );
01324 
01325       work_indices_.resize(rowLength);
01326       work_data1D_.resize(rowLength);
01327 
01328       int* indPtr = &work_indices_[0];
01329       double* coefPtr = &work_data1D_[0];
01330 
01331       CHK_ERR( copyOutRow(i, rowLength, coefPtr, indPtr) );
01332 
01333       for(int j=0; j<rowLength; ++j) {
01334         ostrm << i << " " << indPtr[j] << " " << coefPtr[j] << FEI_ENDL;
01335       }
01336     }
01337   }
01338 
01339   ostrm.setf(oldf, IOS_FLOATFIELD);
01340 
01341   return(0);
01342 }
01343 
01344 #endif // _fei_Matrix_Impl_hpp_
01345 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends