| 395 | |
| 396 | |
| 397 | void BFGS::BFGSUpdate(IncrementalIntegrator *theIntegrator, |
| 398 | LinearSOE *theSOE, |
| 399 | Vector &du, |
| 400 | Vector &b, |
| 401 | int nBFGS) |
| 402 | { |
| 403 | |
| 404 | static const double eps = 1.0e-16; |
| 405 | |
| 406 | // int systemSize = ( theSOE->getB() ).Size(); |
| 407 | int systemSize = theSOE->getNumEqn( ); |
| 408 | |
| 409 | |
| 410 | //compute z |
| 411 | // theSOE->setB( (*r[nBFGS]) - (*r[nBFGS-1]) ); |
| 412 | // theSOE->setB( (*residNew) - (*residOld) ); |
| 413 | *temp = *residNew; |
| 414 | *temp -= *residOld; |
| 415 | theSOE->setB(*temp); |
| 416 | |
| 417 | |
| 418 | if (theSOE->solve() < 0) { |
| 419 | opserr << "WARNING BFGS::solveCurrentStep() -"; |
| 420 | opserr << "the LinearSysOfEqn failed in solve()\n"; |
| 421 | } |
| 422 | |
| 423 | if ( z[nBFGS] == 0 ) |
| 424 | z[nBFGS] = new Vector(systemSize); |
| 425 | |
| 426 | *z[nBFGS] = theSOE->getX(); |
| 427 | // *z[nBFGS] *= (-1.0); |
| 428 | |
| 429 | int i; |
| 430 | for ( i=1; i<=(nBFGS-1); i++ ) { |
| 431 | |
| 432 | if ( sdotr[i] < eps ) |
| 433 | break; |
| 434 | |
| 435 | double fact1 = 1.0 + ( rdotz[i] / sdotr[i] ); |
| 436 | |
| 437 | fact1 /= sdotr[i]; |
| 438 | |
| 439 | double pdotb = (*s[i]) ^ ( theSOE->getB() ); |
| 440 | |
| 441 | |
| 442 | fact1 *= pdotb; |
| 443 | |
| 444 | // *z[nBFGS] += fact1 * ( *s[i] ); |
| 445 | *temp = *s[i]; |
| 446 | *temp *= fact1; |
| 447 | *z[nBFGS] += *temp; |
| 448 | |
| 449 | |
| 450 | double bdotz = (*z[i]) ^ ( theSOE->getB() ); |
| 451 | |
| 452 | // *z[nBFGS] -= (1.0/sdotr[i]) * |
| 453 | // ( bdotz * (*s[i]) + pdotb * (*z[i]) ); |
| 454 | *temp = *s[i]; |