My title

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.1 $
00022 // $Date: 2001/07/18 16:24:10 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/algorithm/equiSolnAlgo/KrylovNewton.cpp,v $
00024 
00025 // Written: MHS
00026 // Created: June 2001
00027 //
00028 // Description: This file contains the class definition for 
00029 // KrylovNewton.  KrylovNewton is a class which uses a Krylov
00030 // subspace accelerator on the modified Newton method.
00031 // The accelerator is described by Carlson and Miller in
00032 // "Design and Application of a 1D GWMFE Code"
00033 // from SIAM Journal of Scientific Computing (Vol. 19, No. 3,
00034 // pp. 728-765, May 1998)
00035 
00036 #include <KrylovNewton.h>
00037 #include <AnalysisModel.h>
00038 #include <StaticAnalysis.h>
00039 #include <IncrementalIntegrator.h>
00040 #include <LinearSOE.h>
00041 #include <Channel.h>
00042 #include <FEM_ObjectBroker.h>
00043 #include <ConvergenceTest.h>
00044 #include <Matrix.h>
00045 #include <Vector.h>
00046 #include <ID.h>
00047 
00048 #include <fstream.h>
00049 
00050 // Constructor
00051 KrylovNewton::KrylovNewton(int theTangentToUse)
00052 :EquiSolnAlgo(EquiALGORITHM_TAGS_KrylovNewton),
00053  theTest(0), tangent(theTangentToUse),
00054  r(0), v(0), Av(0),
00055  numTests(0), numEqns(0)
00056 {
00057 
00058 }
00059 
00060 KrylovNewton::KrylovNewton(ConvergenceTest &theT, int theTangentToUse)
00061 :EquiSolnAlgo(EquiALGORITHM_TAGS_KrylovNewton),
00062  theTest(&theT), tangent(theTangentToUse),
00063  r(0), v(0), Av(0),
00064  numTests(0), numEqns(0)
00065 {
00066 
00067 }
00068 
00069 // Destructor
00070 KrylovNewton::~KrylovNewton()
00071 {
00072   if (r != 0) {
00073     for (int i = 0; i < numTests; i++)
00074       delete r[i];
00075     delete [] r;
00076   }
00077 
00078   if (v != 0) {
00079     for (int i = 0; i < numTests; i++)
00080       delete v[i];
00081     delete [] v;
00082   }
00083 
00084   if (Av != 0) {
00085     for (int i = 0; i < numTests; i++)
00086       delete Av[i];
00087     delete [] Av;
00088   }
00089 }
00090 
00091 void 
00092 KrylovNewton::setTest(ConvergenceTest &newTest)
00093 {
00094   theTest = &newTest;
00095 }
00096 
00097 int 
00098 KrylovNewton::solveCurrentStep(void)
00099 {
00100   // set up some pointers and check they are valid
00101   // NOTE this could be taken away if we set Ptrs as protecetd in superclass
00102   AnalysisModel   *theAnaModel = this->getAnalysisModelPtr();
00103   IncrementalIntegrator *theIntegrator = this->getIncrementalIntegratorPtr();
00104   LinearSOE  *theSOE = this->getLinearSOEptr();
00105   
00106   if ((theAnaModel == 0) || (theIntegrator == 0) || (theSOE == 0)
00107       || (theTest == 0)){
00108     cerr << "WARNING KrylovNewton::solveCurrentStep() - setLinks() has";
00109     cerr << " not been called - or no ConvergenceTest has been set\n";
00110     return -5;
00111   } 
00112   
00113   numTests = 2 + theTest->getMaxNumTests();
00114   numEqns  = theSOE->getNumEqn();
00115   
00116   if (r == 0) {
00117     r = new Vector*[numTests];
00118     for (int i = 0; i < numTests; i++)
00119       r[i] = new Vector(numEqns);
00120   }
00121 
00122   if (v == 0) {
00123     v = new Vector*[numTests];
00124     for (int i = 0; i < numTests; i++)
00125       v[i] = new Vector(numEqns);
00126   }
00127 
00128   if (Av == 0) {
00129     Av = new Vector*[numTests];
00130     for (int i = 0; i < numTests; i++)
00131       Av[i] = new Vector(numEqns);
00132   }
00133 
00134   // set itself as the ConvergenceTest objects EquiSolnAlgo
00135   theTest->setEquiSolnAlgo(*this);
00136   if (theTest->start() < 0) {
00137     cerr << "KrylovNewton::solveCurrentStep() -";
00138     cerr << "the ConvergenceTest object failed in start()\n";
00139     return -3;
00140     }
00141   
00142   // Evaluate system residual R(y_0)
00143   if (theIntegrator->formUnbalance() < 0) {
00144     cerr << "WARNING KrylovNewton::solveCurrentStep() -";
00145     cerr << "the Integrator failed in formUnbalance()\n"; 
00146     return -2;
00147   }
00148   
00149   // Evaluate system Jacobian J = R'(y)|y_0
00150   if (theIntegrator->formTangent(tangent) < 0){
00151     cerr << "WARNING KrylovNewton::solveCurrentStep() -";
00152     cerr << "the Integrator failed in formTangent()\n";
00153     return -1;
00154   }    
00155   
00156   // Solve for residual f(y_0) = J^{-1} R(y_0)
00157   if (theSOE->solve() < 0) {
00158     cerr << "WARNING KrylovNewton::solveCurrentStep() -";
00159     cerr << "the LinearSysOfEqn failed in solve()\n"; 
00160     return -3;
00161   }
00162 
00163   // Update y ... y_1 = y_0 + v_1 --> (y_0 is zero) --> y_1 = v_1
00164   // --> (v_1 = f(y_0)) --> y_1 = f(y_0)
00165   *(r[0]) = theSOE->getX();
00166 
00167   // Store first update v_1 = r_0
00168   *(v[1]) = *(r[0]);
00169 
00170   int result = -1;
00171   int k = 1;
00172 
00173   do {
00174 
00175     // Update system with v_k
00176     if (theIntegrator->update(*(v[k])) < 0) {
00177       cerr << "WARNING KrylovNewton::solveCurrentStep() -";
00178       cerr << "the Integrator failed in update()\n"; 
00179       return -4;
00180     } 
00181 
00182     // Evaluate system residual R(y_k)
00183     if (theIntegrator->formUnbalance() < 0) {
00184       cerr << "WARNING KrylovNewton::solveCurrentStep() -";
00185       cerr << "the Integrator failed in formUnbalance()\n"; 
00186       return -2;
00187     }
00188     
00189     // Solve for residual f(y_k) = J^{-1} R(y_k)
00190     if (theSOE->solve() < 0) {
00191       cerr << "WARNING KrylovNewton::solveCurrentStep() -";
00192       cerr << "the LinearSysOfEqn failed in solve()\n"; 
00193       return -3;
00194     }
00195 
00196     // Store residual r_k \equiv f(y_k)
00197     *(r[k]) = theSOE->getX();
00198 
00199     // Compute Av_k = f(y_{k-1}) - f(y_k) = r_{k-1} - r_k
00200     // NOTE: Not using Av[0] so that notation follows paper
00201     *(Av[k]) = *(r[k-1]);
00202     Av[k]->addVector(1.0, *(r[k]), -1.0);
00203 
00204     // Solve least squares A w_{k+1} = r_k
00205     if (this->leastSquares(k) < 0) {
00206       cerr << "WARNING KrylovNewton::solveCurrentStep() -";
00207       cerr << "the Integrator failed in leastSquares()\n";
00208       return -1;
00209     }      
00210 
00211     result = theTest->test();
00212     this->record(k++);
00213 
00214   } while (result == -1);
00215   
00216   if (result == -2) {
00217     cerr << "KrylovNewton::solveCurrentStep() -";
00218     cerr << "the ConvergenceTest object failed in test()\n";
00219     return -3;
00220   }
00221   
00222   // note - if postive result we are returning what the convergence
00223   // test returned which should be the number of iterations
00224   return result;
00225 }
00226 
00227 ConvergenceTest *
00228 KrylovNewton::getTest(void)
00229 {
00230   return theTest;
00231 }
00232 
00233 int
00234 KrylovNewton::sendSelf(int cTag, Channel &theChannel)
00235 {
00236   return -1;
00237 }
00238 
00239 int
00240 KrylovNewton::recvSelf(int cTag, 
00241          Channel &theChannel, 
00242          FEM_ObjectBroker &theBroker)
00243 {
00244   return -1;
00245 }
00246 
00247 void
00248 KrylovNewton::Print(ostream &s, int flag)
00249 {
00250   if (flag == 0)
00251     s << "KrylovNewton\n";
00252 }
00253 
00254 int
00255 KrylovNewton::leastSquares(int k)
00256 {
00257   // NOTE: Not using Av[0] so that notation follows paper
00258   const int maxK = 30;
00259 
00260   if (k > maxK) {
00261     cerr << "WARNING KrylovNewton::leastSquares() -";
00262     cerr << "Number of iterations exceeds matrix limits";
00263     return -2;
00264   }
00265 
00266   static double workA[maxK*maxK];
00267   static double workb[maxK];
00268   static double workc[maxK];
00269 
00270   Matrix A(workA,k,k);
00271   Vector b(workb,k);
00272   Vector c(workc,k);
00273 
00274   int i,j;
00275 
00276   // Form normal equations (Av_i,Av_j) c_j = (Av_i,r_k)
00277   for (i = 1; i <= k; i++) {
00278     for (j = 1; j <= k; j++)
00279       A(i-1,j-1) = *(Av[i]) ^ *(Av[j]);
00280     b(i-1) = *(Av[i]) ^ *(r[k]);
00281   }
00282 
00283   // Solve normal equations
00284   if (A.Solve(b,c) < 0) {
00285     cerr << "WARNING KrylovNewton::leastSquares() -";
00286     cerr << "solution of normal equations failed in Matrix::Solve()\n";
00287     return -1;
00288   }
00289 
00290   // v_{k+1} = w_{k+1} + q_{k+1}
00291   *(v[k+1]) = *(r[k]);
00292 
00293   double cj;
00294   for (j = 1; j <= k; j++) {
00295 
00296     cj = c(j-1);
00297 
00298     // Compute w_{k+1} = c_1 v_1 + ... + c_k v_k
00299     v[k+1]->addVector(1.0, *(v[j]), cj);
00300 
00301     // Compute least squares residual q_{k+1} = r_k - c_1 Av_1 - ... - c_k Av_k
00302     v[k+1]->addVector(1.0, *(Av[j]), -cj);
00303   }
00304 
00305   return 0;
00306 }