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