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