FEI Version of the Day
fei_LinSysCoreFilter.cpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 #include "fei_sstream.hpp"
00045 
00046 #include <limits>
00047 #include <cmath>
00048 #include <assert.h>
00049 
00050 #include "fei_utils.hpp"
00051 #include <fei_impl_utils.hpp>
00052 
00053 #include "fei_defs.h"
00054 
00055 #include <fei_ostream_ops.hpp>
00056 #include "fei_CommUtils.hpp"
00057 #include "fei_TemplateUtils.hpp"
00058 #include "snl_fei_Constraint.hpp"
00059 typedef snl_fei::Constraint<GlobalID> ConstraintType;
00060 
00061 #include "fei_DirichletBCManager.hpp"
00062 #include "fei_FillableMat.hpp"
00063 #include "fei_CSVec.hpp"
00064 #include "FEI_Implementation.hpp"
00065 #include "fei_EqnCommMgr.hpp"
00066 #include "fei_ConnectivityTable.hpp"
00067 #include "fei_NodeDescriptor.hpp"
00068 #include "fei_NodeDatabase.hpp"
00069 #include "fei_BlockDescriptor.hpp"
00070 #include "SNL_FEI_Structure.hpp"
00071 #include "snl_fei_Utils.hpp"
00072 #include "fei_Data.hpp"
00073 #include "fei_LinearSystemCore.hpp"
00074 #include "fei_LinSysCore_flexible.hpp"
00075 
00076 #include "fei_LinSysCoreFilter.hpp"
00077 
00078 #undef fei_file
00079 #define fei_file "fei_LinSysCoreFilter.cpp"
00080 #include "fei_ErrMacros.hpp"
00081 
00082 #define ASSEMBLE_PUT 0
00083 #define ASSEMBLE_SUM 1
00084 
00085 
00086 //------------------------------------------------------------------------------
00087 LinSysCoreFilter::LinSysCoreFilter(FEI_Implementation* owner,
00088                                    MPI_Comm comm,
00089                                    SNL_FEI_Structure* probStruct,
00090                                    LinearSystemCore* lsc,
00091                                    int masterRank)
00092  : Filter(probStruct),
00093    timesInitializeCalled_(0),
00094    lsc_(lsc),
00095    useLookup_(true),
00096    newMatrixData_(false),
00097    newVectorData_(false),
00098    newConstraintData_(false),
00099    newBCData_(false),
00100    connectivitiesInitialized_(false),
00101    firstRemEqnExchange_(true),
00102    needToCallMatrixLoadComplete_(false),
00103    resolveConflictRequested_(false),
00104    localStartRow_(0),             //
00105    localEndRow_(0),               //Initialize all private variables here,
00106    numLocalEqns_(0),              //in the order that they are declared.
00107    numGlobalEqns_(0),             //(Don't bother initializing those that will
00108    numLocallyOwnedNodes_(0),      //be initialized in the body of the 
00109    numGlobalNodes_(0),            //constructor below.)
00110    firstLocalNodeNumber_(-1),     //
00111    blockMatrix_(false),
00112    tooLateToChooseBlock_(false),
00113    numLocalEqnBlks_(0),
00114    localReducedBlkOffset_(0),
00115    numLocalReducedEqnBlks_(0),
00116    iterations_(0),
00117    currentRHS_(0),
00118    rhsIDs_(),
00119    outputLevel_(0),
00120    comm_(comm),
00121    masterRank_(masterRank),
00122    problemStructure_(probStruct),
00123    matrixAllocated_(false),
00124    rowIndices_(),
00125    rowColOffsets_(0),
00126    colIndices_(0),
00127    nodeIDType_(0),
00128    eqnCommMgr_(NULL),
00129    eqnCommMgr_put_(NULL),
00130    maxElemRows_(0),
00131    scatterIndices_(),
00132    blkScatterIndices_(),
00133    iworkSpace_(),
00134    iworkSpace2_(),
00135    dworkSpace_(),
00136    dworkSpace2_(),
00137    eStiff_(NULL),
00138    eStiff1D_(NULL),
00139    eLoad_(NULL)
00140 {
00141     localRank_ = fei::localProc(comm_);
00142     numProcs_ = fei::numProcs(comm_);
00143 
00144     internalFei_ = 0;
00145 
00146     numRHSs_ = 1;
00147     rhsIDs_.resize(numRHSs_);
00148     rhsIDs_[0] = 0;
00149 
00150     bcManager_ = new fei::DirichletBCManager(probStruct);
00151     eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
00152     int err = createEqnCommMgr_put();
00153 
00154    //We need to get the parameters from the owning FEI_Implementation...
00155    int numParams = -1;
00156    char** paramStrings = NULL;
00157    err = owner->getParameters(numParams, paramStrings);
00158 
00159    //Now let's pass them into our own parameter-handling mechanism.
00160    err = parameters(numParams, paramStrings);
00161    if (err != 0) {
00162      fei::console_out() << "LinSysCoreFilter::LinSysCoreFilter ERROR, parameters failed." << FEI_ENDL;
00163      MPI_Abort(comm_, -1);
00164    }
00165 
00166    Kid_ = new fei::FillableMat;
00167    Kdi_ = new fei::FillableMat;
00168    Kdd_ = new fei::FillableMat;
00169    reducedEqnCounter_ = 0;
00170    reducedRHSCounter_ = 0;
00171 
00172    return;
00173 }
00174 
00175 //------------------------------------------------------------------------------
00176 LinSysCoreFilter::~LinSysCoreFilter() {
00177 //
00178 //  Destructor function. Free allocated memory, etc.
00179 //
00180   numRHSs_ = 0;
00181 
00182   delete bcManager_;
00183   delete eqnCommMgr_;
00184   delete eqnCommMgr_put_;
00185 
00186   delete [] eStiff_;
00187   delete [] eStiff1D_;
00188   delete [] eLoad_;
00189 
00190   delete Kid_;
00191   delete Kdi_;
00192   delete Kdd_;
00193 }
00194 
00195 //------------------------------------------------------------------------------
00196 int LinSysCoreFilter::initialize()
00197 {
00198 // Determine final sparsity pattern for setting the structure of the
00199 // underlying sparse matrix.
00200 //
00201   debugOutput("#LinSysCoreFilter::initialize");
00202 
00203   numLocalEqns_ = problemStructure_->getNumLocalEqns();
00204   numLocalEqnBlks_ = problemStructure_->getNumLocalEqnBlks();
00205   numLocalReducedEqnBlks_ = problemStructure_->getNumLocalReducedEqnBlks();
00206 
00207   // now, obtain the global equation info, such as how many equations there
00208   // are globally, and what the local starting and ending row-numbers are.
00209 
00210   // let's also get the number of global nodes, and a first-local-node-number.
00211   // node-number is a globally 0-based number we are assigning to nodes.
00212   // node-numbers are contiguous on a processor -- i.e., a processor owns a
00213   // contiguous block of node-numbers. This provides an easier-to-work-with
00214   // node numbering than the application-supplied nodeIDs which may not be
00215   // assumed to be contiguous or 0-based, or anything else.
00216 
00217   std::vector<int>& eqnOffsets = problemStructure_->getGlobalEqnOffsets();
00218   localStartRow_ = eqnOffsets[localRank_];
00219   localEndRow_ = eqnOffsets[localRank_+1]-1;
00220   numGlobalEqns_ = eqnOffsets[numProcs_];
00221 
00222   std::vector<int>& nodeOffsets = problemStructure_->getGlobalNodeOffsets();
00223   firstLocalNodeNumber_ =  nodeOffsets[localRank_];
00224   numGlobalNodes_ = nodeOffsets[numProcs_];
00225 
00226   //--------------------------------------------------------------------------
00227   //  ----- end active equation calculations -----
00228 
00229   if (eqnCommMgr_ != NULL) delete eqnCommMgr_;
00230   eqnCommMgr_ = NULL;
00231   if (eqnCommMgr_put_ != NULL) delete eqnCommMgr_put_;
00232   eqnCommMgr_put_ = NULL;
00233 
00234   eqnCommMgr_ = problemStructure_->getEqnCommMgr().deepCopy();
00235   if (eqnCommMgr_ == NULL) ERReturn(-1);
00236 
00237   int err = createEqnCommMgr_put();
00238   if (err != 0) ERReturn(err);
00239 
00240   //(we need to set the number of RHSs in the eqn comm manager)
00241   eqnCommMgr_->setNumRHSs(numRHSs_);
00242 
00243   //let's let the underlying linear system know what the global offsets are.
00244   //While we're dealing with global offsets, we'll also calculate the starting
00245   //and ending 'reduced' rows, etc.
00246   CHK_ERR( initLinSysCore() );
00247 
00248   setLinSysCoreCREqns();
00249 
00250   if (timesInitializeCalled_ == 0) {
00251     //
00252     // let's prepare some arrays for handing the matrix structure to
00253     // the linear system.
00254 
00255     std::vector<int> rowLengths;
00256     CHK_ERR( problemStructure_->getMatrixRowLengths(rowLengths) );
00257 
00258     int numReducedEqns = problemStructure_->getNumLocalReducedEqns();
00259     int maxBlkSize = problemStructure_->getGlobalMaxBlkSize();
00260     std::vector<int> blkSizes(numLocalReducedEqnBlks_, 1);
00261 
00262     int numNonzeros = 0;
00263     for(size_t ii=0; ii<rowLengths.size(); ++ii) {
00264       numNonzeros += rowLengths[ii];
00265     }
00266 
00267     std::vector<int> indices_1D(numNonzeros);
00268     std::vector<int*> indices(numReducedEqns);
00269 
00270     int offset = 0;
00271     for(size_t ii=0; ii<rowLengths.size(); ++ii) {
00272       indices[ii] = &(indices_1D[offset]);
00273       offset += rowLengths[ii];
00274     }
00275 
00276     if (maxBlkSize == 0) ERReturn(-1);
00277 
00278     if (maxBlkSize == 1) {
00279       CHK_ERR( problemStructure_->getMatrixStructure(&indices[0], rowLengths) );
00280 
00281       debugOutput("#LinSysCoreFilter calling point lsc_->setMatrixStructure");
00282       CHK_ERR( lsc_->setMatrixStructure(&indices[0], &rowLengths[0],
00283                                         &indices[0], &rowLengths[0], &blkSizes[0]) );
00284     }
00285     else {
00286       std::vector<int> blkRowLengths;
00287       int* blkIndices_1D = numNonzeros > 0 ? new int[numNonzeros] : NULL;
00288       int** blkIndices = numLocalReducedEqnBlks_ ?
00289         new int*[numLocalReducedEqnBlks_] : NULL;
00290       if (blkIndices == NULL && numLocalReducedEqnBlks_ != 0) ERReturn(-1);
00291 
00292       CHK_ERR( problemStructure_->getMatrixStructure(&indices[0], rowLengths,
00293                                                      blkIndices, blkIndices_1D,
00294                                                      blkRowLengths, blkSizes) );
00295 
00296       offset = 0;
00297       for(int ii=0; ii<numLocalReducedEqnBlks_; ++ii) {
00298         blkIndices[ii] = &(blkIndices_1D[offset]);
00299         offset += blkRowLengths[ii];
00300       }
00301 
00302       debugOutput("#LinSysCoreFilter calling block lsc_->setMatrixStructure");
00303       CHK_ERR( lsc_->setMatrixStructure(&indices[0], &rowLengths[0],
00304                                         blkIndices, &blkRowLengths[0], &blkSizes[0]) );
00305 
00306       if (numLocalReducedEqnBlks_ != 0) delete [] blkIndices;
00307       if (numNonzeros > 0) delete [] blkIndices_1D;
00308     }
00309   }
00310 
00311   matrixAllocated_ = true;
00312 
00313   debugOutput("#leaving LinSysCoreFilter::initialize");
00314   ++timesInitializeCalled_;
00315   return(FEI_SUCCESS);
00316 }
00317 
00318 //------------------------------------------------------------------------------
00319 int LinSysCoreFilter::createEqnCommMgr_put()
00320 {
00321   if (eqnCommMgr_put_ != NULL) return(0);
00322 
00323   eqnCommMgr_put_  = eqnCommMgr_->deepCopy();
00324   if (eqnCommMgr_put_ == NULL) ERReturn(-1);
00325 
00326   eqnCommMgr_put_->resetCoefs();
00327   eqnCommMgr_put_->accumulate_ = false;
00328   return(0);
00329 }
00330 
00331 //==============================================================================
00332 int LinSysCoreFilter::initLinSysCore()
00333 {
00334   debugOutput("#LinSysCoreFilter calling lsc_->setLookup");
00335   int err = lsc_->setLookup(*problemStructure_);
00336 
00337   if (err != 0) {
00338     useLookup_ = false;
00339   }
00340 
00341   std::vector<int>& globalNodeOffsets = problemStructure_->getGlobalNodeOffsets();
00342   std::vector<int>& globalEqnOffsets = problemStructure_->getGlobalEqnOffsets();
00343   std::vector<int>& globalBlkEqnOffsets =
00344     problemStructure_->getGlobalBlkEqnOffsets();
00345 
00346   int startRow = localStartRow_;
00347   while(problemStructure_->isSlaveEqn(startRow)) startRow++;
00348   int endRow = localEndRow_;
00349   while(problemStructure_->isSlaveEqn(endRow)) endRow--;
00350 
00351   problemStructure_->translateToReducedEqn(startRow, reducedStartRow_);
00352   problemStructure_->translateToReducedEqn(endRow, reducedEndRow_);
00353   numReducedRows_ = reducedEndRow_ - reducedStartRow_ + 1;
00354 
00355   std::vector<int> reducedEqnOffsets(globalEqnOffsets.size());
00356   std::vector<int> reducedBlkEqnOffsets(globalBlkEqnOffsets.size());
00357   std::vector<int> reducedNodeOffsets(globalNodeOffsets.size());
00358 
00359   int numSlaveEqns = problemStructure_->numSlaveEquations();
00360 
00361   int reducedNodeNum =
00362     problemStructure_->getAssociatedNodeNumber(reducedStartRow_);
00363 
00364   std::vector<int> tmpSend(2), tmpRecv, tmpRecvLengths;
00365   tmpSend[0] = reducedStartRow_;
00366   tmpSend[1] = reducedNodeNum;
00367   CHK_ERR( fei::Allgatherv(comm_, tmpSend, tmpRecvLengths, tmpRecv) );
00368 
00369   for(size_t ii=0; ii<globalEqnOffsets.size()-1; ++ii) {
00370     reducedEqnOffsets[ii] = tmpRecv[ii*2];
00371     reducedNodeOffsets[ii] = tmpRecv[ii*2+1];
00372 
00373     //Major assumption: we're assuming that if there are slave equations, then
00374     //there are no lagrange-multiplier constraints. Because if there are no
00375     //lagrange-multiplier constraints, then blkEqn == nodeNum.
00376     if (numSlaveEqns > 0) {
00377       reducedBlkEqnOffsets[ii] = reducedNodeOffsets[ii];
00378     }
00379     else {
00380       reducedBlkEqnOffsets[ii] = globalBlkEqnOffsets[ii];
00381     }
00382   }
00383 
00384   if (localRank_ == numProcs_-1) {
00385     reducedNodeNum = problemStructure_->
00386       translateToReducedNodeNumber(globalNodeOffsets[numProcs_]-1, localRank_);
00387     int blkEqn = globalBlkEqnOffsets[numProcs_]-1;
00388     if (numSlaveEqns > 0) {
00389       blkEqn = reducedNodeNum;
00390     }
00391 
00392     tmpSend.resize(3);
00393     tmpSend[0] = reducedEndRow_;
00394     tmpSend[1] = reducedNodeNum;
00395     tmpSend[2] = blkEqn;
00396   }
00397   else {
00398     tmpSend.resize(3);
00399   }
00400 
00401   CHK_ERR( fei::Bcast(comm_, tmpSend, numProcs_-1) );
00402   reducedEqnOffsets[numProcs_] = tmpSend[0]+1;
00403   reducedNodeOffsets[numProcs_] = tmpSend[1]+1;
00404   reducedBlkEqnOffsets[numProcs_] = tmpSend[2]+1;
00405 
00406   debugOutput("#LinSysCoreFilter calling lsc_->setGlobalOffsets");
00407   CHK_ERR( lsc_->setGlobalOffsets(numProcs_+1,
00408                                   &reducedNodeOffsets[0],
00409                                   &reducedEqnOffsets[0],
00410                                   &reducedBlkEqnOffsets[0]) );
00411 
00412   if (connectivitiesInitialized_) return(0);
00413 
00414   int numBlocks = problemStructure_->getNumElemBlocks();
00415   NodeDatabase& nodeDB     = problemStructure_->getNodeDatabase();
00416   NodeCommMgr& nodeCommMgr = problemStructure_->getNodeCommMgr();
00417 
00418   int numNodes = nodeDB.getNumNodeDescriptors();
00419   int numRemoteNodes = nodeCommMgr.getSharedNodeIDs().size() -
00420                          nodeCommMgr.getLocalNodeIDs().size();
00421   numNodes -= numRemoteNodes;
00422 
00423   std::vector<int> numElemsPerBlock(numBlocks);
00424   std::vector<int> numNodesPerElem(numBlocks);
00425   std::vector<int> numDofPerNode(0,1);
00426 
00427   for(int blk=0; blk<numBlocks; blk++) {
00428     BlockDescriptor* block = NULL;
00429     CHK_ERR( problemStructure_->getBlockDescriptor_index(blk, block) );
00430 
00431     numElemsPerBlock[blk] = block->getNumElements();
00432     numNodesPerElem[blk]  = block->getNumNodesPerElement();
00433 
00434     int* fieldsPerNode = block->fieldsPerNodePtr();
00435     int** fieldIDsTable = block->fieldIDsTablePtr();
00436 
00437     for(int nn=0; nn<numNodesPerElem[blk]; nn++) {
00438       if (fieldsPerNode[nn] <= 0) ERReturn(-1);
00439       numDofPerNode.push_back(0);
00440       int indx = numDofPerNode.size()-1;
00441       
00442       for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
00443         numDofPerNode[indx] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
00444       }
00445     }
00446   }
00447 
00448   for(int i=0; i<numBlocks; i++) {
00449     BlockDescriptor* block = NULL;
00450     CHK_ERR( problemStructure_->getBlockDescriptor_index(i, block) );
00451 
00452     if (block->getNumElements() == 0) continue;
00453 
00454     ConnectivityTable& ctbl =
00455       problemStructure_->getBlockConnectivity(block->getGlobalBlockID());
00456 
00457     std::vector<int> cNodeList(block->getNumNodesPerElement());
00458 
00459     int* fieldsPerNode = block->fieldsPerNodePtr();
00460     int** fieldIDsTable = block->fieldIDsTablePtr();
00461 
00462     numDofPerNode.resize(0);
00463     for(int nn=0; nn<numNodesPerElem[i]; nn++) {
00464       if (fieldsPerNode[nn] <= 0) ERReturn(-1);
00465       numDofPerNode.push_back(0);
00466       int indx = numDofPerNode.size()-1;
00467 
00468       for(int nf=0; nf<fieldsPerNode[nn]; nf++) {
00469         numDofPerNode[indx] += problemStructure_->getFieldSize(fieldIDsTable[nn][nf]);
00470       }
00471     }
00472 
00473     int nodesPerElement = block->getNumNodesPerElement();
00474     NodeDescriptor** elemNodePtrs = &((*ctbl.elem_conn_ptrs)[0]);
00475     int offset = 0;
00476     for(int j=0; j<block->getNumElements(); j++) {
00477 
00478       for(int k=0; k<nodesPerElement; k++) {
00479         NodeDescriptor* node = elemNodePtrs[offset++];
00480         cNodeList[k] = node->getNodeNumber();
00481       }
00482 
00483       int elemID = ctbl.elemNumbers[j];
00484       int* nodeNumbers = &cNodeList[0];
00485       debugOutput("#LinSysCoreFilter calling lsc->setConnectivities");
00486       CHK_ERR( lsc_->setConnectivities(i, 1, nodesPerElement,
00487                                        &elemID, &nodeNumbers) );
00488     }
00489   }
00490 
00491   connectivitiesInitialized_ = true;
00492 
00493   return(FEI_SUCCESS);
00494 }
00495 
00496 //==============================================================================
00497 void LinSysCoreFilter::setLinSysCoreCREqns()
00498 {
00499    int err, i=0;
00500 
00501    std::vector<int> iwork;
00502 
00503    std::map<GlobalID,ConstraintType*>::const_iterator
00504      cr_iter = problemStructure_->getMultConstRecords().begin(),
00505      cr_end = problemStructure_->getMultConstRecords().end();
00506 
00507    while(cr_iter != cr_end) {
00508       ConstraintType& constraint = *((*cr_iter).second);
00509       int numNodesPerCR = constraint.getMasters().size();
00510       int meqn = constraint.getEqnNumber();
00511 
00512       std::vector<GlobalID>& nodeID_vec = constraint.getMasters();
00513       GlobalID* nodeIDPtr = &nodeID_vec[0];
00514 
00515       if ((int)iwork.size() < 2*numNodesPerCR) {
00516         iwork.resize(2*numNodesPerCR);
00517       }
00518 
00519       int* nodeList = &(iwork[0]);
00520 
00521       int* eqnList = nodeList+numNodesPerCR;
00522 
00523       std::vector<int>& fieldIDs_vec = constraint.getMasterFieldIDs();
00524       int* fieldIDs = &fieldIDs_vec[0];
00525 
00526       for(int k=0; k<numNodesPerCR; k++) {
00527         const NodeDescriptor *node = Filter::findNode(nodeIDPtr[k]);
00528         if(node == NULL)
00529         {
00530           nodeList[k] = -1; // Indicates that the node wasn't found
00531         }
00532         else
00533         {
00534           nodeList[k] = node->getNodeNumber();
00535         }
00536 
00537         int eqn = -1; // Indicates that the equation wasn't found.
00538         if ( node ) {
00539           node->getFieldEqnNumber(fieldIDs[k], eqn);
00540         }
00541         eqnList[k] = eqn;
00542       }
00543 
00544       int crMultID = constraint.getConstraintID() + i++;
00545       if (Filter::logStream() != NULL) {
00546         FEI_OSTREAM& os = *logStream();
00547         os << "#LinSysCoreFilter calling lsc_->setMultCREqns"<<FEI_ENDL;
00548         os << "#  multiplier eqn: " << meqn << ", columns: ";
00549         for(int j=0; j<numNodesPerCR; ++j) os << eqnList[j] << " ";
00550         os << FEI_ENDL;
00551       }
00552 
00553       err = lsc_->setMultCREqns(crMultID, 1, numNodesPerCR,
00554                                 &nodeList, &eqnList,
00555                                 fieldIDs, &meqn);
00556       if (err) voidERReturn;
00557       ++cr_iter;
00558    }
00559 
00560    LinSysCore_flexible* lscf = dynamic_cast<LinSysCore_flexible*>(lsc_);
00561    if (lscf != NULL) {
00562      debugOutput("LinSysCoreFilter calling lscf->setMultCRComplete");
00563      err = lscf->setMultCRComplete();
00564      if (err != 0) {
00565        fei::console_out() << "LinSysCoreFilter::setLinSysCoreCREqns ERROR returned from "
00566             << "lscf->setMultCRComplete()" << FEI_ENDL;
00567      }
00568    }
00569 
00570    //
00571    //Now the penalty CRs...
00572    //
00573    cr_iter = problemStructure_->getPenConstRecords().begin();
00574    cr_end = problemStructure_->getPenConstRecords().end();
00575 
00576    while(cr_iter != cr_end) {
00577       ConstraintType& crset = *((*cr_iter).second);
00578       int numNodesPerCR = crset.getMasters().size();
00579 
00580       std::vector<GlobalID>& nodeIDs_vec = crset.getMasters();
00581       GlobalID* nodeIDsPtr = &nodeIDs_vec[0];
00582 
00583       if ((int)iwork.size() < 2*numNodesPerCR) {
00584         iwork.resize(2*numNodesPerCR);
00585       }
00586 
00587       int* nodeList = &(iwork[0]);
00588 
00589       int* eqnList = nodeList+numNodesPerCR;
00590 
00591       std::vector<int>& fieldIDs_vec = crset.getMasterFieldIDs();
00592       int* fieldIDs = &fieldIDs_vec[0];
00593 
00594       for(int k=0; k<numNodesPerCR; k++) {
00595         const NodeDescriptor& node = Filter::findNodeDescriptor(nodeIDsPtr[k]);
00596         nodeList[k] = node.getNodeNumber();
00597 
00598         int eqn = -1;
00599         if (!node.getFieldEqnNumber(fieldIDs[k], eqn)) voidERReturn;
00600 
00601         eqnList[k] = eqn;
00602       }
00603 
00604       int crPenID = crset.getConstraintID() + i;
00605       err = lsc_->setPenCREqns(crPenID, 1, numNodesPerCR,
00606                                &nodeList, &eqnList,
00607                                fieldIDs);
00608       if (err) voidERReturn;
00609       ++cr_iter;
00610    }
00611 }
00612 
00613 //==============================================================================
00614 int LinSysCoreFilter::storeNodalColumnCoefs(int eqn, const NodeDescriptor& node,
00615                                             int fieldID, int fieldSize,
00616                                             double* coefs)
00617 {
00618   //
00619   //This function stores the coeficients for 'node' at 'fieldID' at the correct
00620   //column indices in row 'eqn' of the system matrix.
00621   //
00622   if ((localStartRow_ > eqn) || (eqn > localEndRow_)) return(0);
00623 
00624   int eqnNumber = -1;
00625   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) ERReturn(FEI_ID_NOT_FOUND);
00626 
00627   int numParams = fieldSize;
00628   if (numParams < 1) {
00629     return(0);
00630   }
00631 
00632   if ((int)iworkSpace2_.size() < numParams) {
00633     iworkSpace2_.resize(numParams);
00634   }
00635 
00636   int* cols = &iworkSpace2_[0];
00637 
00638   for(int j=0; j<numParams; j++) {
00639     cols[j] = eqnNumber + j;
00640   }
00641 
00642   CHK_ERR( giveToMatrix(1, &eqn, numParams, cols, &coefs, ASSEMBLE_SUM) );
00643 
00644   return(FEI_SUCCESS);
00645 }
00646 
00647 //==============================================================================
00648 int LinSysCoreFilter::storeNodalRowCoefs(const NodeDescriptor& node,
00649                                          int fieldID, int fieldSize,
00650                                          double* coefs, int eqn)
00651 {
00652   //
00653   //This function stores coeficients in the equations for 'node', 'fieldID' at
00654   //column index 'eqn' of the system matrix.
00655   //
00656   int eqnNumber = -1;
00657   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) ERReturn(FEI_ID_NOT_FOUND);
00658 
00659   if ((localStartRow_ > eqnNumber) || (eqnNumber > localEndRow_)) return(0);
00660 
00661   int numParams = fieldSize;
00662 
00663   if (numParams < 1) {
00664     return(0);
00665   }
00666 
00667   if ((int)iworkSpace2_.size() < numParams) {
00668     iworkSpace2_.resize(numParams);
00669   }
00670 
00671   int* ptRows = &iworkSpace2_[0];
00672 
00673   if ((int)dworkSpace2_.size() < numParams) {
00674     dworkSpace2_.resize(numParams);
00675   }
00676 
00677   const double* * values = &dworkSpace2_[0];
00678 
00679   for(int j=0; j<numParams; j++) {
00680     ptRows[j] = eqnNumber + j;
00681     values[j] = &(coefs[j]);
00682   }
00683 
00684   CHK_ERR( giveToMatrix(numParams, ptRows, 1, &eqn, values, ASSEMBLE_SUM) );
00685 
00686   return(FEI_SUCCESS);
00687 }
00688 
00689 //==============================================================================
00690 void LinSysCoreFilter::storeNodalSendEqn(const NodeDescriptor& node,
00691                                          int fieldID, int col,
00692                                          double* coefs)
00693 {
00694   //
00695   //This is a private LinSysCoreFilter function. We can safely assume that
00696   //it will only be called with a node that is not locally owned.
00697   //
00698   //This function tells the eqn comm mgr that we'll be sending contributions
00699   //to column 'col' for the equations associated with 'fieldID', on 'node', on
00700   //node's owning processor.
00701   //
00702   int proc = node.getOwnerProc();
00703 
00704   int eqnNumber = -1;
00705   if (!node.getFieldEqnNumber(fieldID, eqnNumber)) voidERReturn;
00706 
00707   int numEqns = problemStructure_->getFieldSize(fieldID);
00708 
00709   for(int i=0; i<numEqns; i++) {
00710     eqnCommMgr_->addRemoteEqn(eqnNumber+i, proc, &coefs[i], &col, 1);
00711   }
00712 }
00713 
00714 //==============================================================================
00715 void LinSysCoreFilter::storePenNodeSendData(const NodeDescriptor& iNode,
00716                                             int iField, int iFieldSize,
00717                                             double* iCoefs,
00718                                             const NodeDescriptor& jNode,
00719                                             int jField, int jFieldSize,
00720                                             double* jCoefs,
00721                                             double penValue, double CRValue)
00722 {
00723 //
00724 //This function will register with the eqn comm mgr the equations associated
00725 //with iNode, field 'iField' having column indices that are the equations
00726 //associated with jNode, field 'jField', to be sent to the owner of iNode.
00727 //
00728    int proc = iNode.getOwnerProc();
00729 
00730    int iEqn = -1, jEqn = -1;
00731    if (!iNode.getFieldEqnNumber(iField, iEqn)) voidERReturn;
00732    if (!jNode.getFieldEqnNumber(jField, jEqn)) voidERReturn;
00733 
00734    int iNumParams = iFieldSize;
00735    int jNumParams = jFieldSize;
00736    if (iNumParams < 1 || jNumParams < 1) {
00737      fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
00738           << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
00739           << jNumParams<<FEI_ENDL;
00740      voidERReturn;
00741    }
00742 
00743    if ((int)dworkSpace_.size() < jNumParams) {
00744      dworkSpace_.resize(jNumParams);
00745    }
00746 
00747    double* coefs = &dworkSpace_[0];
00748 
00749    if ((int)iworkSpace2_.size() < jNumParams) {
00750      iworkSpace2_.resize(jNumParams);
00751    }
00752 
00753    int* cols = &iworkSpace2_[0];
00754 
00755    for(int i=0; i<iNumParams; i++) {
00756       for(int j=0; j<jNumParams; j++) {
00757          cols[j] = jEqn + j;
00758          coefs[j] = penValue*iCoefs[i]*jCoefs[j];
00759       }
00760 
00761       int row = iEqn + i;
00762       eqnCommMgr_->addRemoteEqn(row, proc, coefs, cols, jNumParams);
00763 
00764       double rhsValue = penValue*iCoefs[i]*CRValue;
00765       eqnCommMgr_->addRemoteRHS(row, proc, currentRHS_, rhsValue);
00766    }
00767 }
00768 
00769 //==============================================================================
00770 int LinSysCoreFilter::storePenNodeData(const NodeDescriptor& iNode,
00771                                        int iField, int iFieldSize,
00772                                        double* iCoefs,
00773                                        const NodeDescriptor& jNode,
00774                                        int jField, int jFieldSize,
00775                                        double* jCoefs,
00776                                        double penValue, double CRValue){
00777 //
00778 //This function will add to the local matrix the penalty constraint equations
00779 //associated with iNode at iField, having column indices that are the
00780 //equations associated with jNode at jField.
00781 //Also, add the penalty contribution to the RHS vector.
00782 //
00783    int iEqn = -1, jEqn = -1;
00784    if (!iNode.getFieldEqnNumber(iField, iEqn)) ERReturn(FEI_ID_NOT_FOUND);
00785    if (!jNode.getFieldEqnNumber(jField, jEqn)) ERReturn(FEI_ID_NOT_FOUND);
00786 
00787    int iNumParams = iFieldSize;
00788    int jNumParams = jFieldSize;
00789    if (iNumParams < 1 || jNumParams < 1) {
00790      fei::console_out() << "FEI ERROR, attempt to store indices for field with non-positive size"
00791           << " field "<<iField<<", size "<<iNumParams<<", field "<<jField<<", size "
00792           << jNumParams<<FEI_ENDL;
00793      ERReturn(-1);
00794    }
00795 
00796    if ((int)dworkSpace2_.size() < iNumParams) {
00797      dworkSpace2_.resize(iNumParams);
00798    }
00799 
00800    const double* * coefs = &dworkSpace2_[0];
00801 
00802    if ((int)dworkSpace_.size() < iNumParams * jNumParams) {
00803      dworkSpace_.resize(iNumParams * jNumParams);
00804    }
00805 
00806    double* coefList = &dworkSpace_[0];
00807 
00808    if ((int)iworkSpace2_.size() < jNumParams+iNumParams) {
00809      iworkSpace2_.resize(jNumParams+iNumParams);
00810    }
00811 
00812    int* cols = &iworkSpace2_[0];
00813    int* rows = &iworkSpace2_[jNumParams];
00814 
00815 
00816    for(int i=0; i<iNumParams; i++) {
00817       double* coefPtr = coefList + i*jNumParams;
00818       coefs[i] = coefPtr;
00819 
00820       rows[i] = iEqn + i;
00821 
00822       for(int j=0; j<jNumParams; j++) {
00823          if (i==0) cols[j] = jEqn + j;
00824          coefPtr[j] = penValue*iCoefs[i]*jCoefs[j];
00825       }
00826 
00827       double rhsValue = penValue*iCoefs[i]*CRValue;
00828       CHK_ERR( giveToRHS(1, &rhsValue, &rows[i], ASSEMBLE_SUM) );
00829    }
00830 
00831    CHK_ERR( giveToMatrix(iNumParams, rows,
00832                          jNumParams, cols, coefs, ASSEMBLE_SUM) );
00833 
00834    return(FEI_SUCCESS);
00835 }
00836 
00837 //------------------------------------------------------------------------------
00838 int LinSysCoreFilter::resetSystem(double s)
00839 {
00840   //
00841   //  This puts the value s throughout both the matrix and the vector.
00842   //
00843   if (Filter::logStream() != NULL) {
00844     (*logStream()) << "FEI: resetSystem" << FEI_ENDL << s << FEI_ENDL;
00845   }
00846 
00847   CHK_ERR( resetTheMatrix(s) );
00848   CHK_ERR( resetTheRHSVector(s) );
00849 
00850   debugOutput("#LinSysCoreFilter  leaving resetSystem");
00851 
00852   return(FEI_SUCCESS);
00853 }
00854 
00855 //------------------------------------------------------------------------------
00856 int LinSysCoreFilter::deleteMultCRs()
00857 {
00858   debugOutput("#LinSysCoreFilter::deleteMultCRs");
00859 
00860   LinSysCore_flexible* lscf = dynamic_cast<LinSysCore_flexible*>(lsc_);
00861   if (lscf == NULL) {
00862 //    fei::console_out() << "FEI::LinSysCoreFilter: ERROR deleteMultCRs requested, but "
00863 //         << "underlying solver doesn't support this operation." << FEI_ENDL;
00864     return(-1);
00865   }
00866 
00867   int err = lscf->resetConstraints(0.0);
00868 
00869   debugOutput("#LinSysCoreFilter leaving deleteMultCRs");
00870 
00871   return(err);
00872 }
00873 
00874 //------------------------------------------------------------------------------
00875 int LinSysCoreFilter::resetTheMatrix(double s)
00876 {
00877   CHK_ERR( lsc_->resetMatrix(s) );
00878 
00879   return(FEI_SUCCESS);
00880 }
00881 
00882 //------------------------------------------------------------------------------
00883 int LinSysCoreFilter::resetTheRHSVector(double s)
00884 {
00885   CHK_ERR( lsc_->resetRHSVector(s) );
00886 
00887   return(FEI_SUCCESS);
00888 }
00889 
00890 //------------------------------------------------------------------------------
00891 int LinSysCoreFilter::resetMatrix(double s) {
00892 //
00893 //  This puts the value s throughout both the matrix and the vector.
00894 //
00895 
00896     debugOutput("FEI: resetMatrix");
00897 
00898     CHK_ERR( resetTheMatrix(s) )
00899 
00900    //clear away any boundary condition data.
00901    bcManager_->clearAllBCs();
00902 
00903     eqnCommMgr_->resetCoefs();
00904 
00905     debugOutput("#LinSysCoreFilter leaving resetMatrix");
00906 
00907     return(FEI_SUCCESS);
00908 }
00909 
00910 //------------------------------------------------------------------------------
00911 int LinSysCoreFilter::resetRHSVector(double s) {
00912 //
00913 //  This puts the value s throughout the rhs vector.
00914 //
00915 
00916     debugOutput("FEI: resetRHSVector");
00917 
00918     CHK_ERR( resetTheRHSVector(s) )
00919 
00920     //clear away any boundary condition data.
00921     bcManager_->clearAllBCs();
00922 
00923     eqnCommMgr_->resetCoefs();
00924 
00925     debugOutput("# LinSysCoreFilter leaving resetRHSVector");
00926 
00927     return(FEI_SUCCESS);
00928 }
00929 
00930 //------------------------------------------------------------------------------
00931 int LinSysCoreFilter::resetInitialGuess(double s) {
00932 //
00933 //  This puts the value s throughout the initial guess (solution) vector.
00934 //
00935   if (Filter::logStream() != NULL) {
00936     FEI_OSTREAM& os = *logStream();
00937     os << "FEI: resetInitialGuess" << FEI_ENDL;
00938     os << "#value to which initial guess is to be set" << FEI_ENDL;
00939     os << s << FEI_ENDL;
00940   }
00941 
00942   int* eqns = new int[numReducedRows_];
00943   double* values = new double[numReducedRows_];
00944   if (eqns == NULL || values == NULL) return(-1);
00945 
00946   for(int i=0; i<numReducedRows_; ++i) {
00947     eqns[i] = reducedStartRow_ + i;
00948     values[i] = s;
00949   }
00950 
00951   CHK_ERR( lsc_->putInitialGuess(eqns, values, numReducedRows_) );
00952 
00953   delete [] eqns;
00954   delete [] values;
00955 
00956   debugOutput("# LinSysCoreFilter leaving resetInitialGuess");
00957 
00958   return(FEI_SUCCESS);
00959 }
00960 
00961 //------------------------------------------------------------------------------
00962 int LinSysCoreFilter::loadNodeBCs(int numNodes,
00963                                   const GlobalID *nodeIDs,
00964                                   int fieldID,
00965                                   const int* offsetsIntoField,
00966                                   const double* prescribedValues)
00967 {
00968   //
00969   //  load boundary condition information for a given set of nodes
00970   //
00971   int size = problemStructure_->getFieldSize(fieldID);
00972   if (size < 1) {
00973     fei::console_out() << "FEI Warning: loadNodeBCs called for fieldID "<<fieldID
00974          <<", which was defined with size "<<size<<" (should be positive)."<<FEI_ENDL;
00975     return(0);
00976   }
00977 
00978   if (Filter::logStream() != NULL) {
00979     (*logStream())<<"FEI: loadNodeBCs"<<FEI_ENDL
00980                      <<"#num-nodes"<<FEI_ENDL<<numNodes<<FEI_ENDL
00981                      <<"#fieldID"<<FEI_ENDL<<fieldID<<FEI_ENDL
00982                      <<"#field-size"<<FEI_ENDL<<size<<FEI_ENDL;
00983     (*logStream())<<"#following lines: nodeID offsetIntoField value "<<FEI_ENDL;
00984 
00985     for(int j=0; j<numNodes; j++) {
00986       GlobalID nodeID = nodeIDs[j];
00987       (*logStream())<<static_cast<int>(nodeID)<<"  ";
00988       (*logStream())<<offsetsIntoField[j]<<"  "<<prescribedValues[j];
00989       (*logStream())<<FEI_ENDL;
00990     }
00991   }
00992 
00993   if (numNodes > 0) newBCData_ = true;
00994 
00995   bcManager_->addBCRecords(numNodes, nodeIDType_, fieldID, nodeIDs,
00996                            offsetsIntoField, prescribedValues);
00997 
00998   return(FEI_SUCCESS);
00999 }
01000 
01001 //------------------------------------------------------------------------------
01002 int LinSysCoreFilter::loadElemBCs(int numElems,
01003                                   const GlobalID *elemIDs,
01004                                   int fieldID,
01005                                   const double *const *alpha,
01006                                   const double *const *beta,
01007                                   const double *const *gamma)
01008 {
01009    return(-1);
01010 }
01011 
01012 //------------------------------------------------------------------------------
01013 void LinSysCoreFilter::allocElemStuff()
01014 {
01015    int nb = problemStructure_->getNumElemBlocks();
01016 
01017    for(int i=0; i<nb; i++) {
01018      BlockDescriptor* block = NULL;
01019      int err = problemStructure_->getBlockDescriptor_index(i, block);
01020      if (err) voidERReturn;
01021 
01022       int numEqns = block->getNumEqnsPerElement();
01023       if (maxElemRows_ < numEqns) maxElemRows_ = numEqns;
01024    }
01025 
01026    scatterIndices_.resize(maxElemRows_);
01027 
01028    if (eStiff_ != NULL) delete [] eStiff_;
01029    if (eStiff1D_ != NULL) delete [] eStiff1D_;
01030    if (eLoad_ != NULL) delete [] eLoad_;
01031 
01032    eStiff_ = new double*[maxElemRows_];
01033    eStiff1D_ = new double[maxElemRows_*maxElemRows_];
01034 
01035    if (eStiff_ == NULL || eStiff1D_ == NULL) voidERReturn
01036 
01037    for(int r=0; r<maxElemRows_; r++) {
01038       eStiff_[r] = eStiff1D_ + r*maxElemRows_;
01039    }
01040 
01041    eLoad_ = new double[maxElemRows_];
01042 
01043    if (eLoad_ == NULL) voidERReturn
01044 }
01045 
01046 //------------------------------------------------------------------------------
01047 int LinSysCoreFilter::sumInElem(GlobalID elemBlockID,
01048                                 GlobalID elemID,
01049                                 const GlobalID* elemConn,
01050                                 const double* const* elemStiffness,
01051                                 const double* elemLoad,
01052                                 int elemFormat)
01053 {
01054   if (Filter::logStream() != NULL && outputLevel_ > 2) {
01055     (*logStream()) << "FEI: sumInElem" << FEI_ENDL <<"#blkID" << FEI_ENDL
01056                       << static_cast<int>(elemBlockID) << FEI_ENDL
01057                       << "#elID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
01058     BlockDescriptor* block = NULL;
01059     CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01060     int numNodes = block->getNumNodesPerElement();
01061     (*logStream()) << "#n-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
01062     (*logStream()) << "#nodes" << FEI_ENDL;
01063     for(int i=0; i<numNodes; ++i) {
01064       GlobalID nodeID = elemConn[i];
01065       (*logStream())<<static_cast<int>(nodeID)<<" ";
01066     }
01067     (*logStream())<<FEI_ENDL;
01068   }
01069 
01070   return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
01071                           elemLoad, elemFormat));
01072 }
01073 
01074 //------------------------------------------------------------------------------
01075 int LinSysCoreFilter::sumInElemMatrix(GlobalID elemBlockID,
01076                                       GlobalID elemID,
01077                                       const GlobalID* elemConn,
01078                                       const double* const* elemStiffness,
01079                                       int elemFormat)
01080 {
01081   if (Filter::logStream() != NULL && outputLevel_ > 2) {
01082     (*logStream()) << "FEI: sumInElemMatrix"<<FEI_ENDL
01083                       << "#blkID" << FEI_ENDL << static_cast<int>(elemBlockID) << FEI_ENDL
01084                       << "#elID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
01085     BlockDescriptor* block = NULL;
01086     CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01087     int numNodes = block->getNumNodesPerElement();
01088     (*logStream()) << "#n-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
01089     (*logStream()) << "#nodes" << FEI_ENDL;
01090     for(int i=0; i<numNodes; ++i) {
01091       GlobalID nodeID = elemConn[i];
01092       (*logStream())<<static_cast<int>(nodeID)<<" ";
01093     }
01094     (*logStream())<<FEI_ENDL;
01095   }
01096 
01097   return(generalElemInput(elemBlockID, elemID, elemConn, elemStiffness,
01098                           NULL, elemFormat));
01099 }
01100 
01101 //------------------------------------------------------------------------------
01102 int LinSysCoreFilter::sumInElemRHS(GlobalID elemBlockID,
01103                            GlobalID elemID,
01104                            const GlobalID* elemConn,
01105                            const double* elemLoad)
01106 {
01107   if (Filter::logStream() != NULL && outputLevel_ > 2) {
01108     (*logStream()) << "FEI: sumInElemRHS"<<FEI_ENDL<<"#blID" << FEI_ENDL
01109                       <<(int)elemBlockID << FEI_ENDL
01110                       << "#elID" << FEI_ENDL << static_cast<int>(elemID) << FEI_ENDL;
01111     BlockDescriptor* block = NULL;
01112     CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01113     int numNodes = block->getNumNodesPerElement();
01114     (*logStream()) << "#n-nodes" << FEI_ENDL << numNodes << FEI_ENDL;
01115     (*logStream()) << "#nodes" << FEI_ENDL;
01116     for(int i=0; i<numNodes; ++i) {
01117       GlobalID nodeID = elemConn[i];
01118       (*logStream())<<static_cast<int>(nodeID)<<" ";
01119     }
01120     (*logStream())<<FEI_ENDL;
01121   }
01122 
01123   return(generalElemInput(elemBlockID, elemID, elemConn, NULL,
01124                           elemLoad, -1));
01125 }
01126 
01127 //------------------------------------------------------------------------------
01128 int LinSysCoreFilter::generalElemInput(GlobalID elemBlockID,
01129                                        GlobalID elemID,
01130                                        const GlobalID* elemConn,
01131                                        const double* const* elemStiffness,
01132                                        const double* elemLoad,
01133                                        int elemFormat)
01134 {
01135   (void)elemConn;
01136   return(generalElemInput(elemBlockID, elemID, elemStiffness, elemLoad,
01137                           elemFormat) );
01138 }
01139 
01140 //------------------------------------------------------------------------------
01141 int LinSysCoreFilter::generalElemInput(GlobalID elemBlockID,
01142                                        GlobalID elemID,
01143                                        const double* const* elemStiffness,
01144                                        const double* elemLoad,
01145                                        int elemFormat)
01146 {
01147   //first get the block-descriptor for this elemBlockID...
01148 
01149   BlockDescriptor* block = NULL;
01150   CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
01151 
01152   //now allocate our local stiffness/load copy if we haven't already.
01153 
01154   if (maxElemRows_ <= 0) allocElemStuff();
01155 
01156   int numElemRows = block->getNumEqnsPerElement();
01157   int numBlkElemRows = block->getNumBlkEqnsPerElement();
01158   int interleave = block->getInterleaveStrategy();
01159 
01160   //an std::vector.resize operation is free if the size is either shrinking or
01161   //staying the same.
01162 
01163   scatterIndices_.resize(numElemRows);
01164   blkScatterIndices_.resize(numBlkElemRows*2);
01165 
01166   const double* const* stiff = NULL;
01167   if (elemStiffness != NULL) stiff = elemStiffness;
01168 
01169   const double* load = NULL;
01170   if (elemLoad != NULL) load = elemLoad;
01171 
01172   //we'll make a local dense copy of the element stiffness array
01173   //if the stiffness array was passed in using elemFormat != FEI_DENSE_ROW
01174   //AND if the stiffness array is non-null.
01175   //Note that this is not optimal if elemFormat is one of the
01176   //diagonal or block-diagonal formats.
01177   if (elemFormat != FEI_DENSE_ROW && stiff != NULL) {
01178     if (elemFormat == FEI_BLOCK_DIAGONAL_ROW ||
01179         elemFormat == FEI_BLOCK_DIAGONAL_COL) {
01180       fei::console_out() << "LinSysCoreFilter::generalElemInput ERROR, elemFormat="
01181                << elemFormat << " not supported."<<FEI_ENDL;
01182       ERReturn(-1);
01183     }
01184 
01185     Filter::copyStiffness(stiff, numElemRows, elemFormat, eStiff_);
01186     stiff = eStiff_;
01187   }
01188 
01189   if (stiff != NULL) newMatrixData_ = true;
01190   if (load != NULL) newVectorData_ = true;
01191 
01192   if (Filter::logStream() != NULL && outputLevel_ > 2) {
01193     FEI_OSTREAM& os = *logStream();
01194 
01195     if (stiff != NULL || load != NULL) {
01196       os << "#numRows"<< FEI_ENDL << numElemRows << FEI_ENDL;
01197     }
01198 
01199     if (stiff != NULL) {
01200       os << "#elem-stiff (after copy into dense-row format)" << FEI_ENDL;
01201       for(int i=0; i<numElemRows; i++) {
01202         const double* stiff_i = stiff[i];
01203         for(int j=0; j<numElemRows; j++) {
01204           os << stiff_i[j] << " ";
01205         }
01206         os << FEI_ENDL;
01207       }
01208     }
01209 
01210     if (load != NULL) {
01211       os << "#elem-load" << FEI_ENDL;
01212       for(int i=0; i<numElemRows; i++) {
01213         os << load[i] << " ";
01214       }
01215       os<<FEI_ENDL;
01216     }
01217 
01218     if (stiff != NULL) {
01219       os << "#elemformat" << FEI_ENDL << elemFormat << FEI_ENDL;
01220     }
01221   }
01222 
01223   //now let's obtain the scatter indices for assembling the equations
01224   //into their appropriate places in the global matrix and rhs vectors
01225 
01226   int* indPtr = &scatterIndices_[0];
01227   int* blkIndPtr = &blkScatterIndices_[0];
01228   int* blkSizesPtr = blkIndPtr + numBlkElemRows;
01229 
01230   bool useBlkEqns = false;
01231   if (interleave == 0) {
01232     //interleave==0 is node-major, so we'll get the block-indices too.
01233     problemStructure_->getScatterIndices_ID(elemBlockID, elemID,
01234                                             interleave, indPtr,
01235                                             blkIndPtr, blkSizesPtr);
01236     int sumBlkSizes = 0;
01237     for(int ii=0; ii<numBlkElemRows; ++ii) {
01238       sumBlkSizes += blkSizesPtr[ii];
01239     }
01240     if (sumBlkSizes == numElemRows) {
01241       useBlkEqns = true;
01242     }
01243   }
01244   else {
01245     //interleave!=0 is field-major, and we'll only bother with point-indices.
01246     problemStructure_->getScatterIndices_ID(elemBlockID, elemID,
01247                                             interleave, indPtr);
01248   }
01249 
01250   if (stiff != NULL) {
01251     if (problemStructure_->numSlaveEquations() == 0) {
01252       //I'm not checking the return-value (error-code) on this call, because
01253       //I wasn't even calling it until recently, and I'm not sure if all
01254       //LinearSystemCore implementations even have it implemented.
01255       lsc_->setStiffnessMatrices(elemBlockID, 1, &elemID,
01256                                  &stiff, scatterIndices_.size(),
01257                                  &indPtr);
01258     }
01259 
01260     int len = scatterIndices_.size();
01261     if (interleave == 0) {
01262       CHK_ERR( assembleEqns( len, len, indPtr, indPtr, stiff, true,
01263                              numBlkElemRows, blkIndPtr,
01264                              blkSizesPtr, useBlkEqns, ASSEMBLE_SUM ) );
01265     }
01266     else {
01267       CHK_ERR( assembleEqns( len, len, indPtr, indPtr, stiff, true,
01268                              numBlkElemRows, blkIndPtr,
01269                              blkSizesPtr, false, ASSEMBLE_SUM ) );
01270     }
01271   }
01272 
01273   if (load != NULL) {
01274     if (problemStructure_->numSlaveEquations() == 0) {
01275       //I'm not checking the return-value (error-code) on this call, because
01276       //I wasn't even calling it until recently, and I'm not sure if all
01277       //LinearSystemCore implementations even have it implemented.
01278       lsc_->setLoadVectors(elemBlockID, 1, &elemID,
01279                            &load, scatterIndices_.size(),
01280                            &indPtr);
01281     }
01282 
01283     CHK_ERR( assembleRHS(scatterIndices_.size(), indPtr, load, ASSEMBLE_SUM) );
01284   }
01285 
01286   return(FEI_SUCCESS);
01287 }
01288 
01289 //------------------------------------------------------------------------------
01290 int LinSysCoreFilter::putIntoRHS(int IDType,
01291                           int fieldID,
01292                           int numIDs,
01293                           const GlobalID* IDs,
01294                           const double* rhsEntries)
01295 {
01296   int fieldSize = problemStructure_->getFieldSize(fieldID);
01297 
01298   rowIndices_.resize(fieldSize*numIDs);
01299   int checkNumEqns;
01300 
01301   CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
01302                                             checkNumEqns,
01303                                             &rowIndices_[0]));
01304   if (checkNumEqns != numIDs*fieldSize) {
01305     ERReturn(-1);
01306   }
01307 
01308   CHK_ERR( exchangeRemoteEquations() );
01309 
01310   CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_PUT));
01311 
01312   return(0);
01313 }
01314 
01315 //------------------------------------------------------------------------------
01316 int LinSysCoreFilter::sumIntoRHS(int IDType,
01317                           int fieldID,
01318                           int numIDs,
01319                           const GlobalID* IDs,
01320                           const double* rhsEntries)
01321 {
01322   int fieldSize = problemStructure_->getFieldSize(fieldID);
01323 
01324   rowIndices_.resize(fieldSize*numIDs);
01325   int checkNumEqns;
01326 
01327   CHK_ERR( problemStructure_->getEqnNumbers(numIDs, IDs, IDType, fieldID,
01328                                             checkNumEqns,
01329                                             &rowIndices_[0]));
01330   if (checkNumEqns != numIDs*fieldSize) {
01331     ERReturn(-1);
01332   }
01333 
01334   CHK_ERR(assembleRHS(rowIndices_.size(), &rowIndices_[0], rhsEntries, ASSEMBLE_SUM));
01335 
01336   return(0);
01337 }
01338 
01339 //------------------------------------------------------------------------------
01340 int LinSysCoreFilter::implementAllBCs() {
01341 //
01342 // This function will handle the modifications to the stiffness and load
01343 // necessary to enforce nodal boundary conditions.
01344 //
01345    debugOutput("# implementAllBCs");
01346 
01347    std::vector<int> essEqns;
01348    std::vector<double> essGamma;
01349 
01350    EqnBuffer bcEqns;
01351 
01352    CHK_ERR( bcManager_->finalizeBCEqns(bcEqns) );
01353 
01354    if (resolveConflictRequested_) {
01355      CHK_ERR( resolveConflictingCRs(bcEqns) );
01356    }
01357 
01358    CHK_ERR( eqnCommMgr_->gatherSharedBCs(bcEqns) );
01359 
01360    //now separate the boundary-condition equations into arrays
01361    fei::FillableMat bcEqns_mat(bcEqns);
01362    fei::impl_utils::separate_BC_eqns(bcEqns_mat, essEqns, essGamma);
01363 
01364    std::vector<double> essAlpha(essEqns.size(), 1);
01365 
01366    exchangeRemoteBCs(essEqns, essAlpha, essGamma);
01367 
01368    if (essEqns.size() > 0) {
01369       CHK_ERR( enforceEssentialBCs(&essEqns[0],
01370                                    &essAlpha[0],
01371                                    &essGamma[0], essEqns.size()) );
01372    }
01373 
01374    debugOutput("#LinSysCoreFilter leaving implementAllBCs");
01375    return(FEI_SUCCESS);
01376 }
01377 
01378 //------------------------------------------------------------------------------
01379 int LinSysCoreFilter::enforceEssentialBCs(const int* eqns,
01380                                           const double* alpha,
01381                                           const double* gamma,
01382                                           int numEqns)
01383 {
01384   int* cc_eqns = const_cast<int*>(eqns);
01385   double* cc_alpha = const_cast<double*>(alpha);
01386   double* cc_gamma = const_cast<double*>(gamma);
01387 
01388   if (problemStructure_->numSlaveEquations() == 0) {
01389     CHK_ERR( lsc_->enforceEssentialBC(cc_eqns,
01390                                       cc_alpha, cc_gamma,
01391                                       numEqns) );
01392   }
01393   else {
01394     std::vector<int> reducedEqns(numEqns);
01395     for(int i=0; i<numEqns; i++) {
01396       problemStructure_->translateToReducedEqn(eqns[i], reducedEqns[i]);
01397     }
01398 
01399     CHK_ERR( lsc_->enforceEssentialBC(&reducedEqns[0], cc_alpha, cc_gamma, numEqns) );
01400   }
01401 
01402   return(FEI_SUCCESS);
01403 }
01404 
01405 //------------------------------------------------------------------------------
01406 int LinSysCoreFilter::enforceRemoteEssBCs(int numEqns, const int* eqns, 
01407                            const int* const* colIndices, const int* colIndLens,
01408                            const double* const* BCcoefs)
01409 {
01410   //These eqn-numbers were reduced to the slave-eqn space by
01411   //LinSysCore::exchangeRemoteBCs, which is the only function that calls THIS
01412   //function, so we can simply pass them straight on in to LinearSystemCore.
01413   //
01414 
01415   int* cc_eqns = const_cast<int*>(eqns);
01416   int** cc_colIndices = const_cast<int**>(colIndices);
01417   int* cc_colIndLens = const_cast<int*>(colIndLens);
01418   double** cc_BCcoefs = const_cast<double**>(BCcoefs);
01419 
01420   CHK_ERR( lsc_->enforceRemoteEssBCs(numEqns, cc_eqns, cc_colIndices,
01421                                      cc_colIndLens, cc_BCcoefs) );
01422 
01423   return(FEI_SUCCESS);
01424 }
01425 
01426 //------------------------------------------------------------------------------
01427 int LinSysCoreFilter::resolveConflictingCRs(EqnBuffer& bcEqns)
01428 {
01429   int numMultCRs = problemStructure_->getNumMultConstRecords();
01430   if (numMultCRs < 1) {
01431     return(0);
01432   }
01433 
01434   std::map<GlobalID,ConstraintType*>::const_iterator
01435     cr_iter = problemStructure_->getMultConstRecords().begin(),
01436     cr_end  = problemStructure_->getMultConstRecords().end();
01437 
01438   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
01439 
01440   std::vector<int>& bcEqnNumbers = bcEqns.eqnNumbers();
01441 
01442   double coefs[3];
01443   int indices[3];
01444   indices[0] = 0;
01445   indices[1] = 1;
01446   indices[2] = 2;
01447 
01448   double fei_eps = 1.e-49;
01449 
01450   while(cr_iter != cr_end) {
01451     ConstraintType& multCR = *((*cr_iter).second);
01452 
01453     int lenList = multCR.getMasters().size();
01454 
01455     std::vector<GlobalID>& CRNode_vec = multCR.getMasters();
01456     GlobalID *CRNodePtr = &CRNode_vec[0];
01457     std::vector<int>& CRField_vec = multCR.getMasterFieldIDs();
01458     int* CRFieldPtr = &CRField_vec[0];
01459     std::vector<double>& weights_vec = multCR.getMasterWeights();
01460     double* weights = &weights_vec[0];
01461 
01462     int offset = 0;
01463     for(int j=0; j<lenList; ++j) {
01464       int fieldSize = problemStructure_->getFieldSize(CRFieldPtr[j]);
01465       for(int k=0; k<fieldSize; ++k) {
01466         if (std::abs(weights[offset++] + 1.0) < fei_eps) {
01467           const NodeDescriptor* node = NULL;
01468           CHK_ERR( nodeDB.getNodeWithID(CRNodePtr[j], node) );
01469           int eqn = 0;
01470           node->getFieldEqnNumber(CRFieldPtr[j], eqn);
01471           eqn += k;
01472 
01473           if (fei::binarySearch(eqn, bcEqnNumbers) >= 0) {
01474             coefs[0] = 1.0;
01475             coefs[1] = 0.0;
01476             coefs[2] = 1.0;
01477             int crEqn = multCR.getEqnNumber();
01478             CHK_ERR( bcEqns.addEqn(crEqn, coefs, indices, 3, false) );
01479           }
01480         }
01481       }
01482     }
01483     ++cr_iter;
01484   }
01485 
01486   return(0);
01487 }
01488 
01489 //------------------------------------------------------------------------------
01490 int LinSysCoreFilter::exchangeRemoteEquations()
01491 {
01492   //
01493   // This function is where processors send local contributions to remote
01494   // equations to the owners of those equations, and receive remote
01495   // contributions to local equations.
01496   //
01497   debugOutput("#LinSysCoreFilter::exchangeRemoteEquations");
01498 
01499   if (reducedEqnCounter_ > 0) CHK_ERR( assembleReducedEqns() );
01500 
01501   if (reducedRHSCounter_ > 0) CHK_ERR( assembleReducedRHS() );
01502 
01503   int len = 4;
01504   std::vector<int> flags(len), globalFlags(len);
01505   flags[0] = newMatrixData_ ? 1 : 0;
01506   flags[1] = newVectorData_ ? 1 : 0;
01507   flags[2] = newConstraintData_ ? 1 : 0;
01508   flags[3] = newBCData_ ? 1 : 0;
01509 
01510   CHK_ERR( fei::GlobalMax(comm_, flags, globalFlags) );
01511 
01512   newMatrixData_     = globalFlags[0] > 0 ? true : false;
01513   newVectorData_     = globalFlags[1] > 0 ? true : false;
01514   newConstraintData_ = globalFlags[2] > 0 ? true : false;
01515   newBCData_         = globalFlags[3] > 0 ? true : false;
01516 
01517   if (newMatrixData_ || newVectorData_ || newConstraintData_) { 
01518 
01519     CHK_ERR( eqnCommMgr_->exchangeEqns(logStream()) );
01520 
01521     needToCallMatrixLoadComplete_ = true;
01522 
01523     //so now the remote contributions should be available, let's get them out
01524     //of the eqn comm mgr and put them into our local matrix structure.
01525 
01526     debugOutput("#   putting remote contributions into linear system...");
01527 
01528     CHK_ERR( unpackRemoteContributions(*eqnCommMgr_, ASSEMBLE_SUM) );
01529 
01530     eqnCommMgr_->resetCoefs();
01531 
01532     newMatrixData_ = false;
01533     newVectorData_ = false;
01534     newConstraintData_ = false;
01535   }
01536 
01537   firstRemEqnExchange_ = false;
01538 
01539   debugOutput("#LinSysCoreFilter leaving exchangeRemoteEquations");
01540 
01541   return(FEI_SUCCESS);
01542 }
01543 
01544 //------------------------------------------------------------------------------
01545 int LinSysCoreFilter::unpackRemoteContributions(EqnCommMgr& eqnCommMgr,
01546                                                 int assemblyMode)
01547 {
01548   int numRecvEqns = eqnCommMgr.getNumLocalEqns();
01549   std::vector<int>& recvEqnNumbers = eqnCommMgr.localEqnNumbers();
01550   std::vector<fei::CSVec*>& recvEqns = eqnCommMgr.localEqns();
01551   std::vector<std::vector<double>*>& recvRHSs = *(eqnCommMgr.localRHSsPtr());
01552 
01553   bool newCoefs = eqnCommMgr.newCoefData();
01554   bool newRHSs = eqnCommMgr.newRHSData();
01555 
01556   int i;
01557   double** coefs = new double*[numRecvEqns];
01558 
01559   for(i=0; i<numRecvEqns; i++) {
01560     coefs[i] = &(recvEqns[i]->coefs()[0]);
01561   }
01562 
01563   for(i=0; i<numRecvEqns; i++) {
01564 
01565     int eqn = recvEqnNumbers[i];
01566     if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
01567       fei::console_out() << "LinSysCoreFilter::unpackRemoteContributions: ERROR, recvEqn "
01568            << eqn << " out of range. (localStartRow_: " << reducedStartRow_
01569            << ", localEndRow_: " << reducedEndRow_ << ", localRank_: "
01570            << localRank_ << ")" << FEI_ENDL;
01571       MPI_Abort(comm_, -1);
01572     }
01573 
01574     for(size_t ii=0; ii<recvEqns[i]->size(); ii++) {
01575       if (coefs[i][ii] > 1.e+200) {
01576         fei::console_out() << localRank_ << ": LinSysCoreFilter::unpackRemoteContributions: "
01577              << "WARNING, coefs["<<i<<"]["<<ii<<"]: " << coefs[i][ii]
01578              << FEI_ENDL;
01579         MPI_Abort(comm_, -1);
01580       }
01581     }
01582 
01583     if (recvEqns[i]->size() > 0 && newCoefs) {
01584       //contribute this equation to the matrix,
01585       CHK_ERR( giveToLocalReducedMatrix(1, &(recvEqnNumbers[i]),
01586                                         recvEqns[i]->size(),
01587                                         &(recvEqns[i]->indices()[0]),
01588                                         &(coefs[i]), assemblyMode ) );
01589     }
01590 
01591     //and now the RHS contributions.
01592     if (newRHSs) {
01593       for(int j=0; j<numRHSs_; j++) {
01594         CHK_ERR( giveToLocalReducedRHS(1, &( (*(recvRHSs[i]))[j] ),
01595                                        &eqn, assemblyMode) );
01596       }
01597     }
01598   }
01599 
01600   delete [] coefs;
01601 
01602   return(0);
01603 }
01604 
01605 //------------------------------------------------------------------------------
01606 int LinSysCoreFilter::exchangeRemoteBCs(std::vector<int>& essEqns,
01607                                         std::vector<double>& essAlpha,
01608                                         std::vector<double>& essGamma)
01609 {
01610   //we need to make sure that the right thing happens for essential
01611   //boundary conditions that get applied to nodes on elements that touch
01612   //a processor boundary. (Note that for this case, the BC node itself doesn't
01613   //touch the processor boundary.) For an essential boundary condition, the row
01614   //and column of the corresponding equation must be diagonalized. If there is
01615   //a processor boundary on any side of the element that contains the node,
01616   //then there are column contributions to the matrix on the other processor.
01617   //That other processor must be notified and told to make the adjustments
01618   //necessary to enforce the boundary condition.
01619 
01620   std::vector<int>* eqns = &essEqns;
01621 
01622   if (problemStructure_->numSlaveEquations() > 0) {
01623     int numEqns = essEqns.size();
01624     eqns = new std::vector<int>(numEqns);
01625     int* eqnsPtr = &(*eqns)[0];
01626 
01627     for(int ii=0; ii<numEqns; ++ii) {
01628       problemStructure_->translateToReducedEqn(essEqns[ii], eqnsPtr[ii]);
01629     }
01630   }
01631 
01632   FEI_OSTREAM* dbgOut = NULL;
01633   if (Filter::logStream() != NULL) {
01634     dbgOut = logStream();
01635   }
01636 
01637   eqnCommMgr_->exchangeRemEssBCs(&(*eqns)[0], eqns->size(),
01638                                  &essAlpha[0], &essGamma[0],
01639                                  comm_, dbgOut);
01640 
01641   int numRemoteEssBCEqns = eqnCommMgr_->getNumRemEssBCEqns();
01642   if (numRemoteEssBCEqns > 0) {
01643     std::vector<int>& remEssBCEqnNumbers = eqnCommMgr_->remEssBCEqnNumbersPtr();
01644     fei::CSVec** remEssBCEqns = &(eqnCommMgr_->remEssBCEqns()[0]);
01645     std::vector<int> remEssBCEqnLengths(remEssBCEqnNumbers.size());
01646 
01647     int** indices = new int*[numRemoteEssBCEqns];
01648     double** coefs = new double*[numRemoteEssBCEqns];
01649 
01650     for(int i=0; i<numRemoteEssBCEqns; i++) {
01651       coefs[i] = &(remEssBCEqns[i]->coefs()[0]);
01652       indices[i] = &(remEssBCEqns[i]->indices()[0]);
01653       remEssBCEqnLengths[i] = remEssBCEqns[i]->size();
01654     }
01655 
01656     CHK_ERR( enforceRemoteEssBCs(numRemoteEssBCEqns,
01657                                  &remEssBCEqnNumbers[0], indices,
01658                                  &remEssBCEqnLengths[0], coefs));
01659 
01660     delete [] indices;
01661     delete [] coefs;
01662   }
01663 
01664   if (problemStructure_->numSlaveEquations() > 0) {
01665     delete eqns;
01666   }
01667 
01668   debugOutput("#LinSysCoreFilter exchangeRemoteBCs");
01669 
01670   return(FEI_SUCCESS);
01671 }
01672 
01673 //------------------------------------------------------------------------------
01674 int LinSysCoreFilter::loadCRMult(int CRID,
01675                                  int numCRNodes,
01676                                  const GlobalID* CRNodes,
01677                                  const int* CRFields,
01678                                  const double* CRWeights,
01679                                  double CRValue)
01680 {
01681 //
01682 // Load Lagrange multiplier constraint relation data
01683 //
01684 //   Question: do we really need to pass CRNodes again?  Here, I'm going
01685 //            to ignore it for now (i.e., not store it, but just check it), 
01686 //            as it got passed during the initialization phase, so all we'll 
01687 //            do here is check for errors...
01688 //
01689   if (Filter::logStream() != NULL) {
01690     FEI_OSTREAM& os = *logStream();
01691     os<<"FEI: loadCRMult"<<FEI_ENDL;
01692     os<<"#num-nodes"<<FEI_ENDL<<numCRNodes<<FEI_ENDL;
01693     os<<"#CRNodes:"<<FEI_ENDL;
01694     int i;
01695     for(i=0; i<numCRNodes; ++i) {
01696       GlobalID nodeID = CRNodes[i];
01697       os << static_cast<int>(nodeID) << " ";
01698     }
01699     os << FEI_ENDL << "#fields:"<<FEI_ENDL;
01700     for(i=0; i<numCRNodes; ++i) os << CRFields[i] << " ";
01701     os << FEI_ENDL << "#field-sizes:"<<FEI_ENDL;
01702     for(i=0; i<numCRNodes; ++i) {
01703       int size = problemStructure_->getFieldSize(CRFields[i]);
01704       os << size << " ";
01705     }
01706     os << FEI_ENDL<<"#weights:"<<FEI_ENDL;
01707     int offset = 0;
01708     for(i=0; i<numCRNodes; ++i) {
01709       int size = problemStructure_->getFieldSize(CRFields[i]);
01710       for(int j=0; j<size; ++j) {
01711         os << CRWeights[offset++] << " ";
01712       }
01713     }
01714     os << FEI_ENDL<<"#CRValue:"<<FEI_ENDL<<CRValue<<FEI_ENDL;
01715   }
01716 
01717   ConstraintType* multCR = NULL;
01718   CHK_ERR( problemStructure_->getMultConstRecord(CRID, multCR) );
01719 
01720   int i;
01721 
01722   int lenList = multCR->getMasters().size();
01723   if (lenList < 1) {
01724     fei::console_out() << "ERROR in FEI, constraint with ID="<<CRID<<" appears to have"
01725          <<" a constrained-node list of length "<<lenList<<", should be > 0."<<FEI_ENDL;
01726     ERReturn(-1);
01727   }
01728 
01729   //  recall the data stored earlier and ensure that the passed data (here, the
01730   //  node list) agrees with the initialization data
01731 
01732   std::vector<GlobalID>& CRNode_vec = multCR->getMasters();
01733   GlobalID *CRNodePtr = &CRNode_vec[0];
01734 
01735   for(i=0; i<lenList; i++) {
01736     if (CRNodePtr[i] != CRNodes[i]) {
01737       fei::console_out() << "ERROR in FEI, constraint with ID="<<CRID<<" had different node-list"
01738            << " in initCRMult than it has in loadCRMult."<<FEI_ENDL;
01739       ERReturn(-1);
01740     }
01741   }
01742 
01743   std::vector<int>& CRField_vec = multCR->getMasterFieldIDs();
01744   int *CRFieldPtr = &CRField_vec[0];
01745 
01746   for (i = 0; i < lenList; i++) {
01747     if (CRFieldPtr[i] != CRFields[i]) {
01748       fei::console_out() <<"ERROR in FEI, constraint with CRID="<<CRID<<" had different field-list"
01749            <<" in initCRMult than it has in loadCRMult."<<FEI_ENDL;
01750       ERReturn(-1);
01751     }
01752   }
01753 
01754   newConstraintData_ = true;
01755 
01756   if ((int)iworkSpace_.size() < lenList) {
01757     iworkSpace_.resize(lenList);
01758   }
01759 
01760   int* fieldSizes = &iworkSpace_[0];
01761 
01762   for (i = 0; i < lenList; i++) {
01763     int numSolnParams = problemStructure_->getFieldSize(CRFields[i]);
01764     assert(numSolnParams >= 0);
01765     fieldSizes[i] = numSolnParams;
01766   }
01767 
01768   std::vector<double>& CRWeightArray = multCR->getMasterWeights();
01769 
01770   int offset = 0;
01771 
01772   try {
01773 
01774   for(i = 0; i < lenList; i++) {
01775     for(int j = 0; j < fieldSizes[i]; j++) {
01776       CRWeightArray.push_back(CRWeights[offset++]);
01777     }
01778   }
01779 
01780   }
01781   catch(std::runtime_error& exc) {
01782     fei::console_out() << exc.what() << FEI_ENDL;
01783     ERReturn(-1);
01784   }
01785 
01786   multCR->setRHSValue(CRValue);
01787   double* CRWeightsPtr = &CRWeightArray[0];
01788 
01789 //  next, perform assembly of the various terms into the system arrays
01790 //  (this is a good candidate for a separate function...)
01791 
01792   int irow = multCR->getEqnNumber();
01793   double zero = 0.0;
01794   double* zeroPtr = &zero;
01795   CHK_ERR( giveToMatrix(1, &irow, 1, &irow, &zeroPtr, ASSEMBLE_PUT) );
01796 
01797   CHK_ERR( giveToRHS(1, &(CRValue), &irow, ASSEMBLE_PUT));
01798 
01799   offset = 0;
01800   for(int j = 0; j < lenList; j++) {
01801     int myFieldID = CRFields[j];
01802 
01803     const NodeDescriptor& node = Filter::findNodeDescriptor(CRNodePtr[j]);
01804 
01805     //first, store the column coeficients for equation irow, the
01806     //constraint's equation.
01807     storeNodalColumnCoefs(irow, node, myFieldID, fieldSizes[j],
01808                           &(CRWeightsPtr[offset]));
01809 
01810 
01811     //next, store store the transpose of the above. i.e., column irow,
01812     //in equations associated with 'node' at 'myFieldID'.
01813 
01814     if (node.getOwnerProc() == localRank_) {
01815 
01816       storeNodalRowCoefs(node, myFieldID, fieldSizes[j],
01817                          &(CRWeightsPtr[offset]), irow);
01818     }
01819     else {
01820 
01821       storeNodalSendEqn(node, myFieldID, irow, &(CRWeightsPtr[offset]));
01822     }
01823 
01824     offset += fieldSizes[j];
01825   }
01826     
01827   return(FEI_SUCCESS);
01828 }
01829 
01830 //------------------------------------------------------------------------------
01831 int LinSysCoreFilter::loadCRPen(int CRID, 
01832                                 int numCRNodes,
01833                                 const GlobalID* CRNodes,
01834                                 const int* CRFields,
01835                                 const double* CRWeights,
01836                                 double CRValue,
01837                                 double penValue)
01838 {
01839   //
01840   // Load penalty constraint relation data
01841   //
01842 
01843   debugOutput("FEI: loadCRPen");
01844 
01845   ConstraintType* penCR = NULL;
01846   CHK_ERR( problemStructure_->getPenConstRecord(CRID, penCR) );
01847 
01848   int i;
01849   int lenList = penCR->getMasters().size();
01850   if (lenList < 1) {
01851     fei::console_out() << "ERROR in FEI, constraint with ID="<<CRID<<" appears to have"
01852          <<" a constrained-node list of length "<<lenList<<", should be > 0."<<FEI_ENDL;
01853     ERReturn(-1);
01854   }
01855 
01856   // recall the data stored earlier and ensure that the passed data (here,
01857   // the node list) agrees with the initialization data
01858 
01859   std::vector<GlobalID>& CRNode_vec = penCR->getMasters();
01860   GlobalID* CRNodePtr = &CRNode_vec[0];
01861                                   
01862   for(int j = 0; j < lenList; j++) {
01863     if (CRNodePtr[j] != CRNodes[j]) {
01864       fei::console_out() << "ERROR in FEI, constraint with ID="<<CRID<<" had different node-list"
01865            << " in initCRPen than it has in loadCRPen."<<FEI_ENDL;
01866       ERReturn(-1);
01867     }
01868   }
01869 
01870   newConstraintData_ = true;
01871 
01872   //  store the weights and rhs-value in the constraint records.
01873 
01874   if ((int)iworkSpace_.size() < lenList) {
01875     iworkSpace_.resize(lenList);
01876   }
01877 
01878   int* fieldSizes = &iworkSpace_[0];
01879 
01880   for (i = 0; i < lenList; i++) {
01881     int numSolnParams = problemStructure_->getFieldSize(CRFields[i]);
01882     assert(numSolnParams >= 0);
01883     fieldSizes[i] = numSolnParams;
01884   }
01885 
01886   std::vector<double>& CRWeightArray = penCR->getMasterWeights();
01887 
01888   try {
01889 
01890   int offset = 0;
01891   for (i = 0; i < lenList; i++) {
01892     for (int j = 0; j < fieldSizes[i]; j++) {
01893       CRWeightArray.push_back(CRWeights[offset++]);
01894     }
01895   }
01896 
01897   }
01898   catch(std::runtime_error& exc) {
01899     fei::console_out() << exc.what() << FEI_ENDL;
01900     ERReturn(-1);
01901   }
01902 
01903   penCR->setRHSValue(CRValue);
01904 
01905   double* CRWeightPtr = &CRWeightArray[0];
01906 
01907   int ioffset = 0, joffset = 0;
01908   for(i = 0; i < lenList; i++) {
01909     GlobalID iNodeID = CRNodePtr[i];
01910     int iField = CRFields[i];
01911 
01912     const NodeDescriptor& iNode = Filter::findNodeDescriptor(iNodeID);
01913     double* iweights = &(CRWeightPtr[ioffset]);
01914     ioffset += fieldSizes[i];
01915 
01916     joffset = 0;
01917     for (int j = 0; j < lenList; j++) {
01918       GlobalID jNodeID = CRNodePtr[j];
01919       int jField = CRFields[j];
01920 
01921       const NodeDescriptor& jNode = Filter::findNodeDescriptor(jNodeID);
01922       double* jweights = &(CRWeightPtr[joffset]);
01923       joffset += fieldSizes[j];
01924 
01925       double rhsValue = CRValue;
01926       if (j < lenList-1) {
01927         rhsValue = 0.0;
01928       }
01929 
01930       if (iNode.getOwnerProc() == localRank_) {
01931 
01932         storePenNodeData(iNode, iField, fieldSizes[i], iweights,
01933                          jNode, jField, fieldSizes[j], jweights,
01934                          penValue, rhsValue);
01935       }
01936       else {
01937         storePenNodeSendData(iNode, iField, fieldSizes[i], iweights,
01938                              jNode, jField, fieldSizes[j], jweights,
01939                              penValue, rhsValue);
01940       }
01941     }
01942   }
01943 
01944   return(FEI_SUCCESS);
01945 }
01946 
01947 //------------------------------------------------------------------------------
01948 int LinSysCoreFilter::parameters(int numParams, const char *const* paramStrings) {
01949 //
01950 // this function takes parameters for setting internal things like solver
01951 // and preconditioner choice, etc.
01952 //
01953    if (numParams == 0 || paramStrings == NULL) {
01954       debugOutput("#LinSysCoreFilter::parameters--- no parameters.");
01955    }
01956    else {
01957       const char* param1 = snl_fei::getParamValue("AZ_matrix_type",
01958                                                          numParams,
01959                                                          paramStrings);
01960 
01961       if (param1 != NULL) {
01962         if (!strcmp(param1, "AZ_VBR_MATRIX") ||
01963             !strcmp(param1, "blockMatrix")) {
01964           blockMatrix_ = true;
01965         }        
01966       }
01967       else {
01968         param1 = snl_fei::getParamValue("matrixType",
01969                                                numParams, paramStrings);
01970         if (param1 != NULL) {
01971           if (!strcmp(param1, "AZ_VBR_MATRIX") ||
01972               !strcmp(param1, "blockMatrix")) {
01973             blockMatrix_ = true;
01974           }        
01975         }
01976       }
01977 
01978       param1 = snl_fei::getParamValue("outputLevel",
01979                                              numParams,paramStrings);
01980       if ( param1 != NULL){
01981         std::string str(param1);
01982         FEI_ISTRINGSTREAM isstr(str);
01983         isstr >> outputLevel_;
01984       }
01985 
01986       param1 = snl_fei::getParam("resolveConflict",numParams,paramStrings);
01987       if ( param1 != NULL){
01988          resolveConflictRequested_ = true;
01989       }
01990 
01991       param1 = snl_fei::getParamValue("internalFei", numParams,paramStrings);
01992       if ( param1 != NULL ){
01993         std::string str(param1);
01994         FEI_ISTRINGSTREAM isstr(str);
01995         isstr >> internalFei_;
01996       }
01997 
01998       if (Filter::logStream() != NULL) {
01999 
02000         (*logStream())<<"#LinSysCoreFilter::parameters"<<FEI_ENDL
02001                          <<"# --- numParams: "<< numParams<<FEI_ENDL;
02002          for(int i=0; i<numParams; i++){
02003            (*logStream())<<"#------ paramStrings["<<i<<"]: "
02004                             <<paramStrings[i];
02005            if (paramStrings[i][strlen(paramStrings[i])-1] != '\n') {
02006              (*logStream())<<FEI_ENDL;
02007            }
02008          }
02009       }
02010    }
02011 
02012    CHK_ERR( Filter::parameters(numParams, paramStrings) );
02013 
02014    debugOutput("#LinSysCoreFilter leaving parameters function");
02015 
02016    return(FEI_SUCCESS);
02017 }
02018 
02019 //------------------------------------------------------------------------------
02020 int LinSysCoreFilter::loadComplete()
02021 {
02022   int len = 4;
02023   std::vector<int> flags(len), globalFlags(len);
02024   flags[0] = newMatrixData_ ? 1 : 0;
02025   flags[1] = newVectorData_ ? 1 : 0;
02026   flags[2] = newConstraintData_ ? 1 : 0;
02027   flags[3] = newBCData_ ? 1 : 0;
02028 
02029   CHK_ERR( fei::GlobalMax(comm_, flags, globalFlags) );
02030 
02031   newMatrixData_     = globalFlags[0] > 0 ? true : false;
02032   newVectorData_     = globalFlags[1] > 0 ? true : false;
02033   newConstraintData_ = globalFlags[2] > 0 ? true : false;
02034   newBCData_         = globalFlags[3] > 0 ? true : false;
02035 
02036   bool called_exchange = false;
02037   if (newMatrixData_ || newVectorData_ || newConstraintData_) {
02038     CHK_ERR( exchangeRemoteEquations() );
02039     called_exchange = true;
02040   }
02041 
02042   bool called_implbcs = false;
02043   if (newBCData_) {
02044     CHK_ERR( implementAllBCs() );
02045     called_implbcs = true;
02046   }
02047 
02048   if (called_exchange || called_implbcs ||needToCallMatrixLoadComplete_) {
02049     debugOutput("#LinSysCoreFilter calling LinSysCore matrixLoadComplete");
02050 
02051     CHK_ERR( lsc_->matrixLoadComplete() );
02052     needToCallMatrixLoadComplete_ = false;
02053   }
02054 
02055   newMatrixData_ = false;
02056   newVectorData_ = false;
02057   newConstraintData_ = false;
02058   newBCData_ = false;
02059 
02060   return(0);
02061 }
02062 
02063 //------------------------------------------------------------------------------
02064 int LinSysCoreFilter::residualNorm(int whichNorm, int numFields,
02065                            int* fieldIDs, double* norms, double& residTime)
02066 {
02067   //
02068   //This function can do 3 kinds of norms: infinity-norm (denoted
02069   //by whichNorm==0), 1-norm and 2-norm.
02070   //
02071   debugOutput("FEI: residualNorm");
02072 
02073   if (whichNorm < 0 || whichNorm > 2) return(-1);
02074 
02075   CHK_ERR( loadComplete() );
02076 
02077   std::vector<double> residValues(numReducedRows_, 0.0);
02078 
02079   double start = fei::utils::cpu_time();
02080 
02081   CHK_ERR( formResidual(&(residValues[0]), numReducedRows_) );
02082 
02083   residTime = fei::utils::cpu_time() - start;
02084 
02085   CHK_ERR( Filter::calculateResidualNorms(whichNorm, numFields,
02086                                           fieldIDs, norms, residValues) );
02087 
02088   return(FEI_SUCCESS);
02089 }
02090 
02091 //------------------------------------------------------------------------------
02092 int LinSysCoreFilter::formResidual(double* residValues, int numLocalEqns)
02093 {
02094   CHK_ERR( lsc_->formResidual(residValues, numLocalEqns))
02095 
02096   return(FEI_SUCCESS);
02097 }
02098 
02099 //------------------------------------------------------------------------------
02100 int LinSysCoreFilter::solve(int& status, double& sTime) {
02101 
02102    debugOutput("FEI: solve");
02103 
02104    CHK_ERR( loadComplete() );
02105 
02106    debugOutput("#LinSysCoreFilter in solve, calling launchSolver...");
02107  
02108    double start = fei::utils::cpu_time();
02109 
02110    CHK_ERR( lsc_->launchSolver(status, iterations_) );
02111 
02112    sTime = fei::utils::cpu_time() - start;
02113 
02114    debugOutput("#LinSysCoreFilter... back from solver");
02115  
02116    //now unpack the locally-owned shared entries of the solution vector into
02117    //the eqn-comm-mgr data structures.
02118    CHK_ERR( unpackSolution() );
02119 
02120    debugOutput("#LinSysCoreFilter leaving solve");
02121 
02122    if (status != 0) return(1);
02123    else return(FEI_SUCCESS);
02124 }
02125 
02126 //------------------------------------------------------------------------------
02127 int LinSysCoreFilter::setNumRHSVectors(int numRHSs, int* rhsIDs){
02128 
02129    if (numRHSs < 0) {
02130       fei::console_out() << "LinSysCoreFilter::setNumRHSVectors: ERROR, numRHSs < 0." << FEI_ENDL;
02131       ERReturn(-1);
02132    }
02133 
02134    numRHSs_ = numRHSs;
02135 
02136    rhsIDs_.resize(numRHSs_);
02137    for(int i=0; i<numRHSs_; i++) rhsIDs_[i] = rhsIDs[i];
02138 
02139   //(we need to set the number of RHSs in the eqn comm manager)
02140   eqnCommMgr_->setNumRHSs(numRHSs_);
02141 
02142    return(FEI_SUCCESS);
02143 }
02144 
02145 //------------------------------------------------------------------------------
02146 int LinSysCoreFilter::setCurrentRHS(int rhsID)
02147 {
02148    std::vector<int>::iterator iter =
02149       std::find( rhsIDs_.begin(), rhsIDs_.end(), rhsID);
02150 
02151    if (iter == rhsIDs_.end()) ERReturn(-1)
02152  
02153    int index = iter - rhsIDs_.begin();
02154    currentRHS_ = index;
02155 
02156    lsc_->setRHSID(rhsID);
02157 
02158    return(FEI_SUCCESS);
02159 }
02160 
02161 //------------------------------------------------------------------------------
02162 int LinSysCoreFilter::giveToMatrix_symm_noSlaves(int numPtRows,
02163                                                  const int* ptRowNumbers,
02164                                                  const double* const* coefs,
02165                                                  int mode)
02166 {
02167   for(int i=0; i<numPtRows; i++) {
02168     int row = ptRowNumbers[i];
02169     const double* valptr = coefs[i];
02170     if (row < localStartRow_ || row > localEndRow_) {
02171       eqnCommMgr_->addRemoteEqn(row, valptr, ptRowNumbers, numPtRows);
02172       continue;
02173     }
02174 
02175     if (mode == ASSEMBLE_SUM) {
02176       if (Filter::logStream() != NULL && 0) {
02177         FEI_OSTREAM& os = *logStream();
02178         os << "#  calling sumIntoSystemMatrix, row: " << ptRowNumbers[i]
02179            << ", columns: ";
02180         for(int j=0; j<numPtRows; ++j) os << ptRowNumbers[j] << " ";
02181         os << FEI_ENDL;
02182       }
02183 
02184       CHK_ERR( lsc_->sumIntoSystemMatrix(1, &(ptRowNumbers[i]),
02185                                          numPtRows, ptRowNumbers,
02186                                          &valptr) );
02187     }
02188     else {
02189       CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRowNumbers[i]),
02190                                          numPtRows, ptRowNumbers,
02191                                          &valptr) );
02192     }
02193   }
02194 
02195   return(0);
02196 }
02197 
02198 //------------------------------------------------------------------------------
02199 int LinSysCoreFilter::giveToBlkMatrix_symm_noSlaves(int numPtRows,
02200                                                     const int* ptRowNumbers,
02201                                                     int numBlkRows,
02202                                                     const int* blkRowNumbers,
02203                                                     const int* blkRowSizes,
02204                                                     const double* const* coefs,
02205                                                     int mode)
02206 {
02207   int i;
02208   if ((int)dworkSpace2_.size() < numPtRows) {
02209     dworkSpace2_.resize(numPtRows);
02210   }
02211   const double* * valptr = &dworkSpace2_[0];
02212   for(i=0; i<numPtRows; i++) {
02213     int row = ptRowNumbers[i];
02214     valptr[i] = coefs[i];
02215     if (row < localStartRow_ || row > localEndRow_) {
02216       eqnCommMgr_->addRemoteEqn(row, valptr[i], ptRowNumbers, numPtRows);
02217       continue;
02218     }
02219 
02220     if (mode == ASSEMBLE_PUT) {
02221        CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRowNumbers[i]),
02222                                          numPtRows, ptRowNumbers,
02223                                          &(valptr[i])) );
02224    }
02225   }
02226 
02227   int offset = 0;
02228   for(i=0; i<numBlkRows; i++) {
02229     int row = ptRowNumbers[offset];
02230     if (row < localStartRow_ || row > localEndRow_) {
02231       offset += blkRowSizes[i];
02232       continue;
02233     }
02234 
02235     if (mode == ASSEMBLE_SUM) {
02236       if (Filter::logStream() != NULL && 0) {
02237         FEI_OSTREAM& os = *logStream();
02238         os << "#  calling sumIntoSystemMatrix, row: " << ptRowNumbers[i]
02239            << ", columns: ";
02240         for(int j=0; j<numPtRows; ++j) os << ptRowNumbers[j] << " ";
02241         os << FEI_ENDL;
02242       }
02243       
02244       CHK_ERR(lsc_->sumIntoSystemMatrix(blkRowSizes[i],&(ptRowNumbers[offset]),
02245                                         numPtRows, ptRowNumbers,
02246                                         1, &(blkRowNumbers[i]),
02247                                         numBlkRows, blkRowNumbers,
02248                                         &(valptr[offset])) );
02249     }
02250     
02251     offset += blkRowSizes[i];
02252   }
02253 
02254   return(0);
02255 }
02256 
02257 //------------------------------------------------------------------------------
02258 int LinSysCoreFilter::giveToMatrix(int numPtRows, const int* ptRows,
02259                            int numPtCols, const int* ptCols,
02260                            const double* const* values, int mode)
02261 {
02262   try {
02263 
02264   if (problemStructure_->numSlaveEquations() == 0) {
02265     for(int i=0; i<numPtRows; i++) {
02266       if (ptRows[i] < localStartRow_ || ptRows[i] > localEndRow_) {
02267         eqnCommMgr_->addRemoteEqn(ptRows[i], values[i], ptCols, numPtCols);
02268         continue;
02269       }
02270 
02271       if (mode == ASSEMBLE_SUM) {
02272         if (Filter::logStream() != NULL && 0) {
02273           FEI_OSTREAM& os = *logStream();
02274           os << "#  calling sumIntoSystemMatrix, row: " << ptRows[i]
02275              << ", columns: ";
02276           for(int j=0; j<numPtCols; ++j) os << ptCols[j] << " ";
02277           os << FEI_ENDL;
02278         }
02279 
02280         CHK_ERR( lsc_->sumIntoSystemMatrix(1, &(ptRows[i]),
02281                                            numPtCols, ptCols,
02282                                            &(values[i])) );
02283       }
02284       else {
02285         CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRows[i]),
02286                                            numPtCols, ptCols,
02287                                            &(values[i])) );
02288       }
02289     }
02290   }
02291   else {
02292     iworkSpace_.resize(numPtCols);
02293     iworkSpace2_.resize(numPtCols);
02294     int* iworkPtr = &iworkSpace_[0];
02295     int* iworkPtr2= &iworkSpace2_[0];
02296     int offset = 0;
02297     for(int ii=0; ii<numPtCols; ii++) {
02298       int reducedEqn = -1;
02299       bool isSlave = problemStructure_->translateToReducedEqn(ptCols[ii],
02300                                                               reducedEqn);
02301       if (isSlave) {
02302         reducedEqn = -1;
02303         iworkPtr[ii] = reducedEqn;
02304       }
02305       else {
02306         iworkPtr[ii] = reducedEqn;
02307         iworkPtr2[offset++] = reducedEqn;
02308       }
02309     }
02310     iworkSpace2_.resize(offset);
02311 
02312     for(int i=0; i<numPtRows; i++) {
02313       int row = ptRows[i];
02314 
02315       int reducedRow;
02316       bool isSlave = problemStructure_->translateToReducedEqn(row, reducedRow);
02317       if (isSlave) continue;
02318 
02319       if (reducedStartRow_ > reducedRow || reducedRow > reducedEndRow_) {
02320 
02321         dworkSpace_.resize(0);
02322         for(int j=0; j<numPtCols; j++) {
02323           if (iworkSpace_[j]>=0) {
02324             if (Filter::logStream() != NULL) {
02325               (*logStream())<<"#  giveToMatrix remote("<<reducedRow<<","
02326                             <<iworkSpace_[j]<<","<<values[i][j]<<")"<<FEI_ENDL;
02327             }
02328 
02329             dworkSpace_.push_back(values[i][j]);
02330           }
02331         }
02332 
02333         if (mode == ASSEMBLE_SUM) {
02334           if (Filter::logStream() != NULL) {
02335             (*logStream())<<"sum"<<FEI_ENDL;
02336           }
02337 
02338           eqnCommMgr_->addRemoteEqn(reducedRow,
02339                                     &dworkSpace_[0],
02340                                     &iworkSpace2_[0],
02341                                     iworkSpace2_.size());
02342         }
02343         else {
02344           if (Filter::logStream() != NULL) {
02345             (*logStream())<<"put"<<FEI_ENDL;
02346           }
02347 
02348           eqnCommMgr_put_->addRemoteEqn(reducedRow,
02349                                         &dworkSpace_[0],
02350                                         &iworkSpace2_[0],
02351                                         iworkSpace2_.size());
02352         }
02353 
02354         continue;
02355       }
02356 
02357       for(int j=0; j<numPtCols; j++) {
02358 
02359         int reducedCol = iworkPtr[j];
02360         if (reducedCol<0) continue;
02361 
02362         double* tmpCoef = const_cast<double*>(&(values[i][j]));
02363 
02364         if (Filter::logStream() != NULL) {
02365           (*logStream())<< "#  giveToMatrix local("<<reducedRow
02366                         <<","<<reducedCol<<","<<*tmpCoef<<")"<<FEI_ENDL;
02367         }
02368 
02369         if (mode == ASSEMBLE_SUM) {
02370           if (Filter::logStream() != NULL && 0) {
02371             FEI_OSTREAM& os = *logStream();
02372             os << "#  calling sumIntoSystemMatrix, row: " << reducedRow
02373                << ", columns: " << reducedCol << FEI_ENDL;
02374           }
02375 
02376           CHK_ERR( lsc_->sumIntoSystemMatrix(1, &reducedRow, 1, &reducedCol,
02377                                              &tmpCoef ) );
02378         }
02379         else {
02380           CHK_ERR( lsc_->putIntoSystemMatrix(1, &reducedRow, 1, &reducedCol,
02381                                              &tmpCoef ) );
02382         }
02383       }
02384     }
02385   }
02386 
02387   }
02388   catch(std::runtime_error& exc) {
02389     fei::console_out() << exc.what() << FEI_ENDL;
02390     ERReturn(-1);
02391   }
02392 
02393   return(FEI_SUCCESS);
02394 }
02395 
02396 //------------------------------------------------------------------------------
02397 int LinSysCoreFilter::giveToLocalReducedMatrix(int numPtRows, const int* ptRows,
02398                                        int numPtCols, const int* ptCols,
02399                                        const double* const* values, int mode)
02400 {
02401   bool specialCase = (!firstRemEqnExchange_ && newConstraintData_
02402                       && !newMatrixData_) ? true : false;
02403 
02404   double fei_min = std::numeric_limits<double>::min();
02405 
02406   for(int i=0; i<numPtRows; i++) {
02407 
02408     if (mode == ASSEMBLE_SUM) {
02409       const double* values_i = values[i];
02410 
02411       for(int j=0; j<numPtCols; ++j) {
02412         if (specialCase && std::abs(values_i[j]) < fei_min) continue;
02413 
02414         const double* valPtr = &(values_i[j]);
02415         CHK_ERR( lsc_->sumIntoSystemMatrix(1, &(ptRows[i]), 1, &(ptCols[j]),
02416                                            &valPtr) );
02417       }
02418     }
02419     else {
02420       CHK_ERR( lsc_->putIntoSystemMatrix(1, &(ptRows[i]), numPtCols, ptCols,
02421                                          &(values[i])) );
02422     }
02423   }
02424 
02425   return(FEI_SUCCESS);
02426 }
02427 
02428 //------------------------------------------------------------------------------
02429 int LinSysCoreFilter::sumIntoMatrix(fei::CSRMat& mat)
02430 {
02431   const std::vector<int>& rowNumbers = mat.getGraph().rowNumbers;
02432   const std::vector<int>& rowOffsets = mat.getGraph().rowOffsets;
02433   const std::vector<int>& pckColInds = mat.getGraph().packedColumnIndices;
02434   const std::vector<double>& pckCoefs = mat.getPackedCoefs();
02435 
02436   for(size_t i=0; i<rowNumbers.size(); ++i) {
02437     int row = rowNumbers[i];
02438     int offset = rowOffsets[i];
02439     int rowlen = rowOffsets[i+1]-offset;
02440     const int* indices = &pckColInds[offset];
02441     const double* coefs = &pckCoefs[offset];
02442 
02443     if (giveToMatrix(1, &row, rowlen, indices, &coefs, ASSEMBLE_SUM) != 0) {
02444       throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.sumIn.");
02445     }
02446   }
02447 
02448   return(FEI_SUCCESS);
02449 }
02450 
02451 //------------------------------------------------------------------------------
02452 int LinSysCoreFilter::getFromMatrix(int numPtRows, const int* ptRows,
02453                             const int* rowColOffsets, const int* ptCols,
02454                             int numColsPerRow, double** values)
02455 {
02456   //This function may be attempting to retrieve matrix rows that are not
02457   //locally owned. If those rows correspond to finite-element nodes that we
02458   //share, AND if the owning processor is also making this function call, then
02459   //we can communicate with that processor and obtain those matrix rows.
02460   //
02461 
02462   ProcEqns remoteProcEqns;
02463 
02464   //Let's populate this ProcEqns object with the remote equations and the procs
02465   //that we need to receive the remote equations from.
02466   for(int re=0; re<numPtRows; re++) {
02467     int eqn = ptRows[re];
02468     int owner = problemStructure_->getOwnerProcForEqn(eqn);
02469     if (owner == localRank_) continue;
02470 
02471     remoteProcEqns.addEqn(eqn, owner);
02472   }
02473 
02474   //so now we know which of the requested equations are remotely owned, and we
02475   //know which processors own them.
02476   //Next we're going to need to know which locally-owned equations are needed
02477   //by other processors.
02478   ProcEqns localProcEqns;
02479   CHK_ERR( eqnCommMgr_->mirrorProcEqns(remoteProcEqns, localProcEqns) )
02480 
02481   //ok, now we know which local equations we'll need to send, so let's extract
02482   //those from the matrix
02483   EqnBuffer localEqns;
02484   CHK_ERR( getEqnsFromMatrix(localProcEqns, localEqns) )
02485 
02486   //now we can set the lengths in localProcEqns.
02487   std::vector<int>& eqnNumbers = localEqns.eqnNumbers();
02488   fei::CSVec** localEqnsPtr = (localEqns.eqns().size() ? &(localEqns.eqns()[0]) : 0);
02489   std::vector<int> eqnLengths(eqnNumbers.size());
02490   for(size_t i=0; i<eqnNumbers.size(); ++i) {
02491     eqnLengths[i] = localEqnsPtr[i]->size();
02492   }
02493 
02494   localProcEqns.setProcEqnLengths(&eqnNumbers[0], &eqnLengths[0],
02495                                   eqnNumbers.size());
02496 
02497   //now mirror those lengths in the remoteProcEqns objects to get ready for the
02498   //all-to-all exchange of equation data.
02499   CHK_ERR( eqnCommMgr_->mirrorProcEqnLengths(localProcEqns, remoteProcEqns) );
02500 
02501   EqnBuffer remoteEqns;
02502   //we're now ready to do the exchange.
02503   CHK_ERR( EqnCommMgr::exchangeEqnBuffers(comm_, &localProcEqns, &localEqns,
02504                                           &remoteProcEqns, &remoteEqns, false));
02505 
02506   std::vector<int>& remEqnNumbers = remoteEqns.eqnNumbers();
02507   fei::CSVec** remEqnsPtr = (remoteEqns.eqns().size() ? &(remoteEqns.eqns()[0]) : 0);
02508   std::vector<fei::CSVec*>& remEqns   = remoteEqns.eqns();
02509 
02510   //now we're ready to fill the values array with the remote coefficients.
02511   for(int i=0; i<numPtRows; i++) {
02512     int row = ptRows[i];
02513 
02514     int eqnIndex = fei::binarySearch(row, remEqnNumbers);
02515 
02516     //if eqnIndex < 0, this is a local equation, so skip to the next loop iter.
02517     if (eqnIndex < 0) continue;
02518 
02519     //the equation is remote, so stride across it copying out the coefs.
02520     //if ptCols is NULL, then we're going to copy all coefficients (the whole
02521     //row) into 'values'.
02522     if (ptCols == NULL) {
02523       for(size_t j=0; j<remEqnsPtr[eqnIndex]->size(); j++) {
02524         values[i][j] = remEqns[eqnIndex]->coefs()[j];
02525       }
02526       continue;
02527     }
02528 
02529     for(int j=0; j<numColsPerRow; j++) {
02530       int offset = rowColOffsets[i] + j;
02531       int colIndex = fei::binarySearch(ptCols[offset], remEqns[eqnIndex]->indices());
02532       if (colIndex < 0) ERReturn(-1);
02533 
02534       values[i][j] = remEqns[eqnIndex]->coefs()[colIndex];
02535     }
02536   }
02537 
02538   //and now, get the local stuff out of the matrix.
02539   for(int i=0; i<numPtRows; i++) {
02540     int row = ptRows[i];
02541     if (row < localStartRow_ || localEndRow_ < row) continue;
02542 
02543     int rowLen = 0, checkRowLen;
02544     CHK_ERR( lsc_->getMatrixRowLength(row, rowLen) )
02545       if (rowLen <= 0) ERReturn(-1);
02546 
02547     //for each local row, establish some temp arrays and get the row from
02548     //the matrix.
02549 
02550     std::vector<double> coefs(rowLen);
02551     std::vector<int> indices(rowLen);
02552 
02553     CHK_ERR( lsc_->getMatrixRow(row, &coefs[0], &indices[0], rowLen, checkRowLen) );
02554     if (rowLen != checkRowLen) ERReturn(-1);
02555 
02556     //now stride across the list of requested column-indices, and find the
02557     //corresponding location in the matrix row. Copy that location into the
02558     //values array.
02559 
02560     //again, if ptCols is NULL, then we're going to copy all coefficients 
02561     //(the whole row) into 'values'.
02562     if (ptCols == NULL) {
02563       for(int j=0; j<rowLen; j++) {
02564         values[i][j] = coefs[j];
02565       }
02566       continue;
02567     }
02568 
02569     for(int j=0; j<numColsPerRow; j++) {
02570       std::vector<int>::iterator iter =
02571           std::find(indices.begin(), indices.end(), ptCols[rowColOffsets[i]+j]);
02572       if (iter == indices.end()) {
02573         ERReturn(-1);
02574       }
02575 
02576       int index = iter - indices.begin();
02577       values[i][j] = coefs[index];
02578     }
02579   }
02580 
02581   return(FEI_SUCCESS);
02582 }
02583 
02584 //------------------------------------------------------------------------------
02585 int LinSysCoreFilter::getEqnsFromMatrix(ProcEqns& procEqns, EqnBuffer& eqnData)
02586 {
02587   //Given a ProcEqns object containing lists of equation-numbers, get the data
02588   //for those equations from the local portion of the matrix and store that data
02589   //in the eqnData object.
02590 
02591   std::vector<std::vector<int>*>& eqnNumbers = procEqns.procEqnNumbersPtr();
02592 
02593   for(unsigned p=0; p<eqnNumbers.size(); p++) {
02594     for(unsigned i=0; i<eqnNumbers[p]->size(); i++) {
02595       int eqn = (*(eqnNumbers[p]))[i];
02596 
02597       if (localStartRow_ > eqn || localEndRow_ < eqn) continue;
02598 
02599       //if this equation is already in eqnData, then don't put it in again...
02600       std::vector<int>& eqnDataEqns = eqnData.eqnNumbers();
02601       if (fei::binarySearch(eqn, eqnDataEqns) >= 0) continue;
02602 
02603       int len = 0;
02604       CHK_ERR( lsc_->getMatrixRowLength(eqn, len) );
02605 
02606       if (len <= 0) continue;
02607       std::vector<double> coefs(len);
02608       std::vector<int> indices(len);
02609       int outputLen = 0;
02610 
02611       CHK_ERR( lsc_->getMatrixRow(eqn, &coefs[0], &indices[0],
02612                                   len, outputLen) );
02613       if (outputLen != len) ERReturn(-1);
02614 
02615       CHK_ERR( eqnData.addEqn(eqn, &coefs[0], &indices[0], len, false) );
02616     }
02617   }
02618   return(FEI_SUCCESS);
02619 }
02620 
02621 //------------------------------------------------------------------------------
02622 int LinSysCoreFilter::getEqnsFromRHS(ProcEqns& procEqns, EqnBuffer& eqnData)
02623 {
02624   //Given a ProcEqns object containing lists of equation-numbers, get the data
02625   //for those equations from the local portion of the RHS vector and store that
02626   // data in the eqnData object. We're only storing rhs coefs in an EqnBuffer
02627   //that was designed for also storing equations with column-indices. So we'll
02628   //put a bogus column-index in with each equation just to make sure the
02629   //EqnBuffer does the right stuff internally...
02630 
02631   int numSendProcs = procEqns.getNumProcs();
02632   std::vector<int>& eqnsPerProc = procEqns.eqnsPerProcPtr();
02633   std::vector<std::vector<int>*>& eqnNumbers = procEqns.procEqnNumbersPtr();
02634 
02635   eqnData.setNumRHSs(1);
02636 
02637   for(int p=0; p<numSendProcs; p++) {
02638     for(int i=0; i<eqnsPerProc[p]; i++) {
02639       int reducedEqn;
02640       problemStructure_->translateToReducedEqn((*(eqnNumbers[p]))[i], reducedEqn);
02641 
02642       if (reducedStartRow_ > reducedEqn || reducedEndRow_ < reducedEqn) continue;
02643 
02644       //if this equation is already in eqnData, then don't put it in again...
02645       std::vector<int>& eqnDataEqns = eqnData.eqnNumbers();
02646       if (fei::binarySearch(reducedEqn, eqnDataEqns) >= 0) continue;
02647 
02648       double rhsValue;
02649 
02650       CHK_ERR( lsc_->getFromRHSVector(1, &rhsValue, &reducedEqn) );
02651 
02652       int bogusIndex = 19;
02653       CHK_ERR( eqnData.addIndices(reducedEqn, &bogusIndex, 1) );
02654       CHK_ERR( eqnData.addRHS(reducedEqn, 0, rhsValue) );
02655     }
02656   }
02657   return(FEI_SUCCESS);
02658 }
02659 
02660 //------------------------------------------------------------------------------
02661 int LinSysCoreFilter::giveToRHS(int num, const double* values,
02662                         const int* indices, int mode)
02663 {
02664   if (problemStructure_->numSlaveEquations() == 0) {
02665     for(int i=0; i<num; i++) {
02666       if (indices[i] < localStartRow_ || indices[i] > localEndRow_) {
02667         if (mode == ASSEMBLE_SUM) {
02668           eqnCommMgr_->addRemoteRHS(indices[i], currentRHS_, values[i]);
02669         }
02670         else {
02671           eqnCommMgr_put_->addRemoteRHS(indices[i], currentRHS_, values[i]);
02672         }
02673 
02674         continue;
02675       }
02676 
02677       if (mode == ASSEMBLE_SUM) {
02678         CHK_ERR( lsc_->sumIntoRHSVector(1, &(values[i]), &(indices[i])) );
02679       }
02680       else {
02681         CHK_ERR( lsc_->putIntoRHSVector(1, &(values[i]), &(indices[i])) );
02682       }
02683     }
02684   }
02685   else {
02686     for(int i=0; i<num; i++) {
02687       int reducedEqn;
02688       bool isSlave = problemStructure_->
02689         translateToReducedEqn(indices[i], reducedEqn);
02690       if (isSlave) continue;
02691 
02692       if (reducedEqn < reducedStartRow_ || reducedEqn > reducedEndRow_) {
02693         if (mode == ASSEMBLE_SUM) {
02694           eqnCommMgr_->addRemoteRHS(reducedEqn, currentRHS_, values[i]);
02695         }
02696         else {
02697           eqnCommMgr_put_->addRemoteRHS(reducedEqn, currentRHS_, values[i]);
02698         }
02699 
02700         continue;
02701       }
02702 
02703       if (mode == ASSEMBLE_SUM) {
02704         CHK_ERR( lsc_->sumIntoRHSVector(1, &(values[i]), &reducedEqn) );
02705       }
02706       else {
02707         CHK_ERR( lsc_->putIntoRHSVector(1, &(values[i]), &reducedEqn) );
02708       }
02709     }
02710   }
02711 
02712   return(FEI_SUCCESS);
02713 }
02714 
02715 //------------------------------------------------------------------------------
02716 int LinSysCoreFilter::giveToLocalReducedRHS(int num, const double* values,
02717                                     const int* indices, int mode)
02718 {
02719   for(int i=0; i<num; i++) {
02720 
02721     if (mode == ASSEMBLE_SUM) {
02722       CHK_ERR( lsc_->sumIntoRHSVector(1, &(values[i]), &(indices[i])) );
02723     }
02724     else {
02725       CHK_ERR( lsc_->putIntoRHSVector(1, &(values[i]), &(indices[i])) );
02726     }
02727   }
02728 
02729   return(FEI_SUCCESS);
02730 }
02731 
02732 //------------------------------------------------------------------------------
02733 int LinSysCoreFilter::sumIntoRHS(fei::CSVec& vec)
02734 {
02735   std::vector<int>& indices = vec.indices();
02736   std::vector<double>& coefs = vec.coefs();
02737 
02738   CHK_ERR( giveToRHS(indices.size(), &coefs[0], &indices[0], ASSEMBLE_SUM) );
02739 
02740   return(FEI_SUCCESS);
02741 }
02742 
02743 //------------------------------------------------------------------------------
02744 int LinSysCoreFilter::getFromRHS(int num, double* values, const int* indices)
02745 {
02746   //We need to do similar things here as we do in getFromMatrix, with respect to
02747   //communications to obtain values for equations that are remotely owned.
02748 
02749   ProcEqns remoteProcEqns;
02750 
02751   //Let's populate this ProcEqns object with the remote equations and the procs
02752   //that we need to receive the remote equations from.
02753   for(int re=0; re<num; re++) {
02754     int eqn = indices[re];
02755     int owner = problemStructure_->getOwnerProcForEqn(eqn);
02756     if (owner==localRank_) continue;
02757 
02758     remoteProcEqns.addEqn(eqn, owner);
02759   }
02760 
02761   //so now we know which of the requested equations are remotely owned, and we
02762   //know which processors own them.
02763   //Next we're going to need to know which locally-owned equations are needed
02764   //by other processors.
02765   ProcEqns localProcEqns;
02766   CHK_ERR( eqnCommMgr_->mirrorProcEqns(remoteProcEqns, localProcEqns) );
02767 
02768   //ok, now we know which equations we'll need to send, so let's extract
02769   //them from the rhs vector.
02770   EqnBuffer localEqns;
02771   CHK_ERR( getEqnsFromRHS(localProcEqns, localEqns) );
02772 
02773   //now we can set the lengths in localProcEqns.
02774   std::vector<int>& eqnNumbers = localEqns.eqnNumbers();
02775   fei::CSVec** localEqnsPtr = &(localEqns.eqns()[0]);
02776   std::vector<int> eqnLengths(eqnNumbers.size());
02777   for(size_t i=0; i<eqnNumbers.size(); ++i) {
02778     eqnLengths[i] = localEqnsPtr[i]->size();
02779   }
02780 
02781   localProcEqns.setProcEqnLengths(&eqnNumbers[0], &eqnLengths[0],
02782                                   eqnNumbers.size());
02783 
02784   //now mirror those lengths in the remoteProcEqns objects to get ready for the
02785   //all-to-all exchange of equation data.
02786   CHK_ERR( eqnCommMgr_->mirrorProcEqnLengths(localProcEqns, remoteProcEqns) );
02787 
02788   EqnBuffer remoteEqns;
02789   //we're now ready to do the exchange.
02790   CHK_ERR( EqnCommMgr::exchangeEqnBuffers(comm_, &localProcEqns, &localEqns,
02791                                            &remoteProcEqns, &remoteEqns, false))
02792 
02793   //now we're ready to get the rhs data we've received from other processors.
02794   std::vector<int>& remEqnNumbers = remoteEqns.eqnNumbers();
02795   std::vector<std::vector<double>*>& remRhsCoefs = *(remoteEqns.rhsCoefsPtr());
02796 
02797   for(int i=0; i<num; i++) {
02798     int row = indices[i];
02799 
02800     int eqnIndex = fei::binarySearch(row, remEqnNumbers);
02801     if (eqnIndex < 0) continue;
02802 
02803     values[i] = (*(remRhsCoefs[eqnIndex]))[0];
02804   }
02805 
02806   //and now get the local stuff.
02807   for(int i=0; i<num; i++) {
02808     if (indices[i] < localStartRow_ || indices[i] > localEndRow_) continue;
02809 
02810     CHK_ERR( lsc_->getFromRHSVector(num, values, indices) );
02811   }
02812 
02813   return(FEI_SUCCESS);
02814 }
02815 
02816 //------------------------------------------------------------------------------
02817 int LinSysCoreFilter::getEqnSolnEntry(int eqnNumber, double& solnValue)
02818 {
02819   //This function's task is to retrieve the solution-value for a global
02820   //equation-number. eqnNumber may or may not be a slave-equation, and may or
02821   //may not be locally owned. If it is not locally owned, it should at least
02822   //be shared.
02823   //return 0 if the solution is successfully retrieved, otherwise return 1.
02824   //
02825 
02826   //
02827   //First and probably most common case: there are no slave equations.
02828   //
02829   if (problemStructure_->numSlaveEquations() == 0) {
02830     if (localStartRow_ > eqnNumber || eqnNumber > localEndRow_) {
02831       //Dig into the eqn-comm-mgr for the shared-remote solution value.
02832       CHK_ERR( getSharedRemoteSolnEntry(eqnNumber, solnValue) );
02833     }
02834     else {
02835       //It's local, simply get the solution from the assembled linear system.
02836       CHK_ERR( getReducedSolnEntry( eqnNumber, solnValue ) );
02837     }
02838     return(0);
02839   }
02840 
02841   //
02842   //If we reach this point, there are slave equations to account for.
02843   //So first translate this equation into 'assembled-linear-system'
02844   //equation-numbers.
02845   //
02846   int reducedEqn;
02847   bool isSlave = problemStructure_->translateToReducedEqn(eqnNumber,reducedEqn);
02848 
02849   if (isSlave) {
02850     //This is a slave-equation, so construct its solution-value as the linear-
02851     //combination of the master-equations it is defined in terms of.
02852 
02853     std::vector<int>* masterEqns = NULL;
02854     std::vector<double>* masterCoefs = NULL;
02855     CHK_ERR( problemStructure_->getMasterEqnNumbers(eqnNumber, masterEqns) );
02856     CHK_ERR( problemStructure_->getMasterEqnCoefs(eqnNumber, masterCoefs) );
02857 
02858     int len = masterEqns->size();
02859     solnValue = 0.0;
02860     CHK_ERR( problemStructure_->getMasterEqnRHS(eqnNumber, solnValue) );
02861 
02862     double coef = 0.0;
02863     for(int i=0; i<len; i++) {
02864       int mEqn = (*masterEqns)[i];
02865       int mReducedeqn;
02866       problemStructure_->translateToReducedEqn(mEqn, mReducedeqn);
02867 
02868       if (reducedStartRow_ > mReducedeqn || mReducedeqn > reducedEndRow_) {
02869         CHK_ERR( getSharedRemoteSolnEntry(mReducedeqn, coef) );
02870       }
02871       else {
02872         CHK_ERR( getReducedSolnEntry(mReducedeqn, coef) );
02873       }
02874       solnValue += coef * (*masterCoefs)[i];
02875     }
02876   }
02877   else {
02878     //This is not a slave-equation, so retrieve the solution from either the
02879     //assembled linear system or the shared-remote data structures.
02880 
02881     if (reducedStartRow_ > reducedEqn || reducedEqn > reducedEndRow_) {
02882       CHK_ERR( getSharedRemoteSolnEntry(reducedEqn, solnValue) );
02883     }
02884     else {
02885       CHK_ERR( getReducedSolnEntry(reducedEqn, solnValue) );
02886     }
02887   }
02888 
02889   return(0);
02890 }
02891 
02892 //------------------------------------------------------------------------------
02893 int LinSysCoreFilter::getSharedRemoteSolnEntry(int eqnNumber, double& solnValue)
02894 {
02895   std::vector<int>& remoteEqnNumbers = eqnCommMgr_->sendEqnNumbersPtr();
02896   double* remoteSoln = eqnCommMgr_->sendEqnSolnPtr();
02897 
02898   int index = fei::binarySearch(eqnNumber, remoteEqnNumbers);
02899   if (index < 0) {
02900     fei::console_out() << "LinSysCoreFilter::getSharedRemoteSolnEntry: ERROR, eqn "
02901          << eqnNumber << " not found." << FEI_ENDL;
02902     ERReturn(-1);
02903   }
02904   solnValue = remoteSoln[index];
02905   return(0);
02906 }
02907 
02908 //------------------------------------------------------------------------------
02909 int LinSysCoreFilter::getReducedSolnEntry(int eqnNumber, double& solnValue)
02910 {
02911   //We may safely assume that this function is called with 'eqnNumber' that is
02912   //local in the underlying assembled linear system. i.e., it isn't a slave-
02913   //equation, it isn't remotely owned, etc.
02914   //
02915   CHK_ERR( lsc_->getSolnEntry(eqnNumber, solnValue) );
02916 
02917   return(FEI_SUCCESS);
02918 }
02919 
02920 //------------------------------------------------------------------------------
02921 int LinSysCoreFilter::unpackSolution()
02922 {
02923   //
02924   //This function should be called after the solver has returned,
02925   //and we know that there is a solution in the underlying vector.
02926   //This function ensures that any locally-owned shared solution values are
02927   //available on the sharing processors.
02928   //
02929   if (Filter::logStream() != NULL) {
02930     (*logStream())<< "#  entering unpackSolution, outputLevel: "
02931                      <<outputLevel_<<FEI_ENDL;
02932   }
02933 
02934   //what we need to do is as follows.
02935   //The eqn comm mgr has a list of what it calls 'recv eqns'. These are
02936   //equations that we own, for which we received contributions from other
02937   //processors. The solution values corresponding to these equations need
02938   //to be made available to those remote contributing processors.
02939 
02940   int numRecvEqns = eqnCommMgr_->getNumLocalEqns();
02941   std::vector<int>& recvEqnNumbers = eqnCommMgr_->localEqnNumbers();
02942 
02943   for(int i=0; i<numRecvEqns; i++) {
02944     int eqn = recvEqnNumbers[i];
02945 
02946     if ((reducedStartRow_ > eqn) || (reducedEndRow_ < eqn)) {
02947       fei::console_out() << "LinSysCoreFilter::unpackSolution: ERROR, 'recv' eqn (" << eqn
02948            << ") out of local range." << FEI_ENDL;
02949       MPI_Abort(comm_, -1);
02950     }
02951 
02952     double solnValue = 0.0;
02953 
02954     CHK_ERR( getReducedSolnEntry(eqn, solnValue) );
02955 
02956     eqnCommMgr_->addSolnValues(&eqn, &solnValue, 1);
02957   }
02958 
02959   eqnCommMgr_->exchangeSoln();
02960 
02961   debugOutput("#LinSysCoreFilter leaving unpackSolution");
02962   return(FEI_SUCCESS);
02963 }
02964              
02965 //------------------------------------------------------------------------------
02966 void LinSysCoreFilter::setEqnCommMgr(EqnCommMgr* eqnCommMgr)
02967 {
02968   delete eqnCommMgr_;
02969   eqnCommMgr_ = eqnCommMgr;
02970 }
02971 
02972 //------------------------------------------------------------------------------
02973 int LinSysCoreFilter::getBlockNodeSolution(GlobalID elemBlockID,  
02974                                    int numNodes, 
02975                                    const GlobalID *nodeIDs, 
02976                                    int *offsets,
02977                                    double *results) {
02978         
02979    debugOutput("FEI: getBlockNodeSolution");
02980 
02981    int numActiveNodes = problemStructure_->getNumActiveNodes();
02982    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
02983 
02984    if (numActiveNodes <= 0) return(0);
02985 
02986    int numSolnParams = 0;
02987 
02988    BlockDescriptor* block = NULL;
02989    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
02990 
02991    //Traverse the node list, checking if nodes are associated with this block.
02992    //If so, put its 'answers' in the results list.
02993 
02994    int offset = 0;
02995    for(int i=0; i<numActiveNodes; i++) {
02996      const NodeDescriptor* node_i = NULL;
02997      nodeDB.getNodeAtIndex(i, node_i);
02998 
02999       if (offset == numNodes) break;
03000 
03001       GlobalID nodeID = nodeIDs[offset];
03002 
03003       //first let's set the offset at which this node's solution coefs start.
03004       offsets[offset++] = numSolnParams;
03005 
03006       const NodeDescriptor* node = NULL;
03007       int err = 0;
03008       //Obtain the NodeDescriptor of nodeID in the activeNodes list...
03009       //Don't call the getActiveNodeDesc_ID function unless we have to.
03010 
03011       if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
03012         node = node_i;
03013       }
03014       else {
03015          err = nodeDB.getNodeWithID(nodeID, node);
03016       }
03017 
03018       //ok. If err is not 0, meaning nodeID is NOT in the
03019       //activeNodes list, then skip to the next loop iteration.
03020 
03021       if (err != 0) {
03022         continue;
03023       }
03024 
03025       int numFields = node->getNumFields();
03026       const int* fieldIDs = node->getFieldIDList();
03027 
03028       for(int j=0; j<numFields; j++) {
03029         if (block->containsField(fieldIDs[j])) {
03030           int size = problemStructure_->getFieldSize(fieldIDs[j]);
03031           assert(size >= 0);
03032 
03033           int thisEqn = -1;
03034           node->getFieldEqnNumber(fieldIDs[j], thisEqn);
03035 
03036           double answer;
03037           for(int k=0; k<size; k++) {
03038             CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
03039               results[numSolnParams++] = answer;
03040           }
03041         }
03042       }//for(j<numFields)loop
03043    }
03044 
03045    offsets[numNodes] = numSolnParams;
03046 
03047    return(FEI_SUCCESS);
03048 }
03049 
03050 //------------------------------------------------------------------------------
03051 int LinSysCoreFilter::getNodalSolution(int numNodes, 
03052                                        const GlobalID *nodeIDs, 
03053                                        int *offsets,
03054                                        double *results)
03055 {
03056   debugOutput("FEI: getNodalSolution");
03057 
03058   int numActiveNodes = problemStructure_->getNumActiveNodes();
03059   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03060 
03061   if (numActiveNodes <= 0) return(0);
03062 
03063   int numSolnParams = 0;
03064 
03065   //Traverse the node list, checking if nodes are local.
03066   //If so, put 'answers' in the results list.
03067 
03068   int offset = 0;
03069   for(int i=0; i<numActiveNodes; i++) {
03070     const NodeDescriptor* node_i = NULL;
03071     nodeDB.getNodeAtIndex(i, node_i);
03072 
03073     if (offset == numNodes) break;
03074 
03075     GlobalID nodeID = nodeIDs[offset];
03076 
03077     //first let's set the offset at which this node's solution coefs start.
03078     offsets[offset++] = numSolnParams;
03079 
03080     const NodeDescriptor* node = NULL;
03081     int err = 0;
03082     //Obtain the NodeDescriptor of nodeID in the activeNodes list...
03083     //Don't call the getNodeWithID function unless we have to.
03084 
03085     if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
03086       node = node_i;
03087     }
03088     else {
03089       err = nodeDB.getNodeWithID(nodeID, node);
03090     }
03091 
03092     //ok. If err is not 0, meaning nodeID is NOT in the
03093     //activeNodes list, then skip to the next loop iteration.
03094 
03095     if (err != 0) {
03096       continue;
03097     }
03098 
03099     int numFields = node->getNumFields();
03100     const int* fieldIDs = node->getFieldIDList();
03101 
03102     for(int j=0; j<numFields; j++) {
03103       int size = problemStructure_->getFieldSize(fieldIDs[j]);
03104       assert(size >= 0);
03105 
03106       int thisEqn = -1;
03107       node->getFieldEqnNumber(fieldIDs[j], thisEqn);
03108 
03109       double answer;
03110       for(int k=0; k<size; k++) {
03111         CHK_ERR( getEqnSolnEntry(thisEqn+k, answer) )
03112           results[numSolnParams++] = answer;
03113       }
03114     }//for(j<numFields)loop
03115   }
03116 
03117   offsets[numNodes] = numSolnParams;
03118 
03119   return(FEI_SUCCESS);
03120 }
03121 
03122 //------------------------------------------------------------------------------
03123 int LinSysCoreFilter::getBlockFieldNodeSolution(GlobalID elemBlockID,
03124                                         int fieldID,
03125                                         int numNodes, 
03126                                         const GlobalID *nodeIDs, 
03127                                         double *results)
03128 {
03129   //Note: if the user-supplied nodeIDs list containts nodes which are not in
03130   //the specified element-block, then the corresponding positions in the
03131   //results array are simply not referenced. This is dangerous behavior that
03132   //hasn't gotten me into trouble yet.
03133 
03134   debugOutput("FEI: getBlockFieldNodeSolution");
03135 
03136   int numActiveNodes = problemStructure_->getNumActiveNodes();
03137   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03138 
03139   if (numActiveNodes <= 0) return(0);
03140 
03141   BlockDescriptor* block = NULL;
03142   CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
03143 
03144   int fieldSize = problemStructure_->getFieldSize(fieldID);
03145   if (fieldSize <= 0) ERReturn(-1);
03146 
03147   if (!block->containsField(fieldID)) {
03148     fei::console_out() << "LinSysCoreFilter::getBlockFieldNodeSolution WARNING: fieldID " << fieldID
03149          << " not contained in element-block " << (int)elemBlockID << FEI_ENDL;
03150     return(1);
03151   }
03152 
03153    //Traverse the node list, checking if nodes are associated with this block.
03154    //If so, put the answers in the results list.
03155 
03156    for(int i=0; i<numNodes; i++) {
03157      const NodeDescriptor* node_i = NULL;
03158      nodeDB.getNodeAtIndex(i, node_i);
03159 
03160      GlobalID nodeID = nodeIDs[i];
03161 
03162      const NodeDescriptor* node = NULL;
03163      int err = 0;
03164      //Obtain the NodeDescriptor of nodeID in the activeNodes list...
03165      //Don't call the getNodeWithID function unless we have to. (getNodeWithID
03166      //does a binary-search over all local nodeIDs, while getNodeAtIndex is
03167      //a direct lookup.) Often the user supplies a nodeIDs list that is in the
03168      //"natural" order, so we don't need to call getNodeWithID at all.
03169 
03170      if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
03171        node = node_i;
03172      }
03173      else {
03174        err = nodeDB.getNodeWithID(nodeID, node);
03175      }
03176 
03177      //ok. If err is not 0, meaning nodeID is NOT in the
03178      //activeNodes list, then skip to the next loop iteration.
03179 
03180      if (err != 0) {
03181        continue;
03182      }
03183 
03184      int eqnNumber = -1;
03185      bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
03186      if (!hasField) continue;
03187 
03188      int offset = fieldSize*i;
03189      for(int j=0; j<fieldSize; j++) {
03190        double answer = 0.0;
03191        CHK_ERR( getEqnSolnEntry(eqnNumber+j, answer) );
03192        results[offset+j] = answer;
03193      }
03194    }
03195 
03196    return(FEI_SUCCESS);
03197 }
03198 
03199 //------------------------------------------------------------------------------
03200 int LinSysCoreFilter::getNodalFieldSolution(int fieldID,
03201                                     int numNodes, 
03202                                     const GlobalID *nodeIDs, 
03203                                     double *results)
03204 {
03205   debugOutput("FEI: getNodalFieldSolution");
03206 
03207   int numActiveNodes = problemStructure_->getNumActiveNodes();
03208   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03209 
03210   if (numActiveNodes <= 0) return(0);
03211 
03212   int fieldSize = problemStructure_->getFieldSize(fieldID);
03213   if (fieldSize <= 0) ERReturn(-1);
03214 
03215   //Traverse the node list, checking if nodes have the specified field.
03216   //If so, put the answers in the results list.
03217 
03218   for(int i=0; i<numNodes; i++) {
03219     const NodeDescriptor* node_i = NULL;
03220     nodeDB.getNodeAtIndex(i, node_i);
03221 
03222     GlobalID nodeID = nodeIDs[i];
03223 
03224     const NodeDescriptor* node = NULL;
03225     int err = 0;
03226     //Obtain the NodeDescriptor of nodeID in the activeNodes list...
03227     //Don't call the getNodeWithID function unless we have to.
03228 
03229     if (node_i!=NULL && nodeID == node_i->getGlobalNodeID()) {
03230       node = node_i;
03231     }
03232     else {
03233       err = nodeDB.getNodeWithID(nodeID, node);
03234     }
03235 
03236     //ok. If err is not 0, meaning nodeID is NOT in the
03237     //activeNodes list, then skip to the next loop iteration.
03238 
03239     if (err != 0) {
03240       continue;
03241     }
03242 
03243     int eqnNumber = -1;
03244     bool hasField = node->getFieldEqnNumber(fieldID, eqnNumber);
03245 
03246     //If this node doesn't have the specified field, then skip to the
03247     //next loop iteration.
03248     if (!hasField) continue;
03249 
03250     int offset = fieldSize*i;
03251     for(int j=0; j<fieldSize; j++) {
03252       double answer = 0.0;
03253       CHK_ERR( getEqnSolnEntry(eqnNumber+j, answer) );
03254       results[offset+j] = answer;
03255     }
03256   }
03257 
03258   return(FEI_SUCCESS);
03259 }
03260 
03261 //------------------------------------------------------------------------------
03262 int LinSysCoreFilter::putBlockNodeSolution(GlobalID elemBlockID,
03263                                    int numNodes, 
03264                                    const GlobalID *nodeIDs, 
03265                                    const int *offsets,
03266                                    const double *estimates) {
03267         
03268    debugOutput("FEI: putBlockNodeSolution");
03269 
03270    int numActiveNodes = problemStructure_->getNumActiveNodes();
03271 
03272    if (numActiveNodes <= 0) return(0);
03273 
03274    BlockDescriptor* block = NULL;
03275    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
03276 
03277    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03278 
03279    //traverse the node list, checking for nodes associated with this block
03280    //when an associated node is found, put its 'answers' into the linear system.
03281 
03282    unsigned blk_idx = problemStructure_->getIndexOfBlock(elemBlockID);
03283 
03284    for(int i=0; i<numNodes; i++) {
03285      const NodeDescriptor* node = NULL;
03286      int err = nodeDB.getNodeWithID(nodeIDs[i], node);
03287 
03288       if (err != 0) continue;
03289    
03290       if (!node->hasBlockIndex(blk_idx)) continue;
03291 
03292       if (node->getOwnerProc() != localRank_) continue;
03293 
03294       int numFields = node->getNumFields();
03295       const int* fieldIDs = node->getFieldIDList();
03296       const int* fieldEqnNumbers = node->getFieldEqnNumbers();
03297 
03298       if (fieldEqnNumbers[0] < localStartRow_ ||
03299           fieldEqnNumbers[0] > localEndRow_) continue;
03300 
03301       int offs = offsets[i];
03302 
03303       for(int j=0; j<numFields; j++) {
03304          int size = problemStructure_->getFieldSize(fieldIDs[j]);
03305 
03306          if (block->containsField(fieldIDs[j])) {
03307             for(int k=0; k<size; k++) {
03308                int reducedEqn;
03309                problemStructure_->
03310                  translateToReducedEqn(fieldEqnNumbers[j]+k, reducedEqn);
03311 
03312                CHK_ERR( lsc_->putInitialGuess(&reducedEqn,
03313                                               &estimates[offs+k], 1) );
03314             }
03315          }
03316          offs += size;
03317       }//for(j<numFields)loop
03318    }
03319 
03320    return(FEI_SUCCESS);
03321 }
03322 
03323 //------------------------------------------------------------------------------
03324 int LinSysCoreFilter::putBlockFieldNodeSolution(GlobalID elemBlockID, 
03325                                         int fieldID, 
03326                                         int numNodes, 
03327                                         const GlobalID *nodeIDs, 
03328                                         const double *estimates)
03329 {
03330    int fieldSize = problemStructure_->getFieldSize(fieldID);
03331 
03332    if (Filter::logStream() != NULL) {
03333      FEI_OSTREAM& os = *logStream();
03334      os << "FEI: putBlockFieldNodeSolution" << FEI_ENDL;
03335      os << "#blkID" << FEI_ENDL << (int)elemBlockID << FEI_ENDL
03336         << "#fieldID"<<FEI_ENDL << fieldID << FEI_ENDL
03337         << "#fieldSize"<<FEI_ENDL << fieldSize << FEI_ENDL
03338         << "#numNodes"<<FEI_ENDL << numNodes << FEI_ENDL
03339         << "#nodeIDs" << FEI_ENDL;
03340      int i;
03341      for(i=0; i<numNodes; ++i) os << (int)nodeIDs[i] << FEI_ENDL;
03342      os << "#estimates" << FEI_ENDL;
03343      for(i=0; i<numNodes*fieldSize; ++i) os << estimates[i] << FEI_ENDL;
03344    }
03345 
03346    BlockDescriptor* block = NULL;
03347    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) );
03348    if (!block->containsField(fieldID)) return(1);
03349 
03350    NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03351 
03352    //if we have a negative fieldID, we'll need a list of length numNodes,
03353    //in which to put nodeNumbers for passing to the solver... 
03354 
03355    std::vector<int> numbers(numNodes);
03356 
03357    //if we have a fieldID >= 0, then our numbers list will hold equation numbers
03358    //and we'll need fieldSize*numNodes of them.
03359 
03360    std::vector<double> data;
03361 
03362    if (fieldID >= 0) {
03363       assert(fieldSize >= 0);
03364       numbers.resize(numNodes*fieldSize);
03365       data.resize(numNodes*fieldSize);
03366    }
03367 
03368    int count = 0;
03369 
03370    for(int i=0; i<numNodes; i++) {
03371      const NodeDescriptor* node = NULL;
03372      CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
03373 
03374       if (fieldID < 0) numbers[count++] = node->getNodeNumber();
03375       else {
03376          int eqn = -1;
03377          if (node->getFieldEqnNumber(fieldID, eqn)) {
03378            if (eqn >= localStartRow_ && eqn <= localEndRow_) {
03379              for(int j=0; j<fieldSize; j++) { 
03380                data[count] = estimates[i*fieldSize + j];
03381                problemStructure_->translateToReducedEqn(eqn+j, numbers[count++]);
03382              }
03383            }
03384          }
03385       }
03386    }
03387 
03388    if (fieldID < 0) {
03389      CHK_ERR( lsc_->putNodalFieldData(fieldID, fieldSize, 
03390                                       &numbers[0], numNodes, estimates));
03391    }
03392    else {
03393      CHK_ERR(lsc_->putInitialGuess(&numbers[0], &data[0], count));
03394    }
03395 
03396    return(FEI_SUCCESS);
03397 }
03398 
03399 //------------------------------------------------------------------------------
03400 int LinSysCoreFilter::getBlockElemSolution(GlobalID elemBlockID,
03401                                    int numElems, 
03402                                    const GlobalID *elemIDs,
03403                                    int& numElemDOFPerElement,
03404                                    double *results)
03405 {
03406 //
03407 //  return the elemental solution parameters associated with a
03408 //  particular block of elements
03409 //
03410    debugOutput("FEI: getBlockElemSolution");
03411 
03412    BlockDescriptor* block = NULL;
03413    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
03414 
03415    numElemDOFPerElement = block->getNumElemDOFPerElement();
03416    if (numElemDOFPerElement <= 0) return(0);
03417 
03418    ConnectivityTable& ctable =
03419      problemStructure_->getBlockConnectivity(elemBlockID);
03420    std::map<GlobalID,int>& elemIDList = ctable.elemIDs;
03421 
03422    std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
03423    double answer;
03424 
03425    for(int i=0; i<numElems; i++) {
03426       std::map<GlobalID,int>::const_iterator
03427         iter = elemIDList.find(elemIDs[i]);
03428       if (iter == elemIDList.end()) continue;
03429       int index = iter->second;
03430 
03431       int offset = i*numElemDOFPerElement;
03432 
03433       for(int j=0; j<numElemDOFPerElement; j++) {
03434          int eqn = elemDOFEqnNumbers[index] + j;
03435 
03436          CHK_ERR( getEqnSolnEntry(eqn, answer) )
03437 
03438          results[offset+j] = answer;
03439       }
03440    }
03441 
03442    return(FEI_SUCCESS);
03443 }
03444 
03445 //------------------------------------------------------------------------------
03446 int LinSysCoreFilter::putBlockElemSolution(GlobalID elemBlockID,
03447                                    int numElems,
03448                                    const GlobalID *elemIDs,
03449                                    int dofPerElem,
03450                                    const double *estimates)
03451 {
03452    debugOutput("FEI: putBlockElemSolution");
03453 
03454    BlockDescriptor* block = NULL;
03455    CHK_ERR( problemStructure_->getBlockDescriptor(elemBlockID, block) )
03456 
03457    int DOFPerElement = block->getNumElemDOFPerElement();
03458    assert(DOFPerElement == dofPerElem);
03459    if (DOFPerElement <= 0) return(0);
03460 
03461    ConnectivityTable& ctable =
03462      problemStructure_->getBlockConnectivity(elemBlockID);
03463    std::map<GlobalID,int>& elemIDList = ctable.elemIDs;
03464 
03465    std::vector<int>& elemDOFEqnNumbers = block->elemDOFEqnNumbers();
03466 
03467 
03468    for(int i=0; i<numElems; i++) {
03469      std::map<GlobalID,int>::const_iterator
03470        iter = elemIDList.find(elemIDs[i]);
03471      if (iter == elemIDList.end()) continue;
03472 
03473      int index = iter->second;
03474 
03475       for(int j=0; j<DOFPerElement; j++) {
03476          int reducedEqn;
03477          problemStructure_->
03478            translateToReducedEqn(elemDOFEqnNumbers[index] + j, reducedEqn);
03479          double soln = estimates[i*DOFPerElement + j];
03480 
03481          CHK_ERR( lsc_->putInitialGuess(&reducedEqn, &soln, 1) );
03482       }
03483    }
03484 
03485    return(FEI_SUCCESS);
03486 }
03487 
03488 //------------------------------------------------------------------------------
03489 int LinSysCoreFilter::getCRMultipliers(int numCRs,
03490                                        const int* CRIDs,
03491                                        double* multipliers)
03492 {
03493   int multCRsLen = problemStructure_->getNumMultConstRecords();
03494   if (numCRs > multCRsLen) {
03495     return(-1);
03496   }
03497 
03498   std::map<GlobalID, ConstraintType*>::const_iterator
03499     cr_iter = problemStructure_->getMultConstRecords().begin(),
03500     cr_end  = problemStructure_->getMultConstRecords().end();
03501 
03502   int i = 0;
03503   while(cr_iter != cr_end && i < numCRs) {
03504     GlobalID CRID = (*cr_iter).first;
03505     ConstraintType* multCR = (*cr_iter).second;
03506     if (CRID != CRIDs[i]) {
03507       CHK_ERR( problemStructure_->getMultConstRecord(CRIDs[i], multCR) );
03508     }
03509 
03510     int eqn = multCR->getEqnNumber();
03511 
03512     CHK_ERR( getEqnSolnEntry(eqn, multipliers[i]) );
03513     ++cr_iter;
03514     ++i;
03515   }
03516 
03517   return(FEI_SUCCESS);
03518 }
03519 
03520 //------------------------------------------------------------------------------
03521 int LinSysCoreFilter::putCRMultipliers(int numMultCRs,
03522                                const int* CRIDs,
03523                                const double *multEstimates)
03524 {
03525   debugOutput("FEI: putCRMultipliers");
03526 
03527   for(int j = 0; j < numMultCRs; j++) {
03528     ConstraintType* multCR = NULL;
03529     CHK_ERR( problemStructure_->getMultConstRecord(CRIDs[j], multCR) );
03530 
03531     int eqnNumber = multCR->getEqnNumber();
03532     if (eqnNumber < localStartRow_ || eqnNumber > localEndRow_) continue;
03533     CHK_ERR( lsc_->putInitialGuess(&eqnNumber, &(multEstimates[j]), 1));
03534   }
03535 
03536   return(FEI_SUCCESS);
03537 }
03538 
03539 //------------------------------------------------------------------------------
03540 int LinSysCoreFilter::putNodalFieldData(int fieldID,
03541                                         int numNodes,
03542                                         const GlobalID* nodeIDs,
03543                                         const double* nodeData)
03544 {
03545   debugOutput("FEI: putNodalFieldData");
03546 
03547   if (fieldID > -1) {
03548     return(putNodalFieldSolution(fieldID, numNodes, nodeIDs, nodeData));
03549   }
03550 
03551   int fieldSize = problemStructure_->getFieldSize(fieldID);
03552   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03553 
03554   std::vector<int> nodeNumbers(numNodes);
03555 
03556   for(int i=0; i<numNodes; i++) {
03557     const NodeDescriptor* node = NULL;
03558     CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
03559 
03560     int nodeNumber = node->getNodeNumber();
03561     if (nodeNumber < 0) {
03562       fei::console_out() << "LinSysCoreFilter::putNodalFieldData ERROR, node with ID " 
03563            << (int)nodeIDs[i] << " doesn't have an associated nodeNumber "
03564            << "assigned. putNodalFieldData shouldn't be called until after the "
03565            << "initComplete method has been called." << FEI_ENDL;
03566       ERReturn(-1);
03567     }
03568 
03569     nodeNumbers[i] = nodeNumber;
03570   }
03571 
03572   CHK_ERR( lsc_->putNodalFieldData(fieldID, fieldSize,
03573                                    &nodeNumbers[0], numNodes, nodeData) );
03574 
03575   return(0);
03576 }
03577 
03578 //------------------------------------------------------------------------------
03579 int LinSysCoreFilter::putNodalFieldSolution(int fieldID,
03580                                 int numNodes,
03581                                 const GlobalID* nodeIDs,
03582                                 const double* nodeData)
03583 {
03584   debugOutput("FEI: putNodalFieldSolution");
03585 
03586   if (fieldID < 0) {
03587     return(putNodalFieldData(fieldID, numNodes, nodeIDs, nodeData));
03588   }
03589 
03590   int fieldSize = problemStructure_->getFieldSize(fieldID);
03591   NodeDatabase& nodeDB = problemStructure_->getNodeDatabase();
03592 
03593   std::vector<int> eqnNumbers(fieldSize);
03594 
03595   for(int i=0; i<numNodes; i++) {
03596     const NodeDescriptor* node = NULL;
03597     CHK_ERR( nodeDB.getNodeWithID(nodeIDs[i], node) );
03598 
03599     int eqn = -1;
03600     bool hasField = node->getFieldEqnNumber(fieldID, eqn);
03601     if (!hasField) continue;
03602 
03603     int reducedEqn = -1;
03604     bool isSlave = problemStructure_->translateToReducedEqn(eqn, reducedEqn);
03605     if (isSlave) continue;
03606 
03607     if (reducedStartRow_ > reducedEqn || reducedEndRow_ < reducedEqn) continue;
03608 
03609     int localLen = fieldSize;
03610     for(int j=0; j<fieldSize; j++) {
03611       int thisEqn = reducedEqn+j;
03612       if (reducedStartRow_ > thisEqn || reducedEndRow_ <thisEqn) {
03613         localLen = j;
03614       }
03615 
03616       eqnNumbers[j] = reducedEqn+j;
03617     }
03618 
03619     int offset = i*fieldSize;
03620     CHK_ERR( lsc_->putInitialGuess(&eqnNumbers[0],
03621                                    &nodeData[offset], localLen) );
03622   }
03623 
03624   return(0);
03625 }
03626 
03627 //------------------------------------------------------------------------------
03628 int LinSysCoreFilter::assembleEqns(int numPtRows, 
03629                                    int numPtCols,
03630                                    const int* rowNumbers,
03631                                    const int* colIndices,
03632                                    const double* const* coefs,
03633                                    bool structurallySymmetric,
03634                                    int numBlkEqns, int* blkEqns,
03635                                    int* blkSizes, bool useBlkEqns,
03636                                    int mode)
03637 {
03638   if (numPtRows == 0) return(FEI_SUCCESS);
03639 
03640   bool anySlaves = false;
03641   int numSlaveEqns = problemStructure_->numSlaveEquations();
03642   if (numSlaveEqns > 0) {
03643     rSlave_.resize(numPtRows);
03644     cSlave_.resize(0);
03645     const int* indPtr = colIndices;
03646     for(int r=0; r<numPtRows; r++) {
03647       rSlave_[r] = problemStructure_->isSlaveEqn(rowNumbers[r]) ? 1 : 0;
03648       if (rSlave_[r] == 1) anySlaves = true;
03649 
03650       for(int j=0; j<numPtCols; j++) {
03651         int isSlave = problemStructure_->isSlaveEqn(indPtr[j]) ? 1 : 0;
03652         cSlave_.push_back(isSlave);
03653         if (isSlave == 1) anySlaves = true;
03654       }
03655 
03656       if (!structurallySymmetric) indPtr += numPtCols;
03657     }
03658   }
03659 
03660   if (numSlaveEqns == 0 || !anySlaves) {
03661     if (numSlaveEqns == 0 && structurallySymmetric) {
03662       if (useBlkEqns) {
03663         CHK_ERR( giveToBlkMatrix_symm_noSlaves(numPtRows, rowNumbers,
03664                                                numBlkEqns, blkEqns, blkSizes,
03665                                                coefs, mode) );
03666       }
03667       else {
03668         CHK_ERR( giveToMatrix_symm_noSlaves(numPtRows, rowNumbers, coefs, mode) );
03669       }
03670     }
03671     else {
03672       if ((int)dworkSpace2_.size() < numPtRows) {
03673         dworkSpace2_.resize(numPtRows);
03674       }
03675       const double* * coefPtr = &dworkSpace2_[0];
03676       for(int i=0; i<numPtRows; i++) {
03677         coefPtr[i] = coefs[i];
03678       }
03679 
03680       if (structurallySymmetric) {
03681         CHK_ERR( giveToMatrix(numPtRows, rowNumbers, numPtRows, rowNumbers,
03682                               coefPtr, mode) );
03683       }
03684       else {
03685         const int* indPtr = colIndices;
03686         for(int i=0; i<numPtRows; i++) {
03687           int row = rowNumbers[i];
03688 
03689           const double* coefPtr1 = coefs[i];
03690 
03691           CHK_ERR(giveToMatrix(1, &row, numPtCols, indPtr, &coefPtr1, mode));
03692           indPtr += numPtCols;
03693         }
03694       }
03695     }
03696   }
03697   else {
03698     int offset = 0;
03699     const int* indicesPtr = colIndices;
03700     for(int i=0; i<numPtRows; i++) {
03701       int row = rowNumbers[i];
03702 
03703       const double* coefPtr = coefs[i];
03704       int* colSlave = &cSlave_[offset];
03705       offset += numPtCols;
03706 
03707       if (rSlave_[i] == 1) {
03708         //Since this is a slave equation, the non-slave columns of this row go
03709         //into 'Kdi_', and the slave columns go into 'Kdd_'.
03710         for(int jj=0; jj<numPtCols; jj++) {
03711           int col = indicesPtr[jj];
03712           if (colSlave[jj]) {
03713             Kdd_->sumInCoef(row, col, coefPtr[jj]);
03714           }
03715           else {
03716             Kdi_->sumInCoef(row, col, coefPtr[jj]);
03717           }
03718         }
03719 
03720         //We also need to put the non-slave rows of column 'row' into 'K_id'.
03721         const int* ii_indicesPtr = colIndices;
03722         for(int ii=0; ii<numPtRows; ii++) {
03723           int rowi = rowNumbers[ii];
03724           if (rSlave_[ii] == 1) continue;
03725 
03726           int index = fei::binarySearch(row, ii_indicesPtr, numPtCols);
03727           if (index < 0) continue;
03728 
03729           const double* coefs_ii = coefs[ii];
03730 
03731           Kid_->sumInCoef(rowi, row, coefs_ii[index]);
03732 
03733           if (!structurallySymmetric) ii_indicesPtr += numPtCols;
03734         }
03735 
03736         reducedEqnCounter_++;
03737 
03738         continue;
03739       }
03740       else {//row is not a slave eqn...
03741 
03742         //put all non-slave columns from this row into the assembled matrix.
03743 
03744         CHK_ERR( giveToMatrix(1, &row, numPtCols, indicesPtr, &coefPtr, mode) );
03745       }
03746 
03747       if (!structurallySymmetric) indicesPtr += numPtCols;
03748     }
03749 
03750     if (reducedEqnCounter_ > 300) CHK_ERR( assembleReducedEqns() );
03751   }
03752 
03753   return(FEI_SUCCESS);
03754 }
03755 
03756 //------------------------------------------------------------------------------
03757 int LinSysCoreFilter::assembleReducedEqns()
03758 {
03759   fei::FillableMat* D = problemStructure_->getSlaveDependencies();
03760 
03761   csrD = *D;
03762   csrKid = *Kid_;
03763   csrKdi = *Kdi_;
03764   csrKdd = *Kdd_;
03765 
03766   //form tmpMat1_ = Kid_*D
03767   fei::multiply_CSRMat_CSRMat(csrKid, csrD, tmpMat1_);
03768 
03769   //form tmpMat2_ = D^T*Kdi_
03770   fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdi, tmpMat2_);
03771 
03772   if (Filter::logStream() != NULL) {
03773     FEI_OSTREAM& os = *Filter::logStream();
03774     os << "#  tmpMat1_"<<FEI_ENDL << tmpMat1_ << FEI_ENDL;
03775     os << "#  tmpMat2_"<<FEI_ENDL << tmpMat2_ << FEI_ENDL;
03776   }
03777 
03778   //accumulate the above two results into the global system matrix.
03779   CHK_ERR( sumIntoMatrix(tmpMat1_) );
03780   CHK_ERR( sumIntoMatrix(tmpMat2_) );
03781 
03782   //form tmpMat1_ = D^T*Kdd_
03783   fei::multiply_trans_CSRMat_CSRMat(csrD, csrKdd, tmpMat1_);
03784 
03785   //form tmpMat2_ = tmpMat1_*D = D^T*Kdd_*D
03786   fei::multiply_CSRMat_CSRMat(tmpMat1_, csrD, tmpMat2_);
03787 
03788   if (Filter::logStream() != NULL) {
03789     FEI_OSTREAM& os = *Filter::logStream();
03790     os << "#  tmpMat2_"<<FEI_ENDL << tmpMat2_ << FEI_ENDL;
03791   }
03792 
03793   //finally, accumulate tmpMat2_ = D^T*Kdd_*D into the global system matrix.
03794   CHK_ERR( sumIntoMatrix(tmpMat2_) );
03795 
03796   Kdi_->clear();
03797   Kid_->clear();
03798   Kdd_->clear();
03799   reducedEqnCounter_ = 0;
03800 
03801   return(0);
03802 }
03803 
03804 //------------------------------------------------------------------------------
03805 int LinSysCoreFilter::assembleRHS(int numValues,
03806                                   const int* indices,
03807                                   const double* coefs,
03808                                   int mode) {
03809 //
03810 //This function hands the data off to the routine that finally
03811 //sticks it into the RHS vector.
03812 //
03813 
03814   if (problemStructure_->numSlaveEquations() == 0) {
03815     CHK_ERR( giveToRHS(numValues, coefs, indices, mode) );
03816     return(FEI_SUCCESS);
03817   }
03818 
03819   for(int i = 0; i < numValues; i++) {
03820     int eqn = indices[i];
03821 
03822     if (problemStructure_->isSlaveEqn(eqn)) {
03823       fei::add_entry( fd_, eqn, coefs[i]);
03824       reducedRHSCounter_++;
03825       continue;
03826     }
03827 
03828     CHK_ERR( giveToRHS(1, &(coefs[i]), &eqn, mode ) );
03829   }
03830 
03831   if (reducedRHSCounter_ > 300) CHK_ERR( assembleReducedRHS() );
03832 
03833   return(FEI_SUCCESS);
03834 }
03835 
03836 //------------------------------------------------------------------------------
03837 int LinSysCoreFilter::assembleReducedRHS()
03838 {
03839   fei::FillableMat* D = problemStructure_->getSlaveDependencies();
03840 
03841   csrD = *D;
03842 
03843   //now form tmpVec1_ = D^T*fd_.
03844   fei::multiply_trans_CSRMat_CSVec(csrD, fd_, tmpVec1_);
03845 
03846   CHK_ERR( sumIntoRHS(tmpVec1_) );
03847 
03848   fd_.clear();
03849   reducedRHSCounter_ = 0;
03850 
03851   return(0);
03852 }
03853 
03854 //==============================================================================
03855 void LinSysCoreFilter::debugOutput(const char* mesg) {
03856    if (Filter::logStream() != NULL) {
03857      (*logStream())<<mesg<<FEI_ENDL;
03858    }
03859 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends