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