FEI Version of the Day
snl_fei_Broker_LinSysCore.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 _snl_fei_Broker_LinSysCore_hpp_
00010 #define _snl_fei_Broker_LinSysCore_hpp_
00011 
00012 #include <fei_macros.hpp>
00013 #include <fei_mpi.h>
00014 #include <fei_CommUtils.hpp>
00015 #include <snl_fei_Broker.hpp>
00016 #include <fei_LinearSystemCore.hpp>
00017 #include <fei_VectorSpace.hpp>
00018 #include <fei_Lookup_Impl.hpp>
00019 #include <fei_MatrixGraph.hpp>
00020 #include <fei_SparseRowGraph.hpp>
00021 #include <fei_Vector_Impl.hpp>
00022 #include <fei_Matrix_Impl.hpp>
00023 #include <fei_MatrixReducer.hpp>
00024 #include <fei_VectorReducer.hpp>
00025 #include <fei_Reducer.hpp>
00026 #include <snl_fei_LinearSystem_General.hpp>
00027 
00028 #undef fei_file
00029 #define fei_file "snl_fei_Broker_LinSysCore.hpp"
00030 #include <fei_ErrMacros.hpp>
00031 
00032 namespace snl_fei {
00033 
00037   class Broker_LinSysCore : public snl_fei::Broker {
00038   public:
00040     Broker_LinSysCore(fei::SharedPtr<LinearSystemCore> lsc,
00041                       fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00042                       fei::SharedPtr<fei::Reducer> reducer,
00043                       bool blockMatrix);
00044 
00046     virtual ~Broker_LinSysCore();
00047 
00057     fei::SharedPtr<fei::Vector> createVector(bool isSolutionVector=false)
00058       {
00059         fei::SharedPtr<fei::Vector> vptr;
00060         if (matrixGraph_.get() == NULL) return(vptr);
00061         if (setGlobalOffsets() != 0) return(vptr);
00062 
00063         fei::SharedPtr<fei::Vector> tmpvec;
00064         tmpvec.reset(new fei::Vector_Impl<LinearSystemCore >(matrixGraph_->getRowSpace(),
00065                                                        linsyscore_.get(),
00066                                                        numLocalEqns_,
00067                                                        isSolutionVector));
00068 
00069         fei::SharedPtr<fei::Vector> vec;
00070         if (reducer_.get() != NULL) {
00071           vec.reset(new fei::VectorReducer(reducer_, tmpvec, isSolutionVector));
00072         }
00073         else {
00074           vec = tmpvec;
00075         }
00076 
00077         return(vec);
00078       }
00079 
00082     fei::SharedPtr<fei::Matrix> createMatrix()
00083       {
00084         fei::SharedPtr<fei::Matrix> mptr;
00085         if (matrixGraph_.get() == NULL) return(mptr);
00086 
00087         if (setMatrixStructure() != 0) return(mptr);
00088 
00089         bool zeroSharedRows = true;
00090         if (reducer_.get() != NULL) {
00091           zeroSharedRows = false;
00092         }
00093 
00094         fei::SharedPtr<fei::Matrix> tmpmat;
00095         tmpmat.reset(new fei::Matrix_Impl<LinearSystemCore>(linsyscore_, matrixGraph_,
00096                                                            numLocalEqns_, zeroSharedRows));
00097         fei::SharedPtr<fei::Matrix> matptr;
00098         if (reducer_.get() != NULL) {
00099           matptr.reset(new fei::MatrixReducer(reducer_, tmpmat));
00100         }
00101         else {
00102           matptr = tmpmat;
00103         }
00104 
00105         return(matptr);
00106       }
00107 
00110     fei::SharedPtr<fei::LinearSystem> createLinearSystem()
00111       {
00112         fei::SharedPtr<fei::LinearSystem> lsptr;
00113         if (matrixGraph_.get() == NULL) return(lsptr);
00114 
00115         if (setMatrixStructure() != 0) return(lsptr);
00116 
00117         fei::SharedPtr<fei::LinearSystem>
00118           linSys(new snl_fei::LinearSystem_General(matrixGraph_));
00119         return(linSys);
00120       }
00121 
00123     void setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00124       { matrixGraph_ = matrixGraph; }
00125 
00126   private:
00127     int setGlobalOffsets()
00128       {
00129         //only set the global offsets once.
00130         if (setGlobalOffsets_) return(0);
00131 
00132         if (matrixGraph_.get() == NULL) return(-1);
00133 
00134         MPI_Comm comm = matrixGraph_->getRowSpace()->getCommunicator();
00135         int num_procs = fei::numProcs(comm);
00136         int local_proc = fei::localProc(comm);
00137 
00138         std::vector<int> globalOffsets;
00139         std::vector<int> globalBlkOffsets;
00140 
00141         if (reducer_.get() != NULL) {
00142           int localsize = reducer_->getLocalReducedEqns().size();
00143           numLocalEqns_ = localsize;
00144           std::vector<int> lsizes(num_procs, 0);
00145           std::vector<int> gsizes(num_procs, 0);
00146           lsizes[local_proc] = localsize;
00147           fei::GlobalMax(comm, lsizes, gsizes);
00148           globalOffsets.resize(num_procs+1);
00149           int offset = 0;
00150           for(int p=0; p<num_procs; ++p) {
00151             globalOffsets[p] = offset;
00152             offset += gsizes[p];
00153           }
00154           globalOffsets[num_procs] = offset;
00155           globalBlkOffsets = globalOffsets;
00156         }
00157         else {
00158           fei::SharedPtr<fei::VectorSpace> vecSpace = 
00159             matrixGraph_->getRowSpace();
00160 
00161           vecSpace->getGlobalIndexOffsets(globalOffsets);
00162           vecSpace->getGlobalBlkIndexOffsets(globalBlkOffsets);
00163 
00164           numLocalEqns_ = globalOffsets[local_proc+1]-globalOffsets[local_proc];
00165         }
00166 
00167         CHK_ERR(linsyscore_->setGlobalOffsets(num_procs+1,
00168                                               &globalBlkOffsets[0],
00169                                               &globalOffsets[0],
00170                                               &globalBlkOffsets[0]));
00171 
00172         setGlobalOffsets_ = true;
00173         return(0);
00174       }
00175 
00176     int setMatrixStructure()
00177       {
00178         //only set the matrix matrixGraph once.
00179         if (setMatrixStructure_) return(0);
00180 
00181         if (matrixGraph_.get() == NULL) return(-1);
00182 
00183         lookup_ = new fei::Lookup_Impl(matrixGraph_, 0);
00184 
00185         CHK_ERR( linsyscore_->setLookup(*lookup_) );
00186 
00187         CHK_ERR( setGlobalOffsets() );
00188 
00189         MPI_Comm comm = matrixGraph_->getRowSpace()->getCommunicator();
00190 
00191         fei::SharedPtr<fei::SparseRowGraph> localSRGraph =
00192           matrixGraph_->createGraph(blockMatrix_);
00193 
00194         std::vector<int>& rowNumbers = localSRGraph->rowNumbers;
00195         int numLocalRows = rowNumbers.size();
00196         int* rowOffsets = &(localSRGraph->rowOffsets[0]);
00197         int numLocalNonzeros = localSRGraph->packedColumnIndices.size();
00198         int* nonzeros = &(localSRGraph->packedColumnIndices[0]);
00199 
00200         int numGlobalNonzeros = 0;
00201         fei::GlobalSum(comm, numLocalNonzeros, numGlobalNonzeros);
00202 
00203         std::vector<int*> colPtrs(numLocalRows);
00204         std::vector<int> ptRowsPerBlkRow(numLocalRows, 1);
00205         std::vector<int> rowLengths(numLocalRows);
00206         int* rowLengthsPtr = &rowLengths[0];
00207 
00208         for(int i=0; i<numLocalRows; ++i) {
00209           colPtrs[i] = &(nonzeros[rowOffsets[i]]);
00210           rowLengthsPtr[i] = rowOffsets[i+1]-rowOffsets[i];
00211           if (blockMatrix_ == true) {
00212             ptRowsPerBlkRow[i] = lookup_->getBlkEqnSize(rowNumbers[i]);
00213           }
00214         }
00215 
00216         CHK_ERR( linsyscore_->setMatrixStructure(&colPtrs[0],
00217                                                  rowLengthsPtr,
00218                                                  &colPtrs[0],
00219                                                  rowLengthsPtr,
00220                                                  &ptRowsPerBlkRow[0]));
00221 
00222         setMatrixStructure_ = true;
00223 
00224         return(0);
00225       }
00226 
00227 
00228     fei::SharedPtr<LinearSystemCore> linsyscore_;
00229     fei::SharedPtr<fei::MatrixGraph> matrixGraph_;
00230     fei::SharedPtr<fei::Reducer> reducer_;
00231     Lookup* lookup_;
00232 
00233     bool setGlobalOffsets_;
00234     int numLocalEqns_;
00235     bool setMatrixStructure_;
00236     bool blockMatrix_;
00237   };//class Broker_LinSysCore
00238 }//namespace snl_fei
00239 
00240 #endif // _snl_fei_Broker_LinSysCore_hpp_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends