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

Generated on Tue Jul 13 09:27:45 2010 for FEI by  doxygen 1.4.7