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 <BandGenLinSOE.h>
00032 #include <BandGenLinSolver.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 #include <Channel.h>
00040 #include <FEM_ObjectBroker.h>
00041
00042 BandGenLinSOE::BandGenLinSOE(BandGenLinSolver &theSolvr)
00043 :LinearSOE(theSolvr, LinSOE_TAGS_BandGenLinSOE),
00044 size(0), numSuperD(0), numSubD(0), A(0), B(0), X(0),
00045 vectX(0), vectB(0), Asize(0), Bsize(0), factored(false)
00046 {
00047 theSolvr.setLinearSOE(*this);
00048 }
00049
00050 BandGenLinSOE::BandGenLinSOE(BandGenLinSolver &theSolvr, int classTag)
00051 :LinearSOE(theSolvr, classTag),
00052 size(0), numSuperD(0), numSubD(0), A(0), B(0), X(0),
00053 vectX(0), vectB(0), Asize(0), Bsize(0), factored(false)
00054 {
00055 theSolvr.setLinearSOE(*this);
00056 }
00057
00058
00059 BandGenLinSOE::BandGenLinSOE(int N, int numSuperDiag, int numSubDiag,
00060 BandGenLinSolver &theSolvr)
00061 :LinearSOE(theSolvr, LinSOE_TAGS_BandGenLinSOE),
00062 size(N), numSuperD(numSuperDiag), numSubD(numSubDiag), A(0), B(0),
00063 X(0), vectX(0), vectB(0), Asize(0), Bsize(0), factored(false)
00064 {
00065 Asize = N * (2*numSubD + numSuperD +1);
00066 A = new double[Asize];
00067
00068 if (A == 0) {
00069 opserr << "WARNING BandGenLinSOE::BandGenLinSOE :";
00070 opserr << " ran out of memory for A (size,super,sub) (";
00071 opserr << size <<", " << numSuperDiag << ", " << numSubDiag << ") \n";
00072 Asize = 0; size = 0; numSubD = 0; numSuperD = 0;
00073 } else {
00074
00075 for (int i=0; i<Asize; i++)
00076 A[i] = 0;
00077
00078 B = new double[size];
00079 X = new double[size];
00080
00081 if (B == 0 || X == 0) {
00082 opserr << "WARNING BandGenLinSOE::BandGenLinSOE :";
00083 opserr << " ran out of memory for vectors (size) (";
00084 opserr << size << ") \n";
00085 Bsize = 0; size = 0; numSubD = 0; numSuperD = 0;
00086 } else {
00087 Bsize = size;
00088
00089 for (int j=0; j<size; j++) {
00090 B[j] = 0;
00091 X[j] = 0;
00092 }
00093 }
00094 }
00095
00096 vectX = new Vector(X,size);
00097 vectB = new Vector(B,size);
00098
00099 theSolvr.setLinearSOE(*this);
00100
00101 int solverOK = theSolvr.setSize();
00102 if (solverOK < 0) {
00103 opserr << "WARNING BandGenLinSOE::BandGenLinSOE :";
00104 opserr << " solver failed setSize() in constructor\n";
00105 }
00106 }
00107
00108 int
00109 BandGenLinSOE::getNumEqn(void) const
00110 {
00111 return size;
00112 }
00113
00114 BandGenLinSOE::~BandGenLinSOE()
00115 {
00116 if (A != 0) delete [] A;
00117 if (B != 0) delete [] B;
00118 if (X != 0) delete [] X;
00119 if (vectX != 0) delete vectX;
00120 if (vectB != 0) delete vectB;
00121 }
00122
00123
00124
00125 int
00126 BandGenLinSOE::setSize(Graph &theGraph)
00127 {
00128 int result = 0;
00129 int oldSize = size;
00130 size = theGraph.getNumVertex();
00131
00132
00133
00134
00135
00136 numSubD = 0;
00137 numSuperD = 0;
00138
00139 Vertex *vertexPtr;
00140 VertexIter &theVertices = theGraph.getVertices();
00141
00142 while ((vertexPtr = theVertices()) != 0) {
00143 int vertexNum = vertexPtr->getTag();
00144 const ID &theAdjacency = vertexPtr->getAdjacency();
00145 for (int i=0; i<theAdjacency.Size(); i++) {
00146 int otherNum = theAdjacency(i);
00147 int diff = vertexNum - otherNum;
00148 if (diff > 0) {
00149 if (diff > numSuperD)
00150 numSuperD = diff;
00151 } else
00152 if (diff < numSubD)
00153 numSubD = diff;
00154 }
00155 }
00156 numSubD *= -1;
00157
00158 int newSize = size * (2*numSubD + numSuperD +1);
00159 if (newSize > Asize) {
00160
00161 if (A != 0)
00162 delete [] A;
00163
00164 A = new double[newSize];
00165
00166 if (A == 0) {
00167 opserr << "WARNING BandGenLinSOE::BandGenLinSOE :";
00168 opserr << " ran out of memory for A (size,super,sub) (";
00169 opserr << size <<", " << numSuperD << ", " << numSubD << ") \n";
00170 Asize = 0; size = 0; numSubD = 0; numSuperD = 0;
00171 result= -1;
00172 }
00173 else
00174 Asize = newSize;
00175 }
00176
00177
00178 for (int i=0; i<Asize; i++)
00179 A[i] = 0;
00180
00181 factored = false;
00182
00183 if (size > Bsize) {
00184
00185
00186 if (B != 0) delete [] B;
00187 if (X != 0) delete [] X;
00188
00189
00190 B = new double[size];
00191 X = new double[size];
00192
00193 if (B == 0 || X == 0) {
00194 opserr << "WARNING BandGenLinSOE::BandGenLinSOE :";
00195 opserr << " ran out of memory for vectors (size) (";
00196 opserr << size << ") \n";
00197 Bsize = 0; size = 0; numSubD = 0; numSuperD = 0;
00198 result = -1;
00199 }
00200 else
00201 Bsize = size;
00202 }
00203
00204
00205 for (int j=0; j<size; j++) {
00206 B[j] = 0;
00207 X[j] = 0;
00208 }
00209
00210
00211 if (oldSize != size) {
00212 if (vectX != 0)
00213 delete vectX;
00214
00215 if (vectB != 0)
00216 delete vectB;
00217
00218 vectX = new Vector(X,size);
00219 vectB = new Vector(B,size);
00220 }
00221
00222
00223 LinearSOESolver *theSolvr = this->getSolver();
00224 int solverOK = theSolvr->setSize();
00225 if (solverOK < 0) {
00226 opserr << "WARNING:BandGenLinSOE::setSize :";
00227 opserr << " solver failed setSize()\n";
00228 return solverOK;
00229 }
00230
00231 return result;
00232 }
00233
00234 int
00235 BandGenLinSOE::addA(const Matrix &m, const ID &id, double fact)
00236 {
00237
00238
00239
00240 if (fact == 0.0) return 0;
00241
00242
00243
00244 int idSize = id.Size();
00245 if (idSize != m.noRows() && idSize != m.noCols()) {
00246 opserr << "BandGenLinSOE::addA() - Matrix and ID not of similar sizes\n";
00247 return -1;
00248 }
00249
00250 int ldA = 2*numSubD + numSuperD + 1;
00251
00252
00253 if (fact == 1.0) {
00254 for (int i=0; i<idSize; i++) {
00255 int col = id(i);
00256 if (col < size && col >= 0) {
00257 double *coliiPtr = A + col*ldA + numSubD + numSuperD;
00258 for (int j=0; j<idSize; j++) {
00259 int row = id(j);
00260 if (row <size && row >= 0) {
00261 int diff = col - row;
00262 if (diff > 0) {
00263 if (diff <= numSuperD) {
00264 double *APtr = coliiPtr - diff;
00265 *APtr += m(j,i);
00266 }
00267
00268 } else {
00269 diff *= -1;
00270 if (diff <= numSubD) {
00271 double *APtr = coliiPtr + diff;
00272 *APtr += m(j,i);
00273 }
00274 }
00275 }
00276 }
00277 }
00278 }
00279 } else {
00280 for (int i=0; i<idSize; i++) {
00281 int col = id(i);
00282 if (col < size && col >= 0) {
00283 double *coliiPtr = A + col*ldA + numSubD + numSuperD;
00284 for (int j=0; j<idSize; j++) {
00285 int row = id(j);
00286 if (row <size && row >= 0) {
00287 int diff = col - row;
00288 if (diff > 0) {
00289 if (diff <= numSuperD) {
00290 double *APtr = coliiPtr - diff;
00291 *APtr += m(j,i) *fact;
00292 }
00293 } else {
00294 diff *= -1;
00295 if (diff <= numSubD) {
00296 double *APtr = coliiPtr + diff;
00297 *APtr += m(j,i) *fact;
00298 }
00299 }
00300 }
00301 }
00302 }
00303 }
00304 }
00305
00306
00307
00308
00309 return 0;
00310 }
00311
00312
00313 int
00314 BandGenLinSOE::addB(const Vector &v, const ID &id, double fact)
00315 {
00316
00317 if (fact == 0.0) return 0;
00318
00319
00320
00321 int idSize = id.Size();
00322 if (idSize != v.Size() ) {
00323 opserr << "BandGenLinSOE::addB() - Vector and ID not of similar sizes\n";
00324 return -1;
00325 }
00326
00327 if (fact == 1.0) {
00328 for (int i=0; i<idSize; i++) {
00329 int pos = id(i);
00330 if (pos <size && pos >= 0)
00331 B[pos] += v(i);
00332 }
00333 } else if (fact == -1.0) {
00334 for (int i=0; i<idSize; i++) {
00335 int pos = id(i);
00336 if (pos <size && pos >= 0)
00337 B[pos] -= v(i);
00338 }
00339 } else {
00340 for (int i=0; i<idSize; i++) {
00341 int pos = id(i);
00342 if (pos <size && pos >= 0)
00343 B[pos] += v(i) * fact;
00344 }
00345 }
00346 return 0;
00347 }
00348
00349
00350 int
00351 BandGenLinSOE::setB(const Vector &v, double fact)
00352 {
00353
00354 if (fact == 0.0) return 0;
00355
00356
00357 if (v.Size() != size) {
00358 opserr << "WARNING BandGenLinSOE::setB() -";
00359 opserr << " incomptable sizes " << size << " and " << v.Size() << endln;
00360 return -1;
00361 }
00362
00363 if (fact == 1.0) {
00364 for (int i=0; i<size; i++) {
00365 B[i] = v(i);
00366 }
00367 } else if (fact == -1.0) {
00368 for (int i=0; i<size; i++) {
00369 B[i] = -v(i);
00370 }
00371 } else {
00372 for (int i=0; i<size; i++) {
00373 B[i] = v(i) * fact;
00374 }
00375 }
00376 return 0;
00377 }
00378
00379
00380 void
00381 BandGenLinSOE::zeroA(void)
00382 {
00383 double *Aptr = A;
00384 int theSize = Asize;
00385 for (int i=0; i<theSize; i++)
00386 *Aptr++ = 0;
00387
00388 factored = false;
00389 }
00390
00391 void
00392 BandGenLinSOE::zeroB(void)
00393 {
00394 double *Bptr = B;
00395 for (int i=0; i<size; i++)
00396 *Bptr++ = 0;
00397 }
00398
00399
00400 const Vector &
00401 BandGenLinSOE::getX(void)
00402 {
00403 if (vectX == 0) {
00404 opserr << "FATAL BandGenLinSOE::getX - vectX == 0!";
00405 exit(-1);
00406 }
00407
00408 return *vectX;
00409 }
00410
00411
00412 const Vector &
00413 BandGenLinSOE::getB(void)
00414 {
00415 if (vectB == 0) {
00416 opserr << "FATAL BandGenLinSOE::getB - vectB == 0!";
00417 exit(-1);
00418 }
00419
00420 return *vectB;
00421 }
00422
00423
00424 double
00425 BandGenLinSOE::normRHS(void)
00426 {
00427 double norm =0.0;
00428 double *Bptr = B;
00429 for (int i=0; i<size; i++) {
00430 double Yi = *Bptr++;
00431 norm += Yi*Yi;
00432 }
00433 return sqrt(norm);
00434 }
00435
00436
00437 void
00438 BandGenLinSOE::setX(int loc, double value)
00439 {
00440 if (loc < size && loc >= 0)
00441 X[loc] = value;
00442 }
00443
00444 void
00445 BandGenLinSOE::setX(const Vector &x)
00446 {
00447 if (x.Size() == size && vectX != 0)
00448 *vectX = x;
00449 }
00450
00451
00452 int
00453 BandGenLinSOE::setBandGenSolver(BandGenLinSolver &newSolver)
00454 {
00455 newSolver.setLinearSOE(*this);
00456
00457 if (size != 0) {
00458 int solverOK = newSolver.setSize();
00459 if (solverOK < 0) {
00460 opserr << "WARNING:BandGenLinSOE::setSolver :";
00461 opserr << "the new solver could not setSeize() - staying with old\n";
00462 return solverOK;
00463 }
00464 }
00465
00466 return this->setSolver(newSolver);
00467 }
00468
00469
00470 int
00471 BandGenLinSOE::sendSelf(int commitTag, Channel &theChannel)
00472 {
00473 return 0;
00474 }
00475
00476
00477 int
00478 BandGenLinSOE::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
00479 {
00480 return 0;
00481 }
00482
00483