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