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 
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() const
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   fei::SharedPtr<fei::VectorSpace> rspace = mgraph->getRowSpace();
00650   rspace->getGlobalIndices(pattern, rowConn, work_indices2_);
00651 
00652   int numRowIndices = work_indices2_.size();
00653   int* rowIndices = &work_indices2_[0];
00654 
00655   if (haveFEMatrix() || haveBlockMatrix()) {
00656     FieldDofMap<int>& fdofmap = rspace->getFieldDofMap();
00657     const std::map<int,int>& connIDs = cblock->getConnectivityIDs();
00658     std::map<int,int>::const_iterator
00659       iter = connIDs.find(connectivityID);
00660     if (iter == connIDs.end()) ERReturn(-1);
00661     int connOffset = iter->second;
00662     int numIDs = pattern->getNumIDs();
00663     const int* indicesPerID = pattern->getNumIndicesPerID();
00664     const int* fieldsPerID = pattern->getNumFieldsPerID();
00665     const int* fieldIDs = pattern->getFieldIDs();
00666 
00667     int numDofs = 0;
00668     for(int ii=0; ii<numIDs; ++ii) {
00669       numDofs += indicesPerID[ii];
00670     }
00671 
00672     const int* numIndicesPerID = pattern->getNumIndicesPerID();
00673     work_indices_.resize(numIDs+numDofs);
00674     int i, *nodeNumbers = &work_indices_[0];
00675     int* dof_ids = nodeNumbers+numIDs;
00676 
00677     int foffset = 0;
00678     int doffset = 0;
00679     for(i=0; i<numIDs; ++i) {
00680       nodeNumbers[i] = rowConn[i]->getNumber();
00681       for(int ii=0; ii<fieldsPerID[i]; ++ii) {
00682         int fieldSize = rspace->getFieldSize(fieldIDs[foffset]);
00683         int dof_id = fdofmap.get_dof_id(fieldIDs[foffset++], 0);
00684         for(int j=0; j<fieldSize; ++j) {
00685           dof_ids[doffset++] = dof_id + j;
00686         }
00687       }
00688     }
00689 
00690     if (haveFEMatrix()) {
00691       CHK_ERR( snl_fei::FEMatrixTraits<T>::sumInElemMatrix(matrix_.get(), blockID, connOffset,
00692                                                   numIDs, nodeNumbers,
00693                                                   numIndicesPerID, dof_ids, values) );
00694       changedSinceMark_ = true;
00695     }
00696 
00697     if (haveBlockMatrix()) {
00698       if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL &&
00699           format != FEI_BLOCK_DIAGONAL_ROW) {
00700         FEI_CERR << "fei::Matrix_Impl::sumIn ERROR, for block-matrix, format must"
00701                  << " be FEI_DENSE_ROW or FEI_DENSE_COL."<<FEI_ENDL;
00702         ERReturn(-1);
00703       }
00704 
00705       int numPtIndices = pattern->getNumIndices();
00706       int* ptIndices = &work_indices2_[0];
00707 
00708       int numPtColIndices = symmetric ? numPtIndices : colpattern->getNumIndices();
00709 
00710       int len = numPtIndices*numPtColIndices;
00711 
00712       if (format == FEI_BLOCK_DIAGONAL_ROW) {
00713         len = 0;
00714         for(i=0; i<numIDs; ++i) {
00715           len += numIndicesPerID[i]*numIndicesPerID[i];
00716         }
00717       }
00718 
00719       double* ccvalues = new double[len];
00720       //Ouch! Data copy! Remember to optimize this later...
00721       if (format == FEI_BLOCK_DIAGONAL_ROW) {
00722         snl_fei::copy2DBlockDiagToColumnContig(numIDs, numIndicesPerID,
00723                                                values, format, ccvalues);
00724       }
00725       else {
00726         snl_fei::copy2DToColumnContig(numPtIndices, numPtColIndices, values, format,
00727                                       ccvalues);
00728       }
00729 
00730       int ptRowOffset = 0;
00731       for(i=0; i<numIDs; ++i) {
00732 
00733         if (rowConn[i]->getOwnerProc() == localProc()) {
00734 
00735           int numColIDs = numIDs;
00736           int* colNodeNums = nodeNumbers;
00737           const int* colDims = numIndicesPerID;
00738           int LDA = numPtColIndices;
00739           if (format == FEI_BLOCK_DIAGONAL_ROW) {
00740             numColIDs = 1;
00741             colNodeNums = &(nodeNumbers[i]);
00742             colDims = &(numIndicesPerID[i]);
00743             LDA = numIndicesPerID[i];
00744           }
00745 
00746           CHK_ERR( snl_fei::BlockMatrixTraits<T>::sumIn(matrix_.get(),
00747                                                nodeNumbers[i],
00748                                                numIndicesPerID[i],
00749                                                numColIDs, colNodeNums,
00750                                                colDims, LDA,
00751                                                &(ccvalues[ptRowOffset])) );
00752           changedSinceMark_ = true;
00753         }
00754         else {
00755           if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
00756             FEI_OSTREAM& os = *output_stream_;
00757             for(int ii=0; ii<numIndicesPerID[i]; ++ii) {
00758               os << "#   remote pt-row " << ptIndices[ptRowOffset]+ii <<" ";
00759               for(int jj=0; jj<numPtIndices; ++jj) {
00760                 os << "("<<ptIndices[jj]<<","<<values[ptRowOffset+ii][jj]<<") ";
00761               }
00762               os << FEI_ENDL;
00763             }
00764           }
00765 
00766           for(int ii=0; ii<numIndicesPerID[i]; ++ii) {
00767             int p=eqnComm_->getOwnerProc(ptIndices[ptRowOffset]+ii);
00768             remote[p]->sumInRow(ptIndices[ptRowOffset]+ii,
00769                                       ptIndices,
00770                                       values[ptRowOffset+ii],
00771                                       numPtIndices);
00772             changedSinceMark_ = true;
00773           }
00774         }
00775 
00776         ptRowOffset += numIndicesPerID[i];
00777       }
00778 
00779       delete [] ccvalues;
00780     }
00781 
00782     return(0);
00783   }
00784 
00785   int numColIndices = symmetric ? numRowIndices : colpattern->getNumIndices();
00786   int* colIndices = rowIndices;
00787   const fei::Record<int>*const* colConn = NULL;
00788 
00789   if (!symmetric) {
00790     colConn = cblock->getColConnectivity(connectivityID);
00791     mgraph->getColSpace()->getGlobalIndices(colpattern,
00792                                             colConn, work_indices_);
00793     colIndices = &work_indices_[0];
00794   }
00795 
00796   if (symmetric) {
00797     if (format == FEI_DENSE_ROW || format == FEI_DENSE_COL) {
00798       CHK_ERR( sumIn(numRowIndices, rowIndices, numRowIndices, rowIndices,
00799                      values, format) );
00800     }
00801     else if (format == FEI_BLOCK_DIAGONAL_ROW) {
00802       int totalnumfields = pattern->getTotalNumFields();
00803       const int* fieldIDs = pattern->getFieldIDs();
00804       fei::SharedPtr<fei::VectorSpace> rowspace = mgraph->getRowSpace();
00805       fei::VectorSpace* rowspaceptr = rowspace.get();
00806       int ioffset = 0;
00807       for(int i=0; i<totalnumfields; ++i) {
00808         int fieldsize = rowspaceptr->getFieldSize(fieldIDs[i]);
00809         if (ioffset+fieldsize > numRowIndices) {
00810           FEI_OSTRINGSTREAM osstr;
00811           osstr<<"snl_fei::sumIn, format=FEI_BLOCK_DIAGONAL_ROW, block-sizes"
00812                << " not consistent with total num-indices."<<FEI_ENDL;
00813           throw std::runtime_error(osstr.str());
00814         }
00815 
00816         CHK_ERR( sumIn(fieldsize, &(rowIndices[ioffset]),
00817                        fieldsize, &(rowIndices[ioffset]),
00818                        &(values[ioffset]), FEI_DENSE_ROW) );
00819         ioffset += fieldsize;
00820       }
00821     }
00822     else {
00823       FEI_OSTRINGSTREAM osstr;
00824       osstr << "fei::Matrix_Impl::sumIn, format="<<format<<" not supported."
00825             << FEI_ENDL;
00826       throw std::runtime_error(osstr.str());
00827     }
00828   }
00829   else {
00830     if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
00831       FEI_OSTRINGSTREAM osstr;
00832       osstr << "fei::Matrix_Impl::sumIn, format="<<format<<" not valid with"
00833             << " un-symmetric matrix contributions."<<FEI_ENDL;
00834       throw std::runtime_error(osstr.str());
00835     }
00836 
00837     CHK_ERR( sumIn(numRowIndices, rowIndices, numColIndices, colIndices,
00838                    values, format) );
00839   }
00840 
00841   return(0);
00842 }
00843 
00844 //----------------------------------------------------------------------------
00845 template<typename T>
00846 int fei::Matrix_Impl<T>::multiply(fei::Vector* x,
00847                                  fei::Vector* y)
00848 {
00849   return( fei::MatrixTraits<T>::matvec(matrix_.get(), x, y) );
00850 }
00851 
00852 //----------------------------------------------------------------------------
00853 template<typename T>
00854 int fei::Matrix_Impl<T>::gatherFromOverlap(bool accumulate)
00855 {
00856   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00857     (*output_stream_) << dbgprefix_<<"gatherFromOverlap"<<FEI_ENDL;
00858   }
00859 
00860   CHK_ERR( Matrix_core::gatherFromOverlap(accumulate) );
00861 
00862   return(0);
00863 }
00864 
00865 //----------------------------------------------------------------------------
00866 template<typename T>
00867 int fei::Matrix_Impl<T>::giveToMatrix(int numRows, const int* rows,
00868                                      int numCols, const int* cols,
00869                                      const double* const* values,
00870                                      bool sumInto,
00871                                      int format)
00872 {
00873   if (numRows == 0 || numCols == 0) {
00874     return(0);
00875   }
00876 
00877   if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
00878     return(-1);
00879   }
00880 
00881   const double** myvalues = const_cast<const double**>(values);
00882   if (format != FEI_DENSE_ROW) {
00883     copyTransposeToWorkArrays(numRows, numCols, values,
00884                               work_data1D_, work_data2D_);
00885     myvalues = &work_data2D_[0];
00886   }
00887 
00888   if ((int)work_ints_.size() < numRows) {
00889     work_ints_.resize(numRows);
00890   }
00891 
00892   if (haveBlockMatrix()) {
00893     return( giveToBlockMatrix(numRows, rows, numCols, cols,
00894                               myvalues, sumInto) );
00895   }
00896 
00897   int i; 
00898   int numRemote = 0;
00899   int* workIntPtr = &work_ints_[0];
00900   for(i=0; i<numRows; ++i) {
00901     int row = rows[i];
00902     if (row < firstLocalOffset() || row > lastLocalOffset()) {
00903       ++numRemote;
00904       workIntPtr[i] = 1;
00905     }
00906     else {
00907       workIntPtr[i] = 0;
00908     }
00909   }
00910 
00911   if (numRemote < 1) {
00912     int err = giveToUnderlyingMatrix(numRows, rows, numCols, cols, myvalues,
00913                                     sumInto, FEI_DENSE_ROW);
00914     if (err != 0) {
00915       FEI_OSTRINGSTREAM osstr;
00916       osstr << "fei::Matrix_Impl::giveToMatrix ERROR: err="<<err
00917         << " returned from giveToUnderlyingMatrix.";
00918       throw std::runtime_error(osstr.str());
00919     }
00920     return(0);
00921   }
00922 
00923   std::vector<FillableMat*>& remote = getRemotelyOwnedMatrix();
00924 
00925   for(i=0; i<numRows; ++i) {
00926     int row = rows[i];
00927     const double*const rowvalues = myvalues[i];
00928 
00929     if (workIntPtr[i] > 0) {
00930       int proc = eqnComm_->getOwnerProc(row);
00931 
00932       if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
00933         FEI_OSTREAM& os = *output_stream_;
00934         os << dbgprefix_<<" remote[" << proc<<"]: ";
00935         for(int jj=0; jj<numCols; ++jj) {
00936           os << "("<<row<<","<<cols[jj]<<","<<rowvalues[jj]<<") ";
00937         }
00938         os << FEI_ENDL;
00939       }
00940 
00941       if (sumInto) {
00942         remote[proc]->sumInRow(row, cols, rowvalues, numCols);
00943         changedSinceMark_ = true;
00944       }
00945       else {
00946         remote[proc]->putRow(row, cols, rowvalues, numCols);
00947         changedSinceMark_ = true;
00948       }
00949 
00950     }
00951     else {
00952       CHK_ERR( giveToUnderlyingMatrix(1, &row, numCols, cols, &rowvalues,
00953                                       sumInto, 0) );
00954     }
00955   }
00956 
00957   return(0);
00958 }
00959 
00960 //----------------------------------------------------------------------------
00961 template<typename T>
00962 int fei::Matrix_Impl<T>::giveToBlockMatrix(int numRows, const int* rows,
00963                                           int numCols, const int* cols,
00964                                           const double* const* values,
00965                                           bool sumInto)
00966 {
00967   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
00968     FEI_OSTREAM& os = *output_stream_;
00969     os << "#  giveToBlockMatrix numRows: " << numRows
00970        << ", numCols: " << numCols << FEI_ENDL;
00971   }
00972 
00973   if (numRows == 0 || numCols == 0) {
00974     return(0);
00975   }
00976 
00977   snl_fei::PointBlockMap* pointBlockMap = vecSpace()->getPointBlockMap();
00978 
00979   if (sumInto && numProcs() == 1) {
00980     std::vector<int> temp((numRows+numCols)*2);
00981     int* tempPtr = &temp[0];
00982     int* blkRows = tempPtr;
00983     int* blkRowOffsets = tempPtr+numRows;
00984     int* blkCols = blkRowOffsets+numRows;
00985     int* blkColOffsets = blkCols+numCols;
00986 
00987     CHK_ERR( convertPtToBlk(numRows, rows, numCols, cols,
00988                             blkRows, blkRowOffsets,
00989                             blkCols, blkColOffsets) );
00990 
00991     std::vector<int> blockRows, blockRowSizes;
00992     std::vector<int> blockCols, blockColSizes;
00993     for(int i=0; i<numRows; ++i) fei::sortedListInsert(blkRows[i], blockRows);
00994     for(int i=0; i<numCols; ++i) fei::sortedListInsert(blkCols[i], blockCols);
00995 
00996     int rowSizeTotal = 0, colSizeTotal = 0;
00997 
00998     for(size_t i=0; i<blockRows.size(); ++i) {
00999       int size = pointBlockMap->getBlkEqnSize(blockRows[i]);
01000       blockRowSizes.push_back(size);
01001       rowSizeTotal += size;
01002     }
01003     for(size_t i=0; i<blockCols.size(); ++i) {
01004       int size = pointBlockMap->getBlkEqnSize(blockCols[i]);
01005       blockColSizes.push_back(size);
01006       colSizeTotal += size;
01007     }
01008     std::vector<double> coefs_1d(rowSizeTotal*colSizeTotal, 0.0);
01009     double* coefs1dPtr = &coefs_1d[0];
01010     std::vector<double*> coefs_2d(blockRows.size()*blockCols.size());
01011     double** coefs2dPtr = &coefs_2d[0];
01012 
01013     int blkCounter = 0;
01014     int offset = 0;
01015     for(size_t i=0; i<blockRows.size(); ++i) {
01016       for(size_t j=0; j<blockCols.size(); ++j) {
01017         coefs2dPtr[blkCounter++] = &(coefs1dPtr[offset]);
01018         offset += blockRowSizes[i]*blockColSizes[j];
01019       }
01020     }
01021 
01022     for(int i=0; i<numRows; ++i) {
01023       int rowind = fei::binarySearch(blkRows[i], blockRows);
01024       int rowsize = blockRowSizes[rowind];
01025 
01026       for(int j=0; j<numCols; ++j) {
01027         int colind = fei::binarySearch(blkCols[j], blockCols);
01028         int pos = blkColOffsets[j]*rowsize + blkRowOffsets[i];
01029 
01030         coefs2dPtr[rowind*blockCols.size()+colind][pos] += values[i][j];
01031       }
01032     }
01033 
01034     for(size_t i=0; i<blockRows.size(); ++i) {
01035       CHK_ERR( giveToUnderlyingBlockMatrix(blockRows[i],
01036                                            blockRowSizes[i],
01037                                            blockCols.size(),
01038                                            &blockCols[0],
01039                                            &blockColSizes[0],
01040                                            &blockColSizes[0],
01041                                            &coefs2dPtr[i*blockCols.size()],
01042                                            true) );
01043     }
01044 
01045     return(0);
01046   }
01047 
01048   std::vector<FillableMat*>& remote = getRemotelyOwnedMatrix();
01049 
01050   int maxBlkEqnSize = pointBlockMap->getMaxBlkEqnSize();
01051   int coefBlkLen = maxBlkEqnSize*maxBlkEqnSize*2;
01052 
01053   for(int i=0; i<numRows; ++i) {
01054     int row = rows[i];
01055 
01056     if (row < firstLocalOffset() || row > lastLocalOffset()) {
01057       int proc = eqnComm_->getOwnerProc(row);
01058       if (sumInto) {
01059         remote[proc]->sumInRow(row, cols, values[i], numCols);
01060         changedSinceMark_ = true;
01061       }
01062       else {
01063         remote[proc]->putRow(row, cols, values[i], numCols);
01064         changedSinceMark_ = true;
01065       }
01066       continue;
01067     }
01068 
01069     int blockRow = pointBlockMap->eqnToBlkEqn(row);
01070     int blockRowSize = pointBlockMap->getBlkEqnSize(blockRow);
01071     int blockRowOffset = pointBlockMap->getBlkEqnOffset(blockRow, row);
01072 
01073     int blockRowLength = 0;
01074     CHK_ERR( snl_fei::BlockMatrixTraits<T>::getRowLength(matrix_.get(), blockRow,
01075                                                 blockRowLength) );
01076 
01077     std::vector<int> blkCols(blockRowLength);
01078     int* blkCols_ptr = &blkCols[0];
01079     std::vector<int> blkColDims(blockRowLength);
01080     int* blkColDims_ptr = &blkColDims[0];
01081     std::vector<double> coefs_1D(blockRowLength*coefBlkLen);
01082     double* coefs_1D_ptr = &coefs_1D[0];
01083     int coefsLen = coefs_1D.size();
01084     std::vector<double*> coefs_2D(blockRowLength);
01085     double** coefs_2D_ptr = &coefs_2D[0];
01086 
01087     std::vector<int> LDAs(blockRowLength);
01088     std::fill(LDAs.begin(), LDAs.end(), blockRowSize);
01089     std::fill(coefs_1D.begin(), coefs_1D.end(), 0.0);
01090 
01091     int checkRowLen = 0;
01092     CHK_ERR( snl_fei::BlockMatrixTraits<T>::copyOutRow(matrix_.get(),
01093                                               blockRow, blockRowLength,
01094                                               blockRowSize,
01095                                               blkCols_ptr,
01096                                               blkColDims_ptr,
01097                                               coefs_1D_ptr,
01098                                               coefsLen,
01099                                               checkRowLen) );
01100     int coefs_1D_offset = 0;
01101     for(int j=0; j<checkRowLen; ++j) {
01102       coefs_2D_ptr[j] = &(coefs_1D_ptr[coefs_1D_offset]);
01103       coefs_1D_offset += blockRowSize*blkColDims_ptr[j];
01104     }
01105 
01106     for(int j=0; j<numCols; ++j) {
01107       int blockCol = pointBlockMap->eqnToBlkEqn(cols[j]);
01108       int blkOffset= pointBlockMap->getBlkEqnOffset(blockCol, cols[j]);
01109 
01110       for(int jj=0; jj<blockRowLength; ++jj) {
01111 
01112         if (blockCol == blkCols_ptr[jj]) {
01113           if (sumInto) {
01114             coefs_2D_ptr[jj][blkOffset*blockRowSize+blockRowOffset] += values[i][j];
01115           }
01116           else {
01117             coefs_2D_ptr[jj][blkOffset*blockRowSize+blockRowOffset] = values[i][j];
01118           }
01119 
01120           break;
01121         }
01122       }
01123     }
01124 
01125     //Now put the block-row back into the matrix
01126     CHK_ERR( giveToUnderlyingBlockMatrix(blockRow, blockRowSize,
01127                                          blockRowLength, blkCols_ptr,
01128                                          &LDAs[0],
01129                                          blkColDims_ptr,
01130                                          coefs_2D_ptr,
01131                                          false) );
01132   }
01133 
01134   return(0);
01135 }
01136 
01137 //----------------------------------------------------------------------------
01138 template<typename T>
01139 int fei::Matrix_Impl<T>::writeToFile(const char* filename,
01140                                     bool matrixMarketFormat)
01141 {
01142   fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
01143   if (mgraph.get() == NULL) {
01144     return(-1);
01145   }
01146 
01147   if (haveFEMatrix()) {
01148     return(0);//temporary no-op...
01149   }
01150 
01151   if (haveBlockMatrix()) {
01152     CHK_ERR( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
01153   }
01154 
01155   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
01156 
01157   fei::SharedPtr<fei::VectorSpace> vspace = mgraph->getRowSpace();
01158 
01159   int globalNNZ = 0;
01160   int globalNumRows = vspace->getGlobalNumIndices();
01161   int localNumRows = vspace->getNumIndices_Owned();
01162 
01163   fei::SharedPtr<fei::VectorSpace> cspace = mgraph->getColSpace();
01164   int globalNumCols = globalNumRows;
01165   if (cspace.get() != NULL) {
01166     globalNumCols = cspace->getGlobalNumIndices();
01167   }
01168 
01169   std::vector<int> indices_owned;
01170   int localNNZ = 0;
01171   CHK_ERR( vspace->getIndices_Owned(indices_owned) );
01172   int* rowsPtr = &indices_owned[0];
01173   for(int i=0; i<localNumRows; ++i) {
01174     int len;
01175     CHK_ERR( getRowLength(rowsPtr[i], len) );
01176     localNNZ += len;
01177   }
01178 
01179   CHK_MPI( GlobalSum(getCommunicator(), localNNZ, globalNNZ) );
01180 
01181   for(int p=0; p<numProcs(); ++p) {
01182     fei::Barrier(getCommunicator());
01183     if (p != localProc()) continue;
01184 
01185     FEI_OFSTREAM* outFile = NULL;
01186     if (p==0) {
01187       outFile = new FEI_OFSTREAM(filename, IOS_OUT);
01188       FEI_OFSTREAM& ofs = *outFile;
01189       if (matrixMarketFormat) {
01190         ofs << mmbanner << FEI_ENDL;
01191         ofs <<globalNumRows << " " <<globalNumCols << " " <<globalNNZ <<FEI_ENDL;
01192       }
01193       else {
01194         ofs <<globalNumRows << " " <<globalNumCols <<FEI_ENDL;
01195       }
01196     }
01197     else outFile = new FEI_OFSTREAM(filename, IOS_APP);
01198 
01199     outFile->setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
01200     outFile->precision(13);
01201     FEI_OFSTREAM& ofs = *outFile;
01202 
01203     int rowLength;
01204 
01205     for(int i=firstLocalOffset(); i<=lastLocalOffset(); ++i) {
01206       CHK_ERR( getRowLength(i, rowLength) );
01207 
01208       work_indices_.resize(rowLength);
01209       work_data1D_.resize(rowLength);
01210 
01211       int* indPtr = &work_indices_[0];
01212       double* coefPtr = &work_data1D_[0];
01213 
01214       CHK_ERR( copyOutRow(i, rowLength, coefPtr, indPtr) );
01215 
01216       fei::insertion_sort_with_companions<double>(rowLength,
01217                                             indPtr, coefPtr);
01218 
01219       for(int j=0; j<rowLength; ++j) {
01220         if (matrixMarketFormat) {
01221           ofs << i+1 <<" "<<indPtr[j]+1<<" "<<coefPtr[j]<<FEI_ENDL;
01222         }
01223         else {
01224           ofs << i <<" "<<indPtr[j]<<" "<<coefPtr[j]<<FEI_ENDL;
01225         }
01226       }
01227     }
01228 
01229     delete outFile;
01230   }
01231 
01232   return(0);
01233 }
01234 
01235 //----------------------------------------------------------------------------
01236 template<typename T>
01237 int fei::Matrix_Impl<T>::writeToStream(FEI_OSTREAM& ostrm,
01238                                       bool matrixMarketFormat)
01239 {
01240   fei::SharedPtr<fei::MatrixGraph> mgraph = getMatrixGraph();
01241   if (mgraph.get() == NULL) {
01242     return(-1);
01243   }
01244 
01245   if (haveFEMatrix()) {
01246     return(0);//temporary no-op...
01247   }
01248 
01249   if (haveBlockMatrix()) {
01250     CHK_ERR( snl_fei::BlockMatrixTraits<T>::globalAssemble(matrix_.get()) );
01251   }
01252 
01253   static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
01254 
01255   fei::SharedPtr<fei::VectorSpace> vspace = mgraph->getRowSpace();
01256 
01257   int globalNNZ = 0;
01258   int localNumRows = vspace->getNumIndices_Owned();
01259   std::vector<int> indices_owned;
01260   int localNNZ = 0;
01261   CHK_ERR( vspace->getIndices_Owned(indices_owned));
01262   int* rowsPtr = &indices_owned[0];
01263   for(int i=0; i<localNumRows; ++i) {
01264     int len;
01265     CHK_ERR( getRowLength(rowsPtr[i], len) );
01266     localNNZ += len;
01267   }
01268 
01269   CHK_MPI( fei::GlobalSum(getCommunicator(), localNNZ, globalNNZ) );
01270 
01271   IOS_FMTFLAGS oldf = ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
01272 
01273   for(int p=0; p<numProcs(); ++p) {
01274     fei::Barrier(getCommunicator());
01275     if (p != localProc()) continue;
01276 
01277     if (p==0) {
01278       int globalSize = globalOffsets()[numProcs()]-1;
01279       if (matrixMarketFormat) {
01280         ostrm << mmbanner << FEI_ENDL;
01281         ostrm << globalSize << " " << globalSize << " " << globalNNZ << FEI_ENDL;
01282       }
01283       else {
01284         ostrm << globalSize << " " << globalSize << FEI_ENDL;
01285       }
01286     }
01287 
01288     int rowLength;
01289 
01290     for(int i=firstLocalOffset(); i<=lastLocalOffset(); ++i) {
01291       CHK_ERR( getRowLength(i, rowLength) );
01292 
01293       work_indices_.resize(rowLength);
01294       work_data1D_.resize(rowLength);
01295 
01296       int* indPtr = &work_indices_[0];
01297       double* coefPtr = &work_data1D_[0];
01298 
01299       CHK_ERR( copyOutRow(i, rowLength, coefPtr, indPtr) );
01300 
01301       for(int j=0; j<rowLength; ++j) {
01302         ostrm << i << " " << indPtr[j] << " " << coefPtr[j] << FEI_ENDL;
01303       }
01304     }
01305   }
01306 
01307   ostrm.setf(oldf, IOS_FLOATFIELD);
01308 
01309   return(0);
01310 }
01311 
01312 #endif // _fei_Matrix_Impl_hpp_
01313 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends