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