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.3 $
00022 // $Date: 2001/07/31 01:06:16 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/analysis/algorithm/equiSolnAlgo/BFGS.cpp,v $
00024                                                                         
00025 // Written: Ed Love
00026 // Created: 06/01
00027 
00028 // What: "@(#)BFGS.cpp, revA"
00029 
00030 #include <BFGS.h>
00031 #include <AnalysisModel.h>
00032 #include <StaticAnalysis.h>
00033 #include <IncrementalIntegrator.h>
00034 #include <LinearSOE.h>
00035 #include <Channel.h>
00036 #include <FEM_ObjectBroker.h>
00037 #include <ConvergenceTest.h>
00038 #include <ID.h>
00039 
00040 
00041 #include <fstream.h>
00042 
00043 
00044 // Constructor
00045 BFGS::BFGS(int theTangentToUse, int n )
00046 :EquiSolnAlgo(EquiALGORITHM_TAGS_BFGS),
00047  theTest(0), tangent(theTangentToUse), numberLoops(n) 
00048 {
00049   s  = new Vector*[numberLoops+3];
00050 
00051   z  = new Vector*[numberLoops+3];
00052 
00053   //  r  = new (Vector*)[numberLoops+3];
00054 
00055   residOld = 0;
00056   residNew = 0;
00057   du = 0;
00058   b  = 0;
00059 
00060   temp = 0;
00061 
00062   rdotz = 0;
00063   sdotr = 0;
00064 
00065   for ( int i =0; i < numberLoops+3; i++ ) {
00066     s[i] = 0;
00067     z[i] = 0;
00068   }
00069  
00070   localTest = 0;
00071 
00072 }
00073 
00074 //Constructor
00075 BFGS::BFGS(ConvergenceTest &theT, int theTangentToUse, int n)
00076 :EquiSolnAlgo(EquiALGORITHM_TAGS_BFGS),
00077  theTest(&theT), tangent(theTangentToUse), numberLoops(n) 
00078 {
00079   s  = new Vector*[numberLoops+3];
00080 
00081   z  = new Vector*[numberLoops+3];
00082 
00083   residOld = 0;
00084   residNew = 0;
00085   du = 0;
00086   b  = 0;
00087   temp = 0;
00088 
00089   rdotz = 0;
00090   sdotr = 0;
00091 
00092   for ( int i =0; i < numberLoops+3; i++ ) {
00093     s[i] = 0;
00094     z[i] = 0;
00095     //r[i] = 0;
00096   }
00097 
00098   localTest = theTest->getCopy( this->numberLoops );
00099 
00100 }
00101 
00102 // Destructor
00103 BFGS::~BFGS()
00104 {
00105 
00106   if (temp != 0) delete temp;
00107   temp = 0;
00108 
00109   if (residOld != 0 ) delete residOld;  
00110   residOld = 0;
00111 
00112   if (residNew != 0) delete residNew;
00113   residNew = 0;
00114 
00115   if (du != 0) delete du;
00116   du = 0;
00117 
00118   if (b != 0 ) delete b;
00119   b = 0;
00120 
00121   if (rdotz != 0 ) delete [] rdotz;
00122   rdotz = 0;
00123   
00124   if (sdotr != 0 ) delete [] sdotr;
00125   sdotr = 0;
00126 
00127   for ( int i =0; i < numberLoops+3; i++ ) {
00128     if ( s[i] != 0 ) delete s[i];
00129     if ( z[i] != 0 ) delete z[i];
00130     //delete r[i];
00131     s[i] = 0;
00132     z[i] = 0;
00133     //r[i] = 0;
00134   } //end for i
00135 
00136   if ( s != 0 ) delete[] s; 
00137   if ( z != 0 ) delete[] z;
00138   s = 0;  
00139   z = 0;
00140 
00141   if ( localTest != 0 )
00142      delete localTest;
00143   localTest = 0;
00144 
00145   
00146 }
00147 
00148 
00149 void 
00150 BFGS::setTest(ConvergenceTest &newTest)
00151 {
00152     theTest = &newTest;
00153 
00154     if ( localTest != 0 )  delete localTest;
00155 
00156     localTest = theTest->getCopy( this->numberLoops );
00157 }
00158 
00159 
00160 
00161 int 
00162 BFGS::solveCurrentStep(void)
00163 {
00164  
00165     // set up some pointers and check they are valid
00166     // NOTE this could be taken away if we set Ptrs as protecetd in superclass
00167 
00168     AnalysisModel   *theAnaModel = this->getAnalysisModelPtr();
00169 
00170     IncrementalIntegrator *theIntegrator = this->getIncrementalIntegratorPtr();
00171 
00172     LinearSOE  *theSOE = this->getLinearSOEptr();
00173 
00174     if ((theAnaModel == 0) || (theIntegrator == 0) || (theSOE == 0)
00175  || (theTest == 0)){
00176  cerr << "WARNING BFGS::solveCurrentStep() - setLinks() has";
00177  cerr << " not been called - or no ConvergenceTest has been set\n";
00178  return -5;
00179     } 
00180 
00181     // set itself as the ConvergenceTest objects EquiSolnAlgo
00182     theTest->setEquiSolnAlgo(*this);
00183     if (theTest->start() < 0) {
00184       cerr << "BFGS::solveCurrentStep() -";
00185       cerr << "the ConvergenceTest object failed in start()\n";
00186       return -3;
00187     }
00188 
00189     localTest->setEquiSolnAlgo(*this);
00190 
00191     if (rdotz == 0)
00192        rdotz = new double[numberLoops+3];
00193 
00194     if (sdotr == 0)
00195  sdotr = new double[numberLoops+3];
00196 
00197 
00198     int result = -1;
00199     int count = 0;
00200     do {
00201 
00202       // cerr << "      BFGS -- Forming New Tangent" << endl;
00203 
00204       //form the initial tangent
00205       if (theIntegrator->formTangent(tangent) < 0){
00206          cerr << "WARNING BFGS::solveCurrentStep() -";
00207          cerr << "the Integrator failed in formTangent()\n";
00208          return -1; 
00209       }
00210 
00211       //form the initial residual 
00212       if (theIntegrator->formUnbalance() < 0) {
00213         cerr << "WARNING BFGS::solveCurrentStep() -";
00214         cerr << "the Integrator failed in formUnbalance()\n"; 
00215       }     
00216 
00217       //solve
00218       if (theSOE->solve() < 0) {
00219    cerr << "WARNING BFGS::solveCurrentStep() -";
00220    cerr << "the LinearSysOfEqn failed in solve()\n"; 
00221    return -3;
00222  }     
00223 
00224       //update
00225       if ( theIntegrator->update(theSOE->getX() ) < 0) {
00226  cerr << "WARNING BFGS::solveCurrentStep() -";
00227  cerr << "the Integrator failed in update()\n"; 
00228  return -4;
00229       }         
00230 
00231 
00232       //    int systemSize = ( theSOE->getB() ).Size();
00233       int systemSize = theSOE->getNumEqn( );
00234 
00235       //temporary vector
00236       if (temp == 0 )
00237  temp = new Vector(systemSize);
00238 
00239       //initial displacement increment
00240       if ( s[1] == 0 ) 
00241  s[1] = new Vector(systemSize);
00242 
00243       *s[1] = theSOE->getX( );
00244 
00245       if ( residOld == 0 ) 
00246  residOld = new Vector(systemSize);
00247 
00248       *residOld = theSOE->getB( ) ;
00249       *residOld *= (-1.0 );
00250 
00251       //form the residual again
00252       if (theIntegrator->formUnbalance() < 0) {
00253         cerr << "WARNING BFGS::solveCurrentStep() -";
00254         cerr << "the Integrator failed in formUnbalance()\n"; 
00255       }     
00256 
00257       if ( residNew == 0 ) 
00258  residNew = new Vector(systemSize);
00259  
00260       if ( du == 0 ) 
00261  du = new Vector(systemSize);
00262 
00263       if ( b == 0 )
00264  b = new Vector(systemSize);
00265 
00266       localTest->start();
00267 
00268       int nBFGS = 1;
00269       do {
00270 
00271         //save residual
00272         *residNew =  theSOE->getB( ); 
00273         *residNew *= (-1.0 );
00274 
00275       
00276         //solve
00277         if (theSOE->solve() < 0) {
00278      cerr << "WARNING BFGS::solveCurrentStep() -";
00279      cerr << "the LinearSysOfEqn failed in solve()\n"; 
00280      return -3;
00281         }     
00282 
00283  //save right hand side
00284         *b = theSOE->getB( );
00285 
00286         //save displacement increment
00287         *du = theSOE->getX( );
00288 
00289         //BFGS modifications to du
00290         BFGSUpdate( theIntegrator, theSOE, *du, *b, nBFGS ) ;
00291 
00292         if ( theIntegrator->update( *du ) < 0 ) {
00293     cerr << "WARNING BFGS::solveCurrentStep() -";
00294     cerr << "the Integrator failed in update()\n"; 
00295     return -4;
00296         }         
00297 
00298  /* cerr << "        BFGS Iteration " << nBFGS 
00299             << " Residual Norm = " 
00300             << sqrt( (*residNew) ^ (*residNew) ) << endl;
00301  */
00302         
00303         //increment broyden counter
00304         nBFGS += 1;
00305 
00306         //save displacement increment
00307         if ( s[nBFGS] == 0 ) 
00308    s[nBFGS] = new Vector(systemSize);
00309 
00310         *s[nBFGS] = *du;
00311 
00312         //swap residuals
00313  *residOld = *residNew;
00314 
00315         //form the residual again
00316         if (theIntegrator->formUnbalance() < 0) {
00317           cerr << "WARNING BFGS::solveCurrentStep() -";
00318           cerr << "the Integrator failed in formUnbalance()\n"; 
00319         }     
00320 
00321         result = localTest->test();
00322  
00323         
00324       } while ( result == -1 && nBFGS <= numberLoops );
00325 
00326 
00327       result = theTest->test();
00328       this->record(count++);
00329 
00330     }  while (result == -1);
00331 
00332 
00333     if (result == -2) {
00334       cerr << "BFGS::solveCurrentStep() -";
00335       cerr << "the ConvergenceTest object failed in test()\n";
00336       return -3;
00337     }
00338 
00339     // note - if postive result we are returning what the convergence test returned
00340     // which should be the number of iterations
00341     return result;
00342 }
00343 
00344 
00345 
00346 void  BFGS::BFGSUpdate(IncrementalIntegrator *theIntegrator, 
00347          LinearSOE *theSOE, 
00348          Vector &du, 
00349          Vector &b,
00350          int nBFGS) 
00351 {
00352 
00353   static const double eps = 1.0e-16;
00354 
00355   //  int systemSize = ( theSOE->getB() ).Size();
00356       int systemSize = theSOE->getNumEqn( );
00357 
00358 
00359   //compute z
00360   //  theSOE->setB( (*r[nBFGS]) - (*r[nBFGS-1]) );
00361   //    theSOE->setB( (*residNew) - (*residOld) );
00362   *temp = *residNew;
00363   *temp -= *residOld;
00364   theSOE->setB(*temp);
00365 
00366 
00367   if (theSOE->solve() < 0) {
00368        cerr << "WARNING BFGS::solveCurrentStep() -";
00369        cerr << "the LinearSysOfEqn failed in solve()\n"; 
00370    }     
00371   
00372   if ( z[nBFGS] == 0 ) 
00373     z[nBFGS] = new Vector(systemSize);
00374 
00375   *z[nBFGS] = theSOE->getX(); 
00376   //  *z[nBFGS] *= (-1.0);
00377 
00378     int i;
00379   for ( i=1; i<=(nBFGS-1); i++ ) {
00380 
00381     if ( sdotr[i] < eps ) 
00382       break; 
00383 
00384     double fact1 = 1.0 + ( rdotz[i] / sdotr[i] );
00385 
00386     fact1 /= sdotr[i];
00387 
00388     double pdotb = (*s[i]) ^ ( theSOE->getB() );
00389 
00390 
00391     fact1 *= pdotb;
00392 
00393     //    *z[nBFGS] +=  fact1 * ( *s[i] );
00394     *temp = *s[i];
00395     *temp *= fact1;
00396     *z[nBFGS] += *temp;
00397 
00398 
00399     double bdotz = (*z[i]) ^ ( theSOE->getB() );  
00400 
00401     //    *z[nBFGS] -= (1.0/sdotr[i]) * 
00402     //             ( bdotz * (*s[i])   +  pdotb * (*z[i]) );   
00403     *temp = *s[i];
00404     *temp *= bdotz;
00405     *temp /= sdotr[i];
00406     *z[nBFGS] -= *temp;
00407 
00408     *temp = *z[i];
00409     *temp *= pdotb;
00410     *temp /= sdotr[i];
00411     *z[nBFGS] -= *temp;
00412  
00413   } //end for i
00414 
00415 
00416   //sdotr[nBFGS] = *s[nBFGS] ^ ( *residNew - *residOld );
00417 
00418   //rdotz[nBFGS] = *z[nBFGS] ^ ( *residNew - *residOld );   
00419 
00420   *temp = *residNew;
00421   *temp -= *residOld;
00422 
00423   sdotr[nBFGS] = *s[nBFGS] ^ (*temp);
00424 
00425   rdotz[nBFGS] = *z[nBFGS] ^ (*temp);
00426 
00427 
00428   //BFGS modifications to du
00429   for ( i=1; i<=nBFGS; i++ ) {
00430 
00431     if ( sdotr[i] < eps )
00432       break;
00433 
00434     double fact1 = 1.0 + ( rdotz[i] / sdotr[i] );
00435 
00436     fact1 /= sdotr[i];
00437 
00438     double sdotb = (*s[i]) ^ b;
00439 
00440     fact1 *= sdotb;
00441 
00442     //du +=  fact1 * ( *s[i] );
00443     *temp = *s[i];
00444     *temp *= fact1;
00445     du += *temp;
00446 
00447 
00448     double bdotz = (*z[i]) ^ b;  
00449 
00450     //du -= (1.0/sdotr[i]) * 
00451     //             ( bdotz * (*s[i])   +  sdotb * (*z[i]) );   
00452     *temp = *s[i];
00453     *temp *= bdotz;
00454     *temp /= sdotr[i];
00455     du -= *temp;
00456 
00457     *temp = *z[i];
00458     *temp *= sdotb;
00459     *temp /= sdotr[i];
00460     du -= *temp;
00461 
00462 
00463   } //end for i
00464 
00465 }
00466 
00467 
00468 ConvergenceTest *
00469 BFGS::getTest(void)
00470 {
00471   return theTest;
00472 }
00473 
00474 int
00475 BFGS::sendSelf(int cTag, Channel &theChannel)
00476 {
00477   return -1;
00478 }
00479 
00480 int
00481 BFGS::recvSelf(int cTag, 
00482         Channel &theChannel, 
00483         FEM_ObjectBroker &theBroker)
00484 {
00485     return -1;
00486 }
00487 
00488 
00489 void
00490 BFGS::Print(ostream &s, int flag)
00491 {
00492     if (flag == 0) {
00493       s << "BFGS" << endl;
00494       s << "  Number of Iterations = " << numberLoops << endl;
00495     }
00496 }