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