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         fei::SharedPtr<fei::Matrix> tmpmat;
00125         tmpmat.reset(new fei::Matrix_Impl<LinearSystemCore>(linsyscore_,
00126                                                          matrixGraph_,
00127                                                            numLocalEqns_));
00128         fei::SharedPtr<fei::Matrix> matptr;
00129         if (reducer_.get() != NULL) {
00130           matptr.reset(new fei::MatrixReducer(reducer_, tmpmat));
00131         }
00132         else {
00133           matptr = tmpmat;
00134         }
00135 
00136         return(matptr);
00137       }
00138 
00141     fei::SharedPtr<fei::LinearSystem> createLinearSystem()
00142       {
00143         fei::SharedPtr<fei::LinearSystem> lsptr;
00144         if (matrixGraph_.get() == NULL) return(lsptr);
00145 
00146         if (setMatrixStructure() != 0) return(lsptr);
00147 
00148         fei::SharedPtr<fei::LinearSystem>
00149           linSys(new snl_fei::LinearSystem_General(matrixGraph_));
00150         return(linSys);
00151       }
00152 
00154     void setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00155       { matrixGraph_ = matrixGraph; }
00156 
00157   private:
00158     int setGlobalOffsets()
00159       {
00160         //only set the global offsets once.
00161         if (setGlobalOffsets_) return(0);
00162 
00163         if (matrixGraph_.get() == NULL) return(-1);
00164 
00165         MPI_Comm comm = matrixGraph_->getRowSpace()->getCommunicator();
00166         int num_procs = fei::numProcs(comm);
00167         int local_proc = fei::localProc(comm);
00168 
00169         std::vector<int> globalOffsets;
00170         std::vector<int> globalBlkOffsets;
00171 
00172         if (reducer_.get() != NULL) {
00173           int localsize = reducer_->getLocalReducedEqns().size();
00174           numLocalEqns_ = localsize;
00175           std::vector<int> lsizes(num_procs, 0);
00176           std::vector<int> gsizes(num_procs, 0);
00177           lsizes[local_proc] = localsize;
00178           fei::GlobalMax(comm, lsizes, gsizes);
00179           globalOffsets.resize(num_procs+1);
00180           int offset = 0;
00181           for(int p=0; p<num_procs; ++p) {
00182             globalOffsets[p] = offset;
00183             offset += gsizes[p];
00184           }
00185           globalOffsets[num_procs] = offset;
00186           globalBlkOffsets = globalOffsets;
00187         }
00188         else {
00189           fei::SharedPtr<fei::VectorSpace> vecSpace = 
00190             matrixGraph_->getRowSpace();
00191 
00192           vecSpace->getGlobalIndexOffsets(globalOffsets);
00193           vecSpace->getGlobalBlkIndexOffsets(globalBlkOffsets);
00194 
00195           numLocalEqns_ = globalOffsets[local_proc+1]-globalOffsets[local_proc];
00196         }
00197 
00198         CHK_ERR(linsyscore_->setGlobalOffsets(num_procs+1,
00199                                               &globalBlkOffsets[0],
00200                                               &globalOffsets[0],
00201                                               &globalBlkOffsets[0]));
00202 
00203         setGlobalOffsets_ = true;
00204         return(0);
00205       }
00206 
00207     int setMatrixStructure()
00208       {
00209         //only set the matrix matrixGraph once.
00210         if (setMatrixStructure_) return(0);
00211 
00212         if (matrixGraph_.get() == NULL) return(-1);
00213 
00214         lookup_ = new fei::Lookup_Impl(matrixGraph_, 0);
00215 
00216         CHK_ERR( linsyscore_->setLookup(*lookup_) );
00217 
00218         CHK_ERR( setGlobalOffsets() );
00219 
00220         MPI_Comm comm = matrixGraph_->getRowSpace()->getCommunicator();
00221 
00222         fei::SharedPtr<fei::SparseRowGraph> localSRGraph =
00223           matrixGraph_->createGraph(blockMatrix_);
00224 
00225         std::vector<int>& rowNumbers = localSRGraph->rowNumbers;
00226         int numLocalRows = rowNumbers.size();
00227         int* rowOffsets = &(localSRGraph->rowOffsets[0]);
00228         int numLocalNonzeros = localSRGraph->packedColumnIndices.size();
00229         int* nonzeros = &(localSRGraph->packedColumnIndices[0]);
00230 
00231         int numGlobalNonzeros = 0;
00232         fei::GlobalSum(comm, numLocalNonzeros, numGlobalNonzeros);
00233 
00234         std::vector<int*> colPtrs(numLocalRows);
00235         std::vector<int> ptRowsPerBlkRow(numLocalRows, 1);
00236         std::vector<int> rowLengths(numLocalRows);
00237         int* rowLengthsPtr = &rowLengths[0];
00238 
00239         for(int i=0; i<numLocalRows; ++i) {
00240           colPtrs[i] = &(nonzeros[rowOffsets[i]]);
00241           rowLengthsPtr[i] = rowOffsets[i+1]-rowOffsets[i];
00242           if (blockMatrix_ == true) {
00243             ptRowsPerBlkRow[i] = lookup_->getBlkEqnSize(rowNumbers[i]);
00244           }
00245         }
00246 
00247         CHK_ERR( linsyscore_->setMatrixStructure(&colPtrs[0],
00248                                                  rowLengthsPtr,
00249                                                  &colPtrs[0],
00250                                                  rowLengthsPtr,
00251                                                  &ptRowsPerBlkRow[0]));
00252 
00253         setMatrixStructure_ = true;
00254 
00255         return(0);
00256       }
00257 
00258 
00259     fei::SharedPtr<LinearSystemCore> linsyscore_;
00260     fei::SharedPtr<fei::MatrixGraph> matrixGraph_;
00261     fei::SharedPtr<fei::Reducer> reducer_;
00262     Lookup* lookup_;
00263 
00264     bool setGlobalOffsets_;
00265     int numLocalEqns_;
00266     bool setMatrixStructure_;
00267     bool blockMatrix_;
00268   };//class Broker_LinSysCore
00269 }//namespace snl_fei
00270 
00271 #endif // _snl_fei_Broker_LinSysCore_hpp_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends