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

Generated on Tue Jul 13 09:27:45 2010 for FEI by  doxygen 1.4.7