badPetscSOE.cppGo to the documentation of this file.00001 /* ****************************************************************** ** 00002 ** OpenSees - Open System for Earthquake Engineering Simulation ** 00003 ** Pacific Earthquake Engineering Research Center ** 00004 ** ** 00005 ** ** 00006 ** (C) Copyright 1999, The Regents of the University of California ** 00007 ** All Rights Reserved. ** 00008 ** ** 00009 ** Commercial use of this program without express permission of the ** 00010 ** University of California, Berkeley, is strictly prohibited. See ** 00011 ** file 'COPYRIGHT' in main directory for information on usage and ** 00012 ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** 00013 ** ** 00014 ** Developed by: ** 00015 ** Frank McKenna (fmckenna@ce.berkeley.edu) ** 00016 ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** 00017 ** Filip C. Filippou (filippou@ce.berkeley.edu) ** 00018 ** ** 00019 ** ****************************************************************** */ 00020 00021 // $Revision: 1.2 $ 00022 // $Date: 2003/02/14 23:02:02 $ 00023 // $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/petsc/badPetscSOE.cpp,v $ 00024 00025 00026 // File: ~/system_of_eqn/linearSOE/petsc/PetscSOE.C 00027 // 00028 // Written: fmk & om 00029 // Created: 7/98 00030 // Revision: A 00031 // 00032 // Description: This file contains the implementation for BandGenLinSOE 00033 00034 00035 #include "PetscSOE.h" 00036 #include "PetscSolver.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 #include <Channel.h> 00044 #include <FEM_ObjectBroker.h> 00045 00046 PetscSOE::PetscSOE(PetscSolver &theSOESolver) 00047 :LinearSOE(theSOESolver, LinSOE_TAGS_PetscSOE), 00048 isFactored(0),size(0), numProcessors(0), B(0), X(0), 00049 vectX(0), vectB(0), A(0), x(0), b(0) 00050 { 00051 theSOESolver.setLinearSOE(*this); 00052 MPI_Comm_size(PETSC_COMM_WORLD, &numProcessors); 00053 } 00054 00055 00056 int 00057 PetscSOE::getNumEqn(void) const 00058 { 00059 return size; 00060 } 00061 00062 PetscSOE::~PetscSOE() 00063 { 00064 /* 00065 if (vectX != 0) delete vectX; 00066 if (vectB != 0) delete vectB; 00067 if (B != 0) delete [] B; 00068 if (X != 0) delete [] X; 00069 00070 // invoke the petsc destructors 00071 if (A != 0) MatDestroy(A); 00072 if (b != 0) VecDestroy(b); 00073 if (x != 0) VecDestroy(x); 00074 */ 00075 } 00076 00077 00078 int 00079 PetscSOE::setSize(Graph &theGraph) 00080 { 00081 int result = 0; 00082 int oldSize = size; 00083 size = theGraph.getNumVertex(); 00084 isFactored = 0; 00085 00086 // delete the old 00087 if (B != 0) delete [] B; 00088 if (X != 0) delete [] X; 00089 00090 00091 // invoke the petsc destructors 00092 if (A != 0) MatDestroy(A); 00093 if (b != 0) VecDestroy(b); 00094 if (x != 0) VecDestroy(x); 00095 00096 // create the new 00097 B = new double[size]; 00098 X = new double[size]; 00099 00100 if (B == 0 || X == 0) { 00101 opserr << "WARNING PetscSOE::PetscSOE :"; 00102 opserr << " ran out of memory for vectors (size) ("; 00103 opserr << size << ") \n"; 00104 size = 0; 00105 result = -1; 00106 } 00107 00108 // zero the vectors 00109 for (int j=0; j<size; j++) { 00110 B[j] = 0; 00111 X[j] = 0; 00112 } 00113 00114 if (vectX != 0) 00115 delete vectX; 00116 00117 if (vectB != 0) 00118 delete vectB; 00119 00120 vectX = new Vector(X,size); 00121 vectB = new Vector(B,size); 00122 00123 // invoke the petsc constructors 00124 int ierr = OptionsGetInt(PETSC_NULL,"-n",&size,&flg); CHKERRA(ierr); 00125 ierr = MatCreate(PETSC_COMM_WORLD,size,size,&A); CHKERRA(ierr); 00126 ierr = VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,size,&x); CHKERRA(ierr); 00127 ierr = VecDuplicate(x,&b); CHKERRA(ierr); 00128 // invoke setSize() on the Solver 00129 int solverOK = theSolver->setSize(); 00130 if (solverOK < 0) { 00131 opserr << "WARNING:PetscSOE::setSize :"; 00132 opserr << " solver failed setSize()\n"; 00133 return solverOK; 00134 } 00135 00136 return result; 00137 } 00138 00139 00140 int 00141 PetscSOE::setSizeDirectly(int newSize) 00142 { 00143 int result = 0; 00144 int oldSize = size; 00145 size = newSize; 00146 isFactored = 0; 00147 00148 // delete the old 00149 if (B != 0) delete [] B; 00150 if (X != 0) delete [] X; 00151 00152 // invoke the petsc destructors 00153 if (A != 0) MatDestroy(A); 00154 if (b != 0) VecDestroy(b); 00155 if (x != 0) VecDestroy(x); 00156 00157 // create the new 00158 B = new double[size]; 00159 X = new double[size]; 00160 00161 if (B == 0 || X == 0) { 00162 opserr << "WARNING PetscSOE::PetscSOE :"; 00163 opserr << " ran out of memory for vectors (size) ("; 00164 opserr << size << ") \n"; 00165 size = 0; 00166 result = -1; 00167 } 00168 00169 // zero the vectors 00170 for (int j=0; j<size; j++) { 00171 B[j] = 0; 00172 X[j] = 0; 00173 } 00174 00175 if (vectX != 0) 00176 delete vectX; 00177 00178 if (vectB != 0) 00179 delete vectB; 00180 00181 vectX = new Vector(X,size); 00182 vectB = new Vector(B,size); 00183 00184 // invoke the petsc constructors 00185 int ierr = OptionsGetInt(PETSC_NULL,"-n",&size,&flg); CHKERRA(ierr); 00186 ierr = MatCreate(PETSC_COMM_WORLD,size,size,&A); CHKERRA(ierr); 00187 ierr = VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,size,&x); CHKERRA(ierr); 00188 ierr = VecDuplicate(x,&b); CHKERRA(ierr); 00189 // invoke setSize() on the Solver 00190 int solverOK = theSolver->setSize(); 00191 if (solverOK < 0) { 00192 opserr << "WARNING:PetscSOE::setSize :"; 00193 opserr << " solver failed setSize()\n"; 00194 return solverOK; 00195 } 00196 00197 return result; 00198 } 00199 00200 00201 int 00202 PetscSOE::setSizeParallel(int newSize, int N, 00203 int d_nz, int *d_nnz, 00204 int o_nz, int *o_nnz) 00205 { 00206 int result = 0; 00207 int oldSize = size; 00208 size = newSize; 00209 isFactored = 0; 00210 00211 // delete the old 00212 if (B != 0) delete [] B; 00213 if (X != 0) delete [] X; 00214 00215 // invoke the petsc destructors 00216 if (A != 0) MatDestroy(A); 00217 if (b != 0) VecDestroy(b); 00218 if (x != 0) VecDestroy(x); 00219 00220 // create the new 00221 B = new double[size]; 00222 X = new double[size]; 00223 00224 if (B == 0 || X == 0) { 00225 opserr << "WARNING PetscSOE::PetscSOE :"; 00226 opserr << " ran out of memory for vectors (size) ("; 00227 opserr << size << ") \n"; 00228 size = 0; 00229 result = -1; 00230 } 00231 00232 // zero the vectors 00233 for (int j=0; j<size; j++) { 00234 B[j] = 0; 00235 X[j] = 0; 00236 } 00237 00238 if (vectX != 0) 00239 delete vectX; 00240 00241 if (vectB != 0) 00242 delete vectB; 00243 00244 vectX = new Vector(X,size); 00245 vectB = new Vector(B,size); 00246 00247 // invoke the petsc constructors 00248 int ierr = OptionsGetInt(PETSC_NULL,"-n",&size,&flg); CHKERRA(ierr); 00249 ierr = MatCreate(PETSC_COMM_WORLD,size,size,&A); CHKERRA(ierr); 00250 ierr = VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,size,&x); CHKERRA(ierr); 00251 ierr = VecDuplicate(x,&b); CHKERRA(ierr); 00252 // invoke setSize() on the Solver 00253 int solverOK = theSolver->setSize(); 00254 if (solverOK < 0) { 00255 opserr << "WARNING:PetscSOE::setSize :"; 00256 opserr << " solver failed setSize()\n"; 00257 return solverOK; 00258 } 00259 00260 return result; 00261 } 00262 00263 00264 00265 int 00266 PetscSOE::addA(const Matrix &m, const ID &id, double fact) 00267 { 00268 isFactored = 0; 00269 00270 // check for a quick return 00271 if (fact == 0.0) return 0; 00272 00273 00274 // check that m and id are of similar size 00275 int idSize = id.Size(); 00276 if (idSize != m.noRows() && idSize != m.noCols()) { 00277 opserr << "PetscSOE::addA() - Matrix and ID not of similar sizes\n"; 00278 return -1; 00279 } 00280 00281 int n = id.Size(); 00282 int row; 00283 int col; 00284 double value; 00285 for (int i=0; i<n; i++) { 00286 row = id(i); 00287 if (row >= 0) { 00288 for (int j=0; j<n; j++) { 00289 col = id(j); 00290 if (col >= 0) { 00291 value = m(i,j)*fact; 00292 int ierr = MatSetValues(A,1,&row,1,&col,&value,ADD_VALUES); CHKERRA(ierr); 00293 } 00294 } 00295 } 00296 } 00297 00298 return 0; 00299 } 00300 00301 00302 int 00303 PetscSOE::addB(const Vector &v, const ID &id, double fact) 00304 { 00305 // check for a quick return 00306 if (fact == 0.0) return 0; 00307 00308 00309 // check that m and id are of similar size 00310 int idSize = id.Size(); 00311 if (idSize != v.Size() ) { 00312 opserr << "PetscSOE::addB() - Vector and ID not of similar sizes\n"; 00313 return -1; 00314 } 00315 00316 int n = id.Size(); 00317 int row; 00318 double value; 00319 for (int i=0; i<n; i++) { 00320 row = id(i); 00321 if (row >= 0) { 00322 value = v(i) * fact; 00323 int ierr = VecSetValues(b,1,&row,&value,ADD_VALUES); CHKERRA(ierr); 00324 } 00325 } 00326 00327 return 0; 00328 } 00329 00330 00331 int 00332 PetscSOE::setB(const Vector &v, double fact) 00333 { 00334 // check for a quick return 00335 if (fact == 0.0) return 0; 00336 00337 00338 if (size != v.Size() ) { 00339 opserr << "PetscSOE::addB() - Vector not of appropriate size\n"; 00340 return -1; 00341 } 00342 00343 int n = size; 00344 int row; 00345 double value; 00346 for (int i=0; i<n; i++) { 00347 value = v(i) * fact; 00348 int ierr = VecSetValues(b,1,&i,&value,INSERT_VALUES); CHKERRA(ierr); 00349 } 00350 int ierr = VecAssemblyBegin(b); CHKERRA(ierr); 00351 ierr = VecAssemblyEnd(b); CHKERRA(ierr); 00352 00353 return 0; 00354 } 00355 00356 00357 void 00358 PetscSOE::zeroA(void) 00359 { 00360 isFactored = 0; 00361 int ierr = MatZeroEntries(A); CHKERRA(ierr); 00362 } 00363 00364 void 00365 PetscSOE::zeroB(void) 00366 { 00367 double zero = 0.0; 00368 VecSet(&zero, b); 00369 } 00370 00371 00372 const Vector & 00373 PetscSOE::getX(void) 00374 { 00375 if (vectX == 0) { 00376 opserr << "FATAL PetscSOE::getX - vectX == 0!"; 00377 exit(-1); 00378 } 00379 double *theData =0; 00380 int ierr = VecGetArray(x, &theData); CHKERRA(ierr); 00381 for (int i=0; i<size; i++) 00382 X[i] = theData[i]; 00383 ierr = VecRestoreArray(x, &theData); CHKERRA(ierr); 00384 00385 00386 return *vectX; 00387 } 00388 00389 00390 const Vector & 00391 PetscSOE::getB(void) 00392 { 00393 if (vectB == 0) { 00394 opserr << "FATAL PetscSOE::getB - vectB == 0!"; 00395 exit(-1); 00396 } 00397 double *theData =0; 00398 int ierr = VecGetArray(b, &theData); CHKERRA(ierr); 00399 for (int i=0; i<size; i++) 00400 B[i] = theData[i]; 00401 ierr = VecRestoreArray(b, &theData); CHKERRA(ierr); 00402 00403 return *vectB; 00404 } 00405 00406 00407 double 00408 PetscSOE::normRHS(void) 00409 { 00410 this->getB(); 00411 double norm =0.0; 00412 double *Bptr = B; 00413 for (int i=0; i<size; i++) { 00414 double Yi = *Bptr++; 00415 norm += Yi*Yi; 00416 } 00417 return sqrt(norm); 00418 } 00419 00420 00421 void 00422 PetscSOE::setX(int loc, double value) 00423 { 00424 int ierr = VecSetValues(x,1,&loc,&value,INSERT_VALUES); CHKERRA(ierr); 00425 } 00426 00427 int 00428 PetscSOE::setSolver(PetscSolver &newSolver) 00429 { 00430 newSolver.setLinearSOE(*this); 00431 00432 if (size != 0) { 00433 int solverOK = newSolver.setSize(); 00434 if (solverOK < 0) { 00435 opserr << "WARNING:PetscSOE::setSolver :"; 00436 opserr << "the new solver could not setSeize() - staying with old\n"; 00437 return solverOK; 00438 } 00439 } 00440 00441 return this->LinearSOE::setSolver(newSolver); 00442 } 00443 00444 00445 int 00446 PetscSOE::sendSelf(Channel &theChannel, FEM_ObjectBroker &theBroker) 00447 { 00448 if (size != 0) 00449 opserr << "WARNING PetscSOE::sendSelf - does not send itself YET\n"; 00450 return 0; 00451 } 00452 00453 00454 int 00455 PetscSOE::recvSelf(Channel &theChannel, FEM_ObjectBroker &theBroker) 00456 { 00457 return 0; 00458 } 00459 00460 00461 00462 00463
Generated on Mon Oct 23 15:05:28 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents