DistributedDiagonalSOE.cpp

Go to the documentation of this file.
00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.3 $
00022 // $Date: 2006/01/25 01:44:00 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/diagonal/DistributedDiagonalSOE.cpp,v $
00024 
00025 // Written: fmk 
00026 // Created: 05/05
00027 //
00028 // Description: This file contains the implementation for DistributedDiagonalSOE
00029 
00030 
00031 #include <DistributedDiagonalSOE.h>
00032 #include <DistributedDiagonalSolver.h>
00033 #include <Matrix.h>
00034 #include <Graph.h>
00035 #include <Vertex.h>
00036 #include <VertexIter.h>
00037 #include <f2c.h>
00038 #include <math.h>
00039 
00040 #include <Channel.h>
00041 #include <FEM_ObjectBroker.h>
00042 
00043 #include <AnalysisModel.h>
00044 #include <DOF_Group.h>
00045 #include <FE_Element.h>
00046 #include <FE_EleIter.h>
00047 #include <DOF_GrpIter.h>
00048 
00049 
00050 DistributedDiagonalSOE::DistributedDiagonalSOE(DistributedDiagonalSolver &the_Solver)
00051 :LinearSOE(the_Solver, LinSOE_TAGS_DistributedDiagonalSOE),
00052  size(0), A(0), B(0), X(0), vectX(0), vectB(0), isAfactored(false),
00053  processID(0), numProcesses(0),
00054  numChannels(0), theChannels(0), localCol(0), 
00055  myDOFs(0,32), myDOFsShared(0,16), numShared(0), dataShared(0), vectShared(0), theModel(0)
00056 {
00057     the_Solver.setLinearSOE(*this);
00058 }
00059 
00060 
00061 DistributedDiagonalSOE::~DistributedDiagonalSOE()
00062 {
00063   if (A != 0) delete [] A;
00064   if (B != 0) delete [] B;
00065   if (X != 0) delete [] X;
00066   if (dataShared != 0) delete [] dataShared;
00067   if (vectX != 0) delete vectX;    
00068   if (vectB != 0) delete vectB;    
00069   if (vectShared != 0) delete vectShared;    
00070 
00071   if (theChannels != 0)
00072     delete [] theChannels;
00073 }
00074 
00075 
00076 int 
00077 DistributedDiagonalSOE::getNumEqn(void) const
00078 {
00079   return size;
00080 }
00081 
00082 
00083 int 
00084 DistributedDiagonalSOE::setSize(Graph &theGraph)
00085 {
00086   int result = 0;
00087   size = theGraph.getNumVertex();
00088 
00089   //
00090   // first we build an ID containing all local DOFs
00091   //
00092 
00093   myDOFs.resize(size);
00094 
00095   int count = 0;
00096   Vertex *theVertex;
00097   VertexIter &theVertices = theGraph.getVertices();
00098   while ((theVertex = theVertices()) != 0) {
00099     int vertexTag = theVertex->getTag();
00100     myDOFs(count) = vertexTag;
00101     count++;
00102   }
00103 
00104   static ID otherSize(1);
00105   ID otherDOFS(0, size/10);
00106 
00107   if (processID != 0) {
00108 
00109     //
00110     // each process send it's local and receives all other processes IDs (remote IDs).
00111     // from the remote IDs the local process figures out which DOF's it shares with others.
00112     // it then sends the shared dofs' out.
00113     //
00114 
00115     Channel *theChannel = theChannels[0];
00116     theChannel->recvID(0, 0, otherSize);
00117     numProcesses = otherSize(0);
00118 
00119     // send local
00120     otherSize(0) = size;
00121     theChannel->sendID(0, 0, otherSize);
00122     if (size != 0)
00123       theChannel->sendID(0, 0, myDOFs);
00124 
00125     for (int i=0; i<numProcesses-1; i++) {
00126       // receive remote & check for shared DOFs
00127       theChannel->recvID(0, 0, otherSize);      
00128       int numOther = otherSize(0);
00129 
00130       if (numOther != 0) {
00131         otherDOFS.resize(numOther);
00132         theChannel->recvID(0, 0, otherDOFS);
00133         
00134         for (int j=0; j<numOther; j++) {
00135           int otherTag = otherDOFS(j);
00136           if (myDOFs.getLocation(otherTag) != -1 && myDOFsShared.getLocation(otherTag) == -1) 
00137             myDOFsShared[numShared++] = otherTag;
00138         }
00139       }
00140     }
00141 
00142     // send my shared DOFs
00143     otherSize(0) = myDOFsShared.Size();
00144     theChannel->sendID(0, 0, otherSize);
00145     if (size != 0)
00146       theChannel->sendID(0, 0, myDOFsShared);
00147 
00148     // recv all shared DOFs
00149     theChannel->recvID(0, 0, otherSize);        
00150     numShared = otherSize(0);
00151 
00152     if (numShared != 0) {
00153       myDOFsShared.resize(numShared);
00154       theChannel->recvID(0, 0, myDOFsShared);
00155     } 
00156 
00157   } else { // in absence of a broadcast we need the following
00158 
00159     //
00160     // each process send it's local and receives all other processes IDs (remote IDs).
00161     // from the remote IDs the local process figures out which DOF's it shares with others.
00162     // it then sends the shared dofs' out. lastly receives from P0 list of all shared.
00163     // NOTE:  if all sent to P0 their ID's; P0 could figure out shared; but scalability issues
00164     //
00165 
00166     for (int j=0; j<numChannels; j++) {
00167       Channel *theChannel = theChannels[j];
00168       otherSize(0) = numChannels+1;
00169       theChannel->sendID(0, 0, otherSize);
00170 
00171       otherSize(0) = size;
00172       theChannel->sendID(0, 0, otherSize);
00173       if (size != 0)
00174         theChannel->sendID(0, 0, myDOFs);
00175     }
00176 
00177     for (int j=0; j<numChannels; j++) {
00178       Channel *theChannel = theChannels[j];
00179 
00180       // receive remote & check for shared DOFs
00181 
00182       theChannel->recvID(0, 0, otherSize);      
00183       int numOther = otherSize(0);
00184       if (numOther != 0) {
00185         otherDOFS.resize(numOther);
00186         theChannel->recvID(0, 0, otherDOFS);
00187       }
00188 
00189       // need to send to all others
00190       for (int k=0; k<numChannels; k++) {
00191         Channel *theChannel = theChannels[k];
00192         if (k != j) {
00193           theChannel->sendID(0, 0, otherSize);  
00194           if (numOther != 0) {
00195             theChannel->sendID(0, 0, otherDOFS);            
00196           }
00197         }
00198       }
00199       
00200       // need to merge with mine 
00201       for (int l=0; l<numOther; l++) {
00202         int otherTag = otherDOFS(l);
00203         if (myDOFs.getLocation(otherTag) != -1 && myDOFsShared.getLocation(otherTag) == -1) { 
00204           myDOFsShared[numShared++] = otherTag;
00205         }
00206       }
00207     }
00208 
00209     // now recv each of the shared DOFs & merge with mine to form the master list
00210     for (int j=0; j<numChannels; j++) {
00211       Channel *theChannel = theChannels[j];
00212       theChannel->recvID(0, 0, otherSize);      
00213       int numOther = otherSize(0);
00214       if (numOther != 0) {
00215         otherDOFS.resize(numOther);
00216         theChannel->recvID(0, 0, otherDOFS);
00217       } 
00218       // need to merge with mine 
00219       for (int k=0; k<numOther; k++) {
00220         int otherTag = otherDOFS(k);
00221         if (myDOFsShared.getLocation(otherTag) == -1) 
00222           myDOFsShared(numShared++) = otherTag;
00223       }
00224     }
00225 
00226     // now send the shared ID
00227     for (int j=0; j<numChannels; j++) {
00228       Channel *theChannel = theChannels[j];
00229       otherSize(0) = numShared;
00230       theChannel->sendID(0, 0, otherSize);      
00231       theChannel->sendID(0, 0, myDOFsShared);
00232     }
00233   }
00234 
00235 
00236   if (A != 0) delete [] A; A = 0;
00237   if (B != 0) delete [] B; B = 0;
00238   if (X != 0) delete [] X; X = 0;
00239   if (dataShared != 0) delete [] dataShared; dataShared = 0;
00240   if (vectX != 0) delete vectX; vectX = 0;
00241   if (vectB != 0) delete vectB; vectB = 0;
00242 
00243   A = new double[size];
00244   B = new double[size];
00245   X = new double[size]; 
00246   if (X != 0) 
00247     vectX = new Vector(X,size);
00248   if (B != 0) 
00249     vectB = new Vector(B,size);
00250   
00251   dataShared = new double[2*numShared];  // 2 times for A & B
00252   if (dataShared != 0)
00253     vectShared = new Vector(dataShared, 2*numShared);
00254 
00255 
00256   if (A == 0 || B == 0 || X == 0 || vectX == 0 || vectB == 0 || dataShared == 0 || vectShared == 0) {
00257     opserr << "ERROR DistributedDiagonalSOE::setSize() - ";
00258     opserr << " ran out of memory for size: " << size << endln;
00259     if (A != 0) delete [] A;
00260     if (B != 0) delete [] B;
00261     if (X != 0) delete [] X;
00262     if (dataShared != 0) delete [] dataShared;
00263     size = 0;
00264     return -1;
00265   }
00266 
00267   // zero the vectors
00268   for (int l=0; l<size; l++) {
00269     A[l] = 0;
00270     B[l] = 0;
00271     X[l] = 0;
00272   }
00273 
00274   //
00275   // now let's redo the mapping of the dof's .. locally numbered 0 through size
00276   //
00277 
00278   
00279   if (theModel == 0) {
00280     opserr << "WARNING DistributedDiagonalSOE::setSize - no AnalysisModel\n";
00281   } 
00282   else {
00283     DOF_GrpIter &theDOFs = theModel->getDOFs();
00284     DOF_Group *dofPtr;
00285     
00286     while ((dofPtr = theDOFs()) != 0) {
00287         const ID &theID = dofPtr->getID();
00288         for (int i=0; i<theID.Size(); i++) {
00289           int dof = theID(i);
00290           if (dof >= 0) {
00291             int newDOF = myDOFs.getLocation(dof);
00292             dofPtr->setID(i, newDOF);
00293           }
00294         }   
00295     }
00296     
00297     // iterate through the FE_Element getting them to set their IDs
00298     FE_EleIter &theEle = theModel->getFEs();
00299     FE_Element *elePtr;
00300     while ((elePtr = theEle()) != 0)
00301         elePtr->setID(); 
00302   }  
00303 
00304   // invoke setSize() on the Solver
00305   DistributedDiagonalSolver *the_Solver = (DistributedDiagonalSolver *)this->getSolver();
00306   the_Solver->setLinearSOE(*this);
00307   int solverOK = the_Solver->setSize();
00308   
00309   if (solverOK < 0) {
00310     opserr << "WARNING DistributedDiagonalSOE::setSize :";
00311     opserr << " solver failed setSize()\n";
00312     return solverOK;
00313   }    
00314 
00315   return result;
00316 }
00317 
00318 int 
00319 DistributedDiagonalSOE::addA(const Matrix &m, const ID &id, double fact)
00320 {
00321   // check for a quick return 
00322   if (fact == 0.0)  return 0;
00323   
00324 #ifdef _G3DEBUG
00325   // check that m and id are of similar size
00326   int idSize = id.Size();    
00327   if (idSize != m.noRows() && idSize != m.noCols()) {
00328     opserr << "FullGenLinSOE::addA()    - Matrix and ID not of similar sizes\n";
00329     return -1;
00330   }
00331 #endif
00332 
00333   if (fact == 1.0) { // do not need to multiply if fact == 1.0
00334     for (int i=0; i<id.Size(); i++) {
00335       int pos = id(i);
00336       if (pos <size && pos >= 0) {
00337         A[pos] += m(i,i);
00338       }
00339     }
00340   } else if (fact == -1.0) { // do not need to multiply if fact == -1.0
00341     for (int i=0; i<id.Size(); i++) {
00342       int pos = id(i);
00343       if (pos <size && pos >= 0) {
00344         A[pos] -= m(i,i);
00345       }
00346     }
00347   } else {
00348     for (int i=0; i<id.Size(); i++) {
00349       int pos = id(i);
00350       if (pos <size && pos >= 0) {
00351         A[pos] += m(i,i) * fact;
00352       }
00353     }
00354   }     
00355 
00356   return 0;
00357 }
00358  
00359     
00360 int 
00361 DistributedDiagonalSOE::addB(const Vector &v, const ID &id, double fact)
00362 {
00363     
00364   // check for a quick return 
00365   if (fact == 0.0)  return 0;
00366   
00367 #ifdef _G3DEBUG
00368   // check that m and id are of similar size
00369   int idSize = id.Size();        
00370   if (idSize != v.Size() ) {
00371     opserr << "DistributedDiagonalSOE::addB() -";
00372     opserr << " Vector and ID not of similar sizes\n";
00373     return -1;
00374   }    
00375 #endif
00376   
00377   if (fact == 1.0) { // do not need to multiply if fact == 1.0
00378     for (int i=0; i<id.Size(); i++) {
00379       int pos = id(i);
00380       if (pos <size && pos >= 0)
00381         B[pos] += v(i);
00382     }
00383   } else if (fact == -1.0) { // do not need to multiply if fact == -1.0
00384     for (int i=0; i<id.Size(); i++) {
00385       int pos = id(i);
00386       if (pos <size && pos >= 0)
00387         B[pos] -= v(i);
00388     }
00389   } else {
00390     for (int i=0; i<id.Size(); i++) {
00391       int pos = id(i);
00392       if (pos <size && pos >= 0)
00393         B[pos] += v(i) * fact;
00394     }
00395   }     
00396   return 0;
00397 }
00398 
00399 
00400 int
00401 DistributedDiagonalSOE::setB(const Vector &v, double fact)
00402 {
00403   // check for a quick return 
00404   if (fact == 0.0)  return 0;
00405   
00406   if (v.Size() != size) {
00407     opserr << "WARNING DistributedDiagonalSOE::setB() -";
00408     opserr << " incomptable sizes " << size << " and " << v.Size() << endln;
00409     return -1;
00410   }
00411   
00412   if (fact == 1.0) { // do not need to multiply if fact == 1.0
00413     for (int i=0; i<size; i++) {
00414       B[i] = v(i);
00415     }
00416   } else if (fact == -1.0) {
00417     for (int i=0; i<size; i++) {
00418       B[i] = -v(i);
00419     }
00420   } else {
00421     for (int i=0; i<size; i++) {
00422       B[i] = v(i) * fact;
00423     }
00424   }     
00425   return 0;
00426 }
00427 
00428 void 
00429 DistributedDiagonalSOE::zeroA(void)
00430 {
00431   double *Aptr = A;
00432   for (int i=0; i<size; i++)
00433     *Aptr++ = 0;
00434   
00435   isAfactored = false;
00436 }
00437 
00438 void 
00439 DistributedDiagonalSOE::zeroB(void)
00440 {
00441   double *Bptr = B;
00442   for (int i=0; i<size; i++)
00443     *Bptr++ = 0;
00444 }
00445 
00446 
00447 void 
00448 DistributedDiagonalSOE::setX(int loc, double value)
00449 {
00450   if (loc < size && loc >=0)
00451     X[loc] = value;
00452 }
00453 
00454 void 
00455 DistributedDiagonalSOE::setX(const Vector &x)
00456 {
00457   if (x.Size() == size && vectX != 0)
00458     *vectX = x;
00459 }
00460 
00461 const Vector &
00462 DistributedDiagonalSOE::getX(void)
00463 {
00464   if (vectX == 0) {
00465     opserr << "FATAL DistributedDiagonalSOE::getX - vectX == 0";
00466     exit(-1);
00467   }    
00468   return *vectX;
00469 }
00470 
00471 const Vector &
00472 DistributedDiagonalSOE::getB(void)
00473 {
00474   if (vectB == 0) {
00475     opserr << "FATAL DistributedDiagonalSOE::getB - vectB == 0";
00476     exit(-1);
00477   }        
00478   return *vectB;
00479 }
00480 
00481 double 
00482 DistributedDiagonalSOE::normRHS(void)
00483 {
00484   double norm =0.0;
00485   for (int i=0; i<size; i++) {
00486     double Yi = B[i];
00487     norm += Yi*Yi;
00488   }
00489   return sqrt(norm);
00490   
00491 }    
00492 
00493 
00494 int
00495 DistributedDiagonalSOE::setDiagonalSolver(DistributedDiagonalSolver &newSolver)
00496 {
00497   newSolver.setLinearSOE(*this);
00498   
00499   if (size != 0) {
00500     int solverOK = newSolver.setSize();
00501     if (solverOK < 0) {
00502       opserr << "WARNING:DistributedDiagonalSOE::setSolver :";
00503       opserr << "the new solver could not setSeize() - staying with old\n";
00504       return -1;
00505     }
00506   }
00507   
00508   return this->setSolver(newSolver);
00509 }
00510 
00511 
00512 int 
00513 DistributedDiagonalSOE::sendSelf(int cTag, Channel &theChannel)
00514 {
00515   processID = 0;
00516   int sendID =0;
00517   
00518   // if P0 check if already sent. If already sent use old processID; if not allocate a new process 
00519   // id for remote part of object, enlarge channel * to hold a channel * for this remote object.
00520   
00521   // if not P0, send current processID
00522   
00523   if (processID == 0) {
00524 
00525     // check if already using this object
00526     bool found = false;
00527     for (int i=0; i<numChannels; i++)
00528       if (theChannels[i] == &theChannel) {
00529         sendID = i+1;
00530         found = true;
00531       }
00532 
00533     // if new object, enlarge Channel pointers to hold new channel * & allocate new ID
00534     if (found == false) {
00535       int nextNumChannels = numChannels + 1;
00536       Channel **nextChannels = new Channel *[nextNumChannels];
00537       if (nextNumChannels == 0) {
00538         opserr << "DistributedBandGenLinSOE::sendSelf() - failed to allocate channel array of size: " << 
00539           nextNumChannels << endln;
00540         return -1;
00541       }
00542       for (int i=0; i<numChannels; i++)
00543         nextChannels[i] = theChannels[i];
00544       nextChannels[numChannels] = &theChannel;
00545       numChannels = nextNumChannels;
00546       
00547       if (theChannels != 0)
00548         delete [] theChannels;
00549       
00550       theChannels = nextChannels;
00551       
00552   if (localCol != 0)
00553         delete [] localCol;
00554       localCol = new ID *[numChannels];
00555       if (localCol == 0) {
00556         opserr << "DistributedBandGenLinSOE::sendSelf() - failed to allocate id array of size: " << 
00557           nextNumChannels << endln;
00558         return -1;
00559       }
00560       for (int i=0; i<numChannels; i++)
00561         localCol[i] = 0;    
00562 
00563       // allocate new processID for remote object
00564       sendID = numChannels;
00565     }
00566 
00567   } else 
00568     sendID = processID;
00569 
00570   // send remotes processID
00571   ID idData(1);
00572   idData(0) = sendID;
00573   
00574   int res = theChannel.sendID(0, cTag, idData);
00575   if (res < 0) {
00576     opserr <<"WARNING DistributedDiagonalSOE::sendSelf() - failed to send data\n";
00577     return -1;
00578   }
00579 
00580   return 0;
00581 }
00582 
00583 
00584 int 
00585 DistributedDiagonalSOE::recvSelf(int cTag, Channel &theChannel, 
00586                    FEM_ObjectBroker &theBroker)
00587 {
00588   ID idData(1);
00589   int res = theChannel.recvID(0, cTag, idData);
00590   if (res < 0) {
00591     opserr <<"WARNING DistributedProfileSPDLinSOE::recvSelf() - failed to send data\n";
00592     return -1;
00593   }           
00594   processID = idData(0);
00595 
00596   numChannels = 1;
00597   theChannels = new Channel *[1];
00598   theChannels[0] = &theChannel;
00599 
00600   localCol = new ID *[numChannels];
00601   for (int i=0; i<numChannels; i++)
00602     localCol[i] = 0;
00603 
00604   return 0;
00605 }
00606 
00607 int
00608 DistributedDiagonalSOE::setChannels(int nChannels, Channel **theC)
00609 {
00610   numChannels = nChannels;
00611 
00612   if (theChannels != 0)
00613     delete [] theChannels;
00614 
00615   theChannels = new Channel *[numChannels];
00616   for (int i=0; i<numChannels; i++)
00617     theChannels[i] = theC[i];
00618 
00619 
00620   localCol = new ID *[nChannels];
00621   for (int i=0; i<numChannels; i++)
00622     localCol[i] = 0;
00623 
00624   return 0;
00625 }
00626 
00627 int 
00628 DistributedDiagonalSOE::setAnalysisModel(AnalysisModel &theAnalysisModel)
00629 {
00630   theModel = &theAnalysisModel;
00631   return 0;
00632 }

Generated on Mon Oct 23 15:05:28 2006 for OpenSees by doxygen 1.5.0