ConjugateGradientSolver.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:01 $ 00023 // $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/cg/ConjugateGradientSolver.cpp,v $ 00024 00025 00026 // File: ~/system_of_eqn/linearSOE/cg/ConjugateGradientSolver.C 00027 // 00028 // Written: fmk 00029 // Created: 06/00 00030 // Revision: A 00031 00032 #include "ConjugateGradientSolver.h" 00033 #include <Vector.h> 00034 #include <OPS_Globals.h> 00035 #include <LinearSOE.h> 00036 00037 ConjugateGradientSolver::ConjugateGradientSolver(int classtag, 00038 LinearSOE *theSOE, 00039 double tol) 00040 :LinearSOESolver(classtag), 00041 r(0),p(0),Ap(0),x(0), 00042 theLinearSOE(theSOE), 00043 tolerance(tol) 00044 { 00045 00046 } 00047 00048 ConjugateGradientSolver::~ConjugateGradientSolver() 00049 { 00050 if (r != 0) 00051 delete r; 00052 if (p != 0) 00053 delete p; 00054 if (Ap != 0) 00055 delete Ap; 00056 if (x != 0) 00057 delete x; 00058 } 00059 00060 00061 int 00062 ConjugateGradientSolver::setSize(void) 00063 { 00064 int n = theLinearSOE->getNumEqn(); 00065 if (n <= 0) { 00066 opserr << "ConjugateGradientSolver::setSize() - n < 0 \n"; 00067 return -1; 00068 } 00069 00070 // if old Vector exist and are of incorrect size destroy them 00071 if (r != 0) { 00072 if (r->Size() != n) { 00073 delete r; 00074 delete p; 00075 delete Ap; 00076 delete x; 00077 r = 0; 00078 p = 0; 00079 Ap = 0; 00080 x = 0; 00081 } 00082 } 00083 00084 // create new vector of correct size 00085 if (r == 0) { 00086 r = new Vector(n); 00087 p = new Vector(n); 00088 Ap = new Vector(n); 00089 x = new Vector(n); 00090 if (r == 0 || p == 0 || Ap == 0 || x == 0) { 00091 opserr << "ConjugateGradientSolver::setSize() - out of memory\n"; 00092 if (r != 0) 00093 delete r; 00094 if (p != 0) 00095 delete p; 00096 if (Ap != 0) 00097 delete Ap; 00098 if (x != 0) 00099 delete x; 00100 r = 0; 00101 p = 0; 00102 Ap = 0; 00103 x = 0; 00104 return -2; 00105 00106 } 00107 } 00108 return 0; 00109 } 00110 00111 00112 00113 int 00114 ConjugateGradientSolver::solve(void) 00115 { 00116 // check for successful setSize 00117 if (r == 0) 00118 return -1; 00119 00120 // initialize 00121 x->Zero(); 00122 *r = theLinearSOE->getB(); 00123 *p = *r; 00124 double rdotr = *r ^ *r; 00125 00126 // lopp till convergence 00127 while (r->Norm() > tolerance) { 00128 this->formAp(*p, *Ap); 00129 00130 double alpha = rdotr/(*p ^ *Ap); 00131 00132 // *x += *p * alpha; 00133 x->addVector(1.0, *p, alpha); 00134 00135 // *r -= *Ap * alpha; 00136 r->addVector(1.0, *Ap, -alpha); 00137 00138 double oldrdotr = rdotr; 00139 00140 rdotr = *r ^ *r; 00141 00142 double beta = rdotr / oldrdotr; 00143 00144 // *p = *r + *p * beta; 00145 p->addVector(beta, *r, 1.0); 00146 } 00147 return 0; 00148 } 00149 00150 00151 00152 00153 00154 00155 00156 00157 00158 00159
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