FEI Version of the Day
fei_MatrixGraph_Impl2.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include <fei_sstream.hpp>
00010 
00011 #include <limits>
00012 #include <cmath>
00013 
00014 #include <fei_MatrixGraph_Impl2.hpp>
00015 
00016 #include <fei_utils.hpp>
00017 #include "fei_TemplateUtils.hpp"
00018 
00019 #include <fei_Pattern.hpp>
00020 #include <fei_LogManager.hpp>
00021 #include <fei_TemplateUtils.hpp>
00022 #include <fei_impl_utils.hpp>
00023 #include <snl_fei_Utils.hpp>
00024 #include <fei_FieldMask.hpp>
00025 #include <fei_Record.hpp>
00026 #include <snl_fei_RecordCollection.hpp>
00027 #include <fei_VectorSpace.hpp>
00028 #include <fei_ParameterSet.hpp>
00029 #include <fei_ostream_ops.hpp>
00030 #include <fei_Reducer.hpp>
00031 #include <fei_GraphReducer.hpp>
00032 #include <fei_ConnectivityBlock.hpp>
00033 #include <snl_fei_BlkSizeMsgHandler.hpp>
00034 #include <fei_Graph_Impl.hpp>
00035 #include <snl_fei_Constraint.hpp>
00036 
00037 #include <fei_EqnBuffer.hpp>
00038 #include <fei_EqnCommMgr.hpp>
00039 #include <SNL_FEI_Structure.hpp>
00040 
00041 #undef fei_file
00042 #define fei_file "fei_MatrixGraph.cpp"
00043 
00044 #include <fei_ErrMacros.hpp>
00045 
00046 namespace snl_fei {
00047 static unsigned getFieldSize(int fieldID,
00048                       fei::VectorSpace* space1,
00049                       fei::VectorSpace* space2)
00050 {
00051   unsigned fieldsize = 0;
00052   bool foundfield = false;
00053   if (space1 != NULL) {
00054     try {
00055       fieldsize = space1->getFieldSize(fieldID);
00056       foundfield = true;
00057     }
00058     catch (...) {
00059       foundfield = false;
00060     }
00061   }
00062 
00063   if (!foundfield && space2 != NULL) {
00064     try {
00065       fieldsize = space2->getFieldSize(fieldID);
00066     }
00067     catch (std::runtime_error& exc) {
00068       std::string msg("snl_fei::getFieldSize: ");
00069       msg += exc.what();
00070       throw std::runtime_error(msg);
00071     }
00072   }
00073 
00074   return(fieldsize);
00075 }
00076 
00077 }//namespace snl_fei
00078 
00079 //------------------------------------------------------------------------------
00080 fei::SharedPtr<fei::MatrixGraph>
00081 fei::MatrixGraph_Impl2::Factory::createMatrixGraph(fei::SharedPtr<fei::VectorSpace> rowSpace,
00082                                              fei::SharedPtr<fei::VectorSpace> columnSpace,
00083                                              const char* name)
00084 {
00085   fei::SharedPtr<fei::MatrixGraph> mgptr(new fei::MatrixGraph_Impl2(rowSpace,
00086                                                               columnSpace,
00087                                                               name));
00088 
00089   return(mgptr);
00090 }
00091 
00092 //=====Constructor==============================================================
00093 fei::MatrixGraph_Impl2::MatrixGraph_Impl2(fei::SharedPtr<fei::VectorSpace> rowSpace,
00094                               fei::SharedPtr<fei::VectorSpace> colSpace,
00095                               const char* name)
00096  : comm_(rowSpace->getCommunicator()),
00097    rowSpace_(rowSpace),
00098    colSpace_(colSpace),
00099    haveRowSpace_(false),
00100    haveColSpace_(false),
00101    symmetric_(false),
00102    remotelyOwnedGraphRows_(NULL),
00103    simpleProblem_(false),
00104    blockEntryGraph_(false),
00105    patterns_(),
00106    connectivityBlocks_(),
00107    arbitraryBlockCounter_(1),
00108    sparseBlocks_(),
00109    lagrangeConstraints_(),
00110    penaltyConstraints_(),
00111    slaveConstraints_(),
00112    ptEqualBlk_(false),
00113    newSlaveData_(false),
00114    localNumSlaves_(0),
00115    globalNumSlaves_(0),
00116    D_(NULL),
00117    g_(),
00118    g_nonzero_(false),
00119    reducer_(),
00120    name_(),
00121    dbgprefix_("MatGrph: "),
00122    tmpIntArray1_(),
00123    tmpIntArray2_()
00124 {
00125   localProc_ = fei::localProc(comm_);
00126   numProcs_ = fei::numProcs(comm_);
00127 
00128   if (rowSpace_.get() == NULL) {
00129     voidERReturn;
00130   }
00131   else haveRowSpace_ = true;
00132 
00133   if (colSpace_.get() != NULL) haveColSpace_ = true;
00134   else colSpace_ = rowSpace_;
00135 
00136   setName(name);
00137 }
00138 
00139 //------------------------------------------------------------------------------
00140 fei::MatrixGraph_Impl2::~MatrixGraph_Impl2()
00141 {
00142   fei::destroyValues(patterns_);
00143   patterns_.clear();
00144 
00145   fei::destroyValues(connectivityBlocks_);
00146   connectivityBlocks_.clear();
00147 
00148   int i, len = sparseBlocks_.size();
00149   for(i=0; i<len; ++i) {
00150     delete sparseBlocks_[i];
00151   }
00152 
00153   fei::destroyValues(lagrangeConstraints_);
00154   lagrangeConstraints_.clear();
00155 
00156   fei::destroyValues(penaltyConstraints_);
00157   penaltyConstraints_.clear();
00158 
00159   fei::destroyValues(slaveConstraints_);
00160   slaveConstraints_.clear();
00161 }
00162 
00163 //------------------------------------------------------------------------------
00164 void fei::MatrixGraph_Impl2::setParameters(const fei::ParameterSet& params)
00165 {
00166   const fei::Param* param = 0;
00167   fei::Param::ParamType ptype = fei::Param::BAD_TYPE;
00168 
00169   param = params.get("FEI_OUTPUT_LEVEL");
00170   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00171   if (ptype == fei::Param::STRING) {
00172     fei::LogManager& log_manager = fei::LogManager::getLogManager();
00173     log_manager.setOutputLevel(param->getStringValue().c_str());
00174     setOutputLevel(fei::utils::string_to_output_level(param->getStringValue()));
00175   }
00176 
00177   param = params.get("FEI_LOG_EQN");
00178   ptype =  param != NULL ? param->getType() : fei::Param::BAD_TYPE;  
00179   if (ptype == fei::Param::INT) {
00180     addLogEqn(param->getIntValue());
00181   }
00182 
00183   param = params.get("FEI_LOG_ID");
00184   ptype =  param != NULL ? param->getType() : fei::Param::BAD_TYPE;  
00185   if (ptype == fei::Param::INT) {
00186     addLogID(param->getIntValue());
00187   }
00188 
00189   param = params.get("name");
00190   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00191   if (ptype == fei::Param::STRING) {
00192     setName(param->getStringValue().c_str());
00193   }
00194 
00195   param = params.get("BLOCK_GRAPH");
00196   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00197   if (ptype != fei::Param::BAD_TYPE) {
00198     blockEntryGraph_ = true;
00199   }
00200 
00201   param = params.get("BLOCK_MATRIX");
00202   ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
00203   if (ptype != fei::Param::BAD_TYPE) {
00204     if (ptype == fei::Param::BOOL) {
00205       blockEntryGraph_ = param->getBoolValue();
00206     }
00207     else {
00208       blockEntryGraph_ = true;
00209     }
00210   }
00211 }
00212 
00213 //----------------------------------------------------------------------------
00214 void fei::MatrixGraph_Impl2::setRowSpace(fei::SharedPtr<fei::VectorSpace> rowSpace)
00215 {
00216   rowSpace_ = rowSpace;
00217   haveRowSpace_ = true;
00218   if (!haveColSpace_) symmetric_ = true;
00219 }
00220 
00221 //----------------------------------------------------------------------------
00222 void fei::MatrixGraph_Impl2::setColumnSpace(fei::SharedPtr<fei::VectorSpace> colSpace)
00223 {
00224   colSpace_ = colSpace;
00225   haveColSpace_ = true;
00226   if (colSpace_ == rowSpace_) symmetric_ = true;
00227   else symmetric_ = false;
00228 }
00229 
00230 //----------------------------------------------------------------------------
00231 int fei::MatrixGraph_Impl2::addPattern(fei::Pattern* pattern)
00232 {
00233   std::map<int,fei::Pattern*>::iterator
00234     p_iter = patterns_.begin(), p_iter_end = patterns_.end();
00235 
00236   bool matches_existing_pattern = false;
00237   int pattern_id = -1;
00238   while(p_iter != p_iter_end && !matches_existing_pattern) {
00239     if (*pattern == *(p_iter->second)) {
00240       matches_existing_pattern = true;
00241       pattern_id = p_iter->first;
00242     }
00243     else ++p_iter;
00244   }
00245 
00246   if (!matches_existing_pattern) {
00247     pattern_id = patterns_.size();
00248     patterns_.insert(std::make_pair(pattern_id, pattern));
00249   }
00250   else {
00251     delete pattern;
00252   }
00253 
00254   return pattern_id;
00255 }
00256 
00257 //----------------------------------------------------------------------------
00258 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
00259                                     int idType)
00260 {
00261   snl_fei::RecordCollection* rec_coll = NULL;
00262   rowSpace_->getRecordCollection(idType, rec_coll);
00263 
00264   fei::Pattern* pattern = new fei::Pattern(numIDs, idType, rec_coll);
00265   return addPattern(pattern);
00266 }
00267 
00268 //----------------------------------------------------------------------------
00269 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
00270                                       int idType,
00271                                       int fieldID)
00272 {
00273   unsigned fieldsize = 0;
00274   try {
00275     fieldsize = snl_fei::getFieldSize(fieldID, rowSpace_.get(), colSpace_.get());
00276   }
00277   catch (std::runtime_error& exc) {
00278     FEI_OSTRINGSTREAM osstr;
00279     osstr << "fei::MatrixGraph_Impl2::definePattern caught error: "<<exc.what();
00280     throw std::runtime_error(osstr.str());
00281   }
00282 
00283   snl_fei::RecordCollection* rec_coll = NULL;
00284   rowSpace_->getRecordCollection(idType, rec_coll);
00285 
00286   fei::Pattern* pattern =
00287     new fei::Pattern(numIDs, idType, rec_coll, fieldID, fieldsize);
00288   return addPattern(pattern);
00289 }
00290 
00291 //----------------------------------------------------------------------------
00292 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
00293                                      int idType,
00294                                      const int* numFieldsPerID,
00295                                      const int* fieldIDs)
00296 { 
00297   std::vector<int> fieldSizes;
00298   try {
00299     int offset = 0;
00300     for(int i=0; i<numIDs; ++i) {
00301       for(int j=0; j<numFieldsPerID[i]; ++j) {
00302         fieldSizes.push_back(snl_fei::getFieldSize(fieldIDs[offset++],
00303                                               rowSpace_.get(), colSpace_.get()));
00304       }
00305     }
00306   }
00307   catch (std::runtime_error& exc) {
00308     FEI_OSTRINGSTREAM osstr;
00309     osstr << "fei::MatrixGraph_Impl2::definePattern caught error: "<<exc.what();
00310     throw std::runtime_error(osstr.str());
00311   }
00312 
00313   snl_fei::RecordCollection* rec_coll = NULL;
00314   rowSpace_->getRecordCollection(idType, rec_coll);
00315 
00316   fei::Pattern* pattern = new fei::Pattern(numIDs, idType, rec_coll,
00317                         numFieldsPerID, fieldIDs, &(fieldSizes[0]));
00318   return addPattern(pattern);
00319 }
00320 
00321 //----------------------------------------------------------------------------
00322 int fei::MatrixGraph_Impl2::definePattern(int numIDs,
00323                                       const int* idTypes,
00324                                       const int* numFieldsPerID,
00325                                       const int* fieldIDs)
00326 {
00327   std::vector<int> fieldSizes;
00328   try {
00329     int offset = 0;
00330     for(int i=0; i<numIDs; ++i) {
00331       for(int j=0; j<numFieldsPerID[i]; ++j) {
00332         fieldSizes.push_back(snl_fei::getFieldSize(fieldIDs[offset++],
00333                                               rowSpace_.get(), colSpace_.get()));
00334       }
00335     }
00336   }
00337   catch (std::runtime_error& exc) {
00338     FEI_OSTRINGSTREAM osstr;
00339     osstr << "fei::MatrixGraph_Impl2::definePattern caught error: "<<exc.what();
00340     throw std::runtime_error(osstr.str());
00341   }
00342 
00343   std::vector<snl_fei::RecordCollection*> recordCollections(numIDs);
00344   for(int i=0; i<numIDs; ++i) {
00345     rowSpace_->getRecordCollection(idTypes[i], recordCollections[i]);
00346   }
00347 
00348   fei::Pattern* pattern = new fei::Pattern(numIDs, idTypes, &recordCollections[0],
00349                         numFieldsPerID, fieldIDs, &(fieldSizes[0]));
00350   return addPattern(pattern);
00351 }
00352 
00353 //------------------------------------------------------------------------------
00354 fei::Pattern* fei::MatrixGraph_Impl2::getPattern(int patternID)
00355 {
00356   std::map<int,fei::Pattern*>::iterator
00357     p_iter = patterns_.find(patternID);
00358 
00359   if (p_iter == patterns_.end()) {
00360     return NULL;
00361   }
00362 
00363   fei::Pattern* pattern = (*p_iter).second;
00364   return pattern;
00365 }
00366 
00367 //------------------------------------------------------------------------------
00368 int fei::MatrixGraph_Impl2::initConnectivityBlock(int blockID,
00369                                             int numConnectivityLists,
00370                                             int patternID,
00371                                             bool diagonal)
00372 {
00373   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00374     (*output_stream_) <<dbgprefix_<< "initConnectivityBlock symm, blkID="
00375      << blockID<<", num="<<numConnectivityLists<<", patternID="
00376      <<patternID<<FEI_ENDL;
00377   }
00378 
00379   if (blockID < 0) {
00380     fei::console_out() << "fei::MatrixGraph_Impl2::initConnectivityBlock: blockID ("
00381          << blockID << ") must be non-negative." << FEI_ENDL;
00382     ERReturn(-1);
00383   }
00384 
00385   std::map<int,fei::Pattern*>::iterator
00386     p_iter = patterns_.find(patternID);
00387 
00388   if (p_iter == patterns_.end()) {
00389     ERReturn(-1);
00390   }
00391 
00392   fei::Pattern* pattern = (*p_iter).second;
00393 
00394   if (getConnectivityBlock(blockID) != NULL) ERReturn(-1);
00395 
00396   fei::ConnectivityBlock* cblock =
00397     new fei::ConnectivityBlock(blockID, pattern, numConnectivityLists);
00398 
00399   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00400   if (diagonal) cblock->setIsDiagonal(diagonal);
00401 
00402   return(0);
00403 }
00404 
00405 //------------------------------------------------------------------------------
00406 int fei::MatrixGraph_Impl2::initConnectivityBlock(int numConnectivityLists,
00407                                                   int patternID,
00408                                                   bool diagonal)
00409 {
00410   int blockID = connectivityBlocks_.size();
00411   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00412     (*output_stream_) <<dbgprefix_<< "initConnectivityBlock symm, blkID="
00413      << blockID<<", num="<<numConnectivityLists<<", patternID="
00414      <<patternID<<FEI_ENDL;
00415   }
00416 
00417   std::map<int,fei::Pattern*>::iterator
00418     p_iter = patterns_.find(patternID);
00419 
00420   if (p_iter == patterns_.end()) {
00421     FEI_OSTRINGSTREAM osstr;
00422     osstr << "fei::MatrixGraph_Impl2::initConnectivityBlock: ERROR, patternID "
00423       << patternID << " not found.";
00424     throw std::runtime_error(osstr.str());
00425   }
00426 
00427   fei::Pattern* pattern = (*p_iter).second;
00428 
00429   fei::ConnectivityBlock* cblock =
00430     new fei::ConnectivityBlock(blockID, pattern, numConnectivityLists);
00431 
00432   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00433   if (diagonal) cblock->setIsDiagonal(diagonal);
00434 
00435   return(blockID);
00436 }
00437 
00438 //------------------------------------------------------------------------------
00439 int fei::MatrixGraph_Impl2::initConnectivityBlock(int blockID,
00440                                             int numConnectivityLists,
00441                                             int rowPatternID,
00442                                             int colPatternID)
00443 {
00444   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00445     (*output_stream_) <<dbgprefix_<< "initConnectivityBlock nonsymm, blkID="
00446      << blockID<<", num="<<numConnectivityLists<<", row-patternID="
00447      <<rowPatternID<<", col-patternID="<<colPatternID<<FEI_ENDL;
00448   }
00449 
00450   if (blockID < 0) {
00451     FEI_OSTRINGSTREAM osstr;
00452     osstr << "fei::MatrixGraph_Impl2::initConnectivityBlock: blockID ("
00453           << blockID << ") must be non-negative." << FEI_ENDL;
00454     throw std::runtime_error(osstr.str());
00455   }
00456 
00457   std::map<int,fei::Pattern*>::iterator
00458     rp_iter = patterns_.find(rowPatternID);
00459   if (rp_iter == patterns_.end()) ERReturn(-1);
00460 
00461   fei::Pattern* rowpattern = (*rp_iter).second;
00462 
00463   std::map<int,fei::Pattern*>::iterator
00464     cp_iter = patterns_.find(colPatternID);
00465   if (cp_iter == patterns_.end()) ERReturn(-1);
00466 
00467   fei::Pattern* colpattern = (*cp_iter).second;
00468 
00469   if (getConnectivityBlock(blockID) != NULL) ERReturn(-1);
00470 
00471   fei::ConnectivityBlock* cblock =
00472     new fei::ConnectivityBlock(blockID, rowpattern,
00473                                    colpattern, numConnectivityLists);
00474 
00475   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00476 
00477   return(0);
00478 }
00479 
00480 //------------------------------------------------------------------------------
00481 int fei::MatrixGraph_Impl2::initConnectivity(int blockID,
00482                                        int connectivityID,
00483                                        const int* connectedIdentifiers)
00484 {
00485   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00486     (*output_stream_) <<dbgprefix_<< "initConnectivity blkID="
00487      <<blockID<<", connID="<<connectivityID<<FEI_ENDL;
00488   }
00489 
00490   //for each connected identifier, get a pointer to its record from the
00491   //solution-space.
00492   //store those pointers in the appropriate connectivity-table, mapped to
00493   //the user-provided connectivityID.
00494 
00495   fei::ConnectivityBlock* connblk = getConnectivityBlock(blockID);
00496   if (connblk == NULL) {
00497     FEI_OSTRINGSTREAM osstr;
00498     osstr << "fei::MatrixGraph_Impl2::initConnectivity ERROR, blockID " << blockID
00499           << " doesn't correspond to an existing ConnectivityBlock.";
00500     throw std::runtime_error(osstr.str());
00501   }
00502 
00503   fei::Pattern* pattern = connblk->getRowPattern();
00504   int numIDs = pattern->getNumIDs();
00505 
00506   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
00507     for(int ii=0; ii<numIDs; ++ii) {
00508       if (isLogID(connectedIdentifiers[ii])) {
00509         FEI_OSTREAM& os = *output_stream_;
00510         os << dbgprefix_<<"initConnectivity blockID " << blockID << ", connID "
00511            << connectivityID << ": ";
00512         for(int jj=0; jj<numIDs; ++jj) {
00513           os << connectedIdentifiers[jj] << " ";
00514         }
00515         os << FEI_ENDL;
00516       }
00517     }
00518   }
00519 
00520   if (rowSpace_.get() == NULL) ERReturn(-1);
00521 
00522   std::map<int,int>& connectivityIDs = connblk->getConnectivityIDs();
00523 
00524   int idOffset = -1;
00525   std::map<int,int>::iterator
00526    iter = connectivityIDs.find(connectivityID);
00527   if (iter == connectivityIDs.end()) {
00528     idOffset = connectivityIDs.size();
00529     connectivityIDs.insert(iter, std::make_pair(connectivityID,idOffset));
00530   }
00531   else {
00532     idOffset = iter->second;
00533   }
00534 
00535   idOffset *= pattern->getNumIDs();
00536 
00537   int* rlist = &(connblk->getRowConnectivities()[idOffset]);
00538 
00539   CHK_ERR( getConnectivityRecords(pattern, rowSpace_.get(),
00540                                   connectedIdentifiers, rlist) );
00541 
00542   for(int i=0; i<numIDs; ++i) {
00543     if (pattern->getNumFieldsPerID()[i] > 0) {
00544       pattern->getRecordCollections()[i]->getRecordWithLocalID(rlist[i])->isInLocalSubdomain_ = true;
00545     }
00546   }
00547 
00548   return(0);
00549 }
00550 
00551 //------------------------------------------------------------------------------
00552 int fei::MatrixGraph_Impl2::initConnectivity(int idType,
00553                                            int numRows,
00554                                            const int* rowIDs,
00555                                            const int* rowOffsets,
00556                                            const int* packedColumnIDs)
00557 {
00558   fei::ConnectivityBlock* block = new fei::ConnectivityBlock(numRows,
00559                                                    rowIDs, rowOffsets);
00560 
00561   int max_row_len = 0;
00562   for(int i=0; i<numRows; ++i) {
00563     int row_len = rowOffsets[i+1]-rowOffsets[i];
00564     if (row_len > max_row_len) max_row_len = row_len;
00565   }
00566 
00567   int patternID = definePattern(max_row_len, idType);
00568   fei::Pattern* pattern = getPattern(patternID);
00569   block->setRowPattern(pattern);
00570 
00571   sparseBlocks_.push_back(block);
00572 
00573   int* row_records = &(block->getRowConnectivities()[0]);
00574 
00575   CHK_ERR( getConnectivityRecords(rowSpace_.get(),
00576                                   idType, numRows, rowIDs, row_records) );
00577 
00578   fei::VectorSpace* colSpace = rowSpace_.get();
00579   if (colSpace_.get() != NULL) {
00580     colSpace = colSpace_.get();
00581   }
00582 
00583   int* col_records = &(block->getColConnectivities()[0]);
00584 
00585   CHK_ERR( getConnectivityRecords(colSpace, idType, rowOffsets[numRows],
00586                                   packedColumnIDs, col_records));
00587 
00588   return(0);
00589 }
00590 
00591 //------------------------------------------------------------------------------
00592 int fei::MatrixGraph_Impl2::initConnectivity(int idType,
00593                                        int numRows,
00594                                        const int* rowIDs,
00595                                        const int* rowLengths,
00596                                        const int*const* columnIDs)
00597 {
00598   fei::ConnectivityBlock* block = new fei::ConnectivityBlock(numRows,
00599                                                    rowIDs, rowLengths, true);
00600 
00601   int max_row_len = 0;
00602   for(int i=0; i<numRows; ++i) {
00603     int row_len = rowLengths[i];
00604     if (row_len > max_row_len) max_row_len = row_len;
00605   }
00606 
00607   int patternID = definePattern(max_row_len, idType);
00608   fei::Pattern* pattern = getPattern(patternID);
00609   block->setRowPattern(pattern);
00610 
00611   sparseBlocks_.push_back(block);
00612 
00613   int* row_records = &(block->getRowConnectivities()[0]);
00614 
00615   CHK_ERR( getConnectivityRecords(rowSpace_.get(),
00616                                   idType, numRows, rowIDs, row_records) );
00617 
00618   fei::VectorSpace* colSpace = rowSpace_.get();
00619   if (colSpace_.get() != NULL) {
00620     colSpace = colSpace_.get();
00621   }
00622 
00623   int* col_records = &(block->getColConnectivities()[0]);
00624 
00625   int offset = 0;
00626   for(int i=0; i<numRows; ++i) {
00627     CHK_ERR( getConnectivityRecords(colSpace, idType, rowLengths[i],
00628                                     columnIDs[i], &(col_records[offset])));
00629     offset += rowLengths[i];
00630   }
00631 
00632   return(0);
00633 }
00634 
00635 //------------------------------------------------------------------------------
00636 int fei::MatrixGraph_Impl2::getConnectivityRecords(fei::VectorSpace* vecSpace,
00637                                                  int idType,
00638                                                  int numIDs,
00639                                                  const int* IDs,
00640                                                  int* records)
00641 {
00642   snl_fei::RecordCollection* collection = NULL;
00643   CHK_ERR( vecSpace->getRecordCollection(idType, collection) );
00644 
00645   for(int i=0; i<numIDs; ++i) {
00646     records[i] = collection->getLocalID(IDs[i]);
00647     if (records[i] == -1) {
00648       CHK_ERR( vecSpace->addDOFs(idType, 1, &IDs[i], &records[i]) );
00649     }
00650   }
00651 
00652   return(0);
00653 }
00654 
00655 //------------------------------------------------------------------------------
00656 int fei::MatrixGraph_Impl2::getConnectivityRecords(fei::VectorSpace* vecSpace,
00657                                                  int idType,
00658                                                  int fieldID,
00659                                                  int numIDs,
00660                                                  const int* IDs,
00661                                                  int* records)
00662 {
00663   snl_fei::RecordCollection* collection = NULL;
00664   CHK_ERR( vecSpace->getRecordCollection(idType, collection) );
00665 
00666   for(int i=0; i<numIDs; ++i) {
00667     records[i] = collection->getLocalID(IDs[i]);
00668     if (records[i] == -1) {
00669       CHK_ERR( vecSpace->addDOFs(fieldID, idType, 1, &IDs[i], &records[i]));
00670     }
00671   }
00672 
00673   return(0);
00674 }
00675 
00676 //------------------------------------------------------------------------------
00677 int fei::MatrixGraph_Impl2::getConnectivityRecords(fei::Pattern* pattern,
00678                                              fei::VectorSpace* vecSpace,
00679                                              const int* connectedIdentifiers,
00680                                              int* recordList)
00681 {
00682   fei::Pattern::PatternType pType = pattern->getPatternType();
00683   int i, numIDs = pattern->getNumIDs();
00684   const int* numFieldsPerID = pattern->getNumFieldsPerID();
00685   const int* fieldIDs = pattern->getFieldIDs();
00686 
00687   if (pType == fei::Pattern::GENERAL) {
00688     const int* idTypes = pattern->getIDTypes();
00689 
00690     int fieldOffset = 0;
00691     for(i=0; i<numIDs; ++i) {
00692       int id = connectedIdentifiers[i];
00693 
00694       for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
00695         CHK_ERR( vecSpace->addDOFs(fieldIDs[fieldOffset++],
00696                                                idTypes[i], 1, &id,
00697                                                &recordList[i]));
00698       }
00699     }
00700   }
00701   else if (pType == fei::Pattern::SINGLE_IDTYPE ||
00702            pType == fei::Pattern::SIMPLE ||
00703            pType == fei::Pattern::NO_FIELD) {
00704 
00705     int idType = pattern->getIDTypes()[0];
00706 
00707     switch(pType) {
00708     case fei::Pattern::SINGLE_IDTYPE:
00709       {
00710         int fieldOffset = 0;
00711         for(i=0; i<numIDs; ++i) {
00712           int id = connectedIdentifiers[i];
00713           for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
00714             CHK_ERR( vecSpace->addDOFs(fieldIDs[fieldOffset++],
00715                                                    idType, 1, &id,
00716                                                    &(recordList[i])));
00717           }
00718         }
00719       }
00720       break;
00721     case fei::Pattern::SIMPLE:
00722       {
00723         CHK_ERR( vecSpace->addDOFs(fieldIDs[0], idType,
00724                                    numIDs, connectedIdentifiers,
00725                                    recordList) );
00726       }
00727       break;
00728     case fei::Pattern::NO_FIELD:
00729       CHK_ERR( vecSpace->addDOFs(idType, numIDs,
00730                                  connectedIdentifiers, recordList));
00731       break;
00732     case fei::Pattern::GENERAL:
00733       //Include a stub for this case to avoid compiler warning...
00734       std::abort(); break;
00735     }
00736   }
00737   else {
00738     ERReturn(-1);
00739   }
00740 
00741   return(0);
00742 }
00743 
00744 //------------------------------------------------------------------------------
00745 int fei::MatrixGraph_Impl2::initConnectivity(int blockID,
00746                                        int connectivityID,
00747                                        const int* rowConnectedIdentifiers,
00748                                        const int* colConnectedIdentifiers)
00749 {
00750   //for each connected identifier, get a pointer to its record from the
00751   //solution-space.
00752   //store those pointers in the appropriate connectivity-table, mapped to
00753   //the user-provided connectivityID.
00754 
00755   fei::ConnectivityBlock* connblk = getConnectivityBlock(blockID);
00756   if (connblk == NULL) {
00757     FEI_OSTRINGSTREAM osstr;
00758     osstr << "fei::MatrixGraph_Impl2::initConnectivity ERROR, blockID " << blockID
00759           << " doesn't correspond to an existing ConnectivityBlock.";
00760     throw std::runtime_error(osstr.str());
00761   }
00762 
00763   fei::Pattern* pattern = connblk->getRowPattern();
00764   int numIDs = pattern->getNumIDs();
00765   fei::Pattern* colPattern = connblk->getColPattern();
00766   int numColIDs = colPattern->getNumIDs();
00767   if (rowSpace_.get() == NULL) {
00768     ERReturn(-1);
00769   }
00770   if (colSpace_.get() == NULL) {
00771     ERReturn(-1);
00772   }
00773 
00774   std::map<int,int>& connectivityIDs = connblk->getConnectivityIDs();
00775 
00776   int i, idOffset = -1;
00777   std::map<int,int>::iterator
00778     iter = connectivityIDs.find(connectivityID);
00779   if (iter == connectivityIDs.end()) {
00780     idOffset = connectivityIDs.size();
00781     connectivityIDs.insert(iter, std::make_pair(connectivityID,idOffset)); 
00782   }
00783   else {
00784     idOffset = iter->second;
00785   }
00786   int* row_rlist =
00787     &(connblk->getRowConnectivities()[idOffset*numIDs]);
00788   int* col_rlist =
00789     &(connblk->getColConnectivities()[idOffset*numColIDs]);
00790 
00791   CHK_ERR( getConnectivityRecords(pattern, rowSpace_.get(),
00792                                   rowConnectedIdentifiers, row_rlist) );
00793 
00794   for(i=0; i<numIDs; ++i)
00795     pattern->getRecordCollections()[i]->getRecordWithLocalID(row_rlist[i])->isInLocalSubdomain_ = true;
00796   CHK_ERR( getConnectivityRecords(colPattern, colSpace_.get(),
00797                                   colConnectedIdentifiers, col_rlist) );
00798 
00799   for(i=0; i<numColIDs; ++i)
00800     colPattern->getRecordCollections()[i]->getRecordWithLocalID(col_rlist[i])->isInLocalSubdomain_ = true;
00801 
00802   return(0);
00803 }
00804 
00805 //------------------------------------------------------------------------------
00806 int fei::MatrixGraph_Impl2::initConnectivity(int patternID,
00807                                        const int* connectedIdentifiers)
00808 {
00809   std::map<int,fei::Pattern*>::iterator
00810     p_iter = patterns_.find(patternID);
00811   if (p_iter == patterns_.end()) ERReturn(-1);
00812 
00813   fei::Pattern* pattern = p_iter->second;
00814 
00815   int blockID = -arbitraryBlockCounter_++;
00816   fei::ConnectivityBlock* cblock = new fei::ConnectivityBlock(blockID, pattern, 1);
00817 
00818   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00819 
00820   CHK_ERR( initConnectivity(blockID, 0, connectedIdentifiers) );
00821   return(0);
00822 }
00823 
00824 //------------------------------------------------------------------------------
00825 int fei::MatrixGraph_Impl2::initConnectivity(int rowPatternID,
00826                                        const int* rowConnectedIdentifiers,
00827                                        int colPatternID,
00828                                        const int* colConnectedIdentifiers)
00829 {
00830   std::map<int,fei::Pattern*>::iterator
00831     rp_iter = patterns_.find(rowPatternID);
00832   if (rp_iter == patterns_.end()) ERReturn(-1);
00833 
00834   fei::Pattern* rowpattern = (*rp_iter).second;
00835 
00836   std::map<int,fei::Pattern*>::iterator
00837     cp_iter = patterns_.find(colPatternID);
00838   if (cp_iter == patterns_.end()) ERReturn(-1);
00839 
00840   fei::Pattern* colpattern = (*cp_iter).second;
00841 
00842   int blockID = -arbitraryBlockCounter_++;
00843   fei::ConnectivityBlock* cblock = new fei::ConnectivityBlock(blockID, rowpattern,
00844                                                     colpattern, 1);
00845 
00846   connectivityBlocks_.insert(std::pair<int,fei::ConnectivityBlock*>(blockID, cblock));
00847 
00848   CHK_ERR( initConnectivity(blockID, 0,
00849                             rowConnectedIdentifiers,
00850                             colConnectedIdentifiers) );
00851   return(0);
00852 }
00853 
00854 //------------------------------------------------------------------------------
00855 int fei::MatrixGraph_Impl2::getPatternNumIndices(int patternID,
00856                                            int& numIndices)
00857 {
00858   std::map<int,fei::Pattern*>::iterator
00859     p_iter = patterns_.find(patternID);
00860   if (p_iter == patterns_.end()) ERReturn(-1);
00861   fei::Pattern* pattern = (*p_iter).second;
00862 
00863   numIndices = pattern->getNumIndices();
00864 
00865   return(0);
00866 }
00867 
00868 //------------------------------------------------------------------------------
00869 int fei::MatrixGraph_Impl2::getPatternIndices(int patternID,
00870                                         const int* IDs,
00871                                         std::vector<int>& indices)
00872 {
00873   std::map<int,fei::Pattern*>::iterator
00874     p_iter = patterns_.find(patternID);
00875   if (p_iter == patterns_.end()) ERReturn(-1);
00876 
00877   fei::Pattern* pattern = (*p_iter).second;
00878 
00879   indices.resize(pattern->getNumIndices());
00880 
00881   int numIDs                = pattern->getNumIDs();
00882   const int* idTypes        = pattern->getIDTypes();
00883   const int* numFieldsPerID = pattern->getNumFieldsPerID();
00884   const int* fieldIDs       = pattern->getFieldIDs();
00885 
00886   int offset = 0, fieldOffset = 0;
00887   for(int i=0; i<numIDs; ++i) {
00888     for(int j=0; j<numFieldsPerID[i]; ++j) {
00889       CHK_ERR( rowSpace_->getGlobalIndices(1, &(IDs[i]), idTypes[i],
00890                                                   fieldIDs[fieldOffset],
00891                                                   &(indices[offset])) );
00892       offset += snl_fei::getFieldSize(fieldIDs[fieldOffset++],
00893                                         rowSpace_.get(), colSpace_.get());
00894     }
00895   }
00896 
00897   return(0);
00898 }
00899 
00900 //------------------------------------------------------------------------------
00901 int fei::MatrixGraph_Impl2::initConnectivity(int idType,
00902                                        int fieldID,
00903                                        int numRows,
00904                                        const int* rowIDs,
00905                                        const int* rowOffsets,
00906                                        const int* packedColumnIDs)
00907 {
00908   fei::ConnectivityBlock* block = new fei::ConnectivityBlock(fieldID, numRows,
00909                                                    rowIDs, rowOffsets);
00910 
00911   sparseBlocks_.push_back(block);
00912 
00913   int* row_records = &(block->getRowConnectivities()[0]);
00914 
00915   CHK_ERR( getConnectivityRecords(rowSpace_.get(), idType, fieldID,
00916                                   numRows, rowIDs, row_records) );
00917 
00918   fei::VectorSpace* colSpace = rowSpace_.get();
00919   if (colSpace_.get() != NULL) {
00920     colSpace = colSpace_.get();
00921   }
00922 
00923   int* col_records = &(block->getColConnectivities()[0]);
00924 
00925   CHK_ERR( getConnectivityRecords(colSpace, idType, fieldID, rowOffsets[numRows],
00926                                   packedColumnIDs, col_records));
00927 
00928   return(0);
00929 }
00930 
00931 //------------------------------------------------------------------------------
00932 int fei::MatrixGraph_Impl2::initLagrangeConstraint(int constraintID,
00933                                                int constraintIDType,
00934                                                int numIDs,
00935                                                const int* idTypes,
00936                                                const int* IDs,
00937                                                const int* fieldIDs)
00938 {
00939   if (rowSpace_.get() == NULL) ERReturn(-1);
00940 
00941   ConstraintType* constraint = getLagrangeConstraint(constraintID);
00942 
00943   CHK_ERR( rowSpace_->addDOFs(constraintIDType, 1, &constraintID) );
00944 
00945   if (haveColSpace_) {
00946     if (colSpace_.get() == NULL) {
00947       ERReturn(-1);
00948     }
00949     CHK_ERR( colSpace_->addDOFs(constraintIDType,
00950                                             1, &constraintID) );
00951   }
00952 
00953   ConstraintType* newconstraint =
00954     new ConstraintType(constraintID, constraintIDType,
00955                       false, //isSlave
00956                       false,  //isPenalty
00957                       numIDs, idTypes, IDs, fieldIDs,
00958                       0, //slaveOffset
00959                       0, //offsetIntoSlaveField
00960                       NULL, //weights
00961                       0.0, //rhsValue
00962                       rowSpace_.get());
00963  
00964   if (constraint != NULL) {
00965     if (*constraint != *newconstraint) {
00966       delete constraint;
00967     }
00968     else {
00969       delete newconstraint;
00970       return(0);
00971     }
00972   }
00973 
00974   lagrangeConstraints_[constraintID] = newconstraint;
00975 
00976   return(0);
00977 }
00978 
00979 //------------------------------------------------------------------------------
00980 int fei::MatrixGraph_Impl2::initPenaltyConstraint(int constraintID,
00981                                             int constraintIDType,
00982                                             int numIDs,
00983                                             const int* idTypes,
00984                                             const int* IDs,
00985                                             const int* fieldIDs)
00986 {
00987   if (rowSpace_.get() == NULL) ERReturn(-1);
00988 
00989   ConstraintType* constraint = getPenaltyConstraint(constraintID);
00990 
00991   ConstraintType* newconstraint =
00992     new ConstraintType(constraintID, constraintIDType,
00993                                               false, //isSlave
00994                                               true,  //isPenalty
00995                                               numIDs, idTypes, IDs, fieldIDs,
00996                                               0, //slaveOffset
00997                                               0, //offsetIntoSlaveField
00998                                               NULL, //weights
00999                                               0.0, //rhsValue
01000                                               rowSpace_.get());
01001 
01002   if (constraint != NULL) {
01003     if (*constraint != *newconstraint) {
01004       delete constraint;
01005     }
01006     else {
01007       delete newconstraint;
01008       return(0);
01009     }
01010   }
01011 
01012   penaltyConstraints_[constraintID] = newconstraint;
01013 
01014   return(0);
01015 }
01016 
01017 //------------------------------------------------------------------------------
01018 bool fei::MatrixGraph_Impl2::hasSlaveDof(int ID, int idType)
01019 {
01020   snl_fei::RecordCollection* collection = NULL;
01021   rowSpace_->getRecordCollection(idType, collection);
01022   if (collection == NULL) {
01023     throw std::runtime_error("fei::MatrixGraph_Impl2::hasSlaveDof: ERROR, unknown idType");
01024   }
01025 
01026   fei::Record<int>* rec = collection->getRecordWithID(ID);
01027 
01028   if (rec == NULL) {
01029     FEI_OSTRINGSTREAM osstr;
01030     osstr << "fei::MatrixGraph_Impl2::hasSlaveDof: ERROR, specified ID ("
01031      << ID << ") not found.";
01032     throw std::runtime_error(osstr.str());
01033   }
01034 
01035   return( rec->hasSlaveDof() );
01036 }
01037 
01038 //------------------------------------------------------------------------------
01039 int fei::MatrixGraph_Impl2::initSlaveConstraint(int numIDs,
01040                                           const int* idTypes,
01041                                           const int* IDs,
01042                                           const int* fieldIDs,
01043                                           int offsetOfSlave,
01044                                           int offsetIntoSlaveField,
01045                                           const double* weights,
01046                                           double rhsValue)
01047 {
01048   if (rowSpace_.get() == NULL) ERReturn(-1);
01049 
01050   FEI_OSTRINGSTREAM idosstr;
01051   idosstr << IDs[offsetOfSlave] << fieldIDs[offsetOfSlave] << offsetIntoSlaveField;
01052   FEI_ISTRINGSTREAM idisstr(idosstr.str());
01053   int crID = IDs[offsetOfSlave];
01054   idisstr >> crID;
01055 
01056   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
01057     FEI_OSTREAM& dbgos = *output_stream_;
01058 
01059     dbgos << dbgprefix_<<"initSlaveConstraint crID=" << crID << ", slaveID="
01060        << IDs[offsetOfSlave];
01061 
01062     if (output_level_ >= fei::FULL_LOGS) {
01063       dbgos << " { ";
01064       for(int ii=0; ii<numIDs; ++ii) {
01065         dbgos << "("<<IDs[ii] << ","<<fieldIDs[ii]<<") ";
01066       }
01067       dbgos << "}"<<FEI_ENDL;
01068     }
01069     else dbgos << FEI_ENDL;
01070   }
01071 
01072   ConstraintType* constraint = getSlaveConstraint(crID);
01073 
01074   ConstraintType* newconstraint =
01075     new ConstraintType(crID, 0,
01076                        true, //isSlave
01077                        false,  //isPenalty
01078                        numIDs, idTypes, IDs, fieldIDs,
01079                        offsetOfSlave,
01080                        offsetIntoSlaveField,
01081                        weights,
01082                        rhsValue,
01083                        rowSpace_.get());
01084  
01085   if (constraint != NULL) {
01086     if (*constraint != *newconstraint) {
01087       if (!constraint->structurallySame(*newconstraint)) {
01088         FEI_OSTRINGSTREAM osstr;
01089         osstr << "fei::MatrixGraph_Impl2::initSlaveConstraint: slave ID "<<IDs[offsetOfSlave]
01090               << " is already constrained, with different connectivity. Changing the"
01091               << " the structure of an existing constraint is not allowed.";
01092         throw std::runtime_error(osstr.str());
01093       }
01094       newSlaveData_ = true;
01095       delete constraint;
01096     }
01097     else {
01098       delete newconstraint;
01099       return(0);
01100     }
01101   }
01102   else {
01103     newSlaveData_ = true;
01104   }
01105 
01106   slaveConstraints_[crID] = newconstraint;
01107 
01108   return(0);
01109 }
01110 
01111 //------------------------------------------------------------------------------
01112 bool
01113 fei::MatrixGraph_Impl2::newSlaveData()
01114 {
01115   bool globalBool = false;
01116   fei::Allreduce(comm_, newSlaveData_, globalBool);
01117   newSlaveData_ = globalBool;
01118 
01119   return(newSlaveData_);
01120 }
01121 
01122 //------------------------------------------------------------------------------
01123 fei::ConstraintType*
01124 fei::MatrixGraph_Impl2::getLagrangeConstraint(int constraintID)
01125 {
01126   std::map<int,ConstraintType*>::iterator
01127     c_iter = lagrangeConstraints_.find(constraintID);
01128   if (c_iter == lagrangeConstraints_.end()) return(NULL);
01129 
01130   return( (*c_iter).second );
01131 }
01132 
01133 //------------------------------------------------------------------------------
01134 fei::ConstraintType*
01135 fei::MatrixGraph_Impl2::getSlaveConstraint(int constraintID)
01136 {
01137   std::map<int,ConstraintType*>::iterator
01138     c_iter = slaveConstraints_.find(constraintID);
01139   if (c_iter == slaveConstraints_.end()) return(NULL);
01140 
01141   return( (*c_iter).second );
01142 }
01143 
01144 //------------------------------------------------------------------------------
01145 fei::ConstraintType*
01146 fei::MatrixGraph_Impl2::getPenaltyConstraint(int constraintID)
01147 {
01148   std::map<int,ConstraintType*>::iterator
01149     c_iter = penaltyConstraints_.find(constraintID);
01150   if (c_iter == penaltyConstraints_.end()) return(NULL);
01151 
01152   return( (*c_iter).second );
01153 }
01154 
01155 //------------------------------------------------------------------------------
01156 int fei::MatrixGraph_Impl2::initComplete()
01157 {
01158   if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
01159     (*output_stream_) <<dbgprefix_<< "initComplete" << FEI_ENDL;
01160   }
01161   if (rowSpace_.get() == NULL) {
01162     ERReturn(-1);
01163   }
01164   else {
01165     CHK_ERR( rowSpace_->initComplete() );
01166   }
01167 
01168   if (haveColSpace_ && colSpace_.get() == NULL) ERReturn(-1);
01169   if (haveColSpace_) {
01170     CHK_ERR( colSpace_->initComplete() );
01171   }
01172 
01173   if (rowSpace_->fieldMasks_.size() == 1 &&
01174       rowSpace_->fieldMasks_[0]->getNumFields() == 1) {
01175     simpleProblem_ = true;
01176   }
01177 
01178   std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
01179   vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
01180 
01181   //If there are slave constraints (on any processor), we need to create
01182   //a slave dependency matrix.
01183   localNumSlaves_ = slaveConstraints_.size();
01184   globalNumSlaves_ = 0;
01185   CHK_ERR( fei::GlobalSum(comm_, localNumSlaves_, globalNumSlaves_) );
01186 
01187   if (globalNumSlaves_ > 0) {
01188     //we need to allocate the slave dependency and reduction matrices too.
01189     CHK_ERR( createSlaveMatrices() );
01190   }
01191 
01192   return(0);
01193 }
01194 
01195 //----------------------------------------------------------------------------
01196 int fei::MatrixGraph_Impl2::createAlgebraicGraph(bool blockEntryGraph,
01197                                                  fei::Graph* graph,
01198                                                  bool gatherFromOverlap)
01199 {
01200   //Now run the connectivity blocks, adding vertices (locations of nonzeros) to
01201   //the graph.
01202 
01203   bool wasAlreadyBlockEntry = blockEntryGraph_;
01204 
01205   blockEntryGraph_ = blockEntryGraph;
01206 
01207   for(unsigned i=0; i<sparseBlocks_.size(); ++i) {
01208     fei::ConnectivityBlock& cblock = *(sparseBlocks_[i]);
01209     CHK_ERR( addBlockToGraph_sparse(graph, &cblock) );
01210   }
01211 
01212   std::map<int,fei::ConnectivityBlock*>::const_iterator
01213     cb_iter = connectivityBlocks_.begin(),
01214     cb_end  = connectivityBlocks_.end();
01215 
01216   while(cb_iter != cb_end) {
01217     fei::ConnectivityBlock& cblock = *((*cb_iter).second);
01218 
01219     bool symmetricBlock = cblock.isSymmetric();
01220 
01221     if (symmetricBlock) {
01222       fei::Pattern* pattern = cblock.getRowPattern();
01223       if (pattern == NULL) {
01224         ERReturn(-1);
01225       }
01226 
01227       fei::Pattern::PatternType pType = pattern->getPatternType();
01228       if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
01229         CHK_ERR( addBlockToGraph_multiField_symmetric(graph, &cblock) );
01230       }
01231       else if (pType == fei::Pattern::SIMPLE) {
01232         CHK_ERR( addBlockToGraph_singleField_symmetric(graph, &cblock) );
01233       }
01234       else if (pType == fei::Pattern::NO_FIELD) {
01235         CHK_ERR( addBlockToGraph_noField_symmetric(graph, &cblock) );
01236       }
01237     }
01238     else {
01239       fei::Pattern* pattern = cblock.getRowPattern();
01240       fei::Pattern* colpattern = cblock.getColPattern();
01241       if (pattern == NULL || colpattern == NULL) {
01242         ERReturn(-1);
01243       }
01244 
01245       fei::Pattern::PatternType pType = pattern->getPatternType();
01246       fei::Pattern::PatternType colPType = colpattern->getPatternType();
01247       if (pType == fei::Pattern::SIMPLE && colPType == fei::Pattern::SIMPLE) {
01248         CHK_ERR( addBlockToGraph_singleField_nonsymmetric(graph, &cblock) );
01249       }
01250       else {
01251         CHK_ERR( addBlockToGraph_multiField_nonsymmetric(graph, &cblock) );
01252       }
01253     }
01254     ++cb_iter;
01255   }
01256 
01257   CHK_ERR( addLagrangeConstraintsToGraph(graph) );
01258 
01259   CHK_ERR( addPenaltyConstraintsToGraph(graph) );
01260 
01261   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01262     FEI_OSTREAM& os = *output_stream_;
01263     os << dbgprefix_<<"# before graph->gatherFromOverlap()" << FEI_ENDL;
01264     CHK_ERR( graph->writeLocalGraph(os) );
01265     CHK_ERR( graph->writeRemoteGraph(os) );
01266   }
01267 
01268   if (gatherFromOverlap) {
01269     CHK_ERR( graph->gatherFromOverlap() );
01270   }
01271 
01272   if (blockEntryGraph) {
01273     CHK_ERR( exchangeBlkEqnSizes(graph) );
01274   }
01275 
01276   if (output_level_ >= fei::FULL_LOGS && fei::numProcs(comm_)>1 &&
01277       output_stream_ != NULL && gatherFromOverlap) {
01278     FEI_OSTREAM& os = *output_stream_;
01279     os << dbgprefix_<<" after graph->gatherFromOverlap()" << FEI_ENDL;
01280     CHK_ERR( graph->writeLocalGraph(*output_stream_) );
01281   }
01282 
01283   if (!wasAlreadyBlockEntry) {
01284     setIndicesMode(POINT_ENTRY_GRAPH);
01285   }
01286 
01287   return(0);
01288 }
01289 
01290 //----------------------------------------------------------------------------
01291 fei::SharedPtr<fei::SparseRowGraph>
01292 fei::MatrixGraph_Impl2::createGraph(bool blockEntryGraph,
01293                                     bool localRowGraph_includeSharedRows)
01294 {
01295   fei::SharedPtr<fei::SparseRowGraph> localRows;
01296 
01297   std::vector<int> globalOffsets;
01298 
01299   if (blockEntryGraph) {
01300     if (reducer_.get() != NULL) {
01301       throw std::runtime_error("fei::MatrixGraph_Impl2::createGraph ERROR, can't specify both block-entry assembly and slave-constraint reduction.");
01302     }
01303 
01304     rowSpace_->getGlobalBlkIndexOffsets(globalOffsets);
01305   }
01306   else {
01307     rowSpace_->getGlobalIndexOffsets(globalOffsets);
01308   }
01309 
01310   if ((int)globalOffsets.size() < localProc_+2) return localRows;
01311 
01312   int firstOffset = globalOffsets[localProc_];
01313   int lastOffset = globalOffsets[localProc_+1] - 1;
01314 
01315   if (reducer_.get() != NULL) {
01316     std::vector<int>& reduced_eqns = reducer_->getLocalReducedEqns();
01317     firstOffset = reduced_eqns[0];
01318     lastOffset = reduced_eqns[reduced_eqns.size()-1];
01319   }
01320 
01321   fei::SharedPtr<fei::Graph> inner_graph(new fei::Graph_Impl(comm_, firstOffset, lastOffset) );
01322   fei::SharedPtr<fei::Graph> graph;
01323 
01324   if (reducer_.get() == NULL) {
01325     graph = inner_graph;
01326   }
01327   else {
01328     graph.reset( new fei::GraphReducer(reducer_, inner_graph) );
01329   }
01330 
01331   bool gatherFromOverlap = !localRowGraph_includeSharedRows;
01332   int err = createAlgebraicGraph(blockEntryGraph, graph.get(),
01333                                  gatherFromOverlap);
01334   if (err != 0) {
01335     return(localRows);
01336   }
01337 
01338   localRows = fei::createSparseRowGraph(*(inner_graph->getLocalGraph()));
01339 
01340   remotelyOwnedGraphRows_ = fei::createSparseRowGraph(inner_graph->getRemoteGraph());
01341 
01342   remotelyOwnedGraphRows_->blockEntries = blockEntryGraph;
01343 
01344   if (localRowGraph_includeSharedRows &&
01345       remotelyOwnedGraphRows_->rowNumbers.size() > 0) {
01346 
01347     fei::SharedPtr<fei::SparseRowGraph> localPlusShared =
01348       snl_fei::mergeSparseRowGraphs(localRows.get(), remotelyOwnedGraphRows_.get());
01349     localRows = localPlusShared;
01350   }
01351 
01352   return(localRows);
01353 }
01354 
01355 //----------------------------------------------------------------------------
01356 int fei::MatrixGraph_Impl2::exchangeBlkEqnSizes(fei::Graph* graph)
01357 {
01358   //If point-equals-block (meaning block-equations are of size 1) or
01359   //if blockEntryGraph_ is false (meaning we haven't constructed a block-entry-
01360   //graph) then there is nothing to do in this function.
01361   if ( rowSpace_->getPointBlockMap()->ptEqualBlk() ) {
01362     return(0);
01363   }
01364 
01365   snl_fei::BlkSizeMsgHandler blkHandler(rowSpace_.get(), graph, comm_);
01366 
01367   CHK_ERR( blkHandler.do_the_exchange() );
01368 
01369   return(0);
01370 }
01371 
01372 //----------------------------------------------------------------------------
01373 int fei::MatrixGraph_Impl2::createSlaveMatrices()
01374 {
01375   //In this function we extract the dependency matrix of equation numbers 
01376   //from the slave-constraints locally on each processor, then gather those
01377   //slave-equations onto each processor (so that each processor holds ALL
01378   //slave-equations).
01379   //Then we remove any couplings that may exist among the slave-equations and
01380   //finally create a FillableMat object to hold the final dependency matrix D_.
01381   //
01382 
01383   if (!newSlaveData()) {
01384     return(0);
01385   }
01386 
01387   std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
01388   vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
01389 
01390   fei::SharedPtr<fei::FillableMat> local_D(new fei::FillableMat);
01391   fei::SharedPtr<fei::CSVec> local_g(new fei::CSVec);
01392 
01393   std::vector<int> masterEqns;
01394   std::vector<double> masterCoefs;
01395 
01396   std::map<int,ConstraintType*>::const_iterator
01397     cr_iter = slaveConstraints_.begin(),
01398     cr_end  = slaveConstraints_.end();
01399 
01400   for(; cr_iter != cr_end; ++cr_iter) {
01401     ConstraintType* cr = (*cr_iter).second;
01402 
01403     fei::Record<int>* slaveRecord = cr->getSlave();
01404     int slaveIDType = cr->getIDType();
01405     int slaveFieldID = cr->getSlaveFieldID();
01406     int offsetIntoSlaveField = cr->getOffsetIntoSlaveField();
01407     int slaveEqn = -1;
01408     CHK_ERR( rowSpace_->getGlobalIndex(slaveIDType, slaveRecord->getID(),
01409                                        slaveFieldID, 0, offsetIntoSlaveField,
01410                                        slaveEqn) );
01411 
01412     std::vector<int>& masterRecords_vec = cr->getMasters();
01413     int* masterRecords = &masterRecords_vec[0];
01414     std::vector<snl_fei::RecordCollection*>& masterRecColls = cr->getMasterRecordCollections();
01415     std::vector<int>& masterIDTypes = cr->getMasterIDTypes();
01416     std::vector<int>& masterFieldIDs = cr->getMasterFieldIDs();
01417     std::vector<double>& masterWeights = cr->getMasterWeights();
01418     double* masterWtPtr = &masterWeights[0];
01419 
01420     masterEqns.resize(masterWeights.size());
01421     masterCoefs.resize(masterWeights.size());
01422 
01423     int* masterEqnsPtr = &(masterEqns[0]);
01424     double* masterCoefsPtr = &(masterCoefs[0]);
01425 
01426     int offset = 0;
01427     for(size_t j=0; j<masterIDTypes.size(); ++j) {
01428       fei::Record<int>* masterRecord = masterRecColls[j]->getRecordWithLocalID(masterRecords[j]);
01429       int* eqnNumbers = vspcEqnPtr_+masterRecord->getOffsetIntoEqnNumbers();
01430       fei::FieldMask* mask = masterRecord->getFieldMask();
01431       int eqnOffset = 0;
01432       if (!simpleProblem_) {
01433         mask->getFieldEqnOffset(masterFieldIDs[j], eqnOffset);
01434       }
01435 
01436       unsigned fieldSize = rowSpace_->getFieldSize(masterFieldIDs[j]);
01437 
01438       int eqn = eqnNumbers[eqnOffset];
01439       for(unsigned k=0; k<fieldSize; ++k) {
01440         masterEqnsPtr[offset++] = eqn+k;
01441       }
01442     }
01443 
01444     double fei_eps = 1.e-49;
01445 
01446     offset = 0;
01447     for(size_t jj=0; jj<masterEqns.size(); ++jj) {
01448       if (std::abs(masterWtPtr[jj]) > fei_eps) {
01449         masterCoefsPtr[offset] = masterWtPtr[jj];
01450         masterEqnsPtr[offset++] = masterEqnsPtr[jj];
01451       }
01452     }
01453 
01454     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
01455       bool log_eqn = isLogEqn(slaveEqn);
01456       if (!log_eqn) {
01457         for(unsigned ii=0; ii<masterEqns.size(); ++ii) {
01458           if (isLogEqn(masterEqnsPtr[ii])) {
01459             log_eqn = true;
01460             break;
01461           }
01462         }
01463       }
01464 
01465       if (log_eqn) {
01466         FEI_OSTREAM& os = *output_stream_;
01467         os << "createSlaveMatrices: " << slaveEqn << " = ";
01468         for(unsigned ii=0; ii<masterEqns.size(); ++ii) {
01469           if (ii!=0) os << "+ ";
01470           os << masterCoefsPtr[ii]<<"*"<<masterEqnsPtr[ii]<<" ";
01471         }
01472         os << FEI_ENDL;
01473       }
01474     }
01475 
01476     local_D->putRow(slaveEqn, masterEqnsPtr, masterCoefsPtr, offset);
01477 
01478     if (std::abs(cr->getRHSValue()) > fei_eps) {
01479       fei::put_entry(*local_g, slaveEqn, cr->getRHSValue());
01480     }
01481   }
01482 
01483   if (D_.get() == NULL) {
01484     D_.reset(new fei::FillableMat);
01485   }
01486 
01487   if (g_.get() == NULL) {
01488     g_.reset(new fei::CSVec);
01489   }
01490 
01491 #ifndef FEI_SER
01492   if (numProcs_ > 1) {
01493     fei::impl_utils::global_union(comm_, *local_D, *D_);
01494     fei::impl_utils::global_union(comm_, *local_g, *g_);
01495   }
01496   else {
01497     *D_ = *local_D;
01498     *g_ = *local_g;
01499   }
01500 #else
01501   *D_ = *local_D;
01502   *g_ = *local_g;
01503 #endif
01504 
01505   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01506     (*output_stream_) << "#  D_ (pre-removeCouplings):"<<FEI_ENDL;
01507     (*output_stream_) << *D_;
01508   }
01509 
01510   int levelsOfCoupling = fei::impl_utils::remove_couplings(*D_);
01511   (void)levelsOfCoupling;
01512 
01513   if (reducer_.get() == NULL) {
01514     reducer_.reset(new fei::Reducer(D_, g_, comm_));
01515 
01516     std::vector<int> indices;
01517     rowSpace_->getIndices_Owned(indices);
01518 
01519     reducer_->setLocalUnreducedEqns(indices);
01520   }
01521   else {
01522     reducer_->initialize();
01523   }
01524 
01525   if (output_level_ >= fei::FULL_LOGS && output_stream_ != NULL) {
01526     (*output_stream_) << "#  D_:"<<FEI_ENDL;
01527     (*output_stream_) << *D_;
01528   }
01529 
01530   double fei_eps = 1.e-49;
01531 
01532   g_nonzero_ = false;
01533   for(size_t j=0; j<g_->size(); ++j) {
01534     double coef = g_->coefs()[j];
01535     if (std::abs(coef) > fei_eps) {
01536       g_nonzero_ = true;
01537     }
01538   }
01539 
01540   newSlaveData_ = false;
01541 
01542   return(0);
01543 }
01544 
01545 //----------------------------------------------------------------------------
01546 fei::SharedPtr<fei::Reducer>
01547 fei::MatrixGraph_Impl2::getReducer()
01548 {
01549   return( reducer_ );
01550 }
01551 
01552 //----------------------------------------------------------------------------
01553 fei::SharedPtr<fei::SparseRowGraph>
01554 fei::MatrixGraph_Impl2::getRemotelyOwnedGraphRows()
01555 {
01556   return( remotelyOwnedGraphRows_ );
01557 }
01558 
01559 //----------------------------------------------------------------------------
01560 void fei::MatrixGraph_Impl2::getConstrainedIndices(std::vector<int>& crindices) const
01561 {
01562   if (constrained_indices_.empty()) {
01563     crindices.clear();
01564     return;
01565   }
01566 
01567   std::set<int>::const_iterator
01568     s_iter = constrained_indices_.begin(),
01569     s_end = constrained_indices_.end();
01570 
01571   crindices.resize(constrained_indices_.size());
01572 
01573   int offset = 0;
01574   for(; s_iter != s_end; ++s_iter) {
01575     crindices[offset++] = *s_iter;
01576   }
01577 }
01578 
01579 //----------------------------------------------------------------------------
01580 int fei::MatrixGraph_Impl2::addLagrangeConstraintsToGraph(fei::Graph* graph)
01581 {
01582   std::vector<int> indices;
01583   std::map<int,ConstraintType*>::const_iterator
01584     cr_iter = lagrangeConstraints_.begin(),
01585     cr_end  = lagrangeConstraints_.end();
01586 
01587   constrained_indices_.clear();
01588 
01589   while(cr_iter != cr_end) {
01590     ConstraintType* cr = (*cr_iter).second;
01591     int crID = cr->getConstraintID();
01592 
01593     CHK_ERR( getConstraintConnectivityIndices(cr, indices) );
01594 
01595     int numIndices = indices.size();
01596     int* indicesPtr = &(indices[0]);
01597 
01598     for(int i=0; i<numIndices; ++i) {
01599       constrained_indices_.insert(indicesPtr[i]);
01600     }
01601 
01602     int crEqnRow = -1, tmp;
01603     if (blockEntryGraph_) {
01604       CHK_ERR( rowSpace_->getGlobalBlkIndex(cr->getIDType(),
01605                                                    crID, crEqnRow) );
01606       cr->setBlkEqnNumber(crEqnRow);
01607       CHK_ERR( rowSpace_->getGlobalIndex(cr->getIDType(),
01608                                                 crID, tmp) );
01609       cr->setEqnNumber(tmp);
01610     }
01611     else {
01612       CHK_ERR( rowSpace_->getGlobalIndex(cr->getIDType(),
01613                                                 crID, crEqnRow) );
01614       cr->setEqnNumber(crEqnRow);
01615       CHK_ERR( rowSpace_->getGlobalBlkIndex(cr->getIDType(),
01616                                                    crID, tmp) );
01617       cr->setBlkEqnNumber(tmp);
01618     }
01619 
01620     //now add the row contribution
01621     CHK_ERR( graph->addIndices(crEqnRow, numIndices, &(indices[0])) );
01622 
01623     //Let's add a diagonal entry to the graph for this constraint-equation,
01624     //just in case we need to fiddle with this equation later (e.g. discard the
01625     //constraint equation and replace it with a dirichlet boundary condition...).
01626     CHK_ERR( graph->addIndices(crEqnRow, 1, &crEqnRow) );
01627 
01628     //and finally, add the column contribution (which is simply the transpose
01629     //of the row contribution)
01630     for(int k=0; k<numIndices; ++k) {
01631       CHK_ERR( graph->addIndices(indicesPtr[k], 1, &crEqnRow) );
01632     }
01633 
01634     ++cr_iter;
01635   }
01636 
01637   return(0);
01638 }
01639 
01640 //----------------------------------------------------------------------------
01641 int fei::MatrixGraph_Impl2::
01642 getConstraintConnectivityIndices(ConstraintType* cr,
01643                                  std::vector<int>& globalIndices)
01644 {
01645   std::vector<int>& fieldSizes = tmpIntArray1_;
01646   std::vector<int>& ones = tmpIntArray2_;
01647 
01648   std::vector<int>& constrainedRecords = cr->getMasters();
01649   std::vector<int>& constrainedFieldIDs = cr->getMasterFieldIDs();
01650   std::vector<snl_fei::RecordCollection*>& recordCollections = cr->getMasterRecordCollections();
01651 
01652   int len = constrainedRecords.size();
01653   fieldSizes.resize(len);
01654 
01655   ones.assign(len, 1);
01656 
01657   int numIndices = 0;
01658   if (blockEntryGraph_) {
01659     numIndices = len;
01660   }
01661   else {
01662     for(int j=0; j<len; ++j) {
01663       unsigned fieldSize = rowSpace_->getFieldSize(constrainedFieldIDs[j]);
01664       fieldSizes[j] = fieldSize;
01665       numIndices += fieldSize;
01666     }
01667   }
01668 
01669   globalIndices.resize(numIndices);
01670 
01671   int checkNum;
01672   CHK_ERR( getConnectivityIndices_multiField(&recordCollections[0],
01673                                              &constrainedRecords[0],
01674                                              len, &ones[0],
01675                                              &constrainedFieldIDs[0],
01676                                              &fieldSizes[0],
01677                                              numIndices, &globalIndices[0],
01678                                              checkNum) );
01679   if (numIndices != checkNum) {
01680     ERReturn(-1);
01681   }
01682 
01683   return(0);
01684 }
01685 
01686 //----------------------------------------------------------------------------
01687 int fei::MatrixGraph_Impl2::addPenaltyConstraintsToGraph(fei::Graph* graph)
01688 {
01689   std::vector<int> indices;
01690   std::map<int,ConstraintType*>::const_iterator
01691     cr_iter = penaltyConstraints_.begin(),
01692     cr_end  = penaltyConstraints_.end();
01693 
01694   while(cr_iter != cr_end) {
01695     ConstraintType* cr = (*cr_iter).second;
01696 
01697     CHK_ERR( getConstraintConnectivityIndices(cr, indices) );
01698 
01699     int numIndices = indices.size();
01700 
01701     //now add the symmetric contributions to the graph
01702     CHK_ERR( graph->addSymmetricIndices(numIndices, &(indices[0])) );
01703 
01704     ++cr_iter;
01705   }
01706 
01707   return(0);
01708 }
01709 
01710 //----------------------------------------------------------------------------
01711 int fei::MatrixGraph_Impl2::compareStructure(const fei::MatrixGraph& matrixGraph,
01712                                        bool& equivalent) const
01713 {
01714   equivalent = false;
01715   int localCode = 1;
01716   int globalCode = 0;
01717 
01718   int myNumBlocks = getNumConnectivityBlocks();
01719   int numBlocks = matrixGraph.getNumConnectivityBlocks();
01720 
01721   if (myNumBlocks != numBlocks) {
01722     CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01723     equivalent = globalCode > 0 ? false : true;
01724     return(0);
01725   }
01726 
01727   if (numBlocks > 0) {
01728     std::vector<int> myBlockIDs;
01729     std::vector<int> blockIDs;
01730 
01731     CHK_ERR( getConnectivityBlockIDs(myBlockIDs) );
01732     CHK_ERR( matrixGraph.getConnectivityBlockIDs(blockIDs) );
01733 
01734     for(int i=0; i<numBlocks; ++i) {
01735       if (myBlockIDs[i] != blockIDs[i]) {
01736         CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01737         equivalent = globalCode > 0 ? false : true;
01738         return(0);
01739       }
01740 
01741       const fei::ConnectivityBlock* mycblock = getConnectivityBlock(myBlockIDs[i]);
01742       const fei::ConnectivityBlock* cblock = matrixGraph.getConnectivityBlock(blockIDs[i]);
01743 
01744       int myNumLists = mycblock->getConnectivityIDs().size();
01745       int numLists = cblock->getConnectivityIDs().size();
01746 
01747       if (myNumLists != numLists ||
01748           mycblock->isSymmetric() != cblock->isSymmetric()) {
01749         CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01750         equivalent = globalCode > 0 ? false : true;
01751         return(0);
01752       }
01753 
01754       int myNumIDsPerList = mycblock->getRowPattern()->getNumIDs();
01755       int numIDsPerList = cblock->getRowPattern()->getNumIDs();
01756 
01757       if (myNumIDsPerList != numIDsPerList) {
01758         CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01759         equivalent = globalCode > 0 ? false : true;
01760         return(0);
01761       }
01762     }
01763   }
01764 
01765   int numMyLagrangeConstraints = getLocalNumLagrangeConstraints();
01766   int numMySlaveConstraints = getGlobalNumSlaveConstraints();
01767   int numLagrangeConstraints = matrixGraph.getLocalNumLagrangeConstraints();
01768   int numSlaveConstraints = matrixGraph.getGlobalNumSlaveConstraints();
01769 
01770   if (numMyLagrangeConstraints != numLagrangeConstraints ||
01771       numMySlaveConstraints != numSlaveConstraints) {
01772     CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01773     equivalent = globalCode > 0 ? false : true;
01774     return(0);
01775   }
01776 
01777   localCode = 0;
01778   CHK_ERR( fei::GlobalMax(comm_, localCode, globalCode) );
01779   equivalent = globalCode > 0 ? false : true;
01780   return(0);
01781 }
01782 
01783 //----------------------------------------------------------------------------
01784 int fei::MatrixGraph_Impl2::getNumConnectivityBlocks() const
01785 {
01786   return(connectivityBlocks_.size());
01787 }
01788 
01789 //----------------------------------------------------------------------------
01790 int fei::MatrixGraph_Impl2::getConnectivityBlockIDs(std::vector<int>& blockIDs) const
01791 {
01792   blockIDs.resize(connectivityBlocks_.size());
01793 
01794   std::map<int,fei::ConnectivityBlock*>::const_iterator
01795     cdb_iter = connectivityBlocks_.begin();
01796 
01797   for(size_t i=0; i<blockIDs.size(); ++i, ++cdb_iter) {
01798     int blockID = (*cdb_iter).first;
01799     blockIDs[i] = blockID;
01800   }
01801 
01802   return(0);
01803 }
01804 
01805 //----------------------------------------------------------------------------
01806 int fei::MatrixGraph_Impl2::getNumIDsPerConnectivityList(int blockID) const
01807 {
01808   const fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01809   if (cblock == NULL) return(-1);
01810 
01811   const fei::Pattern* pattern = cblock->getRowPattern();
01812   return(pattern->getNumIDs());
01813 }
01814 
01815 //----------------------------------------------------------------------------
01816 int fei::MatrixGraph_Impl2::getConnectivityNumIndices(int blockID) const
01817 {
01818   const fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01819   if (cblock == NULL) return(-1);
01820 
01821   const fei::Pattern* pattern = cblock->getRowPattern();
01822   return( blockEntryGraph_ ?
01823     pattern->getNumIDs() : pattern->getNumIndices());
01824 }
01825 
01826 //----------------------------------------------------------------------------
01827 int fei::MatrixGraph_Impl2::getConnectivityNumIndices(int blockID,
01828                                                   int& numRowIndices,
01829                                                     int& numColIndices)
01830 {
01831   fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01832   if (cblock == NULL) return(-1);
01833 
01834   fei::Pattern* pattern = cblock->getRowPattern();
01835   numRowIndices = pattern->getNumIndices();
01836   fei::Pattern* colpattern = cblock->isSymmetric() ?
01837     pattern : cblock->getColPattern();
01838   numColIndices = colpattern->getNumIndices();
01839 
01840   return(0);
01841 }
01842 
01843 //----------------------------------------------------------------------------
01844 const fei::ConnectivityBlock* fei::MatrixGraph_Impl2::getConnectivityBlock(int blockID) const
01845 {
01846   std::map<int,fei::ConnectivityBlock*>::const_iterator
01847     c_iter = connectivityBlocks_.find(blockID);
01848   if (c_iter == connectivityBlocks_.end()) return(0);
01849 
01850   return( (*c_iter).second );
01851 }
01852 
01853 //----------------------------------------------------------------------------
01854 fei::ConnectivityBlock* fei::MatrixGraph_Impl2::getConnectivityBlock(int blockID)
01855 {
01856   std::map<int,fei::ConnectivityBlock*>::iterator
01857     c_iter = connectivityBlocks_.find(blockID);
01858   if (c_iter == connectivityBlocks_.end()) return(0);
01859 
01860   return( (*c_iter).second );
01861 }
01862 
01863 //----------------------------------------------------------------------------
01864 int fei::MatrixGraph_Impl2::getConnectivityIndices(int blockID,
01865                                              int connectivityID,
01866                                              int indicesAllocLen,
01867                                              int* indices,
01868                                              int& numIndices)
01869 {
01870   fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01871   if (cblock == NULL) return(-1);
01872 
01873   std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
01874   vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
01875 
01876   fei::Pattern* pattern = cblock->getRowPattern();
01877   numIndices = pattern->getNumIndices();
01878 
01879   int len = numIndices > indicesAllocLen ? indicesAllocLen : numIndices;
01880 
01881   int* records = cblock->getRowConnectivity(connectivityID);
01882   if (records == NULL) {
01883     ERReturn(-1);
01884   }
01885 
01886   fei::Pattern::PatternType pType = pattern->getPatternType();
01887 
01888   if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
01889     const int* numFieldsPerID = pattern->getNumFieldsPerID();
01890     const int* fieldIDs = pattern->getFieldIDs();
01891     int totalNumFields = pattern->getTotalNumFields();
01892 
01893     std::vector<int> fieldSizes(totalNumFields);
01894 
01895     for(int j=0; j<totalNumFields; ++j) {
01896       fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
01897                                             colSpace_.get());
01898     }
01899 
01900     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
01901                                             records, pattern->getNumIDs(),
01902                                                numFieldsPerID,
01903                                                fieldIDs, &fieldSizes[0],
01904                                                len, indices, numIndices) );
01905   }
01906   else if (pType == fei::Pattern::SIMPLE) {
01907     const int* fieldIDs = pattern->getFieldIDs();
01908 
01909     int fieldID = fieldIDs[0];
01910     unsigned fieldSize = snl_fei::getFieldSize(fieldID,
01911                                                rowSpace_.get(),
01912                                                colSpace_.get());
01913 
01914     CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
01915                                             records, pattern->getNumIDs(),
01916                                                 fieldID, fieldSize,
01917                                                 len, indices, numIndices) );
01918   }
01919   else if (pType == fei::Pattern::NO_FIELD) {
01920     CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(),
01921                                             records, pattern->getNumIDs(),
01922                                             len, indices, numIndices) );
01923   }
01924 
01925   return(0);
01926 }
01927 
01928 //----------------------------------------------------------------------------
01929 int fei::MatrixGraph_Impl2::getConnectivityIndices(int blockID,
01930                                              int connectivityID,
01931                                              int rowIndicesAllocLen,
01932                                              int* rowIndices,
01933                                              int& numRowIndices,
01934                                              int colIndicesAllocLen,
01935                                              int* colIndices,
01936                                              int& numColIndices)
01937 {
01938   fei::ConnectivityBlock* cblock = getConnectivityBlock(blockID);
01939   if (cblock == NULL) return(-1);
01940 
01941   std::vector<int>& eqnNums = rowSpace_->getEqnNumbers();
01942   vspcEqnPtr_ = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
01943 
01944   fei::Pattern* pattern = cblock->getRowPattern();
01945   numRowIndices = pattern->getNumIndices();
01946 
01947   int len = numRowIndices > rowIndicesAllocLen ?
01948     rowIndicesAllocLen : numRowIndices;
01949 
01950   int* records = cblock->getRowConnectivity(connectivityID);
01951   if (records == NULL) {
01952     ERReturn(-1);
01953   }
01954 
01955   fei::Pattern::PatternType pType = pattern->getPatternType();
01956 
01957   if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
01958     const int* numFieldsPerID = pattern->getNumFieldsPerID();
01959     const int* fieldIDs = pattern->getFieldIDs();
01960     int totalNumFields = pattern->getTotalNumFields();
01961 
01962     std::vector<int> fieldSizes(totalNumFields);
01963 
01964     for(int j=0; j<totalNumFields; ++j) {
01965       fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
01966                                             colSpace_.get());
01967     }
01968 
01969     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
01970                                                records, pattern->getNumIDs(),
01971                                                numFieldsPerID,
01972                                                fieldIDs, &fieldSizes[0],
01973                                                len, rowIndices, numRowIndices) );
01974   }
01975   else if (pType == fei::Pattern::SIMPLE) {
01976     const int* fieldIDs = pattern->getFieldIDs();
01977 
01978     int fieldID = fieldIDs[0];
01979     unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
01980                                             colSpace_.get());
01981 
01982     CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
01983                                                 records, pattern->getNumIDs(),
01984                                                 fieldID, fieldSize,
01985                                                 len, rowIndices, numRowIndices) );
01986   }
01987   else if (pType == fei::Pattern::NO_FIELD) {
01988     CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(),
01989                                             records, pattern->getNumIDs(),
01990                                             len, rowIndices, numRowIndices) );
01991   }
01992 
01993   fei::Pattern* colpattern = cblock->isSymmetric() ? 
01994     pattern : cblock->getColPattern();
01995   pType = colpattern->getPatternType();
01996   numColIndices = colpattern->getNumIndices();
01997   len = numColIndices > colIndicesAllocLen ?
01998     colIndicesAllocLen : numColIndices;
01999 
02000   if (!cblock->isSymmetric()) {
02001     records = cblock->getColConnectivity(connectivityID);
02002   }
02003   if (records == NULL) {
02004     return(-1);
02005   }
02006 
02007   if (pType == fei::Pattern::GENERAL || pType == fei::Pattern::SINGLE_IDTYPE) {
02008     const int* numFieldsPerID = colpattern->getNumFieldsPerID();
02009     const int* fieldIDs = colpattern->getFieldIDs();
02010     int totalNumFields = colpattern->getTotalNumFields();
02011 
02012     std::vector<int> fieldSizes(totalNumFields);
02013 
02014     for(int j=0; j<totalNumFields; ++j) {
02015       fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
02016                                             colSpace_.get());
02017     }
02018 
02019     CHK_ERR( getConnectivityIndices_multiField(colpattern->getRecordCollections(),
02020                                                records, colpattern->getNumIDs(),
02021                                                numFieldsPerID,
02022                                                fieldIDs, &fieldSizes[0],
02023                                                len, colIndices, numColIndices) );
02024   }
02025   else if (pType == fei::Pattern::SIMPLE) {
02026     const int* fieldIDs = colpattern->getFieldIDs();
02027 
02028     int fieldID = fieldIDs[0];
02029     unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
02030                                                colSpace_.get());
02031 
02032     CHK_ERR( getConnectivityIndices_singleField(colpattern->getRecordCollections(),
02033                                                 records, colpattern->getNumIDs(),
02034                                                 fieldID, fieldSize,
02035                                                 len, colIndices, numColIndices) );
02036   }
02037   else if (pType == fei::Pattern::NO_FIELD) {
02038     CHK_ERR( getConnectivityIndices_noField(colpattern->getRecordCollections(),
02039                                             records, colpattern->getNumIDs(),
02040                                             len, colIndices, numColIndices) );
02041   }
02042 
02043   return(0);
02044 }
02045 
02046 //----------------------------------------------------------------------------
02047 int fei::MatrixGraph_Impl2::getLocalNumLagrangeConstraints() const
02048 {
02049   return( lagrangeConstraints_.size() );
02050 }
02051 
02052 //----------------------------------------------------------------------------
02053 int fei::MatrixGraph_Impl2::addBlockToGraph_multiField_symmetric(fei::Graph* graph,
02054                                                         fei::ConnectivityBlock* cblock)
02055 {
02056   fei::Pattern* pattern = cblock->getRowPattern();
02057   int j;
02058   int numIDs = pattern->getNumIDs();
02059   int numIndices = blockEntryGraph_ ?
02060     pattern->getNumIDs() : pattern->getNumIndices();
02061   int checkNumIndices = numIndices;
02062   std::vector<int> indices(numIndices);
02063   int* indicesPtr = &indices[0];
02064 
02065   const int* numFieldsPerID = pattern->getNumFieldsPerID();
02066   const int* fieldIDs = pattern->getFieldIDs();
02067   int totalNumFields = pattern->getTotalNumFields();
02068 
02069   std::vector<int> fieldSizes(totalNumFields);
02070 
02071   for(j=0; j<totalNumFields; ++j) {
02072     fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
02073                                           colSpace_.get());
02074   }
02075 
02076   std::map<int,int>& connIDs = cblock->getConnectivityIDs();
02077   std::vector<int>& values = cblock->getRowConnectivities();
02078 
02079   std::map<int,int>::iterator
02080     c_iter = connIDs.begin(),
02081     c_end  = connIDs.end();
02082 
02083   for(; c_iter != c_end; ++c_iter) {
02084     int offset = c_iter->second;
02085     int* records = &values[offset*numIDs];
02086 
02087     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
02088                                                 records, numIDs,
02089                                                  numFieldsPerID,
02090                                                  fieldIDs,
02091                                                  &fieldSizes[0],
02092                                                  numIndices,
02093                                                  indicesPtr,
02094                                                  checkNumIndices) );
02095 
02096     if (checkNumIndices != numIndices) {
02097       ERReturn(-1);
02098     }
02099 
02100     if (output_level_ > fei::BRIEF_LOGS && output_stream_ != NULL) {
02101       FEI_OSTREAM& os = *output_stream_;
02102 
02103       const snl_fei::RecordCollection*const* recordColls = pattern->getRecordCollections();
02104       unsigned thisoffset = 0;
02105       for(int ii=0; ii<numIDs; ++ii) {
02106         const fei::Record<int>* record = recordColls[ii]->getRecordWithLocalID(records[ii]);
02107         int ID = record->getID();
02108         os << dbgprefix_<<"scatterIndices: ID=" <<ID<<": ";
02109         int num = pattern->getNumIndicesPerID()[ii];
02110         for(int jj=0; jj<num; ++jj) {
02111           os << indicesPtr[thisoffset++] << " ";
02112         }
02113         os << FEI_ENDL;
02114       }
02115 
02116       for(int ii=0; ii<numIndices; ++ii) {
02117         if (isLogEqn(indicesPtr[ii])) {
02118           os << "adding Symm inds: ";
02119           for(int jj=0; jj<numIndices; ++jj) {
02120             os << indicesPtr[jj] << " ";
02121           }
02122           os << FEI_ENDL;
02123           break;
02124         }
02125       }
02126     }
02127 
02128     //now we have the indices array, so we're ready to push it into
02129     //the graph container
02130     if (numIndices == numIDs || !cblock->isDiagonal()) {
02131       CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr,
02132                                           cblock->isDiagonal()) );
02133     }
02134     else {
02135       int ioffset = 0;
02136       int* fieldSizesPtr = &fieldSizes[0];
02137       for(int i=0; i<numIDs; ++i) {
02138         CHK_ERR( graph->addSymmetricIndices(fieldSizesPtr[i], &(indicesPtr[ioffset])));
02139         ioffset += fieldSizesPtr[i];
02140       }
02141     }
02142   }
02143 
02144   return(0);
02145 }
02146 
02147 //----------------------------------------------------------------------------
02148 int fei::MatrixGraph_Impl2::addBlockToGraph_multiField_nonsymmetric(fei::Graph* graph,
02149                                                         fei::ConnectivityBlock* cblock)
02150 {
02151   fei::Pattern* pattern = cblock->getRowPattern();
02152   fei::Pattern* colpattern = cblock->getColPattern();
02153   int j;
02154   int numIDs = pattern->getNumIDs();
02155   int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
02156   int checkNumIndices = numIndices;
02157   std::vector<int> indices(numIndices);
02158   int* indicesPtr = &indices[0];
02159 
02160   int numColIDs = colpattern->getNumIDs();
02161   int numColIndices = blockEntryGraph_ ? numColIDs : colpattern->getNumIndices();
02162   int checkNumColIndices = numColIndices;
02163   std::vector<int> colindices(numColIndices);
02164   int* colindicesPtr = &colindices[0];
02165 
02166   const int* numFieldsPerID = pattern->getNumFieldsPerID();
02167   const int* fieldIDs = pattern->getFieldIDs();
02168   int totalNumFields = pattern->getTotalNumFields();
02169 
02170   const int* numFieldsPerColID = colpattern->getNumFieldsPerID();
02171   const int* colfieldIDs = colpattern->getFieldIDs();
02172   int totalNumColFields = colpattern->getTotalNumFields();
02173 
02174   std::vector<int> fieldSizes(totalNumFields);
02175   std::vector<int> colfieldSizes(totalNumColFields);
02176 
02177   for(j=0; j<totalNumFields; ++j) {
02178     fieldSizes[j] = snl_fei::getFieldSize(fieldIDs[j], rowSpace_.get(),
02179                                           colSpace_.get());
02180   }
02181   for(j=0; j<totalNumColFields; ++j) {
02182     colfieldSizes[j] = snl_fei::getFieldSize(colfieldIDs[j], rowSpace_.get(),
02183                                            colSpace_.get());
02184   }
02185 
02186   std::map<int,int>& connIDs = cblock->getConnectivityIDs();
02187   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02188   std::vector<int>& colrecords = cblock->getColConnectivities();
02189 
02190   std::map<int,int>::iterator
02191    c_iter = connIDs.begin(),
02192    c_end  = connIDs.end();
02193 
02194   for(; c_iter != c_end; ++c_iter) {
02195     int offset = c_iter->second;
02196     int* records = &rowrecords[offset*numIDs];
02197 
02198     int* colRecords = &colrecords[offset*numColIDs];
02199 
02200     CHK_ERR( getConnectivityIndices_multiField(pattern->getRecordCollections(),
02201                                                 records, numIDs,
02202                                                  numFieldsPerID,
02203                                                  fieldIDs,
02204                                                  &fieldSizes[0],
02205                                                  numIndices,
02206                                                  indicesPtr,
02207                                                  checkNumIndices) );
02208 
02209     if (checkNumIndices != numIndices) {
02210       ERReturn(-1);
02211     }
02212 
02213       
02214     CHK_ERR( getConnectivityIndices_multiField(colpattern->getRecordCollections(),
02215                                                 colRecords, numColIDs,
02216                                                  numFieldsPerColID,
02217                                                  colfieldIDs,
02218                                                  &colfieldSizes[0],
02219                                                  numColIndices,
02220                                                  colindicesPtr,
02221                                                  checkNumColIndices) );
02222 
02223     if (checkNumColIndices != numColIndices) {
02224       ERReturn(-1);
02225     }
02226 
02227     //now we have the indices arrays, so we're ready to push them into
02228     //the graph container
02229     for(int r=0; r<numIndices; ++r) {
02230       CHK_ERR( graph->addIndices(indicesPtr[r], numColIndices, colindicesPtr));
02231     }
02232   }
02233 
02234   return(0);
02235 }
02236 
02237 //----------------------------------------------------------------------------
02238 int fei::MatrixGraph_Impl2::getConnectivityIndices_multiField(const snl_fei::RecordCollection*const* recordCollections,
02239                                                           int* records,
02240                                                           int numRecords,
02241                                                           const int* numFieldsPerID,
02242                                                           const int* fieldIDs,
02243                                                           const int* fieldSizes,
02244                                                           int indicesAllocLen,
02245                                                           int* indices,
02246                                                           int& numIndices)
02247 {
02248   numIndices = 0;
02249   int fld_offset = 0;
02250 
02251   for(int i=0; i<numRecords; ++i) {
02252     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
02253     if (record==NULL) continue;
02254 
02255     if (blockEntryGraph_) {
02256       indices[numIndices++] = record->getNumber();
02257       continue;
02258     }
02259 
02260     const fei::FieldMask* fieldMask = record->getFieldMask();
02261     int* eqnNumbers = vspcEqnPtr_ + record->getOffsetIntoEqnNumbers();
02262 
02263     for(int nf=0; nf<numFieldsPerID[i]; ++nf) {
02264       int eqnOffset = 0;
02265       if (!simpleProblem_) {
02266         fieldMask->getFieldEqnOffset(fieldIDs[fld_offset], eqnOffset);
02267       }
02268 
02269       for(int fs=0; fs<fieldSizes[fld_offset]; ++fs) {
02270         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
02271       }
02272 
02273       ++fld_offset;
02274     }
02275   }
02276 
02277   return(0);
02278 }
02279 
02280 //----------------------------------------------------------------------------
02281 int fei::MatrixGraph_Impl2::addBlockToGraph_singleField_symmetric(fei::Graph* graph,
02282                                                         fei::ConnectivityBlock* cblock)
02283 {
02284   fei::Pattern* pattern = cblock->getRowPattern();
02285   int numIDs = pattern->getNumIDs();
02286   int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
02287   int checkNumIndices = numIndices;
02288   std::vector<int> indices(numIndices);
02289   int* indicesPtr = &indices[0];
02290 
02291   const int* fieldIDs = pattern->getFieldIDs();
02292 
02293   int fieldID = fieldIDs[0];
02294   unsigned fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
02295                                              colSpace_.get());
02296 
02297   std::map<int,int>& connIDs = cblock->getConnectivityIDs();
02298   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02299 
02300   std::map<int,int>::iterator
02301     c_iter = connIDs.begin(),
02302     c_end  = connIDs.end();
02303 
02304   for(; c_iter != c_end; ++c_iter) {
02305     int offset = c_iter->second;
02306     int* records = &rowrecords[offset*numIDs];
02307 
02308     CHK_ERR( getConnectivityIndices_singleField(pattern->getRecordCollections(),
02309                                                 records, numIDs,
02310                                                   fieldID, fieldSize,
02311                                                   checkNumIndices,
02312                                                   indicesPtr,
02313                                                   numIndices) );
02314 
02315     if (checkNumIndices != numIndices) {
02316       ERReturn(-1);
02317     }
02318 
02319     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
02320         for(int ii=0; ii<numIndices; ++ii) {
02321           if (isLogEqn(indicesPtr[ii])) {
02322             FEI_OSTREAM& os = *output_stream_;
02323             os << "adding Symm inds: ";
02324             for(int jj=0; jj<numIndices; ++jj) {
02325               os << indicesPtr[jj] << " ";
02326             }
02327             os << FEI_ENDL;
02328             break;
02329           }
02330         }
02331     }
02332 
02333     //now we have the indices array, so we're ready to push it into
02334     //the graph container
02335     if (numIndices == numIDs || !cblock->isDiagonal()) {
02336       CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr,
02337                                             cblock->isDiagonal()));
02338     }
02339     else {
02340       int ioffset = 0;
02341       for(int i=0; i<numIDs; ++i) {
02342         CHK_ERR( graph->addSymmetricIndices(fieldSize, &(indicesPtr[ioffset])));
02343         ioffset += fieldSize;
02344       }
02345     }
02346   }
02347 
02348   return(0);
02349 }
02350 
02351 //----------------------------------------------------------------------------
02352 int fei::MatrixGraph_Impl2::addBlockToGraph_singleField_nonsymmetric(fei::Graph* graph,
02353                                                         fei::ConnectivityBlock* cblock)
02354 {
02355   fei::Pattern* pattern = cblock->getRowPattern();
02356   fei::Pattern* colpattern = cblock->getColPattern();
02357   int numIDs = pattern->getNumIDs();
02358   int numIndices = blockEntryGraph_ ? numIDs : pattern->getNumIndices();
02359   int checkNumIndices = numIndices;
02360   std::vector<int> indices(numIndices);
02361   int* indicesPtr = &indices[0];
02362 
02363   int numColIDs = colpattern->getNumIDs();
02364   int numColIndices = blockEntryGraph_ ?
02365     numColIDs : colpattern->getNumIndices();
02366   int checkNumColIndices = numColIndices;
02367   std::vector<int> colindices(numColIndices);
02368   int* colindicesPtr = &colindices[0];
02369 
02370   const int* fieldIDs = pattern->getFieldIDs();
02371 
02372   int rowFieldID = fieldIDs[0];
02373   int rowFieldSize = snl_fei::getFieldSize(rowFieldID, rowSpace_.get(),
02374                                            colSpace_.get());
02375 
02376   std::map<int,int>& connIDs = cblock->getConnectivityIDs();
02377   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02378   std::vector<int>& colrecords = cblock->getColConnectivities();
02379 
02380   int colFieldID = colpattern->getFieldIDs()[0];
02381   int colFieldSize = snl_fei::getFieldSize(colFieldID, rowSpace_.get(),
02382                                            colSpace_.get());
02383 
02384   std::map<int,int>::iterator
02385     c_iter = connIDs.begin(),
02386     c_end  = connIDs.end();
02387 
02388   for(; c_iter != c_end; ++c_iter) {
02389     int offset = c_iter->second;
02390     int* records = &rowrecords[offset*numIDs];
02391 
02392     int* colRecords = &colrecords[offset*numColIDs];
02393 
02394     if (blockEntryGraph_) {
02395       rowSpace_->getGlobalBlkIndicesL(numIDs, pattern->getRecordCollections(),
02396                                       records, checkNumIndices,
02397                                             indicesPtr, numIndices);
02398 
02399       colSpace_->getGlobalBlkIndicesL(numColIDs, colpattern->getRecordCollections(),
02400                                       colRecords, checkNumColIndices,
02401                                             colindicesPtr, numColIndices);
02402     }
02403     else {
02404       rowSpace_->getGlobalIndicesL(numIDs, pattern->getRecordCollections(),
02405                                       records, rowFieldID,
02406                                          rowFieldSize, checkNumIndices,
02407                                          indicesPtr, numIndices);
02408 
02409       colSpace_->getGlobalIndicesL(numColIDs, colpattern->getRecordCollections(),
02410                                      colRecords, colFieldID,
02411                                          colFieldSize, checkNumColIndices,
02412                                          colindicesPtr, numColIndices);
02413     }
02414 
02415     if (checkNumIndices != numIndices || checkNumColIndices != numColIndices) {
02416       ERReturn(-1);
02417     }
02418 
02419     for(int r=0; r<numIndices; ++r) {
02420       CHK_ERR( graph->addIndices(indicesPtr[r], numColIndices, colindicesPtr));
02421     }
02422   }
02423 
02424   return(0);
02425 }
02426 
02427 //----------------------------------------------------------------------------
02428 int fei::MatrixGraph_Impl2::getConnectivityIndices_singleField(const snl_fei::RecordCollection*const* recordCollections,
02429                                                              int* records,
02430                                                              int numRecords,
02431                                                              int fieldID,
02432                                                              int fieldSize,
02433                                                              int indicesAllocLen,
02434                                                              int* indices,
02435                                                              int& numIndices)
02436 {
02437   numIndices = 0;
02438 
02439   for(int i=0; i<numRecords; ++i) {
02440     if (numIndices == indicesAllocLen) break;
02441 
02442     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
02443 
02444     if (blockEntryGraph_) {
02445       indices[numIndices++] = record->getNumber();
02446       continue;
02447     }
02448 
02449     int* eqnNumbers = vspcEqnPtr_+record->getOffsetIntoEqnNumbers();
02450 
02451     int eqnOffset = 0;
02452     if (!simpleProblem_) {
02453       const fei::FieldMask* fieldMask = record->getFieldMask();
02454       fieldMask->getFieldEqnOffset(fieldID, eqnOffset);
02455     }
02456 
02457     indices[numIndices++] = eqnNumbers[eqnOffset];
02458     if (fieldSize > 1) {
02459       for(int fs=1; fs<fieldSize; ++fs) {
02460         indices[numIndices++] = eqnNumbers[eqnOffset+fs];
02461       }
02462     }
02463   }
02464 
02465   return(0);
02466 }
02467 
02468 //----------------------------------------------------------------------------
02469 int fei::MatrixGraph_Impl2::getConnectivityIndices_noField(const snl_fei::RecordCollection*const* recordCollections,
02470                                                          int* records,
02471                                                          int numRecords,
02472                                                          int indicesAllocLen,
02473                                                          int* indices,
02474                                                          int& numIndices)
02475 {
02476   numIndices = 0;
02477 
02478   for(int i=0; i<numRecords; ++i) {
02479 
02480     const fei::Record<int>* record = recordCollections[i]->getRecordWithLocalID(records[i]);
02481     int* eqnNumbers = vspcEqnPtr_+record->getOffsetIntoEqnNumbers();
02482 
02483     if (blockEntryGraph_) {
02484       indices[numIndices++] = record->getNumber();
02485     }
02486     else {
02487       indices[numIndices++] = eqnNumbers[0];
02488     }
02489   }
02490 
02491   return(0);
02492 }
02493 
02494 //----------------------------------------------------------------------------
02495 int fei::MatrixGraph_Impl2::addBlockToGraph_noField_symmetric(fei::Graph* graph,
02496                                                         fei::ConnectivityBlock* cblock)
02497 {
02498   fei::Pattern* pattern = cblock->getRowPattern();
02499   int numIDs = pattern->getNumIDs();
02500   int numIndices = pattern->getNumIndices();
02501   std::vector<int> indices(numIndices);
02502   int* indicesPtr = &indices[0];
02503 
02504   std::map<int,int>& connIDs = cblock->getConnectivityIDs();
02505   std::vector<int>& rowrecords = cblock->getRowConnectivities();
02506 
02507   std::map<int,int>::iterator
02508     c_iter = connIDs.begin(),
02509     c_end  = connIDs.end();
02510 
02511   for(; c_iter != c_end; ++c_iter) {
02512     int offset = c_iter->second;
02513     int* records = &rowrecords[offset*numIDs];
02514 
02515     int checkNumIndices;
02516     CHK_ERR( getConnectivityIndices_noField(pattern->getRecordCollections(), records, numIDs, numIndices,
02517                                             indicesPtr, checkNumIndices) );
02518 
02519     if (checkNumIndices != numIndices) {
02520       ERReturn(-1);
02521     }
02522 
02523     if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
02524       for(int ii=0; ii<numIndices; ++ii) {
02525         if (isLogEqn(indicesPtr[ii])) {
02526           FEI_OSTREAM& os = *output_stream_;
02527           os << "adding Symm inds: ";
02528           for(int jj=0; jj<numIndices; ++jj) {
02529             os << indicesPtr[jj] << " ";
02530           }
02531           os << FEI_ENDL;
02532           break;
02533         }
02534       }
02535     }
02536 
02537     //now we have the indices array, so we're ready to push it into
02538     //the graph container
02539     CHK_ERR( graph->addSymmetricIndices(numIndices, indicesPtr) );
02540   }
02541 
02542   return(0);
02543 }
02544 
02545 //----------------------------------------------------------------------------
02546 int fei::MatrixGraph_Impl2::addBlockToGraph_sparse(fei::Graph* graph,
02547                                              fei::ConnectivityBlock* cblock)
02548 {
02549   std::vector<int> row_indices;
02550   std::vector<int> indices;
02551 
02552   fei::Pattern* pattern = cblock->getRowPattern();
02553   const snl_fei::RecordCollection*const* recordCollections = pattern->getRecordCollections();
02554 
02555   std::map<int,int>& connIDs = cblock->getConnectivityIDs();
02556   std::vector<int>& connOffsets = cblock->getConnectivityOffsets();
02557   int* rowrecords = &(cblock->getRowConnectivities()[0]);
02558   int* colrecords = &(cblock->getColConnectivities()[0]);
02559   bool haveField = cblock->haveFieldID();
02560   int fieldID = cblock->fieldID();
02561   int fieldSize = 1;
02562 
02563   if (haveField) {
02564     fieldSize = snl_fei::getFieldSize(fieldID, rowSpace_.get(),
02565                                       colSpace_.get());
02566   }
02567 
02568   std::map<int,int>::iterator
02569     c_iter = connIDs.begin(),
02570     c_end  = connIDs.end();
02571 
02572   for(; c_iter != c_end; ++c_iter) {
02573     int offset = c_iter->second;
02574     int rowlen = connOffsets[offset+1] - offset;
02575 
02576     int* records = &(rowrecords[offset]);
02577 
02578     int checkNumIndices;
02579     row_indices.resize(fieldSize);
02580 
02581     if (haveField) {
02582       CHK_ERR( getConnectivityIndices_singleField(recordCollections, records, 1,
02583                                                   fieldID, fieldSize,
02584                                                   fieldSize,
02585                                                   &row_indices[0],
02586                                                   checkNumIndices) );
02587     }
02588     else {
02589       CHK_ERR( getConnectivityIndices_noField(recordCollections, records, 1, fieldSize,
02590                                               &row_indices[0],
02591                                               checkNumIndices) );
02592     }
02593 
02594     indices.resize(fieldSize*rowlen);
02595     int* indicesPtr = &indices[0];
02596     int* crecords = &(colrecords[offset]);
02597 
02598     if (haveField) {
02599       CHK_ERR( getConnectivityIndices_singleField(recordCollections, crecords, rowlen,
02600                                                   fieldID, fieldSize,
02601                                                   fieldSize*rowlen,
02602                                                   indicesPtr,
02603                                                   checkNumIndices) );
02604     }
02605     else {
02606       CHK_ERR( getConnectivityIndices_noField(recordCollections, crecords, rowlen, rowlen,
02607                                               indicesPtr, checkNumIndices) );
02608     }
02609     if (checkNumIndices != rowlen) {
02610       ERReturn(-1);
02611     }
02612 
02613     //now we have the indices, so we're ready to push them into
02614     //the graph container
02615     int* row_ind_ptr = &row_indices[0];
02616     for(int i=0; i<fieldSize; ++i) {
02617       CHK_ERR( graph->addIndices(row_ind_ptr[i], fieldSize*rowlen,
02618                                  indicesPtr) );
02619     }
02620   }
02621 
02622   return(0);
02623 }
02624 
02625 //----------------------------------------------------------------------------
02626 void fei::MatrixGraph_Impl2::setName(const char* name)
02627 {
02628   if (name == NULL) return;
02629 
02630   if (name_ == name) return;
02631 
02632   name_ = name;
02633   dbgprefix_ = "MatGraph_"+name_+": ";
02634 }
02635 
02636 //----------------------------------------------------------------------------
02637 void fei::MatrixGraph_Impl2::setIndicesMode(int mode)
02638 {
02639   if (mode == BLOCK_ENTRY_GRAPH) {
02640     blockEntryGraph_ = true;
02641     return;
02642   }
02643 
02644   if (mode == POINT_ENTRY_GRAPH) {
02645     blockEntryGraph_ = false;
02646     return;
02647   }
02648 
02649   voidERReturn;
02650 }
02651 
02652 //----------------------------------------------------------------------------
02653 fei::SharedPtr<fei::FillableMat>
02654 fei::MatrixGraph_Impl2::getSlaveDependencyMatrix()
02655 {
02656   return( D_ );
02657 }
02658 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends