FEI  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
fei_MatrixGraph_Impl2.cpp
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 #include <fei_sstream.hpp>
10 
11 #include <limits>
12 #include <cmath>
13 
14 #include <fei_MatrixGraph_Impl2.hpp>
15 
16 #include <fei_utils.hpp>
17 #include "fei_TemplateUtils.hpp"
18 
19 #include <fei_Pattern.hpp>
20 #include <fei_LogManager.hpp>
21 #include <fei_TemplateUtils.hpp>
22 #include <fei_impl_utils.hpp>
23 #include <snl_fei_Utils.hpp>
24 #include <fei_FieldMask.hpp>
25 #include <fei_Record.hpp>
26 #include <snl_fei_RecordCollection.hpp>
27 #include <fei_VectorSpace.hpp>
28 #include <fei_ParameterSet.hpp>
29 #include <fei_ostream_ops.hpp>
30 #include <fei_Reducer.hpp>
31 #include <fei_GraphReducer.hpp>
32 #include <fei_ConnectivityBlock.hpp>
33 #include <snl_fei_BlkSizeMsgHandler.hpp>
34 #include <fei_Graph_Impl.hpp>
35 #include <snl_fei_Constraint.hpp>
36 
37 #include <fei_EqnBuffer.hpp>
38 #include <fei_EqnCommMgr.hpp>
39 #include <SNL_FEI_Structure.hpp>
40 
41 #undef fei_file
42 #define fei_file "fei_MatrixGraph.cpp"
43 
44 #include <fei_ErrMacros.hpp>
45 
46 namespace snl_fei {
47 static unsigned getFieldSize(int fieldID,
48  fei::VectorSpace* space1,
49  fei::VectorSpace* space2)
50 {
51  unsigned fieldsize = 0;
52  bool foundfield = false;
53  if (space1 != NULL) {
54  try {
55  fieldsize = space1->getFieldSize(fieldID);
56  foundfield = true;
57  }
58  catch (...) {
59  foundfield = false;
60  }
61  }
62 
63  if (!foundfield && space2 != NULL) {
64  try {
65  fieldsize = space2->getFieldSize(fieldID);
66  }
67  catch (std::runtime_error& exc) {
68  std::string msg("snl_fei::getFieldSize: ");
69  msg += exc.what();
70  throw std::runtime_error(msg);
71  }
72  }
73 
74  return(fieldsize);
75 }
76 
77 }//namespace snl_fei
78 
79 //------------------------------------------------------------------------------
80 fei::SharedPtr<fei::MatrixGraph>
81 fei::MatrixGraph_Impl2::Factory::createMatrixGraph(fei::SharedPtr<fei::VectorSpace> rowSpace,
82  fei::SharedPtr<fei::VectorSpace> columnSpace,
83  const char* name)
84 {
85  fei::SharedPtr<fei::MatrixGraph> mgptr(new fei::MatrixGraph_Impl2(rowSpace,
86  columnSpace,
87  name));
88 
89  return(mgptr);
90 }
91 
92 //=====Constructor==============================================================
93 fei::MatrixGraph_Impl2::MatrixGraph_Impl2(fei::SharedPtr<fei::VectorSpace> rowSpace,
94  fei::SharedPtr<fei::VectorSpace> colSpace,
95  const char* name)
96  : comm_(rowSpace->getCommunicator()),
97  rowSpace_(rowSpace),
98  colSpace_(colSpace),
99  haveRowSpace_(false),
100  haveColSpace_(false),
101  symmetric_(false),
102  remotelyOwnedGraphRows_(NULL),
103  simpleProblem_(false),
104  blockEntryGraph_(false),
105  patterns_(),
106  connectivityBlocks_(),
107  arbitraryBlockCounter_(1),
108  sparseBlocks_(),
109  lagrangeConstraints_(),
110  penaltyConstraints_(),
111  slaveConstraints_(),
112  ptEqualBlk_(false),
113  newSlaveData_(false),
114  localNumSlaves_(0),
115  globalNumSlaves_(0),
116  D_(NULL),
117  g_(),
118  g_nonzero_(false),
119  reducer_(),
120  name_(),
121  dbgprefix_("MatGrph: "),
122  tmpIntArray1_(),
123  tmpIntArray2_(),
124  includeAllSlaveConstraints_(false)
125 {
126  localProc_ = fei::localProc(comm_);
127  numProcs_ = fei::numProcs(comm_);
128 
129  if (rowSpace_.get() == NULL) {
130  voidERReturn;
131  }
132  else haveRowSpace_ = true;
133 
134  if (colSpace_.get() != NULL) haveColSpace_ = true;
135  else colSpace_ = rowSpace_;
136 
137  setName(name);
138 }
139 
140 //------------------------------------------------------------------------------
141 fei::MatrixGraph_Impl2::~MatrixGraph_Impl2()
142 {
143  fei::destroyValues(patterns_);
144  patterns_.clear();
145 
146  fei::destroyValues(connectivityBlocks_);
147  connectivityBlocks_.clear();
148 
149  int i, len = sparseBlocks_.size();
150  for(i=0; i<len; ++i) {
151  delete sparseBlocks_[i];
152  }
153 
154  fei::destroyValues(lagrangeConstraints_);
155  lagrangeConstraints_.clear();
156 
157  fei::destroyValues(penaltyConstraints_);
158  penaltyConstraints_.clear();
159 
160  fei::destroyValues(slaveConstraints_);
161  slaveConstraints_.clear();
162 }
163 
164 //------------------------------------------------------------------------------
165 void fei::MatrixGraph_Impl2::setParameters(const fei::ParameterSet& params)
166 {
167  const fei::Param* param = 0;
168  fei::Param::ParamType ptype = fei::Param::BAD_TYPE;
169 
170  param = params.get("FEI_OUTPUT_LEVEL");
171  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
172  if (ptype == fei::Param::STRING) {
173  fei::LogManager& log_manager = fei::LogManager::getLogManager();
174  log_manager.setOutputLevel(param->getStringValue().c_str());
175  setOutputLevel(fei::utils::string_to_output_level(param->getStringValue()));
176  }
177 
178  param = params.get("FEI_LOG_EQN");
179  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
180  if (ptype == fei::Param::INT) {
181  addLogEqn(param->getIntValue());
182  }
183 
184  param = params.get("FEI_LOG_ID");
185  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
186  if (ptype == fei::Param::INT) {
187  addLogID(param->getIntValue());
188  }
189 
190  param = params.get("name");
191  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
192  if (ptype == fei::Param::STRING) {
193  setName(param->getStringValue().c_str());
194  }
195 
196  param = params.get("BLOCK_GRAPH");
197  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
198  if (ptype != fei::Param::BAD_TYPE) {
199  blockEntryGraph_ = true;
200  }
201 
202  param = params.get("BLOCK_MATRIX");
203  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
204  if (ptype != fei::Param::BAD_TYPE) {
205  if (ptype == fei::Param::BOOL) {
206  blockEntryGraph_ = param->getBoolValue();
207  }
208  else {
209  blockEntryGraph_ = true;
210  }
211  }
212 
213  param = params.get("INCLUDE_ALL_SLAVE_CONSTRAINTS");
214  ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
215  if (ptype != fei::Param::BAD_TYPE) {
216  includeAllSlaveConstraints_ = true;
217  }
218 
219 }
220 
221 //----------------------------------------------------------------------------
222 void fei::MatrixGraph_Impl2::setRowSpace(fei::SharedPtr<fei::VectorSpace> rowSpace)
223 {
224  rowSpace_ = rowSpace;
225  haveRowSpace_ = true;
226  if (!haveColSpace_) symmetric_ = true;
227 }
228 
229 //----------------------------------------------------------------------------
230 void fei::MatrixGraph_Impl2::setColumnSpace(fei::SharedPtr<fei::VectorSpace> colSpace)
231 {
232  colSpace_ = colSpace;
233  haveColSpace_ = true;
234  if (colSpace_ == rowSpace_) symmetric_ = true;
235  else symmetric_ = false;
236 }
237 
238 //----------------------------------------------------------------------------
239 int fei::MatrixGraph_Impl2::addPattern(fei::Pattern* pattern)
240 {
241  std::map<int,fei::Pattern*>::iterator
242  p_iter = patterns_.begin(), p_iter_end = patterns_.end();
243 
244  bool matches_existing_pattern = false;
245  int pattern_id = -1;
246  while(p_iter != p_iter_end && !matches_existing_pattern) {
247  if (*pattern == *(p_iter->second)) {
248  matches_existing_pattern = true;
249  pattern_id = p_iter->first;
250  }
251  else ++p_iter;
252  }
253 
254  if (!matches_existing_pattern) {
255  pattern_id = patterns_.size();
256  patterns_.insert(std::make_pair(pattern_id, pattern));
257  }
258  else {
259  delete pattern;
260  }
261 
262  return pattern_id;
263 }
264 
265 //----------------------------------------------------------------------------
266 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
267  int idType)
268 {
269  snl_fei::RecordCollection* rec_coll = NULL;
270  rowSpace_->getRecordCollection(idType, rec_coll);
271 
272  fei::Pattern* pattern = new fei::Pattern(numIDs, idType, rec_coll);
273  return addPattern(pattern);
274 }
275 
276 //----------------------------------------------------------------------------
277 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
278  int idType,
279  int fieldID)
280 {
281  unsigned fieldsize = 0;
282  try {
283  fieldsize = snl_fei::getFieldSize(fieldID, rowSpace_.get(), colSpace_.get());
284  }
285  catch (std::runtime_error& exc) {
286  FEI_OSTRINGSTREAM osstr;
287  osstr << "fei::MatrixGraph_Impl2::definePattern caught error: "<<exc.what();
288  throw std::runtime_error(osstr.str());
289  }
290 
291  snl_fei::RecordCollection* rec_coll = NULL;
292  rowSpace_->getRecordCollection(idType, rec_coll);
293 
294  fei::Pattern* pattern =
295  new fei::Pattern(numIDs, idType, rec_coll, fieldID, fieldsize);
296  return addPattern(pattern);
297 }
298 
299 //----------------------------------------------------------------------------
300 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
301  int idType,
302  const int* numFieldsPerID,
303  const int* fieldIDs)
304 {
305  std::vector<int> fieldSizes;
306  try {
307  int offset = 0;
308  for(int i=0; i<numIDs; ++i) {
309  for(int j=0; j<numFieldsPerID[i]; ++j) {
310  fieldSizes.push_back(snl_fei::getFieldSize(fieldIDs[offset++],
311  rowSpace_.get(), colSpace_.get()));
312  }
313  }
314  }
315  catch (std::runtime_error& exc) {
316  FEI_OSTRINGSTREAM osstr;
317  osstr << "fei::MatrixGraph_Impl2::definePattern caught error: "<<exc.what();
318  throw std::runtime_error(osstr.str());
319  }
320 
321  snl_fei::RecordCollection* rec_coll = NULL;
322  rowSpace_->getRecordCollection(idType, rec_coll);
323 
324  fei::Pattern* pattern = new fei::Pattern(numIDs, idType, rec_coll,
325  numFieldsPerID, fieldIDs, &(fieldSizes[0]));
326  return addPattern(pattern);
327 }
328 
329 //----------------------------------------------------------------------------
330 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
331  const int* idTypes,
332  const int* numFieldsPerID,
333  const int* fieldIDs)
334 {
335  std::vector<int> fieldSizes;
336  try {
337  int offset = 0;
338  for(int i=0; i<numIDs; ++i) {
339  for(int j=0; j<numFieldsPerID[i]; ++j) {
340  fieldSizes.push_back(snl_fei::getFieldSize(fieldIDs[offset++],
341  rowSpace_.get(), colSpace_.get()));
342  }
343  }
344  }
345  catch (std::runtime_error& exc) {
346  FEI_OSTRINGSTREAM osstr;
347  osstr << "fei::MatrixGraph_Impl2::definePattern caught error: "<<exc.what();
348  throw std::runtime_error(osstr.str());
349  }
350 
351  std::vector<snl_fei::RecordCollection*> recordCollections(numIDs);
352  for(int i=0; i<numIDs; ++i) {
353  rowSpace_->getRecordCollection(idTypes[i], recordCollections[i]);
354  }
355 
356  fei::Pattern* pattern = new fei::Pattern(numIDs, idTypes, &recordCollections[0],
357  numFieldsPerID, fieldIDs, &(fieldSizes[0]));
358  return addPattern(pattern);
359 }
360 
361 //------------------------------------------------------------------------------
362 fei::Pattern* fei::MatrixGraph_Impl2::getPattern(int patternID)
363 {
364  std::map<int,fei::Pattern*>::iterator
365  p_iter = patterns_.find(patternID);
366 
367  if (p_iter == patterns_.end()) {
368  return NULL;
369  }
370 
371  fei::Pattern* pattern = (*p_iter).second;
372  return pattern;
373 }
374 
375 //------------------------------------------------------------------------------
376 int fei::MatrixGraph_Impl2::initConnectivityBlock(int blockID,
377  int numConnectivityLists,
378  int patternID,
379  bool diagonal)
380 {
381  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
382  (*output_stream_) <<dbgprefix_<< "initConnectivityBlock symm, blkID="
383  << blockID<<", num="<<numConnectivityLists<<", patternID="
384  <<patternID<<FEI_ENDL;
385  }
386 
387  if (blockID < 0) {
388  fei::console_out() << "fei::MatrixGraph_Impl2::initConnectivityBlock: blockID ("
389  << blockID << ") must be non-negative." << FEI_ENDL;
390  ERReturn(-1);
391  }
392 
393  std::map<int,fei::Pattern*>::iterator
394  p_iter = patterns_.find(patternID);
395 
396  if (p_iter == patterns_.end()) {
397  ERReturn(-1);
398  }
399 
400  fei::Pattern* pattern = (*p_iter).second;
401 
402  if (getConnectivityBlock(blockID) != NULL) ERReturn(-1);
403 
404  fei::ConnectivityBlock* cblock =
405  new fei::ConnectivityBlock(blockID, pattern, numConnectivityLists);
406 
407  connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
408  if (diagonal) cblock->setIsDiagonal(diagonal);
409 
410  return(0);
411 }
412 
413 //------------------------------------------------------------------------------
414 int fei::MatrixGraph_Impl2::initConnectivityBlock(int numConnectivityLists,
415  int patternID,
416  bool diagonal)
417 {
418  int blockID = connectivityBlocks_.size();
419  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
420  (*output_stream_) <<dbgprefix_<< "initConnectivityBlock symm, blkID="
421  << blockID<<", num="<<numConnectivityLists<<", patternID="
422  <<patternID<<FEI_ENDL;
423  }
424 
425  std::map<int,fei::Pattern*>::iterator
426  p_iter = patterns_.find(patternID);
427 
428  if (p_iter == patterns_.end()) {
429  FEI_OSTRINGSTREAM osstr;
430  osstr << "fei::MatrixGraph_Impl2::initConnectivityBlock: ERROR, patternID "
431  << patternID << " not found.";
432  throw std::runtime_error(osstr.str());
433  }
434 
435  fei::Pattern* pattern = (*p_iter).second;
436 
437  fei::ConnectivityBlock* cblock =
438  new fei::ConnectivityBlock(blockID, pattern, numConnectivityLists);
439 
440  connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
441  if (diagonal) cblock->setIsDiagonal(diagonal);
442 
443  return(blockID);
444 }
445 
446 //------------------------------------------------------------------------------
447 int fei::MatrixGraph_Impl2::initConnectivityBlock(int blockID,
448  int numConnectivityLists,
449  int rowPatternID,
450  int colPatternID)
451 {
452  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
453  (*output_stream_) <<dbgprefix_<< "initConnectivityBlock nonsymm, blkID="
454  << blockID<<", num="<<numConnectivityLists<<", row-patternID="
455  <<rowPatternID<<", col-patternID="<<colPatternID<<FEI_ENDL;
456  }
457 
458  if (blockID < 0) {
459  FEI_OSTRINGSTREAM osstr;
460  osstr << "fei::MatrixGraph_Impl2::initConnectivityBlock: blockID ("
461  << blockID << ") must be non-negative." << FEI_ENDL;
462  throw std::runtime_error(osstr.str());
463  }
464 
465  std::map<int,fei::Pattern*>::iterator
466  rp_iter = patterns_.find(rowPatternID);
467  if (rp_iter == patterns_.end()) ERReturn(-1);
468 
469  fei::Pattern* rowpattern = (*rp_iter).second;
470 
471  std::map<int,fei::Pattern*>::iterator
472  cp_iter = patterns_.find(colPatternID);
473  if (cp_iter == patterns_.end()) ERReturn(-1);
474 
475  fei::Pattern* colpattern = (*cp_iter).second;
476 
477  if (getConnectivityBlock(blockID) != NULL) ERReturn(-1);
478 
479  fei::ConnectivityBlock* cblock =
480  new fei::ConnectivityBlock(blockID, rowpattern,
481  colpattern, numConnectivityLists);
482 
483  connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
484 
485  return(0);
486 }
487 
488 //------------------------------------------------------------------------------
489 int fei::MatrixGraph_Impl2::initConnectivity(int blockID,
490  int connectivityID,
491  const int* connectedIdentifiers)
492 {
493  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
494  (*output_stream_) <<dbgprefix_<< "initConnectivity blkID="
495  <<blockID<<", connID="<<connectivityID<<FEI_ENDL;
496  }
497 
498  //for each connected identifier, get a pointer to its record from the
499  //solution-space.
500  //store those pointers in the appropriate connectivity-table, mapped to
501  //the user-provided connectivityID.
502 
503  fei::ConnectivityBlock* connblk = getConnectivityBlock(blockID);
504  if (connblk == NULL) {
505  FEI_OSTRINGSTREAM osstr;
506  osstr << "fei::MatrixGraph_Impl2::initConnectivity ERROR, blockID " << blockID
507  << " doesn't correspond to an existing ConnectivityBlock.";
508  throw std::runtime_error(osstr.str());
509  }
510 
511  fei::Pattern* pattern = connblk->getRowPattern();
512  int numIDs = pattern->getNumIDs();
513 
514  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
515  for(int ii=0; ii<numIDs; ++ii) {
516  if (isLogID(connectedIdentifiers[ii])) {
517  FEI_OSTREAM& os = *output_stream_;
518  os << dbgprefix_<<"initConnectivity blockID " << blockID << ", connID "
519  << connectivityID << ": ";
520  for(int jj=0; jj<numIDs; ++jj) {
521  os << connectedIdentifiers[jj] << " ";
522  }
523  os << FEI_ENDL;
524  }
525  }
526  }
527 
528  if (rowSpace_.get() == NULL) ERReturn(-1);
529 
530  std::map<int,int>& connectivityIDs = connblk->getConnectivityIDs();
531 
532  int idOffset = -1;
533  std::map<int,int>::iterator
534  iter = connectivityIDs.find(connectivityID);
535  if (iter == connectivityIDs.end()) {
536  idOffset = connectivityIDs.size();
537  connectivityIDs.insert(iter, std::make_pair(connectivityID,idOffset));
538  }
539  else {
540  idOffset = iter->second;
541  }
542 
543  idOffset *= pattern->getNumIDs();
544 
545  int* rlist = &(connblk->getRowConnectivities()[idOffset]);
546 
547  CHK_ERR( getConnectivityRecords(pattern, rowSpace_.get(),
548  connectedIdentifiers, rlist) );
549 
550  for(int i=0; i<numIDs; ++i) {
551  if (pattern->getNumFieldsPerID()[i] > 0) {
552  pattern->getRecordCollections()[i]->getRecordWithLocalID(rlist[i])->isInLocalSubdomain_ = true;
553  }
554  }
555 
556  return(0);
557 }
558 
559 //------------------------------------------------------------------------------
560 int fei::MatrixGraph_Impl2::initConnectivity(int idType,
561  int numRows,
562  const int* rowIDs,
563  const int* rowOffsets,
564  const int* packedColumnIDs)
565 {
566  fei::ConnectivityBlock* block = new fei::ConnectivityBlock(numRows,
567  rowIDs, rowOffsets);
568 
569  int max_row_len = 0;
570  for(int i=0; i<numRows; ++i) {
571  int row_len = rowOffsets[i+1]-rowOffsets[i];
572  if (row_len > max_row_len) max_row_len = row_len;
573  }
574 
575  int patternID = definePattern(max_row_len, idType);
576  fei::Pattern* pattern = getPattern(patternID);
577  block->setRowPattern(pattern);
578 
579  sparseBlocks_.push_back(block);
580 
581  int* row_records = &(block->getRowConnectivities()[0]);
582 
583  CHK_ERR( getConnectivityRecords(rowSpace_.get(),
584  idType, numRows, rowIDs, row_records) );
585 
586  fei::VectorSpace* colSpace = rowSpace_.get();
587  if (colSpace_.get() != NULL) {
588  colSpace = colSpace_.get();
589  }
590 
591  int* col_records = &(block->getColConnectivities()[0]);
592 
593  CHK_ERR( getConnectivityRecords(colSpace, idType, rowOffsets[numRows],
594  packedColumnIDs, col_records));
595 
596  return(0);
597 }
598 
599 //------------------------------------------------------------------------------
600 int fei::MatrixGraph_Impl2::initConnectivity(int idType,
601  int numRows,
602  const int* rowIDs,
603  const int* rowLengths,
604  const int*const* columnIDs)
605 {
606  fei::ConnectivityBlock* block = new fei::ConnectivityBlock(numRows,
607  rowIDs, rowLengths, true);
608 
609  int max_row_len = 0;
610  for(int i=0; i<numRows; ++i) {
611  int row_len = rowLengths[i];
612  if (row_len > max_row_len) max_row_len = row_len;
613  }
614 
615  int patternID = definePattern(max_row_len, idType);
616  fei::Pattern* pattern = getPattern(patternID);
617  block->setRowPattern(pattern);
618 
619  sparseBlocks_.push_back(block);
620 
621  int* row_records = &(block->getRowConnectivities()[0]);
622 
623  CHK_ERR( getConnectivityRecords(rowSpace_.get(),
624  idType, numRows, rowIDs, row_records) );
625 
626  fei::VectorSpace* colSpace = rowSpace_.get();
627  if (colSpace_.get() != NULL) {
628  colSpace = colSpace_.get();
629  }
630 
631  int* col_records = &(block->getColConnectivities()[0]);
632 
633  int offset = 0;
634  for(int i=0; i<numRows; ++i) {
635  CHK_ERR( getConnectivityRecords(colSpace, idType, rowLengths[i],
636  columnIDs[i], &(col_records[offset])));
637  offset += rowLengths[i];
638  }
639 
640  return(0);
641 }
642 
643 //------------------------------------------------------------------------------
644 int fei::MatrixGraph_Impl2::getConnectivityRecords(fei::VectorSpace* vecSpace,
645  int idType,
646  int numIDs,
647  const int* IDs,
648  int* records)
649 {
650  snl_fei::RecordCollection* collection = NULL;
651  CHK_ERR( vecSpace->getRecordCollection(idType, collection) );
652 
653  for(int i=0; i<numIDs; ++i) {
654  records[i] = collection->getLocalID(IDs[i]);
655  if (records[i] == -1) {
656  CHK_ERR( vecSpace->addDOFs(idType, 1, &IDs[i], &records[i]) );
657  }
658  }
659 
660  return(0);
661 }
662 
663 //------------------------------------------------------------------------------
664 int fei::MatrixGraph_Impl2::getConnectivityRecords(fei::VectorSpace* vecSpace,
665  int idType,
666  int fieldID,
667  int numIDs,
668  const int* IDs,
669  int* records)
670 {
671  snl_fei::RecordCollection* collection = NULL;
672  CHK_ERR( vecSpace->getRecordCollection(idType, collection) );
673 
674  for(int i=0; i<numIDs; ++i) {
675  records[i] = collection->getLocalID(IDs[i]);
676  if (records[i] == -1) {
677  CHK_ERR( vecSpace->addDOFs(fieldID, idType, 1, &IDs[i], &records[i]));
678  }
679  }
680 
681  return(0);
682 }
683 
684 //------------------------------------------------------------------------------
685 int fei::MatrixGraph_Impl2::getConnectivityRecords(fei::Pattern* pattern,
686  fei::VectorSpace* vecSpace,
687  const int* connectedIdentifiers,
688  int* recordList)
689 {
690  fei::Pattern::PatternType pType = pattern->getPatternType();
691  int i, numIDs = pattern->getNumIDs();
692  const int* numFieldsPerID = pattern->getNumFieldsPerID();
693  const int* fieldIDs = pattern->getFieldIDs();
694 
695  if (pType == fei::Pattern::GENERAL) {
696  const int* idTypes = pattern->getIDTypes();
697 
698  int fieldOffset = 0;
699  for(i=0; i<numIDs; ++i) {
700  int id = connectedIdentifiers[i];
701 
702  for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
703  CHK_ERR( vecSpace->addDOFs(fieldIDs[fieldOffset++],
704  idTypes[i], 1, &id,
705  &recordList[i]));
706  }
707  }
708  }
709  else if (pType == fei::Pattern::SINGLE_IDTYPE ||
710  pType == fei::Pattern::SIMPLE ||
711  pType == fei::Pattern::NO_FIELD) {
712 
713  int idType = pattern->getIDTypes()[0];
714 
715  switch(pType) {
716  case fei::Pattern::SINGLE_IDTYPE:
717  {
718  int fieldOffset = 0;
719  for(i=0; i<numIDs; ++i) {
720  int id = connectedIdentifiers[i];
721  for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
722  CHK_ERR( vecSpace->addDOFs(fieldIDs[fieldOffset++],
723  idType, 1, &id,
724  &(recordList[i])));
725  }
726  }
727  }
728  break;
729  case fei::Pattern::SIMPLE:
730  {
731  CHK_ERR( vecSpace->addDOFs(fieldIDs[0], idType,
732  numIDs, connectedIdentifiers,
733  recordList) );
734  }
735  break;
736  case fei::Pattern::NO_FIELD:
737  CHK_ERR( vecSpace->addDOFs(idType, numIDs,
738  connectedIdentifiers, recordList));
739  break;
740  case fei::Pattern::GENERAL:
741  //Include a stub for this case to avoid compiler warning...
742  std::abort(); break;
743  }
744  }
745  else {
746  ERReturn(-1);
747  }
748 
749  return(0);
750 }
751 
752 //------------------------------------------------------------------------------
753 int fei::MatrixGraph_Impl2::initConnectivity(int blockID,
754  int connectivityID,
755  const int* rowConnectedIdentifiers,
756  const int* colConnectedIdentifiers)
757 {
758  //for each connected identifier, get a pointer to its record from the
759  //solution-space.
760  //store those pointers in the appropriate connectivity-table, mapped to
761  //the user-provided connectivityID.
762 
763  fei::ConnectivityBlock* connblk = getConnectivityBlock(blockID);
764  if (connblk == NULL) {
765  FEI_OSTRINGSTREAM osstr;
766  osstr << "fei::MatrixGraph_Impl2::initConnectivity ERROR, blockID " << blockID
767  << " doesn't correspond to an existing ConnectivityBlock.";
768  throw std::runtime_error(osstr.str());
769  }
770 
771  fei::Pattern* pattern = connblk->getRowPattern();
772  int numIDs = pattern->getNumIDs();
773  fei::Pattern* colPattern = connblk->getColPattern();
774  int numColIDs = colPattern->getNumIDs();
775  if (rowSpace_.get() == NULL) {
776  ERReturn(-1);
777  }
778  if (colSpace_.get() == NULL) {
779  ERReturn(-1);
780  }
781 
782  std::map<int,int>& connectivityIDs = connblk->getConnectivityIDs();
783 
784  int i, idOffset = -1;
785  std::map<int,int>::iterator
786  iter = connectivityIDs.find(connectivityID);
787  if (iter == connectivityIDs.end()) {
788  idOffset = connectivityIDs.size();
789  connectivityIDs.insert(iter, std::make_pair(connectivityID,idOffset));
790  }
791  else {
792  idOffset = iter->second;
793  }
794  int* row_rlist =
795  &(connblk->getRowConnectivities()[idOffset*numIDs]);
796  int* col_rlist =
797  &(connblk->getColConnectivities()[idOffset*numColIDs]);
798 
799  CHK_ERR( getConnectivityRecords(pattern, rowSpace_.get(),
800  rowConnectedIdentifiers, row_rlist) );
801 
802  for(i=0; i<numIDs; ++i)
803  pattern->getRecordCollections()[i]->getRecordWithLocalID(row_rlist[i])->isInLocalSubdomain_ = true;
804  CHK_ERR( getConnectivityRecords(colPattern, colSpace_.get(),
805  colConnectedIdentifiers, col_rlist) );
806 
807  for(i=0; i<numColIDs; ++i)
808  colPattern->getRecordCollections()[i]->getRecordWithLocalID(col_rlist[i])->isInLocalSubdomain_ = true;
809 
810  return(0);
811 }
812 
813 //------------------------------------------------------------------------------
814 int fei::MatrixGraph_Impl2::initConnectivity(int patternID,
815  const int* connectedIdentifiers)
816 {
817  std::map<int,fei::Pattern*>::iterator
818  p_iter = patterns_.find(patternID);
819  if (p_iter == patterns_.end()) ERReturn(-1);
820 
821  fei::Pattern* pattern = p_iter->second;
822 
823  int blockID = -arbitraryBlockCounter_++;
824  fei::ConnectivityBlock* cblock = new fei::ConnectivityBlock(blockID, pattern, 1);
825 
826  connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
827 
828  CHK_ERR( initConnectivity(blockID, 0, connectedIdentifiers) );
829  return(0);
830 }
831 
832 //------------------------------------------------------------------------------
833 int fei::MatrixGraph_Impl2::initConnectivity(int rowPatternID,
834  const int* rowConnectedIdentifiers,
835  int colPatternID,
836  const int* colConnectedIdentifiers)
837 {
838  std::map<int,fei::Pattern*>::iterator
839  rp_iter = patterns_.find(rowPatternID);
840  if (rp_iter == patterns_.end()) ERReturn(-1);
841 
842  fei::Pattern* rowpattern = (*rp_iter).second;
843 
844  std::map<int,fei::Pattern*>::iterator
845  cp_iter = patterns_.find(colPatternID);
846  if (cp_iter == patterns_.end()) ERReturn(-1);
847 
848  fei::Pattern* colpattern = (*cp_iter).second;
849 
850  int blockID = -arbitraryBlockCounter_++;
851  fei::ConnectivityBlock* cblock = new fei::ConnectivityBlock(blockID, rowpattern,
852  colpattern, 1);
853 
854  connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
855 
856  CHK_ERR( initConnectivity(blockID, 0,
857  rowConnectedIdentifiers,
858  colConnectedIdentifiers) );
859  return(0);
860 }
861 
862 //------------------------------------------------------------------------------
863 int fei::MatrixGraph_Impl2::getPatternNumIndices(int patternID,
864  int& numIndices)
865 {
866  std::map<int,fei::Pattern*>::iterator
867  p_iter = patterns_.find(patternID);
868  if (p_iter == patterns_.end()) ERReturn(-1);
869  fei::Pattern* pattern = (*p_iter).second;
870 
871  numIndices = pattern->getNumIndices();
872 
873  return(0);
874 }
875 
876 //------------------------------------------------------------------------------
877 int fei::MatrixGraph_Impl2::getPatternIndices(int patternID,
878  const int* IDs,
879  std::vector<int>& indices)
880 {
881  std::map<int,fei::Pattern*>::iterator
882  p_iter = patterns_.find(patternID);
883  if (p_iter == patterns_.end()) ERReturn(-1);
884 
885  fei::Pattern* pattern = (*p_iter).second;
886 
887  indices.resize(pattern->getNumIndices());
888 
889  int numIDs = pattern->getNumIDs();
890  const int* idTypes = pattern->getIDTypes();
891  const int* numFieldsPerID = pattern->getNumFieldsPerID();
892  const int* fieldIDs = pattern->getFieldIDs();
893 
894  int offset = 0, fieldOffset = 0;
895  for(int i=0; i<numIDs; ++i) {
896  for(int j=0; j<numFieldsPerID[i]; ++j) {
897  CHK_ERR( rowSpace_->getGlobalIndices(1, &(IDs[i]), idTypes[i],
898  fieldIDs[fieldOffset],
899  &(indices[offset])) );
900  offset += snl_fei::getFieldSize(fieldIDs[fieldOffset++],
901  rowSpace_.get(), colSpace_.get());
902  }
903  }
904 
905  return(0);
906 }
907 
908 //------------------------------------------------------------------------------
909 int fei::MatrixGraph_Impl2::initConnectivity(int idType,
910  int fieldID,
911  int numRows,
912  const int* rowIDs,
913  const int* rowOffsets,
914  const int* packedColumnIDs)
915 {
916  fei::ConnectivityBlock* block = new fei::ConnectivityBlock(fieldID, numRows,
917  rowIDs, rowOffsets);
918 
919  sparseBlocks_.push_back(block);
920 
921  int* row_records = &(block->getRowConnectivities()[0]);
922 
923  CHK_ERR( getConnectivityRecords(rowSpace_.get(), idType, fieldID,
924  numRows, rowIDs, row_records) );
925 
926  fei::VectorSpace* colSpace = rowSpace_.get();
927  if (colSpace_.get() != NULL) {
928  colSpace = colSpace_.get();
929  }
930 
931  int* col_records = &(block->getColConnectivities()[0]);
932 
933  CHK_ERR( getConnectivityRecords(colSpace, idType, fieldID, rowOffsets[numRows],
934  packedColumnIDs, col_records));
935 
936  return(0);
937 }
938 
939 //------------------------------------------------------------------------------
940 int fei::MatrixGraph_Impl2::initLagrangeConstraint(int constraintID,
941  int constraintIDType,
942  int numIDs,
943  const int* idTypes,
944  const int* IDs,
945  const int* fieldIDs)
946 {
947  if (rowSpace_.get() == NULL) ERReturn(-1);
948 
949  ConstraintType* constraint = getLagrangeConstraint(constraintID);
950 
951  CHK_ERR( rowSpace_->addDOFs(constraintIDType, 1, &constraintID) );
952 
953  if (haveColSpace_) {
954  if (colSpace_.get() == NULL) {
955  ERReturn(-1);
956  }
957  CHK_ERR( colSpace_->addDOFs(constraintIDType,
958  1, &constraintID) );
959  }
960 
961  ConstraintType* newconstraint =
962  new ConstraintType(constraintID, constraintIDType,
963  false, //isSlave
964  false, //isPenalty
965  numIDs, idTypes, IDs, fieldIDs,
966  0, //slaveOffset
967  0, //offsetIntoSlaveField
968  NULL, //weights
969  0.0, //rhsValue
970  rowSpace_.get());
971 
972  if (constraint != NULL) {
973  if (*constraint != *newconstraint) {
974  delete constraint;
975  }
976  else {
977  delete newconstraint;
978  return(0);
979  }
980  }
981 
982  lagrangeConstraints_[constraintID] = newconstraint;
983 
984  return(0);
985 }
986 
987 //------------------------------------------------------------------------------
988 int fei::MatrixGraph_Impl2::initPenaltyConstraint(int constraintID,
989  int constraintIDType,
990  int numIDs,
991  const int* idTypes,
992  const int* IDs,
993  const int* fieldIDs)
994 {
995  if (rowSpace_.get() == NULL) ERReturn(-1);
996 
997  ConstraintType* constraint = getPenaltyConstraint(constraintID);
998 
999  ConstraintType* newconstraint =
1000  new ConstraintType(constraintID, constraintIDType,
1001  false, //isSlave
1002  true, //isPenalty
1003  numIDs, idTypes, IDs, fieldIDs,
1004  0, //slaveOffset
1005  0, //offsetIntoSlaveField
1006  NULL, //weights
1007  0.0, //rhsValue
1008  rowSpace_.get());
1009 
1010  if (constraint != NULL) {
1011  if (*constraint != *newconstraint) {
1012  delete constraint;
1013  }
1014  else {
1015  delete newconstraint;
1016  return(0);
1017  }
1018  }
1019 
1020  penaltyConstraints_[constraintID] = newconstraint;
1021 
1022  return(0);
1023 }
1024 
1025 //------------------------------------------------------------------------------
1026 bool fei::MatrixGraph_Impl2::hasSlaveDof(int ID, int idType)
1027 {
1028  snl_fei::RecordCollection* collection = NULL;
1029  rowSpace_->getRecordCollection(idType, collection);
1030  if (collection == NULL) {
1031  throw std::runtime_error("fei::MatrixGraph_Impl2::hasSlaveDof: ERROR, unknown idType");
1032  }
1033 
1034  fei::Record<int>* rec = collection->getRecordWithID(ID);
1035 
1036  if (rec == NULL) {
1037  FEI_OSTRINGSTREAM osstr;
1038  osstr << "fei::MatrixGraph_Impl2::hasSlaveDof: ERROR, specified ID ("
1039  << ID << ") not found.";
1040  throw std::runtime_error(osstr.str());
1041  }
1042 
1043  return( rec->hasSlaveDof() );
1044 }
1045 
1046 //------------------------------------------------------------------------------
1047 int fei::MatrixGraph_Impl2::initSlaveConstraint(int numIDs,
1048  const int* idTypes,
1049  const int* IDs,
1050  const int* fieldIDs,
1051  int offsetOfSlave,
1052  int offsetIntoSlaveField,
1053  const double* weights,
1054  double rhsValue)
1055 {
1056  if (rowSpace_.get() == NULL) ERReturn(-1);
1057 
1058  FEI_OSTRINGSTREAM idosstr;
1059  idosstr << IDs[offsetOfSlave] << fieldIDs[offsetOfSlave] << offsetIntoSlaveField;
1060  FEI_ISTRINGSTREAM idisstr(idosstr.str());
1061  int crID = IDs[offsetOfSlave];
1062  idisstr >> crID;
1063 
1064  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
1065  FEI_OSTREAM& dbgos = *output_stream_;
1066 
1067  dbgos << dbgprefix_<<"initSlaveConstraint crID=" << crID << ", slaveID="
1068  << IDs[offsetOfSlave];
1069 
1070  if (output_level_ >= fei::FULL_LOGS) {
1071  dbgos << " { ";
1072  for(int ii=0; ii<numIDs; ++ii) {
1073  dbgos << "("<<IDs[ii] << ","<<fieldIDs[ii]<<") ";
1074  }
1075  dbgos << "}"<<FEI_ENDL;
1076  }
1077  else dbgos << FEI_ENDL;
1078  }
1079 
1080  ConstraintType* constraint = getSlaveConstraint(crID);
1081 
1082  ConstraintType* newconstraint =
1083  new ConstraintType(crID, 0,
1084  true, //isSlave
1085  false, //isPenalty
1086  numIDs, idTypes, IDs, fieldIDs,
1087  offsetOfSlave,
1088  offsetIntoSlaveField,
1089  weights,
1090  rhsValue,
1091  rowSpace_.get());
1092 
1093  if (constraint != NULL) {
1094  if (*constraint != *newconstraint) {
1095  if (!constraint->structurallySame(*newconstraint)) {
1096  FEI_OSTRINGSTREAM osstr;
1097  osstr << "fei::MatrixGraph_Impl2::initSlaveConstraint: slave ID "<<IDs[offsetOfSlave]
1098  << " is already constrained, with different connectivity. Changing the"
1099  << " the structure of an existing constraint is not allowed.";
1100  throw std::runtime_error(osstr.str());
1101  }
1102  newSlaveData_ = true;
1103  delete constraint;
1104  }
1105  else {
1106  delete newconstraint;
1107  return(0);
1108  }
1109  }
1110  else {
1111  newSlaveData_ = true;
1112  }
1113 
1114  slaveConstraints_[crID] = newconstraint;
1115 
1116  return(0);
1117 }
1118 
1119 //------------------------------------------------------------------------------
1120 bool
1121 fei::MatrixGraph_Impl2::newSlaveData()
1122 {
1123  bool globalBool = false;
1124  fei::Allreduce(comm_, newSlaveData_, globalBool);
1125  newSlaveData_ = globalBool;
1126 
1127  return(newSlaveData_);
1128 }
1129 
1130 //------------------------------------------------------------------------------
1131 fei::ConstraintType*
1132 fei::MatrixGraph_Impl2::getLagrangeConstraint(int constraintID)
1133 {
1134  std::map<int,ConstraintType*>::iterator
1135  c_iter = lagrangeConstraints_.find(constraintID);
1136  if (c_iter == lagrangeConstraints_.end()) return(NULL);
1137 
1138  return( (*c_iter).second );
1139 }
1140 
1141 //------------------------------------------------------------------------------
1142 fei::ConstraintType*
1143 fei::MatrixGraph_Impl2::getSlaveConstraint(int constraintID)
1144 {
1145  std::map<int,ConstraintType*>::iterator
1146  c_iter = slaveConstraints_.find(constraintID);
1147  if (c_iter == slaveConstraints_.end()) return(NULL);
1148 
1149  return( (*c_iter).second );
1150 }
1151 
1152 //------------------------------------------------------------------------------
1153 fei::ConstraintType*
1154 fei::MatrixGraph_Impl2::getPenaltyConstraint(int constraintID)
1155 {
1156  std::map<int,ConstraintType*>::iterator
1157  c_iter = penaltyConstraints_.find(constraintID);
1158  if (c_iter == penaltyConstraints_.end()) return(NULL);
1159 
1160  return( (*c_iter).second );
1161 }
1162 
1163 //------------------------------------------------------------------------------
1164 int fei::MatrixGraph_Impl2::initComplete()
1165 {
1166  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
1167  (*output_stream_) <<dbgprefix_<< "initComplete" << FEI_ENDL;
1168  }
1169  if (rowSpace_.get() == NULL) {
1170  ERReturn(-1);
1171  }
1172  else {
1173  CHK_ERR( rowSpace_->initComplete() );
1174  }
1175 
1176  if (haveColSpace_ && colSpace_.get() == NULL) ERReturn(-1);
1177  if (haveColSpace_) {
1178  CHK_ERR( colSpace_->initComplete() );
1179  }
1180 
1181  if (rowSpace_->fieldMasks_.size() == 1 &&
1182  rowSpace_->fieldMasks_[0]->getNumFields() == 1) {
1183  simpleProblem_ = true;
1184  }
1185 
1186  std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
1187  vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
1188 
1189  //If there are slave constraints (on any processor), we need to create
1190  //a slave dependency matrix.
1191  localNumSlaves_ = slaveConstraints_.size();
1192  globalNumSlaves_ = 0;
1193  CHK_ERR( fei::GlobalSum(comm_, localNumSlaves_, globalNumSlaves_) );
1194 
1195  if (globalNumSlaves_ > 0) {
1196  //we need to allocate the slave dependency and reduction matrices too.
1197  CHK_ERR( createSlaveMatrices() );
1198  }
1199 
1200  return(0);
1201 }
1202 
1203 //----------------------------------------------------------------------------
1204 int fei::MatrixGraph_Impl2::createAlgebraicGraph(bool blockEntryGraph,
1205  fei::Graph* graph,
1206  bool gatherFromOverlap)
1207 {
1208  //Now run the connectivity blocks, adding vertices (locations of nonzeros) to
1209  //the graph.
1210 
1211  bool wasAlreadyBlockEntry = blockEntryGraph_;
1212 
1213  blockEntryGraph_ = blockEntryGraph;
1214 
1215  for(unsigned i=0; i<sparseBlocks_.size(); ++i) {
1216  fei::ConnectivityBlock& cblock = *(sparseBlocks_[i]);
1217  CHK_ERR( addBlockToGraph_sparse(graph, &cblock) );
1218  }
1219 
1220  std::map<int,fei::ConnectivityBlock*>::const_iterator
1221  cb_iter = connectivityBlocks_.begin(),
1222  cb_end = connectivityBlocks_.end();
1223 
1224  while(cb_iter != cb_end) {
1225  fei::ConnectivityBlock& cblock = *((*cb_iter).second);
1226 
1227  bool symmetricBlock = cblock.isSymmetric();
1228 
1229  if (symmetricBlock) {
1230  fei::Pattern* pattern = cblock.getRowPattern();
1231  if (pattern == NULL) {
1232  ERReturn(-1);
1233  }
1234 
1235  fei::Pattern::PatternType pType = pattern->getPatternType();
1236  if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
1237  CHK_ERR( addBlockToGraph_multiField_symmetric(graph, &cblock) );
1238  }
1239  else if (pType == fei::Pattern::SIMPLE) {
1240  CHK_ERR( addBlockToGraph_singleField_symmetric(graph, &cblock) );
1241  }
1242  else if (pType == fei::Pattern::NO_FIELD) {
1243  CHK_ERR( addBlockToGraph_noField_symmetric(graph, &cblock) );
1244  }
1245  }
1246  else {
1247  fei::Pattern* pattern = cblock.getRowPattern();
1248  fei::Pattern* colpattern = cblock.getColPattern();
1249  if (pattern == NULL || colpattern == NULL) {
1250  ERReturn(-1);
1251  }
1252 
1253  fei::Pattern::PatternType pType = pattern->getPatternType();
1254  fei::Pattern::PatternType colPType = colpattern->getPatternType();
1255  if (pType == fei::Pattern::SIMPLE && colPType == fei::Pattern::SIMPLE) {
1256  CHK_ERR( addBlockToGraph_singleField_nonsymmetric(graph, &cblock) );
1257  }
1258  else {
1259  CHK_ERR( addBlockToGraph_multiField_nonsymmetric(graph, &cblock) );
1260  }
1261  }
1262  ++cb_iter;
1263  }
1264 
1265  CHK_ERR( addLagrangeConstraintsToGraph(graph) );
1266 
1267  CHK_ERR( addPenaltyConstraintsToGraph(graph) );
1268 
1269  if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
1270  FEI_OSTREAM& os = *output_stream_;
1271  os << dbgprefix_<<"# before graph->gatherFromOverlap()" << FEI_ENDL;
1272  CHK_ERR( graph->writeLocalGraph(os) );
1273  CHK_ERR( graph->writeRemoteGraph(os) );
1274  }
1275 
1276  if (gatherFromOverlap) {
1277  CHK_ERR( graph->gatherFromOverlap() );
1278  }
1279 
1280  if (blockEntryGraph) {
1281  CHK_ERR( exchangeBlkEqnSizes(graph) );
1282  }
1283 
1284  if (output_level_ >= fei::FULL_LOGS && fei::numProcs(comm_)>1 &&
1285  output_stream_ != NULL && gatherFromOverlap) {
1286  FEI_OSTREAM& os = *output_stream_;
1287  os << dbgprefix_<<" after graph->gatherFromOverlap()" << FEI_ENDL;
1288  CHK_ERR( graph->writeLocalGraph(*output_stream_) );
1289  }
1290 
1291  if (!wasAlreadyBlockEntry) {
1292  setIndicesMode(POINT_ENTRY_GRAPH);
1293  }
1294 
1295  return(0);
1296 }
1297 
1298 //----------------------------------------------------------------------------
1299 fei::SharedPtr<fei::SparseRowGraph>
1300 fei::MatrixGraph_Impl2::createGraph(bool blockEntryGraph,
1301  bool localRowGraph_includeSharedRows)
1302 {
1303  fei::SharedPtr<fei::SparseRowGraph> localRows;
1304 
1305  std::vector<int> globalOffsets;
1306 
1307  if (blockEntryGraph) {
1308  if (reducer_.get() != NULL) {
1309  throw std::runtime_error("fei::MatrixGraph_Impl2::createGraph ERROR, can't specify both block-entry assembly and slave-constraint reduction.");
1310  }
1311 
1312  rowSpace_->getGlobalBlkIndexOffsets(globalOffsets);
1313  }
1314  else {
1315  rowSpace_->getGlobalIndexOffsets(globalOffsets);
1316  }
1317 
1318  if ((int)globalOffsets.size() < localProc_+2) return localRows;
1319 
1320  int firstOffset = globalOffsets[localProc_];
1321  int lastOffset = globalOffsets[localProc_+1] - 1;
1322 
1323  if (reducer_.get() != NULL) {
1324  std::vector<int>& reduced_eqns = reducer_->getLocalReducedEqns();
1325  if (!reduced_eqns.empty()) {
1326  firstOffset = reduced_eqns[0];
1327  lastOffset = reduced_eqns[reduced_eqns.size()-1];
1328  }
1329  }
1330 
1331  fei::SharedPtr<fei::Graph> inner_graph(new fei::Graph_Impl(comm_, firstOffset, lastOffset) );
1332  fei::SharedPtr<fei::Graph> graph;
1333 
1334  if (reducer_.get() == NULL) {
1335  graph = inner_graph;
1336  }
1337  else {
1338  graph.reset( new fei::GraphReducer(reducer_, inner_graph) );
1339  }
1340 
1341  bool gatherFromOverlap = !localRowGraph_includeSharedRows;
1342  int err = createAlgebraicGraph(blockEntryGraph, graph.get(),
1343  gatherFromOverlap);
1344  if (err != 0) {
1345  return(localRows);
1346  }
1347 
1348  localRows = fei::createSparseRowGraph(*(inner_graph->getLocalGraph()));
1349 
1350  remotelyOwnedGraphRows_ = fei::createSparseRowGraph(inner_graph->getRemoteGraph());
1351 
1352  remotelyOwnedGraphRows_->blockEntries = blockEntryGraph;
1353 
1354  if (localRowGraph_includeSharedRows &&
1355  remotelyOwnedGraphRows_->rowNumbers.size() > 0) {
1356 
1357  fei::SharedPtr<fei::SparseRowGraph> localPlusShared =
1358  snl_fei::mergeSparseRowGraphs(localRows.get(), remotelyOwnedGraphRows_.get());
1359  localRows = localPlusShared;
1360  }
1361 
1362  return(localRows);
1363 }
1364 
1365 //----------------------------------------------------------------------------
1366 int fei::MatrixGraph_Impl2::exchangeBlkEqnSizes(fei::Graph* graph)
1367 {
1368  //If point-equals-block (meaning block-equations are of size 1) or
1369  //if blockEntryGraph_ is false (meaning we haven't constructed a block-entry-
1370  //graph) then there is nothing to do in this function.
1371  if ( rowSpace_->getPointBlockMap()->ptEqualBlk() ) {
1372  return(0);
1373  }
1374 
1375  snl_fei::BlkSizeMsgHandler blkHandler(rowSpace_.get(), graph, comm_);
1376 
1377  CHK_ERR( blkHandler.do_the_exchange() );
1378 
1379  return(0);
1380 }
1381 
1382 //----------------------------------------------------------------------------
1383 int fei::MatrixGraph_Impl2::createSlaveMatrices()
1384 {
1385  //In this function we extract the dependency matrix of equation numbers
1386  //from the slave-constraints locally on each processor, then gather those
1387  //slave-equations onto each processor (so that each processor holds ALL
1388  //slave-equations).
1389  //Then we remove any couplings that may exist among the slave-equations and
1390  //finally create a FillableMat object to hold the final dependency matrix D_.
1391  //
1392 
1393  if (!newSlaveData()) {
1394  return(0);
1395  }
1396 
1397  std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
1398  vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
1399 
1400  fei::SharedPtr<fei::FillableMat> local_D(new fei::FillableMat);
1401  fei::SharedPtr<fei::CSVec> local_g(new fei::CSVec);
1402 
1403  std::vector<int> masterEqns;
1404  std::vector<double> masterCoefs;
1405 
1406  std::map<int,ConstraintType*>::const_iterator
1407  cr_iter = slaveConstraints_.begin(),
1408  cr_end = slaveConstraints_.end();
1409 
1410  for(; cr_iter != cr_end; ++cr_iter) {
1411  ConstraintType* cr = (*cr_iter).second;
1412 
1413  fei::Record<int>* slaveRecord = cr->getSlave();
1414  int slaveID = slaveRecord->getID();
1415  int slaveIDType = cr->getIDType();
1416  int slaveFieldID = cr->getSlaveFieldID();
1417  int offsetIntoSlaveField = cr->getOffsetIntoSlaveField();
1418  int slaveEqn = -1;
1419  CHK_ERR( rowSpace_->getGlobalIndex(slaveIDType, slaveID,
1420  slaveFieldID, 0, offsetIntoSlaveField,
1421  slaveEqn) );
1422 
1423  std::vector<int>& masterRecords_vec = cr->getMasters();
1424  int* masterRecords = &masterRecords_vec[0];
1425  std::vector<snl_fei::RecordCollection*>& masterRecColls = cr->getMasterRecordCollections();
1426  std::vector<int>& masterIDTypes = cr->getMasterIDTypes();
1427  std::vector<int>& masterFieldIDs = cr->getMasterFieldIDs();
1428  std::vector<double>& masterWeights = cr->getMasterWeights();
1429  double* masterWtPtr = &masterWeights[0];
1430 
1431  masterEqns.resize(masterWeights.size());
1432  masterCoefs.resize(masterWeights.size());
1433 
1434  int* masterEqnsPtr = &(masterEqns[0]);
1435  double* masterCoefsPtr = &(masterCoefs[0]);
1436 
1437  int offset = 0;
1438  for(size_t j=0; j<masterIDTypes.size(); ++j) {
1439  fei::Record<int>* masterRecord = masterRecColls[j]->getRecordWithLocalID(masterRecords[j]);
1440  int* eqnNumbers = vspcEqnPtr_+masterRecord->getOffsetIntoEqnNumbers();
1441  fei::FieldMask* mask = masterRecord->getFieldMask();
1442  int eqnOffset = 0;
1443  if (!simpleProblem_) {
1444  int err = mask->getFieldEqnOffset(masterFieldIDs[j], eqnOffset);
1445  if (err != 0) {
1446  throw std::runtime_error("FEI ERROR, failed to get eqn-offset for constraint master-field.");
1447  }
1448  }
1449 
1450  unsigned fieldSize = rowSpace_->getFieldSize(masterFieldIDs[j]);
1451 
1452  int eqn = eqnNumbers[eqnOffset];
1453  for(unsigned k=0; k<fieldSize; ++k) {
1454  masterEqnsPtr[offset++] = eqn+k;
1455  }
1456  }
1457 
1458  double fei_eps = 1.e-49;
1459 
1460  offset = 0;
1461  for(size_t jj=0; jj<masterEqns.size(); ++jj) {
1462  if ( includeAllSlaveConstraints_ ||
1463  (std::abs(masterWtPtr[jj]) > fei_eps) ) {
1464  masterCoefsPtr[offset] = masterWtPtr[jj];
1465  masterEqnsPtr[offset++] = masterEqnsPtr[jj];
1466  }
1467  }
1468 
1469  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
1470  bool log_eqn = isLogEqn(slaveEqn);
1471  if (!log_eqn) {
1472  for(unsigned ii=0; ii<masterEqns.size(); ++ii) {
1473  if (isLogEqn(masterEqnsPtr[ii])) {
1474  log_eqn = true;
1475  break;
1476  }
1477  }
1478  }
1479 
1480  if (log_eqn) {
1481  FEI_OSTREAM& os = *output_stream_;
1482  os << "createSlaveMatrices: " << slaveEqn << " = ";
1483  for(unsigned ii=0; ii<masterEqns.size(); ++ii) {
1484  if (ii!=0) os << "+ ";
1485  os << masterCoefsPtr[ii]<<"*"<<masterEqnsPtr[ii]<<" ";
1486  }
1487  os << FEI_ENDL;
1488  }
1489  }
1490 
1491  local_D->putRow(slaveEqn, masterEqnsPtr, masterCoefsPtr, offset);
1492 
1493  if ( includeAllSlaveConstraints_ ||
1494  (std::abs(cr->getRHSValue()) > fei_eps) ) {
1495  fei::put_entry(*local_g, slaveEqn, cr->getRHSValue());
1496  }
1497  }
1498 
1499  if (D_.get() == NULL) {
1500  D_.reset(new fei::FillableMat);
1501  }
1502 
1503  if (g_.get() == NULL) {
1504  g_.reset(new fei::CSVec);
1505  }
1506 
1507 #ifndef FEI_SER
1508  if (numProcs_ > 1) {
1509  fei::impl_utils::global_union(comm_, *local_D, *D_);
1510  fei::impl_utils::global_union(comm_, *local_g, *g_);
1511  }
1512  else {
1513  *D_ = *local_D;
1514  *g_ = *local_g;
1515  }
1516 #else
1517  *D_ = *local_D;
1518  *g_ = *local_g;
1519 #endif
1520 
1521  if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
1522  (*output_stream_) << "# D_ (pre-removeCouplings):"<<FEI_ENDL;
1523  (*output_stream_) << *D_;
1524  }
1525 
1526  int levelsOfCoupling = fei::impl_utils::remove_couplings(*D_);
1527  (void)levelsOfCoupling;
1528 
1529  if (reducer_.get() == NULL) {
1530  reducer_.reset(new fei::Reducer(D_, g_, comm_));
1531 
1532  std::vector<int> indices;
1533  rowSpace_->getIndices_Owned(indices);
1534 
1535  reducer_->setLocalUnreducedEqns(indices);
1536  }
1537  else {
1538  reducer_->initialize();
1539  }
1540 
1541  if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
1542  (*output_stream_) << "# D_:"<<FEI_ENDL;
1543  (*output_stream_) << *D_;
1544  }
1545 
1546  double fei_eps = 1.e-49;
1547 
1548  g_nonzero_ = false;
1549  for(size_t j=0; j<g_->size(); ++j) {
1550  double coef = g_->coefs()[j];
1551  if ( includeAllSlaveConstraints_ ||
1552  (std::abs(coef) > fei_eps) ) {
1553  g_nonzero_ = true;
1554  }
1555  }
1556 
1557  newSlaveData_ = false;
1558 
1559  return(0);
1560 }
1561 
1562 //----------------------------------------------------------------------------
1563 fei::SharedPtr<fei::Reducer>
1564 fei::MatrixGraph_Impl2::getReducer()
1565 {
1566  return( reducer_ );
1567 }
1568 
1569 //----------------------------------------------------------------------------
1570 fei::SharedPtr<fei::SparseRowGraph>
1571 fei::MatrixGraph_Impl2::getRemotelyOwnedGraphRows()
1572 {
1573  return( remotelyOwnedGraphRows_ );
1574 }
1575 
1576 //----------------------------------------------------------------------------
1577 void fei::MatrixGraph_Impl2::getConstrainedIndices(std::vector<int>& crindices) const
1578 {
1579  if (constrained_indices_.empty()) {
1580  crindices.clear();
1581  return;
1582  }
1583 
1584  std::set<int>::const_iterator
1585  s_iter = constrained_indices_.begin(),
1586  s_end = constrained_indices_.end();
1587 
1588  crindices.resize(constrained_indices_.size());
1589 
1590  int offset = 0;
1591  for(; s_iter != s_end; ++s_iter) {
1592  crindices[offset++] = *s_iter;
1593  }
1594 }
1595 
1596 //----------------------------------------------------------------------------
1597 int fei::MatrixGraph_Impl2::addLagrangeConstraintsToGraph(fei::Graph* graph)
1598 {
1599  std::vector<int> indices;
1600  std::map<int,ConstraintType*>::const_iterator
1601  cr_iter = lagrangeConstraints_.begin(),
1602  cr_end = lagrangeConstraints_.end();
1603 
1604  constrained_indices_.clear();
1605 
1606  while(cr_iter != cr_end) {
1607  ConstraintType* cr = (*cr_iter).second;
1608  int crID = cr->getConstraintID();
1609 
1610  CHK_ERR( getConstraintConnectivityIndices(cr, indices) );
1611 
1612  int numIndices = indices.size();
1613  int* indicesPtr = &(indices[0]);
1614 
1615  for(int i=0; i<numIndices; ++i) {
1616  constrained_indices_.insert(indicesPtr[i]);
1617  }
1618 
1619  int crEqnRow = -1, tmp;
1620  if (blockEntryGraph_) {
1621  CHK_ERR( rowSpace_->getGlobalBlkIndex(cr->getIDType(),
1622  crID, crEqnRow) );
1623  cr->setBlkEqnNumber(crEqnRow);
1624  CHK_ERR( rowSpace_->getGlobalIndex(cr->getIDType(),
1625  crID, tmp) );
1626  cr->setEqnNumber(tmp);
1627  }
1628  else {
1629  CHK_ERR( rowSpace_->getGlobalIndex(cr->getIDType(),
1630  crID, crEqnRow) );
1631  cr->setEqnNumber(crEqnRow);
1632  CHK_ERR( rowSpace_->getGlobalBlkIndex(cr->getIDType(),
1633  crID, tmp) );
1634  cr->setBlkEqnNumber(tmp);
1635  }
1636 
1637  //now add the row contribution
1638  CHK_ERR( graph->addIndices(crEqnRow, numIndices, &(indices[0])) );
1639 
1640  //Let's add a diagonal entry to the graph for this constraint-equation,
1641  //just in case we need to fiddle with this equation later (e.g. discard the
1642  //constraint equation and replace it with a dirichlet boundary condition...).
1643  CHK_ERR( graph->addIndices(crEqnRow, 1, &crEqnRow) );
1644 
1645  //and finally, add the column contribution (which is simply the transpose
1646  //of the row contribution)
1647  for(int k=0; k<numIndices; ++k) {
1648  CHK_ERR( graph->addIndices(indicesPtr[k], 1, &crEqnRow) );
1649  }
1650 
1651  ++cr_iter;
1652  }
1653 
1654  return(0);
1655 }
1656 
1657 //----------------------------------------------------------------------------
1658 int fei::MatrixGraph_Impl2::
1659 getConstraintConnectivityIndices(ConstraintType* cr,
1660  std::vector<int>& globalIndices)
1661 {
1662  std::vector<int>& fieldSizes = tmpIntArray1_;
1663  std::vector<int>& ones = tmpIntArray2_;
1664 
1665  std::vector<int>& constrainedRecords = cr->getMasters();
1666  std::vector<int>& constrainedFieldIDs = cr->getMasterFieldIDs();
1667  std::vector<snl_fei::RecordCollection*>& recordCollections = cr->getMasterRecordCollections();
1668 
1669  int len = constrainedRecords.size();
1670  fieldSizes.resize(len);
1671 
1672  ones.assign(len, 1);
1673 
1674  int numIndices = 0;
1675  if (blockEntryGraph_) {
1676  numIndices = len;
1677  }
1678  else {
1679  for(int j=0; j<len; ++j) {
1680  unsigned fieldSize = rowSpace_->getFieldSize(constrainedFieldIDs[j]);
1681  fieldSizes[j] = fieldSize;
1682  numIndices += fieldSize;
1683  }
1684  }
1685 
1686  globalIndices.resize(numIndices);
1687 
1688  int checkNum;
1689  CHK_ERR( getConnectivityIndices_multiField(&recordCollections[0],
1690  &constrainedRecords[0],
1691  len, &ones[0],
1692  &constrainedFieldIDs[0],
1693  &fieldSizes[0],
1694  numIndices, &globalIndices[0],
1695  checkNum) );
1696  if (numIndices != checkNum) {
1697  ERReturn(-1);
1698  }
1699 
1700  return(0);
1701 }
1702 
1703 //----------------------------------------------------------------------------
1704 int fei::MatrixGraph_Impl2::addPenaltyConstraintsToGraph(fei::Graph* graph)
1705 {
1706  std::vector<int> indices;
1707  std::map<int,ConstraintType*>::const_iterator
1708  cr_iter = penaltyConstraints_.begin(),
1709  cr_end = penaltyConstraints_.end();
1710 
1711  while(cr_iter != cr_end) {
1712  ConstraintType* cr = (*cr_iter).second;
1713 
1714  CHK_ERR( getConstraintConnectivityIndices(cr, indices) );
1715 
1716  int numIndices = indices.size();
1717 
1718  //now add the symmetric contributions to the graph
1719  CHK_ERR( graph->addSymmetricIndices(numIndices, &(indices[0])) );
1720 
1721  ++cr_iter;
1722  }
1723 
1724  return(0);
1725 }
1726 
1727 //----------------------------------------------------------------------------
1728 int fei::MatrixGraph_Impl2::compareStructure(const fei::MatrixGraph& matrixGraph,
1729  bool& equivalent) const
1730 {
1731  equivalent = false;
1732  int localCode = 1;
1733  int globalCode = 0;
1734 
1735  int myNumBlocks = getNumConnectivityBlocks();
1736  int numBlocks = matrixGraph.getNumConnectivityBlocks();
1737 
1738  if (myNumBlocks != numBlocks) {
1739  CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
1740  equivalent = globalCode > 0 ? false : true;
1741  return(0);
1742  }
1743 
1744  if (numBlocks > 0) {
1745  std::vector<int> myBlockIDs;
1746  std::vector<int> blockIDs;
1747 
1748  CHK_ERR( getConnectivityBlockIDs(myBlockIDs) );
1749  CHK_ERR( matrixGraph.getConnectivityBlockIDs(blockIDs) );
1750 
1751  for(int i=0; i<numBlocks; ++i) {
1752  if (myBlockIDs[i] != blockIDs[i]) {
1753  CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
1754  equivalent = globalCode > 0 ? false : true;
1755  return(0);
1756  }
1757 
1758  const fei::ConnectivityBlock* mycblock = getConnectivityBlock(myBlockIDs[i]);
1759  const fei::ConnectivityBlock* cblock = matrixGraph.getConnectivityBlock(blockIDs[i]);
1760 
1761  int myNumLists = mycblock->getConnectivityIDs().size();
1762  int numLists = cblock->getConnectivityIDs().size();
1763 
1764  if (myNumLists != numLists ||
1765  mycblock->isSymmetric() != cblock->isSymmetric()) {
1766  CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
1767  equivalent = globalCode > 0 ? false : true;
1768  return(0);
1769  }
1770 
1771  int myNumIDsPerList = mycblock->getRowPattern()->getNumIDs();
1772  int numIDsPerList = cblock->getRowPattern()->getNumIDs();
1773 
1774  if (myNumIDsPerList != numIDsPerList) {
1775  CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
1776  equivalent = globalCode > 0 ? false : true;
1777  return(0);
1778  }
1779  }
1780  }
1781 
1782  int numMyLagrangeConstraints = getLocalNumLagrangeConstraints();
1783  int numMySlaveConstraints = getGlobalNumSlaveConstraints();
1784  int numLagrangeConstraints = matrixGraph.getLocalNumLagrangeConstraints();
1785  int numSlaveConstraints = matrixGraph.getGlobalNumSlaveConstraints();
1786 
1787  if (numMyLagrangeConstraints != numLagrangeConstraints ||
1788  numMySlaveConstraints != numSlaveConstraints) {
1789  CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
1790  equivalent = globalCode > 0 ? false : true;
1791  return(0);
1792  }
1793 
1794  localCode = 0;
1795  CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
1796  equivalent = globalCode > 0 ? false : true;
1797  return(0);
1798 }
1799 
1800 //----------------------------------------------------------------------------
1801 int fei::MatrixGraph_Impl2::getNumConnectivityBlocks() const
1802 {
1803  return(connectivityBlocks_.size());
1804 }
1805 
1806 //----------------------------------------------------------------------------
1807 int fei::MatrixGraph_Impl2::getConnectivityBlockIDs(std::vector<int>& blockIDs) const
1808 {
1809  blockIDs.resize(connectivityBlocks_.size());
1810 
1811  std::map<int,fei::ConnectivityBlock*>::const_iterator
1812  cdb_iter = connectivityBlocks_.begin();
1813 
1814  for(size_t i=0; i<blockIDs.size(); ++i, ++cdb_iter) {
1815  int blockID = (*cdb_iter).first;
1816  blockIDs[i] = blockID;
1817  }
1818 
1819  return(0);
1820 }
1821 
1822 //----------------------------------------------------------------------------
1823 int fei::MatrixGraph_Impl2::getNumIDsPerConnectivityList(int blockID) const
1824 {
1825  const fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
1826  if (cblock == NULL) return(-1);
1827 
1828  const fei::Pattern* pattern = cblock->getRowPattern();
1829  return(pattern->getNumIDs());
1830 }
1831 
1832 //----------------------------------------------------------------------------
1833 int fei::MatrixGraph_Impl2::getConnectivityNumIndices(int blockID) const
1834 {
1835  const fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
1836  if (cblock == NULL) return(-1);
1837 
1838  const fei::Pattern* pattern = cblock->getRowPattern();
1839  return( blockEntryGraph_ ?
1840  pattern->getNumIDs() : pattern->getNumIndices());
1841 }
1842 
1843 //----------------------------------------------------------------------------
1844 int fei::MatrixGraph_Impl2::getConnectivityNumIndices(int blockID,
1845  int& numRowIndices,
1846  int& numColIndices)
1847 {
1848  fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
1849  if (cblock == NULL) return(-1);
1850 
1851  fei::Pattern* pattern = cblock->getRowPattern();
1852  numRowIndices = pattern->getNumIndices();
1853  fei::Pattern* colpattern = cblock->isSymmetric() ?
1854  pattern : cblock->getColPattern();
1855  numColIndices = colpattern->getNumIndices();
1856 
1857  return(0);
1858 }
1859 
1860 //----------------------------------------------------------------------------
1861 const fei::ConnectivityBlock* fei::MatrixGraph_Impl2::getConnectivityBlock(int blockID) const
1862 {
1863  std::map<int,fei::ConnectivityBlock*>::const_iterator
1864  c_iter = connectivityBlocks_.find(blockID);
1865  if (c_iter == connectivityBlocks_.end()) return(0);
1866 
1867  return( (*c_iter).second );
1868 }
1869 
1870 //----------------------------------------------------------------------------
1871 fei::ConnectivityBlock* fei::MatrixGraph_Impl2::getConnectivityBlock(int blockID)
1872 {
1873  std::map<int,fei::ConnectivityBlock*>::iterator
1874  c_iter = connectivityBlocks_.find(blockID);
1875  if (c_iter == connectivityBlocks_.end()) return(0);
1876 
1877  return( (*c_iter).second );
1878 }
1879 
1880 //----------------------------------------------------------------------------
1881 int fei::MatrixGraph_Impl2::getConnectivityIndices(int blockID,
1882  int connectivityID,
1883  int indicesAllocLen,
1884  int* indices,
1885  int& numIndices)
1886 {
1887  fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
1888  if (cblock == NULL) return(-1);
1889 
1890  std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
1891  vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
1892 
1893  fei::Pattern* pattern = cblock->getRowPattern();
1894  numIndices = pattern->getNumIndices();
1895 
1896  int len = numIndices > indicesAllocLen ? indicesAllocLen : numIndices;
1897 
1898  int* records = cblock->getRowConnectivity(connectivityID);
1899  if (records == NULL) {
1900  ERReturn(-1);
1901  }
1902 
1903  fei::Pattern::PatternType pType = pattern->getPatternType();
1904 
1905  if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
1906  const int* numFieldsPerID = pattern->getNumFieldsPerID();
1907  const int* fieldIDs = pattern->getFieldIDs();
1908  int totalNumFields = pattern->getTotalNumFields();
1909 
1910  std::vector<int> fieldSizes(totalNumFields);
1911 
1912  for(int j=0; j<totalNumFields; ++j) {
1913  fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
1914  colSpace_.get());
1915  }
1916 
1917  CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
1918  records, pattern->getNumIDs(),
1919  numFieldsPerID,
1920  fieldIDs, &fieldSizes[0],
1921  len, indices, numIndices) );
1922  }
1923  else if (pType == fei::Pattern::SIMPLE) {
1924  const int* fieldIDs = pattern->getFieldIDs();
1925 
1926  int fieldID = fieldIDs[0];
1927  unsigned fieldSize = snl_fei::getFieldSize(fieldID,
1928  rowSpace_.get(),
1929  colSpace_.get());
1930 
1931  CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
1932  records, pattern->getNumIDs(),
1933  fieldID, fieldSize,
1934  len, indices, numIndices) );
1935  }
1936  else if (pType == fei::Pattern::NO_FIELD) {
1937  CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(),
1938  records, pattern->getNumIDs(),
1939  len, indices, numIndices) );
1940  }
1941 
1942  return(0);
1943 }
1944 
1945 //----------------------------------------------------------------------------
1946 int fei::MatrixGraph_Impl2::getConnectivityIndices(int blockID,
1947  int connectivityID,
1948  int rowIndicesAllocLen,
1949  int* rowIndices,
1950  int& numRowIndices,
1951  int colIndicesAllocLen,
1952  int* colIndices,
1953  int& numColIndices)
1954 {
1955  fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
1956  if (cblock == NULL) return(-1);
1957 
1958  std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
1959  vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
1960 
1961  fei::Pattern* pattern = cblock->getRowPattern();
1962  numRowIndices = pattern->getNumIndices();
1963 
1964  int len = numRowIndices > rowIndicesAllocLen ?
1965  rowIndicesAllocLen : numRowIndices;
1966 
1967  int* records = cblock->getRowConnectivity(connectivityID);
1968  if (records == NULL) {
1969  ERReturn(-1);
1970  }
1971 
1972  fei::Pattern::PatternType pType = pattern->getPatternType();
1973 
1974  if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
1975  const int* numFieldsPerID = pattern->getNumFieldsPerID();
1976  const int* fieldIDs = pattern->getFieldIDs();
1977  int totalNumFields = pattern->getTotalNumFields();
1978 
1979  std::vector<int> fieldSizes(totalNumFields);
1980 
1981  for(int j=0; j<totalNumFields; ++j) {
1982  fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
1983  colSpace_.get());
1984  }
1985 
1986  CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
1987  records, pattern->getNumIDs(),
1988  numFieldsPerID,
1989  fieldIDs, &fieldSizes[0],
1990  len, rowIndices, numRowIndices) );
1991  }
1992  else if (pType == fei::Pattern::SIMPLE) {
1993  const int* fieldIDs = pattern->getFieldIDs();
1994 
1995  int fieldID = fieldIDs[0];
1996  unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
1997  colSpace_.get());
1998 
1999  CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
2000  records, pattern->getNumIDs(),
2001  fieldID, fieldSize,
2002  len, rowIndices, numRowIndices) );
2003  }
2004  else if (pType == fei::Pattern::NO_FIELD) {
2005  CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(),
2006  records, pattern->getNumIDs(),
2007  len, rowIndices, numRowIndices) );
2008  }
2009 
2010  fei::Pattern* colpattern = cblock->isSymmetric() ?
2011  pattern : cblock->getColPattern();
2012  pType = colpattern->getPatternType();
2013  numColIndices = colpattern->getNumIndices();
2014  len = numColIndices > colIndicesAllocLen ?
2015  colIndicesAllocLen : numColIndices;
2016 
2017  if (!cblock->isSymmetric()) {
2018  records = cblock->getColConnectivity(connectivityID);
2019  }
2020  if (records == NULL) {
2021  return(-1);
2022  }
2023 
2024  if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
2025  const int* numFieldsPerID = colpattern->getNumFieldsPerID();
2026  const int* fieldIDs = colpattern->getFieldIDs();
2027  int totalNumFields = colpattern->getTotalNumFields();
2028 
2029  std::vector<int> fieldSizes(totalNumFields);
2030 
2031  for(int j=0; j<totalNumFields; ++j) {
2032  fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
2033  colSpace_.get());
2034  }
2035 
2036  CHK_ERR( getConnectivityIndices_multiField(colpattern->getRecordCollections(),
2037  records, colpattern->getNumIDs(),
2038  numFieldsPerID,
2039  fieldIDs, &fieldSizes[0],
2040  len, colIndices, numColIndices) );
2041  }
2042  else if (pType == fei::Pattern::SIMPLE) {
2043  const int* fieldIDs = colpattern->getFieldIDs();
2044 
2045  int fieldID = fieldIDs[0];
2046  unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
2047  colSpace_.get());
2048 
2049  CHK_ERR( getConnectivityIndices_singleField(colpattern->getRecordCollections(),
2050  records, colpattern->getNumIDs(),
2051  fieldID, fieldSize,
2052  len, colIndices, numColIndices) );
2053  }
2054  else if (pType == fei::Pattern::NO_FIELD) {
2055  CHK_ERR( getConnectivityIndices_noField(colpattern->getRecordCollections(),
2056  records, colpattern->getNumIDs(),
2057  len, colIndices, numColIndices) );
2058  }
2059 
2060  return(0);
2061 }
2062 
2063 //----------------------------------------------------------------------------
2064 int fei::MatrixGraph_Impl2::getLocalNumLagrangeConstraints() const
2065 {
2066  return( lagrangeConstraints_.size() );
2067 }
2068 
2069 //----------------------------------------------------------------------------
2070 int fei::MatrixGraph_Impl2::addBlockToGraph_multiField_symmetric(fei::Graph* graph,
2071  fei::ConnectivityBlock* cblock)
2072 {
2073  fei::Pattern* pattern = cblock->getRowPattern();
2074  int j;
2075  int numIDs = pattern->getNumIDs();
2076  int numIndices = blockEntryGraph_ ?
2077  pattern->getNumIDs() : pattern->getNumIndices();
2078  int checkNumIndices = numIndices;
2079  std::vector<int> indices(numIndices);
2080  int* indicesPtr = &indices[0];
2081 
2082  const int* numFieldsPerID = pattern->getNumFieldsPerID();
2083  const int* fieldIDs = pattern->getFieldIDs();
2084  int totalNumFields = pattern->getTotalNumFields();
2085 
2086  std::vector<int> fieldSizes(totalNumFields);
2087 
2088  for(j=0; j<totalNumFields; ++j) {
2089  fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
2090  colSpace_.get());
2091  }
2092 
2093  std::map<int,int>& connIDs = cblock->getConnectivityIDs();
2094  std::vector<int>& values = cblock->getRowConnectivities();
2095 
2096  std::map<int,int>::iterator
2097  c_iter = connIDs.begin(),
2098  c_end = connIDs.end();
2099 
2100  for(; c_iter != c_end; ++c_iter) {
2101  int offset = c_iter->second;
2102  int* records = &values[offset*numIDs];
2103 
2104  CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
2105  records, numIDs,
2106  numFieldsPerID,
2107  fieldIDs,
2108  &fieldSizes[0],
2109  numIndices,
2110  indicesPtr,
2111  checkNumIndices) );
2112 
2113  if (checkNumIndices != numIndices) {
2114  ERReturn(-1);
2115  }
2116 
2117  if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
2118  FEI_OSTREAM& os = *output_stream_;
2119 
2120  const snl_fei::RecordCollection*const* recordColls = pattern->getRecordCollections();
2121  unsigned thisoffset = 0;
2122  for(int ii=0; ii<numIDs; ++ii) {
2123  const fei::Record<int>* record = recordColls[ii]->getRecordWithLocalID(records[ii]);
2124  int ID = record->getID();
2125  os << dbgprefix_<<"scatterIndices: ID=" <<ID<<": ";
2126  int num = pattern->getNumIndicesPerID()[ii];
2127  for(int jj=0; jj<num; ++jj) {
2128  os << indicesPtr[thisoffset++] << " ";
2129  }
2130  os << FEI_ENDL;
2131  }
2132 
2133  for(int ii=0; ii<numIndices; ++ii) {
2134  if (isLogEqn(indicesPtr[ii])) {
2135  os << "adding Symm inds: ";
2136  for(int jj=0; jj<numIndices; ++jj) {
2137  os << indicesPtr[jj] << " ";
2138  }
2139  os << FEI_ENDL;
2140  break;
2141  }
2142  }
2143  }
2144 
2145  //now we have the indices array, so we're ready to push it into
2146  //the graph container
2147  if (numIndices == numIDs || !cblock->isDiagonal()) {
2148  CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr,
2149  cblock->isDiagonal()) );
2150  }
2151  else {
2152  int ioffset = 0;
2153  int* fieldSizesPtr = &fieldSizes[0];
2154  for(int i=0; i<numIDs; ++i) {
2155  CHK_ERR( graph->addSymmetricIndices(fieldSizesPtr[i], &(indicesPtr[ioffset])));
2156  ioffset += fieldSizesPtr[i];
2157  }
2158  }
2159  }
2160 
2161  return(0);
2162 }
2163 
2164 //----------------------------------------------------------------------------
2165 int fei::MatrixGraph_Impl2::addBlockToGraph_multiField_nonsymmetric(fei::Graph* graph,
2166  fei::ConnectivityBlock* cblock)
2167 {
2168  fei::Pattern* pattern = cblock->getRowPattern();
2169  fei::Pattern* colpattern = cblock->getColPattern();
2170  int j;
2171  int numIDs = pattern->getNumIDs();
2172  int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
2173  int checkNumIndices = numIndices;
2174  std::vector<int> indices(numIndices);
2175  int* indicesPtr = &indices[0];
2176 
2177  int numColIDs = colpattern->getNumIDs();
2178  int numColIndices = blockEntryGraph_ ? numColIDs : colpattern->getNumIndices();
2179  int checkNumColIndices = numColIndices;
2180  std::vector<int> colindices(numColIndices);
2181  int* colindicesPtr = &colindices[0];
2182 
2183  const int* numFieldsPerID = pattern->getNumFieldsPerID();
2184  const int* fieldIDs = pattern->getFieldIDs();
2185  int totalNumFields = pattern->getTotalNumFields();
2186 
2187  const int* numFieldsPerColID = colpattern->getNumFieldsPerID();
2188  const int* colfieldIDs = colpattern->getFieldIDs();
2189  int totalNumColFields = colpattern->getTotalNumFields();
2190 
2191  std::vector<int> fieldSizes(totalNumFields);
2192  std::vector<int> colfieldSizes(totalNumColFields);
2193 
2194  for(j=0; j<totalNumFields; ++j) {
2195  fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
2196  colSpace_.get());
2197  }
2198  for(j=0; j<totalNumColFields; ++j) {
2199  colfieldSizes[j] = snl_fei::getFieldSize(colfieldIDs[j], rowSpace_.get(),
2200  colSpace_.get());
2201  }
2202 
2203  std::map<int,int>& connIDs = cblock->getConnectivityIDs();
2204  std::vector<int>& rowrecords = cblock->getRowConnectivities();
2205  std::vector<int>& colrecords = cblock->getColConnectivities();
2206 
2207  std::map<int,int>::iterator
2208  c_iter = connIDs.begin(),
2209  c_end = connIDs.end();
2210 
2211  for(; c_iter != c_end; ++c_iter) {
2212  int offset = c_iter->second;
2213  int* records = &rowrecords[offset*numIDs];
2214 
2215  int* colRecords = &colrecords[offset*numColIDs];
2216 
2217  CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
2218  records, numIDs,
2219  numFieldsPerID,
2220  fieldIDs,
2221  &fieldSizes[0],
2222  numIndices,
2223  indicesPtr,
2224  checkNumIndices) );
2225 
2226  if (checkNumIndices != numIndices) {
2227  ERReturn(-1);
2228  }
2229 
2230 
2231  CHK_ERR( getConnectivityIndices_multiField(colpattern->getRecordCollections(),
2232  colRecords, numColIDs,
2233  numFieldsPerColID,
2234  colfieldIDs,
2235  &colfieldSizes[0],
2236  numColIndices,
2237  colindicesPtr,
2238  checkNumColIndices) );
2239 
2240  if (checkNumColIndices != numColIndices) {
2241  ERReturn(-1);
2242  }
2243 
2244  //now we have the indices arrays, so we're ready to push them into
2245  //the graph container
2246  for(int r=0; r<numIndices; ++r) {
2247  CHK_ERR( graph->addIndices(indicesPtr[r], numColIndices, colindicesPtr));
2248  }
2249  }
2250 
2251  return(0);
2252 }
2253 
2254 //----------------------------------------------------------------------------
2255 int fei::MatrixGraph_Impl2::getConnectivityIndices_multiField(const snl_fei::RecordCollection*const* recordCollections,
2256  int* records,
2257  int numRecords,
2258  const int* numFieldsPerID,
2259  const int* fieldIDs,
2260  const int* fieldSizes,
2261  int indicesAllocLen,
2262  int* indices,
2263  int& numIndices)
2264 {
2265  numIndices = 0;
2266  int fld_offset = 0;
2267 
2268  for(int i=0; i<numRecords; ++i) {
2269  const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
2270  if (record==NULL) continue;
2271 
2272  if (blockEntryGraph_) {
2273  indices[numIndices++] = record->getNumber();
2274  continue;
2275  }
2276 
2277  const fei::FieldMask* fieldMask = record->getFieldMask();
2278  int* eqnNumbers = vspcEqnPtr_ + record->getOffsetIntoEqnNumbers();
2279 
2280  for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
2281  int eqnOffset = 0;
2282  if (!simpleProblem_) {
2283  int err = fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
2284  if (err != 0) {
2285  for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
2286  indices[numIndices++] = -1;
2287  }
2288  continue;
2289  }
2290  }
2291 
2292  for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
2293  indices[numIndices++] = eqnNumbers[eqnOffset+fs];
2294  }
2295 
2296  ++fld_offset;
2297  }
2298  }
2299 
2300  return(0);
2301 }
2302 
2303 //----------------------------------------------------------------------------
2304 int fei::MatrixGraph_Impl2::addBlockToGraph_singleField_symmetric(fei::Graph* graph,
2305  fei::ConnectivityBlock* cblock)
2306 {
2307  fei::Pattern* pattern = cblock->getRowPattern();
2308  int numIDs = pattern->getNumIDs();
2309  int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
2310  int checkNumIndices = numIndices;
2311  std::vector<int> indices(numIndices);
2312  int* indicesPtr = &indices[0];
2313 
2314  const int* fieldIDs = pattern->getFieldIDs();
2315 
2316  int fieldID = fieldIDs[0];
2317  unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
2318  colSpace_.get());
2319 
2320  std::map<int,int>& connIDs = cblock->getConnectivityIDs();
2321  std::vector<int>& rowrecords = cblock->getRowConnectivities();
2322 
2323  std::map<int,int>::iterator
2324  c_iter = connIDs.begin(),
2325  c_end = connIDs.end();
2326 
2327  for(; c_iter != c_end; ++c_iter) {
2328  int offset = c_iter->second;
2329  int* records = &rowrecords[offset*numIDs];
2330 
2331  CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
2332  records, numIDs,
2333  fieldID, fieldSize,
2334  checkNumIndices,
2335  indicesPtr,
2336  numIndices) );
2337 
2338  if (checkNumIndices != numIndices) {
2339  ERReturn(-1);
2340  }
2341 
2342  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
2343  for(int ii=0; ii<numIndices; ++ii) {
2344  if (isLogEqn(indicesPtr[ii])) {
2345  FEI_OSTREAM& os = *output_stream_;
2346  os << "adding Symm inds: ";
2347  for(int jj=0; jj<numIndices; ++jj) {
2348  os << indicesPtr[jj] << " ";
2349  }
2350  os << FEI_ENDL;
2351  break;
2352  }
2353  }
2354  }
2355 
2356  //now we have the indices array, so we're ready to push it into
2357  //the graph container
2358  if (numIndices == numIDs || !cblock->isDiagonal()) {
2359  CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr,
2360  cblock->isDiagonal()));
2361  }
2362  else {
2363  int ioffset = 0;
2364  for(int i=0; i<numIDs; ++i) {
2365  CHK_ERR( graph->addSymmetricIndices(fieldSize, &(indicesPtr[ioffset])));
2366  ioffset += fieldSize;
2367  }
2368  }
2369  }
2370 
2371  return(0);
2372 }
2373 
2374 //----------------------------------------------------------------------------
2375 int fei::MatrixGraph_Impl2::addBlockToGraph_singleField_nonsymmetric(fei::Graph* graph,
2376  fei::ConnectivityBlock* cblock)
2377 {
2378  fei::Pattern* pattern = cblock->getRowPattern();
2379  fei::Pattern* colpattern = cblock->getColPattern();
2380  int numIDs = pattern->getNumIDs();
2381  int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
2382  int checkNumIndices = numIndices;
2383  std::vector<int> indices(numIndices);
2384  int* indicesPtr = &indices[0];
2385 
2386  int numColIDs = colpattern->getNumIDs();
2387  int numColIndices = blockEntryGraph_ ?
2388  numColIDs : colpattern->getNumIndices();
2389  int checkNumColIndices = numColIndices;
2390  std::vector<int> colindices(numColIndices);
2391  int* colindicesPtr = &colindices[0];
2392 
2393  const int* fieldIDs = pattern->getFieldIDs();
2394 
2395  int rowFieldID = fieldIDs[0];
2396  int rowFieldSize = snl_fei::getFieldSize(rowFieldID, rowSpace_.get(),
2397  colSpace_.get());
2398 
2399  std::map<int,int>& connIDs = cblock->getConnectivityIDs();
2400  std::vector<int>& rowrecords = cblock->getRowConnectivities();
2401  std::vector<int>& colrecords = cblock->getColConnectivities();
2402 
2403  int colFieldID = colpattern->getFieldIDs()[0];
2404  int colFieldSize = snl_fei::getFieldSize(colFieldID, rowSpace_.get(),
2405  colSpace_.get());
2406 
2407  std::map<int,int>::iterator
2408  c_iter = connIDs.begin(),
2409  c_end = connIDs.end();
2410 
2411  for(; c_iter != c_end; ++c_iter) {
2412  int offset = c_iter->second;
2413  int* records = &rowrecords[offset*numIDs];
2414 
2415  int* colRecords = &colrecords[offset*numColIDs];
2416 
2417  if (blockEntryGraph_) {
2418  rowSpace_->getGlobalBlkIndicesL(numIDs, pattern->getRecordCollections(),
2419  records, checkNumIndices,
2420  indicesPtr, numIndices);
2421 
2422  colSpace_->getGlobalBlkIndicesL(numColIDs, colpattern->getRecordCollections(),
2423  colRecords, checkNumColIndices,
2424  colindicesPtr, numColIndices);
2425  }
2426  else {
2427  rowSpace_->getGlobalIndicesL(numIDs, pattern->getRecordCollections(),
2428  records, rowFieldID,
2429  rowFieldSize, checkNumIndices,
2430  indicesPtr, numIndices);
2431 
2432  colSpace_->getGlobalIndicesL(numColIDs, colpattern->getRecordCollections(),
2433  colRecords, colFieldID,
2434  colFieldSize, checkNumColIndices,
2435  colindicesPtr, numColIndices);
2436  }
2437 
2438  if (checkNumIndices != numIndices || checkNumColIndices != numColIndices) {
2439  ERReturn(-1);
2440  }
2441 
2442  for(int r=0; r<numIndices; ++r) {
2443  CHK_ERR( graph->addIndices(indicesPtr[r], numColIndices, colindicesPtr));
2444  }
2445  }
2446 
2447  return(0);
2448 }
2449 
2450 //----------------------------------------------------------------------------
2451 int fei::MatrixGraph_Impl2::getConnectivityIndices_singleField(const snl_fei::RecordCollection*const* recordCollections,
2452  int* records,
2453  int numRecords,
2454  int fieldID,
2455  int fieldSize,
2456  int indicesAllocLen,
2457  int* indices,
2458  int& numIndices)
2459 {
2460  numIndices = 0;
2461 
2462  for(int i=0; i<numRecords; ++i) {
2463  if (numIndices == indicesAllocLen) break;
2464 
2465  const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
2466 
2467  if (blockEntryGraph_) {
2468  indices[numIndices++] = record->getNumber();
2469  continue;
2470  }
2471 
2472  int* eqnNumbers = vspcEqnPtr_+record->getOffsetIntoEqnNumbers();
2473 
2474  int eqnOffset = 0;
2475  if (!simpleProblem_) {
2476  const fei::FieldMask* fieldMask = record->getFieldMask();
2477  int err = fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
2478  if (err != 0) {
2479  indices[numIndices++] = -1;
2480  if (fieldSize > 1) {
2481  for(int fs=1; fs<fieldSize; ++fs) {
2482  indices[numIndices++] = -1;
2483  }
2484  }
2485  continue;
2486  }
2487  }
2488 
2489  indices[numIndices++] = eqnNumbers[eqnOffset];
2490  if (fieldSize > 1) {
2491  for(int fs=1; fs<fieldSize; ++fs) {
2492  indices[numIndices++] = eqnOffset >= 0 ? eqnNumbers[eqnOffset+fs] : -1;
2493  }
2494  }
2495  }
2496 
2497  return(0);
2498 }
2499 
2500 //----------------------------------------------------------------------------
2501 int fei::MatrixGraph_Impl2::getConnectivityIndices_noField(const snl_fei::RecordCollection*const* recordCollections,
2502  int* records,
2503  int numRecords,
2504  int indicesAllocLen,
2505  int* indices,
2506  int& numIndices)
2507 {
2508  numIndices = 0;
2509 
2510  for(int i=0; i<numRecords; ++i) {
2511 
2512  const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
2513  int* eqnNumbers = vspcEqnPtr_+record->getOffsetIntoEqnNumbers();
2514 
2515  if (blockEntryGraph_) {
2516  indices[numIndices++] = record->getNumber();
2517  }
2518  else {
2519  indices[numIndices++] = eqnNumbers[0];
2520  }
2521  }
2522 
2523  return(0);
2524 }
2525 
2526 //----------------------------------------------------------------------------
2527 int fei::MatrixGraph_Impl2::addBlockToGraph_noField_symmetric(fei::Graph* graph,
2528  fei::ConnectivityBlock* cblock)
2529 {
2530  fei::Pattern* pattern = cblock->getRowPattern();
2531  int numIDs = pattern->getNumIDs();
2532  int numIndices = pattern->getNumIndices();
2533  std::vector<int> indices(numIndices);
2534  int* indicesPtr = &indices[0];
2535 
2536  std::map<int,int>& connIDs = cblock->getConnectivityIDs();
2537  std::vector<int>& rowrecords = cblock->getRowConnectivities();
2538 
2539  std::map<int,int>::iterator
2540  c_iter = connIDs.begin(),
2541  c_end = connIDs.end();
2542 
2543  for(; c_iter != c_end; ++c_iter) {
2544  int offset = c_iter->second;
2545  int* records = &rowrecords[offset*numIDs];
2546 
2547  int checkNumIndices;
2548  CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(), records, numIDs, numIndices,
2549  indicesPtr, checkNumIndices) );
2550 
2551  if (checkNumIndices != numIndices) {
2552  ERReturn(-1);
2553  }
2554 
2555  if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
2556  for(int ii=0; ii<numIndices; ++ii) {
2557  if (isLogEqn(indicesPtr[ii])) {
2558  FEI_OSTREAM& os = *output_stream_;
2559  os << "adding Symm inds: ";
2560  for(int jj=0; jj<numIndices; ++jj) {
2561  os << indicesPtr[jj] << " ";
2562  }
2563  os << FEI_ENDL;
2564  break;
2565  }
2566  }
2567  }
2568 
2569  //now we have the indices array, so we're ready to push it into
2570  //the graph container
2571  CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr) );
2572  }
2573 
2574  return(0);
2575 }
2576 
2577 //----------------------------------------------------------------------------
2578 int fei::MatrixGraph_Impl2::addBlockToGraph_sparse(fei::Graph* graph,
2579  fei::ConnectivityBlock* cblock)
2580 {
2581  std::vector<int> row_indices;
2582  std::vector<int> indices;
2583 
2584  fei::Pattern* pattern = cblock->getRowPattern();
2585  const snl_fei::RecordCollection*const* recordCollections = pattern->getRecordCollections();
2586 
2587  std::map<int,int>& connIDs = cblock->getConnectivityIDs();
2588  std::vector<int>& connOffsets = cblock->getConnectivityOffsets();
2589  int* rowrecords = &(cblock->getRowConnectivities()[0]);
2590  int* colrecords = &(cblock->getColConnectivities()[0]);
2591  bool haveField = cblock->haveFieldID();
2592  int fieldID = cblock->fieldID();
2593  int fieldSize = 1;
2594 
2595  if (haveField) {
2596  fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
2597  colSpace_.get());
2598  }
2599 
2600  std::map<int,int>::iterator
2601  c_iter = connIDs.begin(),
2602  c_end = connIDs.end();
2603 
2604  for(; c_iter != c_end; ++c_iter) {
2605  int offset = c_iter->second;
2606  int rowlen = connOffsets[offset+1] - offset;
2607 
2608  int* records = &(rowrecords[offset]);
2609 
2610  int checkNumIndices;
2611  row_indices.resize(fieldSize);
2612 
2613  if (haveField) {
2614  CHK_ERR( getConnectivityIndices_singleField(recordCollections, records, 1,
2615  fieldID, fieldSize,
2616  fieldSize,
2617  &row_indices[0],
2618  checkNumIndices) );
2619  }
2620  else {
2621  CHK_ERR( getConnectivityIndices_noField(recordCollections, records, 1, fieldSize,
2622  &row_indices[0],
2623  checkNumIndices) );
2624  }
2625 
2626  indices.resize(fieldSize*rowlen);
2627  int* indicesPtr = &indices[0];
2628  int* crecords = &(colrecords[offset]);
2629 
2630  if (haveField) {
2631  CHK_ERR( getConnectivityIndices_singleField(recordCollections, crecords, rowlen,
2632  fieldID, fieldSize,
2633  fieldSize*rowlen,
2634  indicesPtr,
2635  checkNumIndices) );
2636  }
2637  else {
2638  CHK_ERR( getConnectivityIndices_noField(recordCollections, crecords, rowlen, rowlen,
2639  indicesPtr, checkNumIndices) );
2640  }
2641  if (checkNumIndices != rowlen) {
2642  ERReturn(-1);
2643  }
2644 
2645  //now we have the indices, so we're ready to push them into
2646  //the graph container
2647  int* row_ind_ptr = &row_indices[0];
2648  for(int i=0; i<fieldSize; ++i) {
2649  CHK_ERR( graph->addIndices(row_ind_ptr[i], fieldSize*rowlen,
2650  indicesPtr) );
2651  }
2652  }
2653 
2654  return(0);
2655 }
2656 
2657 //----------------------------------------------------------------------------
2658 void fei::MatrixGraph_Impl2::setName(const char* name)
2659 {
2660  if (name == NULL) return;
2661 
2662  if (name_ == name) return;
2663 
2664  name_ = name;
2665  dbgprefix_ = "MatGraph_"+name_+": ";
2666 }
2667 
2668 //----------------------------------------------------------------------------
2669 void fei::MatrixGraph_Impl2::setIndicesMode(int mode)
2670 {
2671  if (mode == BLOCK_ENTRY_GRAPH) {
2672  blockEntryGraph_ = true;
2673  return;
2674  }
2675 
2676  if (mode == POINT_ENTRY_GRAPH) {
2677  blockEntryGraph_ = false;
2678  return;
2679  }
2680 
2681  voidERReturn;
2682 }
2683 
2684 //----------------------------------------------------------------------------
2685 fei::SharedPtr<fei::FillableMat>
2686 fei::MatrixGraph_Impl2::getSlaveDependencyMatrix()
2687 {
2688  return( D_ );
2689 }
2690 
fei::MatrixGraph_Impl2::initLagrangeConstraint
int initLagrangeConstraint(int constraintID, int constraintIDType, int numIDs, const int *idTypes, const int *IDs, const int *fieldIDs)
Definition: fei_MatrixGraph_Impl2.cpp:940
fei::ConnectivityBlock::isDiagonal
bool isDiagonal() const
fei::GlobalSum
int GlobalSum(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
Definition: fei_CommUtils.hpp:189
fei::Record::getNumber
GlobalIDType getNumber() const
fei::MatrixGraph_Impl2::definePattern
int definePattern(int numIDs, int idType)
Definition: fei_MatrixGraph_Impl2.cpp:266
fei::MatrixGraph_Impl2::getLocalNumLagrangeConstraints
int getLocalNumLagrangeConstraints() const
Definition: fei_MatrixGraph_Impl2.cpp:2064
fei::Pattern::getRecordCollections
snl_fei::RecordCollection *const * getRecordCollections() const
fei::Param::getType
ParamType getType() const
Definition: fei_Param.hpp:98
fei::ConnectivityBlock::getRowConnectivity
const int * getRowConnectivity(int ID) const
Definition: fei_ConnectivityBlock.cpp:175
fei::Pattern::getIDTypes
const int * getIDTypes() const
fei::MatrixGraph::getGlobalNumSlaveConstraints
virtual int getGlobalNumSlaveConstraints() const =0
fei::Graph::writeRemoteGraph
virtual int writeRemoteGraph(FEI_OSTREAM &os)=0
fei::ConnectivityBlock::getConnectivityIDs
const std::map< int, int > & getConnectivityIDs() const
snl_fei::Constraint::getMasterFieldIDs
std::vector< int > & getMasterFieldIDs()
fei::Pattern::getNumIndices
int getNumIndices() const
fei::ParameterSet::get
const Param * get(const char *name) const
Definition: fei_ParameterSet.hpp:260
fei::MatrixGraph_Impl2::initComplete
int initComplete()
Definition: fei_MatrixGraph_Impl2.cpp:1164
snl_fei::Constraint::getRHSValue
double getRHSValue() const
fei::Pattern::getNumFieldsPerID
const int * getNumFieldsPerID() const
fei::MatrixGraph_Impl2::createGraph
fei::SharedPtr< fei::SparseRowGraph > createGraph(bool blockEntryGraph, bool localRowGraph_includeSharedRows=false)
Definition: fei_MatrixGraph_Impl2.cpp:1300
fei::MatrixGraph_Impl2::getPattern
fei::Pattern * getPattern(int patternID)
Definition: fei_MatrixGraph_Impl2.cpp:362
fei::MatrixGraph_Impl2::createSlaveMatrices
int createSlaveMatrices()
Definition: fei_MatrixGraph_Impl2.cpp:1383
fei::MatrixGraph_Impl2::getConnectivityIndices
int getConnectivityIndices(int blockID, int connectivityID, int indicesAllocLen, int *indices, int &numIndices)
Definition: fei_MatrixGraph_Impl2.cpp:1881
fei::utils::string_to_output_level
fei::OutputLevel string_to_output_level(const std::string &str)
Definition: fei_utils.cpp:58
fei::Param::getBoolValue
bool getBoolValue() const
Definition: fei_Param.hpp:122
fei::ConstraintType
snl_fei::Constraint< fei::Record< int > * > ConstraintType
fei::Graph::gatherFromOverlap
virtual int gatherFromOverlap()=0
fei::MatrixGraph::getConnectivityBlock
virtual const fei::ConnectivityBlock * getConnectivityBlock(int blockID) const =0
fei::Pattern::getFieldIDs
const int * getFieldIDs() const
fei::Pattern::getTotalNumFields
int getTotalNumFields() const
fei::MatrixGraph_Impl2::initConnectivityBlock
int initConnectivityBlock(int blockID, int numConnectivityLists, int patternID, bool diagonal=false)
Definition: fei_MatrixGraph_Impl2.cpp:376
fei::ConnectivityBlock::getColConnectivities
std::vector< int > & getColConnectivities()
fei::MatrixGraph_Impl2::getNumConnectivityBlocks
int getNumConnectivityBlocks() const
Definition: fei_MatrixGraph_Impl2.cpp:1801
snl_fei::Constraint::setEqnNumber
void setEqnNumber(int eqn)
fei::VectorSpace::addDOFs
int addDOFs(int fieldID, int idType, int numIDs, const int *IDs)
Definition: fei_VectorSpace.cpp:325
snl_fei::Constraint::structurallySame
bool structurallySame(const Constraint< RecordType > &rhs)
Definition: snl_fei_Constraint.hpp:340
fei::Graph::writeLocalGraph
virtual int writeLocalGraph(FEI_OSTREAM &os, bool debug=false, bool prefixLinesWithPoundSign=true)=0
fei::SharedPtr::reset
void reset(T *p=0)
fei::Param::getIntValue
int getIntValue() const
Definition: fei_Param.hpp:116
snl_fei::Constraint::getSlaveFieldID
int getSlaveFieldID() const
fei::Graph_Impl
Definition: fei_Graph_Impl.hpp:21
fei::ConnectivityBlock::getRowConnectivities
std::vector< int > & getRowConnectivities()
snl_fei::Constraint::setBlkEqnNumber
void setBlkEqnNumber(int blkEqn)
fei::MatrixGraph_Impl2::getConnectivityBlockIDs
int getConnectivityBlockIDs(std::vector< int > &blockIDs) const
Definition: fei_MatrixGraph_Impl2.cpp:1807
snl_fei::BlkSizeMsgHandler
Definition: snl_fei_BlkSizeMsgHandler.hpp:21
fei::Graph::addIndices
virtual int addIndices(int row, int len, const int *indices)=0
snl_fei::Constraint::getMasterWeights
std::vector< double > & getMasterWeights()
fei::MatrixGraph_Impl2::getSlaveDependencyMatrix
fei::SharedPtr< FillableMat > getSlaveDependencyMatrix()
Definition: fei_MatrixGraph_Impl2.cpp:2686
fei::ConnectivityBlock::getColPattern
const fei::Pattern * getColPattern() const
fei::ConnectivityBlock::isSymmetric
bool isSymmetric() const
snl_fei::Constraint::getConstraintID
int getConstraintID() const
fei::MatrixGraph_Impl2::MatrixGraph_Impl2
MatrixGraph_Impl2(fei::SharedPtr< fei::VectorSpace > rowSpace, fei::SharedPtr< fei::VectorSpace > colSpace, const char *name=NULL)
Definition: fei_MatrixGraph_Impl2.cpp:93
fei::MatrixGraph_Impl2::setIndicesMode
void setIndicesMode(int mode)
Definition: fei_MatrixGraph_Impl2.cpp:2669
fei::MatrixGraph_Impl2::getConnectivityNumIndices
int getConnectivityNumIndices(int blockID) const
Definition: fei_MatrixGraph_Impl2.cpp:1833
fei::MatrixGraph_Impl2::setRowSpace
void setRowSpace(fei::SharedPtr< fei::VectorSpace > rowSpace)
Definition: fei_MatrixGraph_Impl2.cpp:222
fei::MatrixGraph_Impl2::getConnectivityBlock
const fei::ConnectivityBlock * getConnectivityBlock(int blockID) const
Definition: fei_MatrixGraph_Impl2.cpp:1861
fei::MatrixGraph_Impl2::initSlaveConstraint
int initSlaveConstraint(int numIDs, const int *idTypes, const int *IDs, const int *fieldIDs, int offsetOfSlave, int offsetIntoSlaveField, const double *weights, double rhsValue)
Definition: fei_MatrixGraph_Impl2.cpp:1047
fei::MatrixGraph_Impl2
Definition: fei_MatrixGraph_Impl2.hpp:33
fei::MatrixGraph_Impl2::getNumIDsPerConnectivityList
int getNumIDsPerConnectivityList(int blockID) const
Definition: fei_MatrixGraph_Impl2.cpp:1823
fei::MatrixGraph_Impl2::Factory::createMatrixGraph
virtual fei::SharedPtr< fei::MatrixGraph > createMatrixGraph(fei::SharedPtr< fei::VectorSpace > rowSpace, fei::SharedPtr< fei::VectorSpace > columnSpace, const char *name)
Definition: fei_MatrixGraph_Impl2.cpp:81
fei::MatrixGraph_Impl2::getPenaltyConstraint
ConstraintType * getPenaltyConstraint(int constraintID)
Definition: fei_MatrixGraph_Impl2.cpp:1154
fei::LogManager::setOutputLevel
void setOutputLevel(OutputLevel olevel)
Definition: fei_LogManager.cpp:35
snl_fei::RecordCollection
Definition: snl_fei_RecordCollection.hpp:32
fei::SharedPtr::get
T * get() const
fei::Pattern::getPatternType
PatternType getPatternType() const
fei::Allreduce
int Allreduce(MPI_Comm comm, bool localBool, bool &globalBool)
Definition: fei_CommUtils.cpp:265
fei::MatrixGraph_Impl2::getSlaveConstraint
ConstraintType * getSlaveConstraint(int constraintID)
Definition: fei_MatrixGraph_Impl2.cpp:1143
fei::MatrixGraph_Impl2::~MatrixGraph_Impl2
virtual ~MatrixGraph_Impl2()
Definition: fei_MatrixGraph_Impl2.cpp:141
fei::MatrixGraph_Impl2::getPatternNumIndices
int getPatternNumIndices(int patternID, int &numIndices)
Definition: fei_MatrixGraph_Impl2.cpp:863
fei::MatrixGraph_Impl2::getConstrainedIndices
void getConstrainedIndices(std::vector< int > &crindices) const
Definition: fei_MatrixGraph_Impl2.cpp:1577
fei::createSparseRowGraph
fei::SharedPtr< fei::SparseRowGraph > createSparseRowGraph(const std::vector< snl_fei::RaggedTable< MAP_TYPE, SET_TYPE > * > &tables)
Definition: fei_TemplateUtils.hpp:217
fei::MatrixGraph_Impl2::initConnectivity
int initConnectivity(int blockID, int connectivityID, const int *connectedIdentifiers)
Definition: fei_MatrixGraph_Impl2.cpp:489
fei::Pattern::getNumIndicesPerID
const int * getNumIndicesPerID() const
snl_fei::Constraint::getOffsetIntoSlaveField
int getOffsetIntoSlaveField() const
fei::MatrixGraph_Impl2::getPatternIndices
int getPatternIndices(int patternID, const int *IDs, std::vector< int > &indices)
Definition: fei_MatrixGraph_Impl2.cpp:877
fei::Graph::addSymmetricIndices
virtual int addSymmetricIndices(int numIndices, int *indices, bool diagonal=false)=0
fei::console_out
std::ostream & console_out()
Definition: fei_console_ostream.cpp:26
fei::MatrixGraph_Impl2::setParameters
void setParameters(const fei::ParameterSet &params)
Definition: fei_MatrixGraph_Impl2.cpp:165
fei::GlobalMax
int GlobalMax(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
Definition: fei_CommUtils.hpp:91
fei::Pattern::getNumIDs
int getNumIDs() const
fei::MatrixGraph_Impl2::compareStructure
int compareStructure(const fei::MatrixGraph &matrixGraph, bool &equivalent) const
Definition: fei_MatrixGraph_Impl2.cpp:1728
fei::Param::getStringValue
const std::string & getStringValue() const
Definition: fei_Param.hpp:104
fei::MatrixGraph_Impl2::getLagrangeConstraint
ConstraintType * getLagrangeConstraint(int constraintID)
Definition: fei_MatrixGraph_Impl2.cpp:1132
fei::MatrixGraph_Impl2::getConstraintConnectivityIndices
int getConstraintConnectivityIndices(ConstraintType *cr, std::vector< int > &globalIndices)
Definition: fei_MatrixGraph_Impl2.cpp:1659
fei::MatrixGraph_Impl2::getReducer
fei::SharedPtr< fei::Reducer > getReducer()
Definition: fei_MatrixGraph_Impl2.cpp:1564
fei::localProc
int localProc(MPI_Comm comm)
Definition: fei_CommUtils.cpp:13
fei::LogManager::getLogManager
static LogManager & getLogManager()
Definition: fei_LogManager.cpp:24
fei::FieldMask::getFieldEqnOffset
int getFieldEqnOffset(int fieldID, int &offset) const
Definition: fei_FieldMask.cpp:56
fei::Record::getID
GlobalIDType getID() const
fei::MatrixGraph_Impl2::initPenaltyConstraint
int initPenaltyConstraint(int constraintID, int constraintIDType, int numIDs, const int *idTypes, const int *IDs, const int *fieldIDs)
Definition: fei_MatrixGraph_Impl2.cpp:988
fei::VectorSpace::getFieldSize
unsigned getFieldSize(int fieldID)
Definition: fei_VectorSpace.cpp:1365
snl_fei::Constraint::getMasterRecordCollections
std::vector< snl_fei::RecordCollection * > & getMasterRecordCollections()
snl_fei::Constraint::getIDType
int getIDType() const
snl_fei::Constraint::getMasters
std::vector< int > & getMasters()
fei::MatrixGraph::getLocalNumLagrangeConstraints
virtual int getLocalNumLagrangeConstraints() const =0
fei::destroyValues
void destroyValues(MAP_TYPE &map_obj)
Definition: fei_TemplateUtils.hpp:134
fei::ConnectivityBlock::getConnectivityOffsets
std::vector< int > & getConnectivityOffsets()
fei::ConnectivityBlock::getColConnectivity
const int * getColConnectivity(int ID) const
Definition: fei_ConnectivityBlock.cpp:203
fei::Record::getOffsetIntoEqnNumbers
int getOffsetIntoEqnNumbers() const
fei::MatrixGraph::getNumConnectivityBlocks
virtual int getNumConnectivityBlocks() const =0
fei::MatrixGraph::getConnectivityBlockIDs
virtual int getConnectivityBlockIDs(std::vector< int > &blockIDs) const =0
fei::Record::getFieldMask
fei::FieldMask * getFieldMask()
fei::MatrixGraph_Impl2::hasSlaveDof
bool hasSlaveDof(int ID, int idType)
Definition: fei_MatrixGraph_Impl2.cpp:1026
fei::VectorSpace::getRecordCollection
int getRecordCollection(int idType, snl_fei::RecordCollection *&records)
Definition: fei_VectorSpace.cpp:1701
snl_fei::RecordCollection::getRecordWithID
fei::Record< int > * getRecordWithID(int ID)
Definition: snl_fei_RecordCollection.cpp:254
fei::ConnectivityBlock::setIsDiagonal
void setIsDiagonal(bool flag)
snl_fei::Constraint
Definition: snl_fei_Constraint.hpp:23
fei::MatrixGraph_Impl2::getRemotelyOwnedGraphRows
fei::SharedPtr< fei::SparseRowGraph > getRemotelyOwnedGraphRows()
Definition: fei_MatrixGraph_Impl2.cpp:1571
snl_fei::Constraint::getSlave
RecordType getSlave()
fei::ConnectivityBlock::getRowPattern
const fei::Pattern * getRowPattern() const
snl_fei::Constraint::getMasterIDTypes
std::vector< int > & getMasterIDTypes()
fei::numProcs
int numProcs(MPI_Comm comm)
Definition: fei_CommUtils.cpp:25
fei::MatrixGraph_Impl2::setColumnSpace
void setColumnSpace(fei::SharedPtr< fei::VectorSpace > columnSpace)
Definition: fei_MatrixGraph_Impl2.cpp:230
Generated on Fri Sep 1 2023 07:58:42 for FEI by   doxygen 1.8.5