Logo Search packages:      
Sourcecode: objcryst-fox version File versions  Download package

template<typename FloatType>
void scitbx::lbfgsb::raw::mainlb ( int const &  n,
int const &  m,
ref1< FloatType > const &  x,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< int > const &  nbd,
FloatType &  f,
ref1< FloatType > const &  g,
FloatType const &  factr,
FloatType const &  pgtol,
ref2< FloatType > const &  ws,
ref2< FloatType > const &  wy,
ref2< FloatType > const &  sy,
ref2< FloatType > const &  ss,
ref2< FloatType > const &  ,
ref2< FloatType > const &  wt,
ref2< FloatType > const &  wn,
ref2< FloatType > const &  snd,
ref1< FloatType > const &  z,
ref1< FloatType > const &  r,
ref1< FloatType > const &  d,
ref1< FloatType > const &  t,
ref1< FloatType > const &  wa,
ref1< FloatType > const &  sg,
ref1< FloatType > const &  ,
ref1< FloatType > const &  yg,
ref1< FloatType > const &  ,
ref1< int > const &  index,
ref1< int > const &  iwhere,
ref1< int > const &  indx2,
std::string &  task,
int const &  iprint,
std::string &  csave,
ref1< bool > const &  lsave,
ref1< int > const &  isave,
ref1< FloatType > const &  dsave 
) [inline]

Solution of bound constrained optimization problems.

This subroutine solves bound constrained optimization problems by using the compact formula of the limited memory BFGS updates.

n is an integer variable. On entry n is the number of variables. On exit n is unchanged.

m is an integer variable. On entry m is the maximum number of variable metric corrections allowed in the limited memory matrix. On exit m is unchanged.

x is a double precision array of dimension n. On entry x is an approximation to the solution. On exit x is the current approximation.

l is a double precision array of dimension n. On entry l is the lower bound of x. On exit l is unchanged.

u is a double precision array of dimension n. On entry u is the upper bound of x. On exit u is unchanged.

nbd is an integer array of dimension n. On entry nbd represents the type of bounds imposed on the variables, and must be specified as follows: nbd(i)=0 if x(i) is unbounded, 1 if x(i) has only a lower bound, 2 if x(i) has both lower and upper bounds, 3 if x(i) has only an upper bound. On exit nbd is unchanged.

f is a double precision variable. On first entry f is unspecified. On final exit f is the value of the function at x.

g is a double precision array of dimension n. On first entry g is unspecified. On final exit g is the value of the gradient at x.

factr is a double precision variable. On entry factr >= 0 is specified by the user. The iteration will stop when

(f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch

where epsmch is the machine precision, which is automatically generated by the code. On exit factr is unchanged.

pgtol is a double precision variable. On entry pgtol >= 0 is specified by the user. The iteration will stop when

max{|proj g_i | i = 1, ..., n} <= pgtol

where pg_i is the ith component of the projected gradient. On exit pgtol is unchanged.

ws, wy, sy, and wt are double precision working arrays used to store the following information defining the limited memory BFGS matrix: ws, of dimension n x m, stores S, the matrix of s-vectors; wy, of dimension n x m, stores Y, the matrix of y-vectors; sy, of dimension m x m, stores S'Y; ss, of dimension m x m, stores S'S; yy, of dimension m x m, stores Y'Y; wt, of dimension m x m, stores the Cholesky factorization of (theta*S'S+LD^(-1)L'); see eq. (2.26) in [3].

wn is a double precision working array of dimension 2m x 2m used to store the LEL^T factorization of the indefinite matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] [L_a -R_z theta*S'AA'S ]

where E = [-I 0] [ 0 I]

snd is a double precision working array of dimension 2m x 2m used to store the lower triangular part of N = [Y' ZZ'Y L_a'+R_z'] [L_a +R_z S'AA'S ]

z(n),r(n),d(n),t(n),wa(8*m) are double precision working arrays. z is used at different times to store the Cauchy point and the Newton point.

sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays.

index is an integer working array of dimension n. In subroutine freev, index is used to store the free and fixed variables at the Generalized Cauchy Point (GCP).

iwhere is an integer working array of dimension n used to record the status of the vector x for GCP computation. iwhere(i)=0 or -3 if x(i) is free and has bounds, 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) -1 if x(i) is always free, i.e., no bounds on it.

