FEI Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fei_Matrix_Local.cpp
Go to the documentation of this file.
1 /*--------------------------------------------------------------------*/
2 /* Copyright 2007 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 #include <fei_ParameterSet.hpp>
10 #include "fei_Matrix_Local.hpp"
11 #include "fei_Matrix_core.hpp"
12 #include "fei_sstream.hpp"
13 #include "fei_fstream.hpp"
14 
15 namespace fei {
16 
19  : matrixGraph_(matrixGraph),
20  sparseRowGraph_(sparseRowGraph),
21  coefs_(sparseRowGraph->packedColumnIndices.size(), 0.0),
22  stateChanged_(false),
23  work_data1D_(),
24  work_data2D_()
25 {
26 }
27 
29 {
30 }
31 
32 //----------------------------------------------------------------------------
35  bool blockEntry)
36 {
38  matrixGraph->createGraph(blockEntry, true);
39  fei::SharedPtr<fei::Matrix> mat(new fei::Matrix_Local(matrixGraph, srg));
40  return(mat);
41 }
42 
43 //----------------------------------------------------------------------------
44 const char*
46 { return( "fei::Matrix_Local" ); }
47 
48 //----------------------------------------------------------------------------
49 int
51 {
52  return(0);
53 }
54 
55 //----------------------------------------------------------------------------
56 int
57 Matrix_Local::parameters(int /*numParams*/, const char* const* /*paramStrings*/)
58 {
59  return(0);
60 }
61 
64 { return( matrixGraph_ ); }
65 
66 void
68 { matrixGraph_ = matrixGraph; }
69 
70 int
72 { return( sparseRowGraph_->rowNumbers.size() ); }
73 
74 int
76 { return( getGlobalNumRows() ); }
77 
78 int
79 Matrix_Local::getRowIndex(int rowNumber) const
80 {
81  int* rows = &(sparseRowGraph_->rowNumbers[0]);
82  int numRows = getLocalNumRows();
83  return( fei::binarySearch(rowNumber, rows, numRows) );
84 }
85 
86 int
87 Matrix_Local::getRowLength(int row, int& length) const
88 {
89  int idx = getRowIndex(row);
90  if (idx < 0) return(idx);
91 
92  length = sparseRowGraph_->rowOffsets[idx+1] -
94  return(0);
95 }
96 
97 int
99 {
100  for(unsigned i=0; i<coefs_.size(); ++i) coefs_[i] = scalar;
101  stateChanged_ = true;
102  return(0);
103 }
104 
105 int
106 Matrix_Local::copyOutRow(int row, int len, double* coefs, int* indices) const
107 {
108  int idx = getRowIndex(row);
109  if (idx < 0) return(idx);
110 
111  int offset = sparseRowGraph_->rowOffsets[idx];
112  int length = sparseRowGraph_->rowOffsets[idx+1]-offset;
113  if (length > len) length = len;
114 
115  for(int i=0; i<length; ++i) {
116  indices[i] = sparseRowGraph_->packedColumnIndices[offset+i];
117  coefs[i] = coefs_[offset+i];
118  }
119 
120  return(0);
121 }
122 
123 int
124 Matrix_Local::giveToMatrix(int numRows, const int* rows,
125  int numCols, const int* cols,
126  const double* const* values,
127  bool sumInto,
128  int format)
129 {
130  if (numRows == 0 || numCols == 0) {
131  return(0);
132  }
133 
134  if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
135  return(-1);
136  }
137 
138  const double** myvalues = const_cast<const double**>(values);
139  if (format != FEI_DENSE_ROW) {
140  fei::Matrix_core::copyTransposeToWorkArrays(numRows, numCols, values,
142  myvalues = &work_data2D_[0];
143  }
144 
145  for(int i=0; i<numRows; ++i) {
146  int idx = getRowIndex(rows[i]);
147  if (idx < 0) {
148  throw std::runtime_error("fei::Matrix_Local::sumIn ERROR, row not found.");
149  }
150 
151  int offset = sparseRowGraph_->rowOffsets[idx];
152  int len = sparseRowGraph_->rowOffsets[idx+1] - offset;
153 
154  int* colInds = &(sparseRowGraph_->packedColumnIndices[offset]);
155  double* coefs = &(coefs_[offset]);
156 
157  for(int j=0; j<numCols; ++j) {
158  int idx2 = fei::binarySearch(cols[j], colInds, len);
159  if (idx2 < 0) {
160  throw std::runtime_error("fei::Matrix_Local::sumIn ERROR, col not found.");
161  }
162 
163  if (sumInto) {
164  coefs[idx2] += myvalues[i][j];
165  }
166  else {
167  coefs[idx2] = myvalues[i][j];
168  }
169  }
170  }
171 
172  stateChanged_ = true;
173  return(0);
174 }
175 
176 int
177 Matrix_Local::sumIn(int numRows, const int* rows,
178  int numCols, const int* cols,
179  const double* const* values,
180  int format)
181 {
182  return( giveToMatrix(numRows, rows, numCols, cols, values,
183  true, format) );
184 }
185 
186 int
187 Matrix_Local::copyIn(int numRows, const int* rows,
188  int numCols, const int* cols,
189  const double* const* values,
190  int format)
191 {
192  return( giveToMatrix(numRows, rows, numCols, cols, values,
193  false, format) );
194 }
195 
196 int
198  int idType,
199  int rowID,
200  int colID,
201  const double* const* data,
202  int format)
203 {
206 
207  int fieldSize = (int)rspace->getFieldSize(fieldID);
208  std::vector<int> indices(2*fieldSize);
209 
210  rspace->getGlobalIndex(idType, rowID, fieldID, indices[0]);
211  for(int i=1; i<fieldSize; ++i) {
212  indices[i] = indices[0]+i;
213  }
214 
215  cspace->getGlobalIndex(idType, colID, fieldID, indices[fieldSize]);
216  for(int i=1; i<fieldSize; ++i) {
217  indices[fieldSize+i] = indices[fieldSize]+i;
218  }
219 
220  return( giveToMatrix(fieldSize, &indices[0], fieldSize, &indices[fieldSize],
221  data, true, format) );
222 }
223 
224 int
226  int idType,
227  int rowID,
228  int colID,
229  const double* data,
230  int format)
231 {
233 
234  int fieldSize = (int)rspace->getFieldSize(fieldID);
235  std::vector<const double*> data2D(fieldSize);
236 
237  int offset = 0;
238  for(int i=0; i<fieldSize; ++i) {
239  data2D[i] = &data[offset];
240  offset += fieldSize;
241  }
242 
243  return( sumInFieldData(fieldID, idType, rowID, colID,
244  &data2D[0], format) );
245 }
246 
247 int
248 Matrix_Local::sumIn(int blockID, int connectivityID,
249  const double* const* values,
250  int format)
251 {
252  int numIndices = matrixGraph_->getConnectivityNumIndices(blockID);
253  std::vector<int> indices(numIndices);
254 
255  matrixGraph_->getConnectivityIndices(blockID, connectivityID,
256  numIndices, &indices[0], numIndices);
257 
258  return( giveToMatrix(numIndices, &indices[0], numIndices, &indices[0],
259  values, true, format) );
260 }
261 
262 int
264 { return(0); }
265 
266 int
268  fei::Vector* y)
269 {
270  FEI_COUT << "fei::Matrix_Local::multiply NOT IMPLEMENTED."<<FEI_ENDL;
271  return(-1);
272 }
273 
274 void
276 {
277 }
278 
279 int
281 {
282  (void)accumulate;
283  return(0);
284 }
285 
286 int
287 Matrix_Local::writeToFile(const char* filename,
288  bool matrixMarketFormat)
289 {
291 
292  MPI_Comm comm = vspace->getCommunicator();
293 
294  int localProc = fei::localProc(comm);
295 
296  FEI_OSTRINGSTREAM osstr;
297  osstr << filename << "." << localProc << ".mtx";
298  std::string fullname = osstr.str();
299 
300  FEI_OFSTREAM ofstr(fullname.c_str(), IOS_OUT);
301 
302  return( writeToStream(ofstr, matrixMarketFormat) );
303 }
304 
305 int
307  bool matrixMarketFormat)
308 {
309  static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
310 
313 
314  int numRows = getLocalNumRows();
315  int numCols = cspace->getEqnNumbers().size();
316  int nnz = coefs_.size();
317 
318  if (matrixMarketFormat) {
319  ostrm << mmbanner << FEI_ENDL;
320  ostrm << numRows << " " << numCols << " " << nnz << FEI_ENDL;
321  }
322  else {
323  ostrm << numRows << " " << numCols << " "<< FEI_ENDL;
324  }
325 
326  std::vector<int>& rowNumbers = sparseRowGraph_->rowNumbers;
327  std::vector<int>& rowOffsets = sparseRowGraph_->rowOffsets;
328  std::vector<int>& colIndices = sparseRowGraph_->packedColumnIndices;
329 
330  ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
331  ostrm.precision(13);
332 
333  int offset = 0;
334  for(unsigned i=0; i<rowNumbers.size(); ++i) {
335  int rowlen = rowOffsets[i+1]-rowOffsets[i];
336 
337  for(int j=0; j<rowlen; ++j) {
338  if (matrixMarketFormat) {
339  ostrm << rowNumbers[i]+1 << " " << colIndices[offset]+1
340  << " " << coefs_[offset] << FEI_ENDL;
341  }
342  else {
343  ostrm << rowNumbers[i] << " " << colIndices[offset]
344  << " " << coefs_[offset] << FEI_ENDL;
345  }
346  ++offset;
347  }
348  }
349 
350  return(0);
351 }
352 
353 bool
355 { return( false ); }
356 
357 void
359 {
360  stateChanged_ = false;
361 }
362 
363 bool
365 { return(stateChanged_); }
366 
367 const std::vector<int>&
369 { return( sparseRowGraph_->rowNumbers ); }
370 
371 const std::vector<int>&
373 { return( sparseRowGraph_->rowOffsets ); }
374 
375 const std::vector<int>&
377 { return( sparseRowGraph_->packedColumnIndices ); }
378 
379 const std::vector<double>&
381 { return( coefs_ ); }
382 
383 }//namespace fei
384 
Matrix_Local(fei::SharedPtr< fei::MatrixGraph > matrixGraph, fei::SharedPtr< fei::SparseRowGraph > sparseRowGraph)
MPI_Comm getCommunicator() const
const std::vector< int > & getRowOffsets() const
int getRowIndex(int rowNumber) const
#define IOS_SCIENTIFIC
int sumInFieldData(int fieldID, int idType, int rowID, int colID, const double *const *data, int format=0)
#define FEI_COUT
int copyIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)
int sumIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)
int getLocalNumRows() const
#define IOS_OUT
std::vector< int > rowNumbers
fei::SharedPtr< fei::SparseRowGraph > sparseRowGraph_
std::vector< double > coefs_
int multiply(fei::Vector *x, fei::Vector *y)
std::vector< const double * > work_data2D_
virtual fei::SharedPtr< fei::SparseRowGraph > createGraph(bool blockEntryGraph, bool localRowGraph_includeSharedRows=false)=0
int giveToMatrix(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, bool sumInto, int format)
#define MPI_Comm
Definition: fei_mpi.h:56
virtual int getConnectivityNumIndices(int blockID) const =0
std::vector< int > packedColumnIndices
int parameters(const fei::ParameterSet &paramset)
std::vector< int > rowOffsets
const std::vector< int > & getRowNumbers() const
const std::vector< int > & getColumnIndices() const
std::vector< double > work_data1D_
int getGlobalIndex(int idType, int ID, int fieldID, int fieldOffset, int whichComponentOfField, int &globalIndex)
int gatherFromOverlap(bool accumulate=true)
virtual fei::SharedPtr< fei::VectorSpace > getRowSpace()=0
int binarySearch(const T &item, const T *list, int len)
const char * typeName()
#define IOS_FLOATFIELD
#define FEI_OSTREAM
Definition: fei_iosfwd.hpp:24
int writeToStream(FEI_OSTREAM &ostrm, bool matrixMarketFormat=true)
#define FEI_OFSTREAM
Definition: fei_fstream.hpp:14
#define FEI_DENSE_COL
Definition: fei_defs.h:80
virtual int getConnectivityIndices(int blockID, int connectivityID, int indicesAllocLen, int *indices, int &numIndices)=0
virtual fei::SharedPtr< fei::VectorSpace > getColSpace()=0
int writeToFile(const char *filename, bool matrixMarketFormat=true)
int getGlobalNumRows() const
int putScalar(double scalar)
const std::vector< double > & getCoefs() const
void setMatrixGraph(fei::SharedPtr< fei::MatrixGraph > matrixGraph)
#define FEI_ENDL
static void copyTransposeToWorkArrays(int numRows, int numCols, const double *const *values, std::vector< double > &work_1D, std::vector< const double * > &work_2D)
fei::SharedPtr< fei::MatrixGraph > matrixGraph_
int localProc(MPI_Comm comm)
unsigned getFieldSize(int fieldID)
int copyOutRow(int row, int len, double *coefs, int *indices) const
static fei::SharedPtr< fei::Matrix > create_Matrix_Local(fei::SharedPtr< fei::MatrixGraph > matrixGraph, bool blockEntry)
#define FEI_DENSE_ROW
Definition: fei_defs.h:77
#define FEI_OSTRINGSTREAM
Definition: fei_sstream.hpp:32
fei::SharedPtr< fei::MatrixGraph > getMatrixGraph() const
std::vector< int > & getEqnNumbers()
int getRowLength(int row, int &length) const