00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
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
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
00111
00112
00113
00114
00115 Channel *theChannel = theChannels[0];
00116 theChannel->recvID(0, 0, otherSize);
00117 numProcesses = otherSize(0);
00118
00119
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
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
00143 otherSize(0) = myDOFsShared.Size();
00144 theChannel->sendID(0, 0, otherSize);
00145 if (size != 0)
00146 theChannel->sendID(0, 0, myDOFsShared);
00147
00148
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 {
00158
00159
00160
00161
00162
00163
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
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
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
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
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
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
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];
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
00268 for (int l=0; l<size; l++) {
00269 A[l] = 0;
00270 B[l] = 0;
00271 X[l] = 0;
00272 }
00273
00274
00275
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
00298 FE_EleIter &theEle = theModel->getFEs();
00299 FE_Element *elePtr;
00300 while ((elePtr = theEle()) != 0)
00301 elePtr->setID();
00302 }
00303
00304
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
00322 if (fact == 0.0) return 0;
00323
00324 #ifdef _G3DEBUG
00325
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) {
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) {
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
00365 if (fact == 0.0) return 0;
00366
00367 #ifdef _G3DEBUG
00368
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) {
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) {
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
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) {
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
00519
00520
00521
00522
00523 if (processID == 0) {
00524
00525
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
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
00564 sendID = numChannels;
00565 }
00566
00567 } else
00568 sendID = processID;
00569
00570
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 }