FEI Version of the Day
snl_fei_Broker_LinSysCore.hpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 #ifndef _snl_fei_Broker_LinSysCore_hpp_
00045 #define _snl_fei_Broker_LinSysCore_hpp_
00046 
00047 #include <fei_macros.hpp>
00048 #include <fei_mpi.h>
00049 #include <fei_CommUtils.hpp>
00050 #include <snl_fei_Broker.hpp>
00051 #include <fei_LinearSystemCore.hpp>
00052 #include <fei_VectorSpace.hpp>
00053 #include <fei_Lookup_Impl.hpp>
00054 #include <fei_MatrixGraph.hpp>
00055 #include <fei_SparseRowGraph.hpp>
00056 #include <fei_Vector_Impl.hpp>
00057 #include <fei_Matrix_Impl.hpp>
00058 #include <fei_MatrixReducer.hpp>
00059 #include <fei_VectorReducer.hpp>
00060 #include <fei_Reducer.hpp>
00061 #include <snl_fei_LinearSystem_General.hpp>
00062 
00063 #undef fei_file
00064 #define fei_file "snl_fei_Broker_LinSysCore.hpp"
00065 #include <fei_ErrMacros.hpp>
00066 
00067 namespace snl_fei {
00068 
00072   class Broker_LinSysCore : public snl_fei::Broker {
00073   public:
00075     Broker_LinSysCore(fei::SharedPtr<LinearSystemCore> lsc,
00076                       fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00077                       fei::SharedPtr<fei::Reducer> reducer,
00078                       bool blockMatrix);
00079 
00081     virtual ~Broker_LinSysCore();
00082 
00092     fei::SharedPtr<fei::Vector> createVector(bool isSolutionVector=false)
00093       {
00094         fei::SharedPtr<fei::Vector> vptr;
00095         if (matrixGraph_.get() == NULL) return(vptr);
00096         if (setGlobalOffsets() != 0) return(vptr);
00097 
00098         fei::SharedPtr<fei::Vector> tmpvec;
00099         tmpvec.reset(new fei::Vector_Impl<LinearSystemCore >(matrixGraph_->getRowSpace(),
00100                                                        linsyscore_.get(),
00101                                                        numLocalEqns_,
00102                                                        isSolutionVector));
00103 
00104         fei::SharedPtr<fei::Vector> vec;
00105         if (reducer_.get() != NULL) {
00106           vec.reset(new fei::VectorReducer(reducer_, tmpvec, isSolutionVector));
00107         }
00108         else {
00109           vec = tmpvec;
00110         }
00111 
00112         return(vec);
00113       }
00114 
00117     fei::SharedPtr<fei::Matrix> createMatrix()
00118       {
00119         fei::SharedPtr<fei::Matrix> mptr;
00120         if (matrixGraph_.get() == NULL) return(mptr);
00121 
00122         if (setMatrixStructure() != 0) return(mptr);
00123 
00124         bool zeroSharedRows = true;
00125         if (reducer_.get() != NULL) {
00126           zeroSharedRows = false;
00127         }
00128 
00129         fei::SharedPtr<fei::Matrix> tmpmat;
00130         tmpmat.reset(new fei::Matrix_Impl<LinearSystemCore>(linsyscore_, matrixGraph_,
00131                                                            numLocalEqns_, zeroSharedRows));
00132         fei::SharedPtr<fei::Matrix> matptr;
00133         if (reducer_.get() != NULL) {
00134           matptr.reset(new fei::MatrixReducer(reducer_, tmpmat));
00135         }
00136         else {
00137           matptr = tmpmat;
00138         }
00139 
00140         return(matptr);
00141       }
00142 
00145     fei::SharedPtr<fei::LinearSystem> createLinearSystem()
00146       {
00147         fei::SharedPtr<fei::LinearSystem> lsptr;
00148         if (matrixGraph_.get() == NULL) return(lsptr);
00149 
00150         if (setMatrixStructure() != 0) return(lsptr);
00151 
00152         fei::SharedPtr<fei::LinearSystem>
00153           linSys(new snl_fei::LinearSystem_General(matrixGraph_));
00154         return(linSys);
00155       }
00156 
00158     void setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00159       { matrixGraph_ = matrixGraph; }
00160 
00161   private:
00162     int setGlobalOffsets()
00163       {
00164         //only set the global offsets once.
00165         if (setGlobalOffsets_) return(0);
00166 
00167         if (matrixGraph_.get() == NULL) return(-1);
00168 
00169         MPI_Comm comm = matrixGraph_->getRowSpace()->getCommunicator();
00170         int num_procs = fei::numProcs(comm);
00171         int local_proc = fei::localProc(comm);
00172 
00173         std::vector<int> globalOffsets;
00174         std::vector<int> globalBlkOffsets;
00175 
00176         if (reducer_.get() != NULL) {
00177           int localsize = reducer_->getLocalReducedEqns().size();
00178           numLocalEqns_ = localsize;
00179           std::vector<int> lsizes(num_procs, 0);
00180           std::vector<int> gsizes(num_procs, 0);
00181           lsizes[local_proc] = localsize;
00182           fei::GlobalMax(comm, lsizes, gsizes);
00183           globalOffsets.resize(num_procs+1);
00184           int offset = 0;
00185           for(int p=0; p<num_procs; ++p) {
00186             globalOffsets[p] = offset;
00187             offset += gsizes[p];
00188           }
00189           globalOffsets[num_procs] = offset;
00190           globalBlkOffsets = globalOffsets;
00191         }
00192         else {
00193           fei::SharedPtr<fei::VectorSpace> vecSpace = 
00194             matrixGraph_->getRowSpace();
00195 
00196           vecSpace->getGlobalIndexOffsets(globalOffsets);
00197           vecSpace->getGlobalBlkIndexOffsets(globalBlkOffsets);
00198 
00199           numLocalEqns_ = globalOffsets[local_proc+1]-globalOffsets[local_proc];
00200         }
00201 
00202         CHK_ERR(linsyscore_->setGlobalOffsets(num_procs+1,
00203                                               &globalBlkOffsets[0],
00204                                               &globalOffsets[0],
00205                                               &globalBlkOffsets[0]));
00206 
00207         setGlobalOffsets_ = true;
00208         return(0);
00209       }
00210 
00211     int setMatrixStructure()
00212       {
00213         //only set the matrix matrixGraph once.
00214         if (setMatrixStructure_) return(0);
00215 
00216         if (matrixGraph_.get() == NULL) return(-1);
00217 
00218         lookup_ = new fei::Lookup_Impl(matrixGraph_, 0);
00219 
00220         CHK_ERR( linsyscore_->setLookup(*lookup_) );
00221 
00222         CHK_ERR( setGlobalOffsets() );
00223 
00224         MPI_Comm comm = matrixGraph_->getRowSpace()->getCommunicator();
00225 
00226         fei::SharedPtr<fei::SparseRowGraph> localSRGraph =
00227           matrixGraph_->createGraph(blockMatrix_);
00228 
00229         std::vector<int>& rowNumbers = localSRGraph->rowNumbers;
00230         int numLocalRows = rowNumbers.size();
00231         int* rowOffsets = &(localSRGraph->rowOffsets[0]);
00232         int numLocalNonzeros = localSRGraph->packedColumnIndices.size();
00233         int* nonzeros = &(localSRGraph->packedColumnIndices[0]);
00234 
00235         int numGlobalNonzeros = 0;
00236         fei::GlobalSum(comm, numLocalNonzeros, numGlobalNonzeros);
00237 
00238         std::vector<int*> colPtrs(numLocalRows);
00239         std::vector<int> ptRowsPerBlkRow(numLocalRows, 1);
00240         std::vector<int> rowLengths(numLocalRows);
00241         int* rowLengthsPtr = &rowLengths[0];
00242 
00243         for(int i=0; i<numLocalRows; ++i) {
00244           colPtrs[i] = &(nonzeros[rowOffsets[i]]);
00245           rowLengthsPtr[i] = rowOffsets[i+1]-rowOffsets[i];
00246           if (blockMatrix_ == true) {
00247             ptRowsPerBlkRow[i] = lookup_->getBlkEqnSize(rowNumbers[i]);
00248           }
00249         }
00250 
00251         CHK_ERR( linsyscore_->setMatrixStructure(&colPtrs[0],
00252                                                  rowLengthsPtr,
00253                                                  &colPtrs[0],
00254                                                  rowLengthsPtr,
00255                                                  &ptRowsPerBlkRow[0]));
00256 
00257         setMatrixStructure_ = true;
00258 
00259         return(0);
00260       }
00261 
00262 
00263     fei::SharedPtr<LinearSystemCore> linsyscore_;
00264     fei::SharedPtr<fei::MatrixGraph> matrixGraph_;
00265     fei::SharedPtr<fei::Reducer> reducer_;
00266     Lookup* lookup_;
00267 
00268     bool setGlobalOffsets_;
00269     int numLocalEqns_;
00270     bool setMatrixStructure_;
00271     bool blockMatrix_;
00272   };//class Broker_LinSysCore
00273 }//namespace snl_fei
00274 
00275 #endif // _snl_fei_Broker_LinSysCore_hpp_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends