Method solveCurrentStep

SRC/analysis/algorithm/equiSolnAlgo/KrylovNewton.cpp:94–233

}

int 
KrylovNewton::solveCurrentStep(void)
{
  // set up some pointers and check they are valid
  // NOTE this could be taken away if we set Ptrs as protecetd in superclass
  AnalysisModel *theAnaModel = this->getAnalysisModelPtr();
  IncrementalIntegrator *theIntegrator = this->getIncrementalIntegratorPtr();
  LinearSOE *theSOE = this->getLinearSOEptr();
  
  if ((theAnaModel == 0) || (theIntegrator == 0) || (theSOE == 0)
      || (theTest == 0)){
    opserr << "WARNING KrylovNewton::solveCurrentStep() - setLinks() has";
    opserr << " not been called - or no ConvergenceTest has been set\n";
    return -5;
  }	

  // Get size information from SOE
  numEqns  = theSOE->getNumEqn();
  if (maxDimension > numEqns)
    maxDimension = numEqns;

  if (v == 0) {
    // Need to allocate an extra vector for "next" update
    v = new Vector*[maxDimension+1];
    for (int i = 0; i < maxDimension+1; i++)
      v[i] = new Vector(numEqns);
  }

  if (Av == 0) {
    Av = new Vector*[maxDimension+1];
    for (int i = 0; i < maxDimension+1; i++)
      Av[i] = new Vector(numEqns);
  }

  if (AvData == 0)
    AvData = new double [maxDimension*numEqns];

  if (rData == 0)
    // The LAPACK least squares subroutine overwrites the RHS vector
    // with the solution vector ... these vectors are not the same
    // size, so we need to use the max size
    rData = new double [(numEqns > maxDimension) ? numEqns : maxDimension];

  // Length of work vector should be >= 2*min(numEqns,maxDimension)
  // See dgels subroutine documentation
  lwork = 2 * ((numEqns < maxDimension) ? numEqns : maxDimension);
  
  if (work == 0)
    work = new double [lwork];

  // Evaluate system residual R(y_0)
  if (theIntegrator->formUnbalance() < 0) {
    opserr << "WARNING KrylovNewton::solveCurrentStep() -";
    opserr << "the Integrator failed in formUnbalance()\n";	
    return -2;
  }