indx2 is an integer working array of dimension n. Within subroutine cauchy, indx2 corresponds to the array iorder. In subroutine freev, a list of variables entering and leaving the free set is stored in indx2, and it is passed on to subroutine formk with this information.

task is a working string of characters of length 60 indicating the current job when entering and leaving this subroutine.

iprint is an INTEGER variable that must be set by the user. It controls the frequency and type of output generated: iprint<0 no output is generated; iprint=0 print only one line at the last iteration; 0<iprint<99 print also f and |proj g| every iprint iterations; iprint=99 print details of every iteration except n-vectors; iprint=100 print also the changes of active set and final x; iprint>100 print details of every iteration including x and g; When iprint > 0, the file iterate.dat will be created to summarize the iteration.

csave is a working string of characters of length 60.

lsave is a logical working array of dimension 4.

isave is an integer working array of dimension 23.

dsave is a double precision working array of dimension 29.

Subprograms called

L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk,

errclb, prn1lb, prn2lb, prn3lb, active, projgr,

freev, cmprlb, matupd, formt.

Minpack2 Library ... timer, dpmeps.

Linpack Library ... dcopy, ddot.

References:

[1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited memory algorithm for bound constrained optimization'', SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.

[2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN Subroutines for Large Scale Bound Constrained Optimization'' Tech. Report, NAM-11, EECS Department, Northwestern University, 1994.

[3] R. Byrd, J. Nocedal and R. Schnabel "Representations of Quasi-Newton Matrices and their use in Limited Memory Methods'', Mathematical Programming 63 (1994), no. 4, pp. 129-156.

(Postscript files of these papers are available via anonymous ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)

* *

NEOS, November 1994. (Latest revision June 1996.) Optimization Technology Center. Argonne National Laboratory and Northwestern University. Written by Ciyou Zhu in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.

Definition at line 3609 of file raw.h.

References active(), scitbx::af::ref< ElementType, AccessorType >::begin(), cauchy(), cmprlb(), dcopy(), dscal(), errclb(), formk(), formt(), freev(), scitbx::math::floating_point_epsilon< FloatType >::get(), scitbx::lbfgsb::raw::ref1< ElementType >::get1(), scitbx::lbfgsb::raw::ref2< ElementType >::get2(), lnsrlb(), matupd(), max3(), prn1lb(), prn2lb(), prn3lb(), projgr(), SCITBX_ASSERT, scitbx::af::const_ref< ElementType, AccessorType >::size(), subsm(), and timer().

Referenced by setulb().

  {
    SCITBX_ASSERT(x.size() == n);
    SCITBX_ASSERT(l.size() == n);
    SCITBX_ASSERT(u.size() == n);
    SCITBX_ASSERT(nbd.size() == n);
    SCITBX_ASSERT(g.size() == n);
#if (SCITBX_LBFGSB_RAW_ASSERTION_FLAG != 0)
    SCITBX_ASSERT(ws.size() == n*m);
    SCITBX_ASSERT(wy.size() == n*m);
    SCITBX_ASSERT(sy.size() == m*m);
    SCITBX_ASSERT(ss.size() == m*m);
    SCITBX_ASSERT(yy.size() == m*m);
    SCITBX_ASSERT(wt.size() == m*m);
    SCITBX_ASSERT(wn.size() == 2*m*2*m);
    SCITBX_ASSERT(snd.size() == 2*m*2*m);
    SCITBX_ASSERT(z.size() == n);
    SCITBX_ASSERT(r.size() == n);
    SCITBX_ASSERT(d.size() == n);
    SCITBX_ASSERT(t.size() == n);
    SCITBX_ASSERT(wa.size() == 8*m);
    SCITBX_ASSERT(sg.size() == m);
    SCITBX_ASSERT(sgo.size() == m);
    SCITBX_ASSERT(yg.size() == m);
    SCITBX_ASSERT(ygo.size() == m);
    SCITBX_ASSERT(index.size() == n);
    SCITBX_ASSERT(iwhere.size() == n);
    SCITBX_ASSERT(indx2.size() == n);
    SCITBX_ASSERT(lsave.size() == 4);
    SCITBX_ASSERT(isave.size() == 23);
    SCITBX_ASSERT(dsave.size() == 29);
#endif
    FloatType zero = 0;
    FloatType one = 1;
    // begin variables in [lid]save arrays
    bool prjctd = false; // uninitialized
    bool cnstnd = false; // uninitialized
    bool boxed = false; // uninitialized
    bool updatd = false; // uninitialized
    int nintol = 0; // uninitialized
    int itfile = 0; // uninitialized
    int iback = 0; // uninitialized
    int nskip = 0; // uninitialized
    int head = 0; // uninitialized
    int col = 0; // uninitialized
    int itail = 0; // uninitialized
    int iter = 0; // uninitialized
    int iupdat = 0; // uninitialized
    int nint = 0; // uninitialized
    int nfgv = 0; // uninitialized
    int info = 0; // uninitialized
    int ifun = 0; // uninitialized
    int iword = 0; // uninitialized
    int nfree = 0; // uninitialized
    int nact = 0; // uninitialized
    int ileave = 0; // uninitialized
    int nenter = 0; // uninitialized
    FloatType theta = 0; // uninitialized
    FloatType fold = 0; // uninitialized
    FloatType tol = 0; // uninitialized
    FloatType dnorm = 0; // uninitialized
    FloatType epsmch = 0; // uninitialized
    FloatType cpu1 = 0; // uninitialized
    FloatType cachyt = 0; // uninitialized
    FloatType sbtime = 0; // uninitialized
    FloatType lnscht = 0; // uninitialized
    FloatType time1 = 0; // uninitialized
    FloatType gd = 0; // uninitialized
    FloatType stpmx = 0; // uninitialized
    FloatType sbgnrm = 0; // uninitialized
    FloatType stp = 0; // uninitialized
    FloatType gdold = 0; // uninitialized
    FloatType dtd = 0; // uninitialized
    // end variables in [lid]save arrays
    bool wrk = false; // uninitialized
    int k = 0; // uninitialized
    FloatType cpu2 = 0; // uninitialized
    FloatType dr = 0; // uninitialized
    FloatType rr = 0; // uninitialized
    FloatType xstep = 0; // uninitialized
    std::string word; // uninitialized
    if (task.substr(0,5) == "START") {
      timer(time1);
      // Generate the current machine precision.
      epsmch = scitbx::math::floating_point_epsilon<FloatType>::get();
      // Initialize counters and scalars when task='START'.
      // for the limited memory BFGS matrices:
      col    = 0;
      head   = 1;
      theta  = one;
      iupdat = 0;
      updatd = false;
      // for operation counts:
      iter   = 0;
      nfgv   = 0;
      nint   = 0;
      nintol = 0;
      nskip  = 0;
      nfree  = n;
      // for stopping tolerance:
      tol = factr*epsmch;
      // for measuring running time:
      cachyt = 0;
      sbtime = 0;
      lnscht = 0;
      // 'word' records the status of subspace solutions.
      word = "---";
      // 'info' records the termination information.
      info = 0;
      if (iprint >= 1) {
        // open a summary file 'iterate.dat'
        //PR open (8, file = 'iterate.dat', status = 'unknown');
        itfile = 8;
      }
      // Check the input arguments for errors.
      errclb(n,m,factr,l,u,nbd,task,info,k);
      if (task.substr(0,5) == "ERROR") {
        prn3lb(n,x,f,task,iprint,info,itfile,
               iter,nfgv,nintol,nskip,nact,sbgnrm,
               zero,nint,word,iback,stp,xstep,k,
               cachyt,sbtime,lnscht);
        return;
      }
      prn1lb(n,m,l,u,x,iprint,itfile,epsmch);
      // Initialize iwhere & project x onto the feasible set.
      active(n,l,u,nbd,x,iwhere,iprint,prjctd,cnstnd,boxed);
    }
    else {
      // restore local variables.
      prjctd = lsave(1);
      cnstnd = lsave(2);
      boxed  = lsave(3);
      updatd = lsave(4);
      nintol = isave(1);
      itfile = isave(3);
      iback  = isave(4);
      nskip  = isave(5);
      head   = isave(6);
      col    = isave(7);
      itail  = isave(8);
      iter   = isave(9);
      iupdat = isave(10);
      nint   = isave(12);
      nfgv   = isave(13);
      info   = isave(14);
      ifun   = isave(15);
      iword  = isave(16);
      nfree  = isave(17);
      nact   = isave(18);
      ileave = isave(19);
      nenter = isave(20);
      theta  = dsave(1);
      fold   = dsave(2);
      tol    = dsave(3);
      dnorm  = dsave(4);
      epsmch = dsave(5);
      cpu1   = dsave(6);
      cachyt = dsave(7);
      sbtime = dsave(8);
      lnscht = dsave(9);
      time1  = dsave(10);
      gd     = dsave(11);
      stpmx  = dsave(12);
      sbgnrm = dsave(13);
      stp    = dsave(14);
      gdold  = dsave(15);
      dtd    = dsave(16);
      // After returning from the driver go to the point where execution
      // is to resume.
      if (task.substr(0,5) == "FG_LN") goto lbl_666;
      if (task.substr(0,5) == "NEW_X") goto lbl_777;
      if (task.substr(0,5) == "FG_ST") goto lbl_111;
      if (task.substr(0,4) == "STOP") {
        if (task.substr(6,3) == "CPU") {
          // restore the previous iterate.
          dcopy(n,t,1,x,1);
          dcopy(n,r,1,g,1);
          f = fold;
        }
        goto lbl_999;
      }
    }
    // Compute f0 and g0.
    task = "FG_START";
    // return to the driver to calculate f and g; reenter at 111.
    goto lbl_1000;
    lbl_111:
    nfgv = 1;
    // Compute the infinity norm of the (-) projected gradient.
    projgr(n,l,u,nbd,x,g,sbgnrm);
    if (iprint >= 1) {
      printf("\nAt iterate%5d    f= %12.5E    |proj g|= %12.5E\n",
             iter,f,sbgnrm);
      printf(
        " %4d %4d     -     -   -     -     -        -    %10.3E %10.3E\n",
        iter,nfgv,sbgnrm,f);
    }
    if (sbgnrm <= pgtol) {
      // terminate the algorithm.
      task = "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL";
      goto lbl_999;
    }
    //----------------- the beginning of the loop --------------------------
    lbl_222:
    if (iprint >= 99) {
      printf("\n\nITERATION %5d\n", iter + 1);
    }
    iword = -1;
    if (!cnstnd && col > 0) {
      // skip the search for GCP.
      dcopy(n,x,1,z,1);
      wrk = updatd;
      nint = 0;
      goto lbl_333;
    }
    //cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    //
    //    Compute the Generalized Cauchy Point (GCP).
    //
    //cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    timer(cpu1);
    cauchy(
           n,x,l,u,nbd,g,indx2,iwhere,t,d,z,
           m,wy.get2(1,1,n,col),ws.get2(1,1,n,col),sy,wt,theta,col,head,
           wa.get1(1, 2*m),
           wa.get1(2*m+1, 2*m),
           wa.get1(4*m+1, 2*m),
           wa.get1(6*m+1, 2*m),
           nint,sg,yg,iprint,sbgnrm,info,epsmch);
    static const char* fmt_1005 =
      "\n"
      " Singular triangular system detected;\n"
      "   refresh the lbfgs memory and restart the iteration.\n";
    if (info != 0) {
      // singular triangular system detected; refresh the lbfgs memory.
      if(iprint >= 1) printf(fmt_1005);
      info   = 0;
      col    = 0;
      head   = 1;
      theta  = one;
      iupdat = 0;
      updatd = false;
      timer(cpu2);
      cachyt = cachyt + cpu2 - cpu1;
      goto lbl_222;
    }
    timer(cpu2);
    cachyt = cachyt + cpu2 - cpu1;
    nintol = nintol + nint;
    // Count the entering and leaving variables for iter > 0;
    // find the index set of free and active variables at the GCP.
    freev(n,nfree,index,nenter,ileave,indx2,
          iwhere,wrk,updatd,cnstnd,iprint,iter);
    nact = n - nfree;
    lbl_333:
    // If there are no free variables or B=theta*I, then
    // skip the subspace minimization.
    if (nfree == 0 || col == 0) goto lbl_555;
    //cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    //
    //    Subspace minimization.
    //
    //cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    timer(cpu1);
    // Form  the LEL^T factorization of the indefinite
    //   matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
    //                 [L_a -R_z           theta*S'AA'S ]
    //   where     E = [-I  0]
    //                 [ 0  I]
    if (wrk) formk(n,nfree,index,nenter,ileave,indx2,iupdat,
                   updatd,wn,snd,m,ws,wy,sy,theta,col,head,info);
    if (info != 0) {
      // nonpositive definiteness in Cholesky factorization;
      // refresh the lbfgs memory and restart the iteration.
      if(iprint >= 1) {
        printf("\n"
          " Nonpositive definiteness in Cholesky factorization in formk;\n"
          "   refresh the lbfgs memory and restart the iteration.\n");
      }
      info   = 0;
      col    = 0;
      head   = 1;
      theta  = one;
      iupdat = 0;
      updatd = false;
      timer(cpu2);
      sbtime = sbtime + cpu2 - cpu1;
      goto lbl_222;
    }
    // compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x)
    // from 'cauchy').
    cmprlb(n,m,x,g,ws,wy,sy,wt,z,r,wa.get1(1,4*m),index,
           theta,col,head,nfree,cnstnd,info);
    if (info == 0) {
      // call the direct method.
      subsm(n,m,nfree,index.get1(1,nfree),l,u,nbd,z,r,ws,wy,theta,
            col,head,iword,wa.get1(1,2*m),wn,iprint,info);
    }
    if (info != 0) {
      // singular triangular system detected;
      // refresh the lbfgs memory and restart the iteration.
      if(iprint >= 1)  printf(fmt_1005);
      info   = 0;
      col    = 0;
      head   = 1;
      theta  = one;
      iupdat = 0;
      updatd = false;
      timer(cpu2);
      sbtime = sbtime + cpu2 - cpu1;
      goto lbl_222;
    }
    timer(cpu2);
    sbtime = sbtime + cpu2 - cpu1;
    lbl_555:
    //cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    //
    //    Line search and optimality tests.
    //
    //cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    // Generate the search direction d:=z-x.
    for(int i=1;i<=n;i++) {
      d(i) = z(i) - x(i);
    }
    timer(cpu1);
    lbl_666:
    lnsrlb(n,l,u,nbd,x,f,fold,gd,gdold,g,d,r,t,z,stp,dnorm,
           dtd,xstep,stpmx,iter,ifun,iback,nfgv,info,task,
           boxed,cnstnd,csave,isave.get1(22,2),dsave.get1(17,13));
    if (info != 0 || iback >= 20) {
      // restore the previous iterate.
      dcopy(n,t,1,x,1);
      dcopy(n,r,1,g,1);
      f = fold;
      if (col == 0) {
        // abnormal termination.
        if (info == 0) {
          info = -9;
          // restore the actual number of f and g evaluations etc.
          nfgv = nfgv - 1;
          ifun = ifun - 1;
          iback = iback - 1;
        }
        task = "ABNORMAL_TERMINATION_IN_LNSRCH";
        iter = iter + 1;
        goto lbl_999;
      }
      else {
        // refresh the lbfgs memory and restart the iteration.
        if(iprint >= 1) {
          printf("\n"
                 " Bad direction in the line search;\n"
                 "   refresh the lbfgs memory and restart the iteration.\n");
        }
        if (info == 0) nfgv = nfgv - 1;
        info   = 0;
        col    = 0;
        head   = 1;
        theta  = one;
        iupdat = 0;
        updatd = false;
        task   = "RESTART_FROM_LNSRCH";
        timer(cpu2);
        lnscht = lnscht + cpu2 - cpu1;
        goto lbl_222;
      }
    }
    else if (task.substr(0,5) == "FG_LN") {
      // return to the driver for calculating f and g; reenter at 666.
      goto lbl_1000;
    }
    else {
      // calculate and print out the quantities related to the new X.
      timer(cpu2);
      lnscht = lnscht + cpu2 - cpu1;
      iter = iter + 1;
      // Compute the infinity norm of the projected (-)gradient.
      projgr(n,l,u,nbd,x,g,sbgnrm);
      // Print iteration information.
      prn2lb(n,x,f,g,iprint,itfile,iter,nfgv,nact,
             sbgnrm,nint,word,iword,iback,stp,xstep);
      goto lbl_1000;
    }
    lbl_777:
    { // scope for variables
      // Test for termination.
      if (sbgnrm <= pgtol) {
        // terminate the algorithm.
        task = "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL";
        goto lbl_999;
      }
      FloatType ddum = max3(fn::absolute(fold), fn::absolute(f), one);
      if ((fold - f) <= tol*ddum) {
        // terminate the algorithm.
        task = "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH";
        if (iback >= 10) info = -5;
        // i.e., to issue a warning if iback>10 in the line search.
        goto lbl_999;
      }
      // Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's.
      for(int i=1;i<=n;i++) {
        r(i) = g(i) - r(i);
      }
      rr = lbfgs::detail::ddot(n,r.begin(),r.begin());
      if (stp == one) {
        dr = gd - gdold;
        ddum = -gdold;
      }
      else {
        dr = (gd - gdold)*stp;
        dscal(n,stp,d,1);
        ddum = -gdold*stp;
      }
      if (dr <= epsmch*ddum) {
        // skip the L-BFGS update.
        nskip = nskip + 1;
        updatd = false;
        if (iprint >= 1) {
          printf("  ys=%10.3E  -gs=%10.3E BFGS update SKIPPED\n",
                 dr, ddum);
        }
        goto lbl_222; // goto lbl_888;
      }
    }
    //cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    //
    //    Update the L-BFGS matrix.
    //
    //cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    updatd = true;
    iupdat = iupdat + 1;
    // Update matrices WS and WY and form the middle matrix in B.
    matupd(n,m,ws,wy,sy,ss,d,r,itail,
           iupdat,col,head,theta,rr,dr,stp,dtd);
    // Form the upper half of the pds T = theta*SS + L*D^(-1)*L';
    // Store T in the upper triangular of the array wt;
    // Cholesky factorize T to J*J' with
    // J' stored in the upper triangular of wt.
    formt(m,wt,sy,ss,col,theta,info);
    if (info != 0) {
      // nonpositive definiteness in Cholesky factorization;
      // refresh the lbfgs memory and restart the iteration.
      if(iprint >= 1) {
        printf("\n"
          " Nonpositive definiteness in Cholesky factorization in formt;\n"
          "   refresh the lbfgs memory and restart the iteration.\n");
      }
      info = 0;
      col = 0;
      head = 1;
      theta = one;
      iupdat = 0;
      updatd = false;
      goto lbl_222;
    }
    // Now the inverse of the middle matrix in B is
    //   [  D^(1/2)      O ] [ -D^(1/2)  D^(-1/2)*L' ]
    //   [ -L*D^(-1/2)   J ] [  0        J'          ]
    // lbl_888:
    //-------------------- the end of the loop -----------------------------
    goto lbl_222;
    lbl_999:
    { // scope for variables
      FloatType time2;
      timer(time2);
      FloatType time = time2 - time1;
      prn3lb(n,x,f,task,iprint,info,itfile,
             iter,nfgv,nintol,nskip,nact,sbgnrm,
             time,nint,word,iback,stp,xstep,k,
             cachyt,sbtime,lnscht);
    }
    lbl_1000:
    // Save local variables.
    lsave(1)  = prjctd;
    lsave(2)  = cnstnd;
    lsave(3)  = boxed;
    lsave(4)  = updatd;
    isave(1)  = nintol;
    isave(3)  = itfile;
    isave(4)  = iback;
    isave(5)  = nskip;
    isave(6)  = head;
    isave(7)  = col;
    isave(8)  = itail;
    isave(9)  = iter;
    isave(10) = iupdat;
    isave(12) = nint;
    isave(13) = nfgv;
    isave(14) = info;
    isave(15) = ifun;
    isave(16) = iword;
    isave(17) = nfree;
    isave(18) = nact;
    isave(19) = ileave;
    isave(20) = nenter;
    dsave(1)  = theta;
    dsave(2)  = fold;
    dsave(3)  = tol;
    dsave(4)  = dnorm;
    dsave(5)  = epsmch;
    dsave(6)  = cpu1;
    dsave(7)  = cachyt;
    dsave(8)  = sbtime;
    dsave(9)  = lnscht;
    dsave(10) = time1;
    dsave(11) = gd;
    dsave(12) = stpmx;
    dsave(13) = sbgnrm;
    dsave(14) = stp;
    dsave(15) = gdold;
    dsave(16) = dtd;
  }


Generated by  Doxygen 1.6.0   Back to index