| 379 | |
| 380 | |
| 381 | void Broyden::BroydenUpdate( IncrementalIntegrator *theIntegrator, |
| 382 | LinearSOE *theSOE, |
| 383 | Vector &du, |
| 384 | int nBroyden ) |
| 385 | { |
| 386 | |
| 387 | static const double eps = 1.0e-16 ; |
| 388 | |
| 389 | // int systemSize = ( theSOE->getB() ).Size(); |
| 390 | int systemSize = theSOE->getNumEqn( ) ; |
| 391 | |
| 392 | |
| 393 | //compute z |
| 394 | // theSOE->setB( (*r[nBroyden]) - (*r[nBroyden-1]) ) ; |
| 395 | // theSOE->setB( (*residNew) - (*residOld) ) ; |
| 396 | *temp = (*residNew) ; |
| 397 | *temp -= (*residOld) ; |
| 398 | theSOE->setB( *temp ) ; |
| 399 | |
| 400 | if (theSOE->solve() < 0) { |
| 401 | opserr << "WARNING Broyden::solveCurrentStep() -"; |
| 402 | opserr << "the LinearSysOfEqn failed in solve()\n"; |
| 403 | } |
| 404 | |
| 405 | if ( z[nBroyden] == 0 ) |
| 406 | z[nBroyden] = new Vector(systemSize) ; |
| 407 | |
| 408 | *z[nBroyden] = theSOE->getX() ; |
| 409 | *z[nBroyden] *= (-1.0) ; |
| 410 | |
| 411 | int i; |
| 412 | for ( i=1; i<=(nBroyden-1); i++ ) { |
| 413 | |
| 414 | double p = - ( (*s[i]) ^ (*z[i]) ) ; |
| 415 | |
| 416 | if ( fabs(p) < eps ) break ; |
| 417 | |
| 418 | double sdotz = (*s[i]) ^ (*z[nBroyden]) ; |
| 419 | |
| 420 | //*z[nBroyden] += (1.0/p) * sdotz * ( *s[i] + *z[i] ) ; |
| 421 | *temp = (*s[i]) ; |
| 422 | *temp += (*z[i]) ; |
| 423 | *temp *= ( (1.0/p) * sdotz ) ; |
| 424 | *z[nBroyden] += (*temp) ; |
| 425 | |
| 426 | |
| 427 | } //end for i |
| 428 | |
| 429 | |
| 430 | //broyden modifications to du |
| 431 | for ( i=1; i<=nBroyden; i++ ) { |
| 432 | |
| 433 | double p = - ( (*s[i]) ^ (*z[i]) ) ; |
| 434 | |
| 435 | if ( fabs(p) < eps ) break ; |
| 436 | |
| 437 | double sdotdu = (*s[i]) ^ du ; |
| 438 | |