FEI  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
snl_fei_Broker_LinSysCore.hpp
1 /*--------------------------------------------------------------------*/
2 /* Copyright 2005 Sandia Corporation. */
3 /* Under the terms of Contract DE-AC04-94AL85000, there is a */
4 /* non-exclusive license for use of this work by or on behalf */
5 /* of the U.S. Government. Export of this program may require */
6 /* a license from the United States Government. */
7 /*--------------------------------------------------------------------*/
8 
9 #ifndef _snl_fei_Broker_LinSysCore_hpp_
10 #define _snl_fei_Broker_LinSysCore_hpp_
11 
12 #include <fei_macros.hpp>
13 #include <fei_mpi.h>
14 #include <fei_CommUtils.hpp>
15 #include <snl_fei_Broker.hpp>
16 #include <fei_LinearSystemCore.hpp>
17 #include <fei_VectorSpace.hpp>
18 #include <fei_Lookup_Impl.hpp>
19 #include <fei_MatrixGraph.hpp>
20 #include <fei_SparseRowGraph.hpp>
21 #include <fei_Vector_Impl.hpp>
22 #include <fei_Matrix_Impl.hpp>
23 #include <fei_MatrixReducer.hpp>
24 #include <fei_VectorReducer.hpp>
25 #include <fei_Reducer.hpp>
26 #include <snl_fei_LinearSystem_General.hpp>
27 
28 #undef fei_file
29 #define fei_file "snl_fei_Broker_LinSysCore.hpp"
30 #include <fei_ErrMacros.hpp>
31 
32 namespace snl_fei {
33 
37  class Broker_LinSysCore : public snl_fei::Broker {
38  public:
43  bool blockMatrix);
44 
46  virtual ~Broker_LinSysCore();
47 
57  fei::SharedPtr<fei::Vector> createVector(bool isSolutionVector=false)
58  {
60  if (matrixGraph_.get() == NULL) return(vptr);
61  if (setGlobalOffsets() != 0) return(vptr);
62 
64  tmpvec.reset(new fei::Vector_Impl<LinearSystemCore >(matrixGraph_->getRowSpace(),
65  linsyscore_.get(),
66  numLocalEqns_,
67  isSolutionVector));
68 
70  if (reducer_.get() != NULL) {
71  vec.reset(new fei::VectorReducer(reducer_, tmpvec, isSolutionVector));
72  }
73  else {
74  vec = tmpvec;
75  }
76 
77  return(vec);
78  }
79 
83  {
85  if (matrixGraph_.get() == NULL) return(mptr);
86 
87  if (setMatrixStructure() != 0) return(mptr);
88 
89  bool zeroSharedRows = true;
90  if (reducer_.get() != NULL) {
91  zeroSharedRows = false;
92  }
93 
95  tmpmat.reset(new fei::Matrix_Impl<LinearSystemCore>(linsyscore_, matrixGraph_,
96  numLocalEqns_, zeroSharedRows));
98  if (reducer_.get() != NULL) {
99  matptr.reset(new fei::MatrixReducer(reducer_, tmpmat));
100  }
101  else {
102  matptr = tmpmat;
103  }
104 
105  return(matptr);
106  }
107 
111  {
113  if (matrixGraph_.get() == NULL) return(lsptr);
114 
115  if (setMatrixStructure() != 0) return(lsptr);
116 
118  linSys(new snl_fei::LinearSystem_General(matrixGraph_));
119  return(linSys);
120  }
121 
124  { matrixGraph_ = matrixGraph; }
125 
126  private:
127  int setGlobalOffsets()
128  {
129  //only set the global offsets once.
130  if (setGlobalOffsets_) return(0);
131 
132  if (matrixGraph_.get() == NULL) return(-1);
133 
134  MPI_Comm comm = matrixGraph_->getRowSpace()->getCommunicator();
135  int num_procs = fei::numProcs(comm);
136  int local_proc = fei::localProc(comm);
137 
138  std::vector<int> globalOffsets;
139  std::vector<int> globalBlkOffsets;
140 
141  if (reducer_.get() != NULL) {
142  int localsize = reducer_->getLocalReducedEqns().size();
143  numLocalEqns_ = localsize;
144  std::vector<int> lsizes(num_procs, 0);
145  std::vector<int> gsizes(num_procs, 0);
146  lsizes[local_proc] = localsize;
147  fei::GlobalMax(comm, lsizes, gsizes);
148  globalOffsets.resize(num_procs+1);
149  int offset = 0;
150  for(int p=0; p<num_procs; ++p) {
151  globalOffsets[p] = offset;
152  offset += gsizes[p];
153  }
154  globalOffsets[num_procs] = offset;
155  globalBlkOffsets = globalOffsets;
156  }
157  else {
159  matrixGraph_->getRowSpace();
160 
161  vecSpace->getGlobalIndexOffsets(globalOffsets);
162  vecSpace->getGlobalBlkIndexOffsets(globalBlkOffsets);
163 
164  numLocalEqns_ = globalOffsets[local_proc+1]-globalOffsets[local_proc];
165  }
166 
167  CHK_ERR(linsyscore_->setGlobalOffsets(num_procs+1,
168  &globalBlkOffsets[0],
169  &globalOffsets[0],
170  &globalBlkOffsets[0]));
171 
172  setGlobalOffsets_ = true;
173  return(0);
174  }
175 
176  int setMatrixStructure()
177  {
178  //only set the matrix matrixGraph once.
179  if (setMatrixStructure_) return(0);
180 
181  if (matrixGraph_.get() == NULL) return(-1);
182 
183  lookup_ = new fei::Lookup_Impl(matrixGraph_, 0);
184 
185  CHK_ERR( linsyscore_->setLookup(*lookup_) );
186 
187  CHK_ERR( setGlobalOffsets() );
188 
189  MPI_Comm comm = matrixGraph_->getRowSpace()->getCommunicator();
190 
192  matrixGraph_->createGraph(blockMatrix_);
193 
194  std::vector<int>& rowNumbers = localSRGraph->rowNumbers;
195  int numLocalRows = rowNumbers.size();
196  int* rowOffsets = &(localSRGraph->rowOffsets[0]);
197  int numLocalNonzeros = localSRGraph->packedColumnIndices.size();
198  int* nonzeros = &(localSRGraph->packedColumnIndices[0]);
199 
200  int numGlobalNonzeros = 0;
201  fei::GlobalSum(comm, numLocalNonzeros, numGlobalNonzeros);
202 
203  std::vector<int*> colPtrs(numLocalRows);
204  std::vector<int> ptRowsPerBlkRow(numLocalRows, 1);
205  std::vector<int> rowLengths(numLocalRows);
206  int* rowLengthsPtr = &rowLengths[0];
207 
208  for(int i=0; i<numLocalRows; ++i) {
209  colPtrs[i] = &(nonzeros[rowOffsets[i]]);
210  rowLengthsPtr[i] = rowOffsets[i+1]-rowOffsets[i];
211  if (blockMatrix_ == true) {
212  ptRowsPerBlkRow[i] = lookup_->getBlkEqnSize(rowNumbers[i]);
213  }
214  }
215 
216  CHK_ERR( linsyscore_->setMatrixStructure(&colPtrs[0],
217  rowLengthsPtr,
218  &colPtrs[0],
219  rowLengthsPtr,
220  &ptRowsPerBlkRow[0]));
221 
222  setMatrixStructure_ = true;
223 
224  return(0);
225  }
226 
227 
231  Lookup* lookup_;
232 
233  bool setGlobalOffsets_;
234  int numLocalEqns_;
235  bool setMatrixStructure_;
236  bool blockMatrix_;
237  };//class Broker_LinSysCore
238 }//namespace snl_fei
239 
240 #endif // _snl_fei_Broker_LinSysCore_hpp_
int GlobalSum(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
MPI_Comm getCommunicator() const
fei::SharedPtr< fei::Matrix > createMatrix()
virtual int setLookup(Lookup &lookup)=0
std::vector< int > rowNumbers
virtual int getBlkEqnSize(int blkEqn)=0
virtual fei::SharedPtr< fei::SparseRowGraph > createGraph(bool blockEntryGraph, bool localRowGraph_includeSharedRows=false)=0
void reset(T *p=0)
std::vector< int > packedColumnIndices
std::vector< int > rowOffsets
fei::SharedPtr< fei::LinearSystem > createLinearSystem()
virtual int setGlobalOffsets(int len, int *nodeOffsets, int *eqnOffsets, int *blkEqnOffsets)=0
virtual fei::SharedPtr< fei::VectorSpace > getRowSpace()=0
virtual int setMatrixStructure(int **ptColIndices, int *ptRrowLengths, int **blkColIndices, int *blkRowLengths, int *ptRowsPerBlkRow)=0
T * get() const
void getGlobalBlkIndexOffsets(std::vector< int > &globalBlkOffsets) const
void setMatrixGraph(fei::SharedPtr< fei::MatrixGraph > matrixGraph)
void getGlobalIndexOffsets(std::vector< int > &globalOffsets) const
int GlobalMax(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
int localProc(MPI_Comm comm)
fei::SharedPtr< fei::Vector > createVector(bool isSolutionVector=false)
int numProcs(MPI_Comm comm)
Broker_LinSysCore(fei::SharedPtr< LinearSystemCore > lsc, fei::SharedPtr< fei::MatrixGraph > matrixGraph, fei::SharedPtr< fei::Reducer > reducer, bool blockMatrix)