00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include <fstream.h>
00016
00017 #include "SymArpackSOE.h"
00018 #include "SymArpackSolver.h"
00019 #include <Matrix.h>
00020 #include <Graph.h>
00021 #include <Vertex.h>
00022 #include <VertexIter.h>
00023 #include <math.h>
00024 #include <Channel.h>
00025 #include <FEM_ObjectBroker.h>
00026 #include <AnalysisModel.h>
00027 #include <Vector.h>
00028
00029 SymArpackSOE::SymArpackSOE(SymArpackSolver &the_Solver, AnalysisModel &aModel,
00030 double theShift)
00031 :EigenSOE(the_Solver, EigenSOE_TAGS_SymArpackSOE),
00032 size(0), nnz(0), colA(0), rowStartA(0),
00033 factored(false), shift(theShift), theModel(&aModel),
00034 nblks(0), xblk(0), invp(0), diag(0), penv(0), rowblks(0),
00035 begblk(0), first(0)
00036 {
00037 the_Solver.setEigenSOE(*this);
00038 }
00039
00040
00041 SymArpackSOE::~SymArpackSOE()
00042 {
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058 }
00059
00060
00061 int
00062 SymArpackSOE::getNumEqn(void) const
00063 {
00064 return size;
00065 }
00066
00067
00068 extern "C" int symFactorization(int *fxadj, int *adjncy, int neq, int LSPARSE,
00069 int **xblkMY, int **invpMY, int **rowblksMY,
00070 OFFDBLK ***begblkMY, OFFDBLK **firstMY,
00071 double ***penvMY, double **diagMY);
00072
00073 int
00074 SymArpackSOE::setSize(Graph &theGraph)
00075 {
00076 int result = 0;
00077 int oldSize = size;
00078 size = theGraph.getNumVertex();
00079
00080
00081 Vertex *theVertex;
00082 int newNNZ = 0;
00083 VertexIter &theVertices = theGraph.getVertices();
00084 while ((theVertex = theVertices()) != 0) {
00085 const ID &theAdjacency = theVertex->getAdjacency();
00086 newNNZ += theAdjacency.Size();
00087 }
00088 nnz = newNNZ;
00089
00090 if (colA != 0)
00091 delete [] colA;
00092
00093 colA = new int[newNNZ];
00094 if (colA == 0) {
00095 cerr << "WARNING SymArpackSOE::SymArpackSOE :";
00096 cerr << " ran out of memory for colA with nnz = ";
00097 cerr << newNNZ << " \n";
00098 size = 0; nnz = 0;
00099 result = -1;
00100 }
00101
00102 factored = false;
00103
00104 if (rowStartA != 0)
00105 delete [] rowStartA;
00106 rowStartA = new int[size+1];
00107 if (rowStartA == 0) {
00108 cerr << "SymArpackSOE::ran out of memory for rowStartA." << endl;
00109 result = -1;
00110 }
00111
00112
00113 if (size != 0) {
00114 rowStartA[0] = 0;
00115 int startLoc = 0;
00116 int lastLoc = 0;
00117
00118 for (int a=0; a<size; a++) {
00119 theVertex = theGraph.getVertexPtr(a);
00120 if (theVertex == 0) {
00121 cerr << "WARNING:SymArpackSOE::setSize :";
00122 cerr << " vertex " << a << " not in graph! - size set to 0\n";
00123 size = 0;
00124 return -1;
00125 }
00126
00127 const ID &theAdjacency = theVertex->getAdjacency();
00128 int idSize = theAdjacency.Size();
00129
00130
00131 for (int i=0; i<idSize; i++) {
00132 int row = theAdjacency(i);
00133 bool foundPlace = false;
00134
00135 for (int j=startLoc; j<lastLoc; j++)
00136 if (colA[j] > row) {
00137
00138
00139 for (int k=lastLoc; k>j; k--)
00140 colA[k] = colA[k-1];
00141 colA[j] = row;
00142 foundPlace = true;
00143 j = lastLoc;
00144 }
00145
00146 if (foundPlace == false)
00147 colA[lastLoc] = row;
00148
00149 lastLoc++;
00150 }
00151 rowStartA[a+1] = lastLoc;;
00152 startLoc = lastLoc;
00153 }
00154 }
00155
00156
00157
00158
00159 int LSPARSE = 1;
00160
00161
00162
00163
00164 nblks = symFactorization(rowStartA, colA, size, LSPARSE,
00165 &xblk, &invp, &rowblks, &begblk, &first, &penv, &diag);
00166
00167
00168 EigenSolver *theSolvr = this->getSolver();
00169 int solverOK = theSolvr->setSize();
00170 if (solverOK < 0) {
00171 cerr << "WARNING:BandArpackSOE::setSize :";
00172 cerr << " solver failed setSize()\n";
00173 return solverOK;
00174 }
00175
00176
00177 return result;
00178 }
00179
00180
00181 int
00182 SymArpackSOE::addA(const Matrix &m, const ID &id, double fact)
00183 {
00184
00185 if (fact == 0.0)
00186 return 0;
00187
00188 int idSize = id.Size();
00189
00190
00191 if (idSize != m.noRows() && idSize != m.noCols()) {
00192 cerr << "SymArpackSOE::addA() ";
00193 cerr << " - Matrix and ID not of similar sizes\n";
00194 return -1;
00195 }
00196
00197 int *newID = new int [idSize];
00198 int *isort = new int[idSize];
00199 if (newID == 0 || isort ==0) {
00200 cerr << "WARNING SymArpackSOE::addA() :";
00201 cerr << " ran out of memory for vectors (newID, isort)";
00202 return -1;
00203 }
00204
00205 int i;
00206 for (i=0; i<idSize; i++) {
00207 newID[i] = id(i);
00208 if (newID[i] >= 0)
00209 newID[i] = invp[newID[i]];
00210 }
00211
00212 long int i_eq, j_eq, iadd;
00213 int j, nee, lnee;
00214 int k, ipos, jpos;
00215 int it, jt;
00216 int iblk;
00217 int jblk;
00218 OFFDBLK *ptr;
00219 OFFDBLK *saveblk;
00220 double *fpt, *iloc, *loc;
00221
00222 nee = idSize;
00223 lnee = nee;
00224
00225
00226 for(i = 0, k = 0; i < lnee ; i++ )
00227 {
00228 if( newID[i] >= 0 ) {
00229 isort[k] = i;
00230 k++;
00231 }
00232 }
00233
00234 lnee = k;
00235
00236 i = k - 1;
00237 do
00238 {
00239 k = 0 ;
00240 for (j = 0 ; j < i ; j++)
00241 {
00242 if ( newID[isort[j]] > newID[isort[j+1]]) {
00243 int temp = isort[j];
00244 isort[j] = isort[j+1];
00245 isort[j+1] = temp;
00246
00247
00248
00249
00250
00251 k = j ;
00252 }
00253 }
00254 i = k ;
00255 } while ( k > 0) ;
00256
00257 i = 0 ;
00258 ipos = isort[i] ;
00259 k = rowblks[newID[ipos]] ;
00260 saveblk = begblk[k] ;
00261
00262 for (i=0; i<lnee; i++)
00263 {
00264 ipos = isort[i] ;
00265 i_eq = newID[ipos] ;
00266 iblk = rowblks[i_eq] ;
00267 iloc = penv[i_eq +1] - i_eq ;
00268 if (k < iblk)
00269 while (saveblk->row != i_eq)
00270 saveblk = saveblk->bnext ;
00271
00272 ptr = saveblk ;
00273 for (j=0; j< i ; j++)
00274 {
00275 jpos = isort[j] ;
00276 j_eq = newID[jpos] ;
00277
00278 if (ipos > jpos) {
00279 jt = ipos;
00280 it = jpos;
00281 } else {
00282 it = ipos;
00283 jt = jpos;
00284 }
00285
00286 if (j_eq >= xblk[iblk])
00287 { loc = iloc + j_eq ;
00288 *loc += m(it, jt) * fact;
00289 } else
00290 { while((j_eq >= (ptr->next)->beg) && ((ptr->next)->row == i_eq))
00291 ptr = ptr->next ;
00292 fpt = ptr->nz ;
00293 fpt[j_eq - ptr->beg] += m(it,jt) * fact;
00294 }
00295 }
00296 diag[i_eq] += m(ipos, ipos) * fact;
00297 }
00298
00299 delete [] newID;
00300 delete [] isort;
00301
00302 return 0;
00303 }
00304
00305
00306 int
00307 SymArpackSOE::addM(const Matrix &m, const ID &id, double fact)
00308 {
00309 return this->addA(m, id, -shift);
00310 }
00311
00312
00313 double
00314 SymArpackSOE::getShift(void)
00315 {
00316 return shift;
00317 }
00318
00319 void
00320 SymArpackSOE::zeroA(void)
00321 {
00322 factored = false;
00323 }
00324
00325 void
00326 SymArpackSOE::zeroM(void)
00327 {
00328
00329 }
00330
00331 int
00332 SymArpackSOE::sendSelf(int cTag, Channel &theChannel)
00333 {
00334 return 0;
00335 }
00336
00337 int
00338 SymArpackSOE::recvSelf(int cTag, Channel &theChannel,
00339 FEM_ObjectBroker &theBroker)
00340 {
00341 return 0;
00342 }
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354