EMAN2
lbfgsb.cpp
Go to the documentation of this file.
00001 #include <cstdio>
00002 #include <cstring>
00003 #include <cmath>
00004 #include <cstdlib>
00005 
00006 #define TRUE_ (1)
00007 #define FALSE_ (0)
00008 #define min_(a,b) ((a) <= (b) ? (a) : (b))
00009 #define max_(a,b) ((a) >= (b) ? (a) : (b))
00010 #define abs_(x) ((x) >= 0 ? (x) : -(x))
00011 
00012 typedef struct
00013 {       long int cierr;
00014         long int ciunit;
00015         long int ciend;
00016         char *cifmt;
00017         long int cirec;
00018 } cilist;
00019 
00020 typedef struct
00021 {       long int oerr;
00022         long int ounit;
00023         char *ofnm;
00024         long int ofnmlen;
00025         char *osta;
00026         char *oacc;
00027         char *ofm;
00028         long int orl;
00029         char *oblnk;
00030 } olist;
00031 
00032 /* Table of constant values */
00033 
00034 static double c_b9 = 0.;
00035 static long int c__1 = 1;
00036 //static long int c__9 = 9;
00037 static long int c__11 = 11;
00038 //static long int c__3 = 3;
00039 static double c_b275 = .001;
00040 static double c_b276 = .9;
00041 static double c_b277 = .1;
00042 //static long int c__5 = 5;
00043 
00044 const int SIXTY=60;
00045 
00046 // function prototype for L-BFGS-B routine
00047 // NOTES: All arguments must be passed as pointers. Must add an underscore
00048 //  to the function name. 'extern "C"' keeps c++ from mangling the function
00049 //  name. Be aware that fortran indexes 1d arrays [1,n] not [0,n-1], and
00050 //  expects 2d arrays to be in col-major order not row-major order.
00051 /*int setulb_(long int *n, long int *m, double *x, double *l, double *u, long int *nbd, double *f, double *g, 
00052                    double *factr, double *pgtol, double *wa, long int *iwa, char *task, long int *iprint, char *csave, long int *lsave,
00053                    long int *isave, double *dsave, long int aaa1, long int aaa2);          */
00054            
00055 //int MAIN__() { abort(); return 1; }
00056 
00057 /*int setulb_(long int *, long int *, double *, double *, double *, long int *, double *, double *, 
00058                    double *, double *, double *, long int *, char *, long int *, char *, long int *,
00059                    long int *, double *, long int , long int);*/
00060 
00061 
00062 
00063 /*int main() {
00064   
00065   //c     This simple driver demonstrates how to call the L-BFGS-B code to
00066   //c       solve a sample problem (the extended Rosenbrock function 
00067   //c       subject to bounds on the variables). The dimension n of this
00068   //c       problem is variable.
00069 
00070 
00071   int nmax=1024, mmax=17;
00072   //c        nmax is the dimension of the largest problem to be solved.
00073   //c        mmax is the maximum number of limited memory corrections.
00074  
00075   //c     Declare the variables needed by the code.
00076   //c       A description of all these variables is given at the end of 
00077   //c       the driver.
00078   char task[SIXTY], csave[SIXTY];
00079   long int lsave[4];
00080   long int n, m, iprint, nbd[nmax], iwa[3*nmax], isave[44];
00081   double f, factr, pgtol, x[nmax], l[nmax], u[nmax], g[nmax], dsave[29], 
00082     wa[2*mmax*nmax+4*nmax+12*mmax*mmax+12*mmax];
00083 
00084   //c     Declare a few additional variables for this sample problem.
00085   double t1, t2;
00086   int i;
00087  
00088   //c     We wish to have output at every iteration.
00089   iprint = -1;
00090 
00091   //c     We specify the tolerances in the stopping criteria.
00092   factr=1.0e1;
00093   pgtol=1.0e-9;
00094 
00095   //c     We specify the dimension n of the sample problem and the number
00096   //c        m of limited memory corrections stored.  (n and m should not
00097   //c        exceed the limits nmax and mmax respectively.)
00098   n=25;
00099   m=5;
00100  
00101   //c     We now provide nbd which defines the bounds on the variables:
00102   //c                    l   specifies the lower bounds,
00103   //c                    u   specifies the upper bounds. 
00104  
00105   //c     First set bounds on the odd-numbered variables. (*even-numbered in c*)
00106   for (i=0; i<n; i+=2){
00107     nbd[i]=2;
00108     l[i]=1.0;
00109     u[i]=100.0;
00110   }
00111 
00112   //c     Next set bounds on the even-numbered variables. (*odd-numbered in c*)
00113   for (i=1; i<n; i+=2){
00114     nbd[i]=2;
00115     l[i]=-100.0;
00116     u[i]=100.0;
00117   }
00118 
00119   //c     We now define the starting point.
00120   for (i=0; i<n; i++)
00121     x[i]=3.0;
00122  
00123   printf("Solving sample problem.\n");
00124   printf(" (f=0.0 at the optimal solution.)\n");
00125 
00126 
00127   //c     We start the iteration by initializing task.
00128   // (**MUST clear remaining chars in task with spaces (else crash)!**)
00129   strcpy(task,"START");
00130   for (i=5;i<SIXTY;i++)
00131     task[i]=' ';
00132 
00133   //c     This is the call to the L-BFGS-B code.
00134   // (* call the L-BFGS-B routine with task='START' once before loop *)
00135   setulb_(&n,&m,x,l,u,nbd,&f,g,&factr,&pgtol,wa,iwa,task,&iprint,
00136          csave,lsave,isave,dsave,(long int)60,(long int)60);
00137 
00138   // (* while routine returns "FG" or "NEW_X" in task, keep calling it *)
00139   while (strncmp(task,"FG",2)==0 || strncmp(task,"NEW_X",5)==0) {
00140 
00141     if (strncmp(task,"FG",2)==0) {
00142       //c   the minimization routine has returned to request the
00143       //c   function f and gradient g values at the current x
00144 
00145       //c        Compute function value f for the sample problem.
00146 
00147       f = 0.25 * (x[0]-1.0)*(x[0]-1.0);
00148       for (i=1; i<n; i++)
00149         f += pow(x[i] - x[i-1]*x[i-1], 2);
00150       f *= 4.0;
00151 
00152       //c        Compute gradient g for the sample problem.
00153       //t1 = (x[1]-x[0])*(x[1]-x[0]);
00154       t1 = x[1]-x[0]*x[0];
00155       g[0] = 2.0*(x[0]-1.0) - 16.0*x[0]*t1;
00156       for (i=1; i<n-1; i++) {
00157         t2=t1;
00158         //t1=(x[i+1]-x[i])*(x[i+1]-x[i]);
00159         t1 = x[i+1]-x[i]*x[i];
00160         g[i]=8.0*t2-16.0*x[i]*t1;
00161       }
00162       g[n-1]=8.0*t1;
00163 
00164     }
00165     else {
00166       // the minimization routine has returned with a new iterate,
00167       // and we have opted to continue the iteration (do nothing here)
00168     }
00169 
00170     //c          go back to the minimization routine.
00171     setulb_(&n,&m,x,l,u,nbd,&f,g,&factr,&pgtol,wa,iwa,task,&iprint,
00172            csave,lsave,isave,dsave,(long int)60,(long int)60);
00173   }
00174 
00175   for (i=0; i<25; i++) { printf("i=%2d     x[i]=%f\n", i, x[i]); }
00176   //c     If task is neither FG nor NEW_X we terminate execution.
00177   // (* the minimization routine has returned with one of
00178   //  {CONV, ABNO, ERROR} *)
00179 
00180   return 0;
00181 }*/
00182 
00183 long int s_cmp(char *str1, const char *const str2, long int l1, long int l2) {
00184         long int l;
00185         int i = 0;
00186         long int flag = 0;
00187         
00188         l = min_(l1,l2);
00189         
00190         while (i<l && flag==0) {
00191                 if (str1[i] != str2[i]) { flag = 1; }
00192                 i++;
00193         }
00194         return flag;
00195 }
00196 
00197 int s_copy(char *str1, const char *const str2, long int l1, long int l2) {
00198         long int l;
00199         int i = 0;
00200         l = min_(l1,l2);
00201         
00202         while (i<l) {
00203                 str1[i]=str2[i];
00204                 i++;
00205         }       
00206         return 0;
00207 }
00208 
00209 /*******************************************************************
00210 c     --------------------------------------------------------------
00211 c             DESCRIPTION OF THE VARIABLES IN L-BFGS-B
00212 c     --------------------------------------------------------------
00213 c
00214 c     n is an long int variable that must be set by the user to the
00215 c       number of variables.  It is not altered by the routine.
00216 c
00217 c     m is an long int variable that must be set by the user to the
00218 c       number of corrections used in the limited memory matrix.
00219 c       It is not altered by the routine.  Values of m < 3  are
00220 c       not recommended, and large values of m can result in excessive
00221 c       computing time. The range  3 <= m <= 20 is recommended. 
00222 c
00223 c     x is a DOUBLE PRECISION array of length n.  On initial entry
00224 c       it must be set by the user to the values of the initial
00225 c       estimate of the solution vector.  Upon successful exit, it
00226 c       contains the values of the variables at the best point
00227 c       found (usually an approximate solution).
00228 c
00229 c     l is a DOUBLE PRECISION array of length n that must be set by
00230 c       the user to the values of the lower bounds on the variables. If
00231 c       the i-th variable has no lower bound, l(i) need not be defined.
00232 c
00233 c     u is a DOUBLE PRECISION array of length n that must be set by
00234 c       the user to the values of the upper bounds on the variables. If
00235 c       the i-th variable has no upper bound, u(i) need not be defined.
00236 c
00237 c     nbd is an long int array of dimension n that must be set by the
00238 c       user to the type of bounds imposed on the variables:
00239 c       nbd(i)=0 if x(i) is unbounded,
00240 c              1 if x(i) has only a lower bound,
00241 c              2 if x(i) has both lower and upper bounds, 
00242 c              3 if x(i) has only an upper bound.
00243 c
00244 c     f is a DOUBLE PRECISION variable.  If the routine setulb returns
00245 c       with task(1:2)= 'FG', then f must be set by the user to
00246 c       contain the value of the function at the point x.
00247 c
00248 c     g is a DOUBLE PRECISION array of length n.  If the routine setulb
00249 c       returns with taskb(1:2)= 'FG', then g must be set by the user to
00250 c       contain the components of the gradient at the point x.
00251 c
00252 c     factr is a DOUBLE PRECISION variable that must be set by the user.
00253 c       It is a tolerance in the termination test for the algorithm.
00254 c       The iteration will stop when
00255 c
00256 c        (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
00257 c
00258 c       where epsmch is the machine precision which is automatically
00259 c       generated by the code. Typical values for factr on a computer
00260 c       with 15 digits of accuracy in double precision are:
00261 c       factr=1.d+12 for low accuracy;
00262 c             1.d+7  for moderate accuracy; 
00263 c             1.d+1  for extremely high accuracy.
00264 c       The user can suppress this termination test by setting factr=0.
00265 c
00266 c     pgtol is a double precision variable.
00267 c       On entry pgtol >= 0 is specified by the user.  The iteration
00268 c         will stop when
00269 c
00270 c                 max{|proj g_i | i = 1, ..., n} <= pgtol
00271 c
00272 c         where pg_i is the ith component of the projected gradient.
00273 c       The user can suppress this termination test by setting pgtol=0.
00274 c
00275 c     wa is a DOUBLE PRECISION  array of length 
00276 c       (2mmax + 4)nmax + 12mmax^2 + 12mmax used as workspace.
00277 c       This array must not be altered by the user.
00278 c
00279 c     iwa is an long int  array of length 3nmax used as
00280 c       workspace. This array must not be altered by the user.
00281 c
00282 c     task is a CHARACTER string of length 60.
00283 c       On first entry, it must be set to 'START'.
00284 c       On a return with task(1:2)='FG', the user must evaluate the
00285 c         function f and gradient g at the returned value of x.
00286 c       On a return with task(1:5)='NEW_X', an iteration of the
00287 c         algorithm has concluded, and f and g contain f(x) and g(x)
00288 c         respectively.  The user can decide whether to continue or stop
00289 c         the iteration. 
00290 c       When
00291 c         task(1:4)='CONV', the termination test in L-BFGS-B has been 
00292 c           satisfied;
00293 c         task(1:4)='ABNO', the routine has terminated abnormally
00294 c           without being able to satisfy the termination conditions,
00295 c           x contains the best approximation found,
00296 c           f and g contain f(x) and g(x) respectively;
00297 c         task(1:5)='ERROR', the routine has detected an error in the
00298 c           input parameters;
00299 c       On exit with task = 'CONV', 'ABNO' or 'ERROR', the variable task
00300 c         contains additional information that the user can print.
00301 c       This array should not be altered unless the user wants to
00302 c          stop the run for some reason.  See driver2 or driver3
00303 c          for a detailed explanation on how to stop the run 
00304 c          by assigning task(1:4)='STOP' in the driver.
00305 c
00306 c     iprint is an long int variable that must be set by the user.
00307 c       It controls the frequency and type of output generated:
00308 c        iprint<0    no output is generated;
00309 c        iprint=0    print only one line at the last iteration;
00310 c        0<iprint<99 print also f and |proj g| every iprint iterations;
00311 c        iprint=99   print details of every iteration except n-vectors;
00312 c        iprint=100  print also the changes of active set and final x;
00313 c        iprint>100  print details of every iteration including x and g;
00314 c       When iprint > 0, the file iterate.dat will be created to
00315 c                        summarize the iteration.
00316 c
00317 c     csave  is a CHARACTER working array of length 60.
00318 c
00319 c     lsave is a long int working array of dimension 4.
00320 c       On exit with task = 'NEW_X', the following information is
00321 c         available:
00322 c       lsave(1) = .true.  the initial x did not satisfy the bounds;
00323 c       lsave(2) = .true.  the problem contains bounds;
00324 c       lsave(3) = .true.  each variable has upper and lower bounds.
00325 c
00326 c     isave is an long int working array of dimension 44.
00327 c       On exit with task = 'NEW_X', it contains information that
00328 c       the user may want to access:
00329 c         isave(30) = the current iteration number;
00330 c         isave(34) = the total number of function and gradient
00331 c                         evaluations;
00332 c         isave(36) = the number of function value or gradient
00333 c                                  evaluations in the current iteration;
00334 c         isave(38) = the number of free variables in the current
00335 c                         iteration;
00336 c         isave(39) = the number of active constraints at the current
00337 c                         iteration;
00338 c
00339 c         see the subroutine setulb.f for a description of other 
00340 c         information contained in isave
00341 c
00342 c     dsave is a DOUBLE PRECISION working array of dimension 29.
00343 c       On exit with task = 'NEW_X', it contains information that
00344 c         the user may want to access:
00345 c         dsave(2) = the value of f at the previous iteration;
00346 c         dsave(5) = the machine precision epsmch generated by the code;
00347 c         dsave(13) = the infinity norm of the projected gradient;
00348 c
00349 c         see the subroutine setulb.f for a description of other 
00350 c         information contained in dsave
00351 c
00352 c     --------------------------------------------------------------
00353 c           END OF THE DESCRIPTION OF THE VARIABLES IN L-BFGS-B
00354 c     --------------------------------------------------------------
00355 c
00356 c     << An example of subroutine 'timer' for AIX Version 3.2 >>
00357 c
00358 c     subroutine timer(ttime)
00359 c     double precision ttime
00360 c     long int itemp, long int mclock
00361 c     itemp = mclock()
00362 c     ttime = dble(itemp)*1.0d-2
00363 c     return
00364 c     end
00365 
00366 ********************************************************************/
00367 
00368 /* Subroutine */ int setulb_(long int *n, long int *m, double *x, double *l, double *u, long int *nbd, double *f, double *g, 
00369                          double *factr, double *pgtol, double *wa, long int *iwa, char *task, long int *iprint, char *csave, long int *lsave,
00370                          long int *isave, double *dsave, long int task_len, long int csave_len)
00371 /*long int *n, *m;
00372 double *x, *l, *u;
00373 long int *nbd;
00374 double *f, *g, *factr, *pgtol, *wa;
00375 long int *iwa;
00376 char *task;
00377 long int *iprint;
00378 char *csave;
00379 long int *lsave;
00380 long int *isave;
00381 double *dsave;
00382 long int task_len;
00383 long int csave_len;*/
00384 {
00385     /* System generated locals */
00386     long int i__1;
00387 
00388     /* Builtin functions */
00389     long int s_cmp(char *, const char *const, long int, long int);
00390 
00391     /* Local variables */
00392     static long int lsnd, lsgo, lygo, l1, l2, l3, ld, lr, lt;
00393     extern /* Subroutine */ int mainlb_(long int *n, long int *m, double *x, double *l, double *u, long int *nbd, double *f, 
00394                          double *g, double *factr, double *pgtol, double *ws, double *wy,
00395                          double *sy, double *ss, double *yy, double *wt, double *wn, double *snd, 
00396                          double *z__, double *r__, double *d__, double *t, double *wa, double *sg, 
00397                          double *sgo, double *yg, double *ygo, long int *index, long int *iwhere, long int *indx2, 
00398                          char *task, long int *iprint, char *csave, long int *lsave, long int *isave, double *dsave, 
00399                          long int task_len, long int csave_len);
00400     static long int lz, lwa, lsg, lyg, lwn, lss, lws, lwt, lsy, lwy, lyy;
00401 
00402 /*     ************ */
00403 
00404 /*     Subroutine setulb */
00405 
00406 /*     This subroutine partitions the working arrays wa and iwa, and */
00407 /*       then uses the limited memory BFGS method to solve the bound */
00408 /*       constrained optimization problem by calling mainlb. */
00409 /*       (The direct method will be used in the subspace minimization.) */
00410 
00411 /*     n is an long int variable. */
00412 /*       On entry n is the dimension of the problem. */
00413 /*       On exit n is unchanged. */
00414 
00415 /*     m is an long int variable. */
00416 /*       On entry m is the maximum number of variable metric corrections */
00417 /*         used to define the limited memory matrix. */
00418 /*       On exit m is unchanged. */
00419 
00420 /*     x is a double precision array of dimension n. */
00421 /*       On entry x is an approximation to the solution. */
00422 /*       On exit x is the current approximation. */
00423 
00424 /*     l is a double precision array of dimension n. */
00425 /*       On entry l is the lower bound on x. */
00426 /*       On exit l is unchanged. */
00427 
00428 /*     u is a double precision array of dimension n. */
00429 /*       On entry u is the upper bound on x. */
00430 /*       On exit u is unchanged. */
00431 
00432 /*     nbd is an long int array of dimension n. */
00433 /*       On entry nbd represents the type of bounds imposed on the */
00434 /*         variables, and must be specified as follows: */
00435 /*         nbd(i)=0 if x(i) is unbounded, */
00436 /*                1 if x(i) has only a lower bound, */
00437 /*                2 if x(i) has both lower and upper bounds, and */
00438 /*                3 if x(i) has only an upper bound. */
00439 /*       On exit nbd is unchanged. */
00440 
00441 /*     f is a double precision variable. */
00442 /*       On first entry f is unspecified. */
00443 /*       On final exit f is the value of the function at x. */
00444 
00445 /*     g is a double precision array of dimension n. */
00446 /*       On first entry g is unspecified. */
00447 /*       On final exit g is the value of the gradient at x. */
00448 
00449 /*     factr is a double precision variable. */
00450 /*       On entry factr >= 0 is specified by the user.  The iteration */
00451 /*         will stop when */
00452 
00453 /*         (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch */
00454 
00455 /*         where epsmch is the machine precision, which is automatically */
00456 /*         generated by the code. Typical values for factr: 1.d+12 for */
00457 /*         low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely */
00458 /*         high accuracy. */
00459 /*       On exit factr is unchanged. */
00460 
00461 /*     pgtol is a double precision variable. */
00462 /*       On entry pgtol >= 0 is specified by the user.  The iteration */
00463 /*         will stop when */
00464 
00465 /*                 max{|proj g_i | i = 1, ..., n} <= pgtol */
00466 
00467 /*         where pg_i is the ith component of the projected gradient. */
00468 /*       On exit pgtol is unchanged. */
00469 
00470 /*     wa is a double precision working array of length */
00471 /*       (2mmax + 4)nmax + 12mmax^2 + 12mmax. */
00472 
00473 /*     iwa is an long int working array of length 3nmax. */
00474 
00475 /*     task is a working string of characters of length 60 indicating */
00476 /*       the current job when entering and quitting this subroutine. */
00477 
00478 /*     iprint is an long int variable that must be set by the user. */
00479 /*       It controls the frequency and type of output generated: */
00480 /*        iprint<0    no output is generated; */
00481 /*        iprint=0    print only one line at the last iteration; */
00482 /*        0<iprint<99 print also f and |proj g| every iprint iterations; */
00483 /*        iprint=99   print details of every iteration except n-vectors; */
00484 /*        iprint=100  print also the changes of active set and final x; */
00485 /*        iprint>100  print details of every iteration including x and g; */
00486 /*       When iprint > 0, the file iterate.dat will be created to */
00487 /*                        summarize the iteration. */
00488 
00489 /*     csave is a working string of characters of length 60. */
00490 
00491 /*     lsave is a long int working array of dimension 4. */
00492 /*       On exit with 'task' = NEW_X, the following information is */
00493 /*                                                             available: */
00494 /*         If lsave(1) = .true.  then  the initial X has been replaced by */
00495 /*                                     its projection in the feasible set; */
00496 /*         If lsave(2) = .true.  then  the problem is constrained; */
00497 /*         If lsave(3) = .true.  then  each variable has upper and lower */
00498 /*                                     bounds; */
00499 
00500 /*     isave is an long int working array of dimension 44. */
00501 /*       On exit with 'task' = NEW_X, the following information is */
00502 /*                                                             available: */
00503 /*         isave(22) = the total number of intervals explored in the */
00504 /*                         search of Cauchy points; */
00505 /*         isave(26) = the total number of skipped BFGS updates before */
00506 /*                         the current iteration; */
00507 /*         isave(30) = the number of current iteration; */
00508 /*         isave(31) = the total number of BFGS updates prior the current */
00509 /*                         iteration; */
00510 /*         isave(33) = the number of intervals explored in the search of */
00511 /*                         Cauchy point in the current iteration; */
00512 /*         isave(34) = the total number of function and gradient */
00513 /*                         evaluations; */
00514 /*         isave(36) = the number of function value or gradient */
00515 /*                                  evaluations in the current iteration; */
00516 /*         if isave(37) = 0  then the subspace argmin is within the box; */
00517 /*         if isave(37) = 1  then the subspace argmin is beyond the box; */
00518 /*         isave(38) = the number of free variables in the current */
00519 /*                         iteration; */
00520 /*         isave(39) = the number of active constraints in the current */
00521 /*                         iteration; */
00522 /*         n + 1 - isave(40) = the number of variables leaving the set of */
00523 /*                           active constraints in the current iteration; */
00524 /*         isave(41) = the number of variables entering the set of active */
00525 /*                         constraints in the current iteration. */
00526 
00527 /*     dsave is a double precision working array of dimension 29. */
00528 /*       On exit with 'task' = NEW_X, the following information is */
00529 /*                                                             available: */
00530 /*         dsave(1) = current 'theta' in the BFGS matrix; */
00531 /*         dsave(2) = f(x) in the previous iteration; */
00532 /*         dsave(3) = factr*epsmch; */
00533 /*         dsave(4) = 2-norm of the line search direction vector; */
00534 /*         dsave(5) = the machine precision epsmch generated by the code; */
00535 /*         dsave(7) = the accumulated time spent on searching for */
00536 /*                                                         Cauchy points; */
00537 /*         dsave(8) = the accumulated time spent on */
00538 /*                                                 subspace minimization; */
00539 /*         dsave(9) = the accumulated time spent on line search; */
00540 /*         dsave(11) = the slope of the line search function at */
00541 /*                                  the current point of line search; */
00542 /*         dsave(12) = the maximum relative step length imposed in */
00543 /*                                                           line search; */
00544 /*         dsave(13) = the infinity norm of the projected gradient; */
00545 /*         dsave(14) = the relative step length in the line search; */
00546 /*         dsave(15) = the slope of the line search function at */
00547 /*                                 the starting point of the line search; */
00548 /*         dsave(16) = the square of the 2-norm of the line search */
00549 /*                                                      direction vector. */
00550 
00551 /*     Subprograms called: */
00552 
00553 /*       L-BFGS-B Library ... mainlb. */
00554 
00555 
00556 /*     References: */
00557 
00558 /*       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
00559 /*       memory algorithm for bound constrained optimization'', */
00560 /*       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
00561 
00562 /*       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a */
00563 /*       limited memory FORTRAN code for solving bound constrained */
00564 /*       optimization problems'', Tech. Report, NAM-11, EECS Department, */
00565 /*       Northwestern University, 1994. */
00566 
00567 /*       (Postscript files of these papers are available via anonymous */
00568 /*        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
00569 
00570 /*                           *  *  * */
00571 
00572 /*     NEOS, November 1994. (Latest revision June 1996.) */
00573 /*     Optimization Technology Center. */
00574 /*     Argonne National Laboratory and Northwestern University. */
00575 /*     Written by */
00576 /*                        Ciyou Zhu */
00577 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
00578 
00579 
00580 /*     ************ */
00581     /* Parameter adjustments */
00582     --iwa;
00583     --g;
00584     --nbd;
00585     --u;
00586     --l;
00587     --x;
00588     --wa;
00589     --lsave;
00590     --isave;
00591     --dsave;
00592 
00593     /* Function Body */
00594     if (s_cmp(task, "START", (long int)60, (long int)5) == 0) {
00595         isave[1] = *m * *n;
00596 /* Computing 2nd power */
00597         i__1 = *m;
00598         isave[2] = i__1 * i__1;
00599 /* Computing 2nd power */
00600         i__1 = *m;
00601         isave[3] = i__1 * i__1 << 2;
00602         isave[4] = 1;
00603         isave[5] = isave[4] + isave[1];
00604         isave[6] = isave[5] + isave[1];
00605         isave[7] = isave[6] + isave[2];
00606         isave[8] = isave[7] + isave[2];
00607         isave[9] = isave[8] + isave[2];
00608         isave[10] = isave[9] + isave[2];
00609         isave[11] = isave[10] + isave[3];
00610         isave[12] = isave[11] + isave[3];
00611         isave[13] = isave[12] + *n;
00612         isave[14] = isave[13] + *n;
00613         isave[15] = isave[14] + *n;
00614         isave[16] = isave[15] + *n;
00615         isave[17] = isave[16] + (*m << 3);
00616         isave[18] = isave[17] + *m;
00617         isave[19] = isave[18] + *m;
00618         isave[20] = isave[19] + *m;
00619     }
00620     l1 = isave[1];
00621     l2 = isave[2];
00622     l3 = isave[3];
00623     lws = isave[4];
00624     lwy = isave[5];
00625     lsy = isave[6];
00626     lss = isave[7];
00627     lyy = isave[8];
00628     lwt = isave[9];
00629     lwn = isave[10];
00630     lsnd = isave[11];
00631     lz = isave[12];
00632     lr = isave[13];
00633     ld = isave[14];
00634     lt = isave[15];
00635     lwa = isave[16];
00636     lsg = isave[17];
00637     lsgo = isave[18];
00638     lyg = isave[19];
00639     lygo = isave[20];
00640     mainlb_(n, m, &x[1], &l[1], &u[1], &nbd[1], f, &g[1], factr, pgtol, &wa[
00641             lws], &wa[lwy], &wa[lsy], &wa[lss], &wa[lyy], &wa[lwt], &wa[lwn], 
00642             &wa[lsnd], &wa[lz], &wa[lr], &wa[ld], &wa[lt], &wa[lwa], &wa[lsg],
00643              &wa[lsgo], &wa[lyg], &wa[lygo], &iwa[1], &iwa[*n + 1], &iwa[(*n 
00644             << 1) + 1], task, iprint, csave, &lsave[1], &isave[22], &dsave[1],
00645              (long int)60, (long int)60);
00646     return 0;
00647 } /* setulb_ */
00648 
00649 /* ======================= The end of setulb ============================= */
00650 /* Subroutine */ int mainlb_(long int *n, long int *m, double *x, double *l, double *u, long int *nbd, double *f, 
00651                          double *g, double *factr, double *pgtol, double *ws, double *wy,
00652                          double *sy, double *ss, double *yy, double *wt, double *wn, double *snd, 
00653                          double *z__, double *r__, double *d__, double *t, double *wa, double *sg, 
00654                          double *sgo, double *yg, double *ygo, long int *index, long int *iwhere, long int *indx2, 
00655                          char *task, long int *iprint, char *csave, long int *lsave, long int *isave, double *dsave, 
00656                          long int task_len, long int csave_len)
00657 /*long int *n, *m;
00658 double *x, *l, *u;
00659 long int *nbd;
00660 double *f, *g, *factr, *pgtol, *ws, *wy, *sy, *ss, *yy, *wt, *wn, *snd, *
00661         z__, *r__, *d__, *t, *wa, *sg, *sgo, *yg, *ygo;
00662 long int *index, *iwhere, *indx2;
00663 char *task;
00664 long int *iprint;
00665 char *csave;
00666 long int *lsave;
00667 long int *isave;
00668 double *dsave;
00669 long int task_len;
00670 long int csave_len;*/
00671 {
00672     /* Format strings *//*
00673     static char fmt_1002[] = "(/,\002At iterate\002,i5,4x,\002f= \002,1p,d12\
00674 .5,4x,\002|proj g|= \002,1p,d12.5)";
00675     static char fmt_1003[] = "(2(1x,i4),5x,\002-\002,5x,\002-\002,3x,\002\
00676 -\002,5x,\002-\002,5x,\002-\002,8x,\002-\002,3x,1p,2(1x,d10.3))";
00677     static char fmt_1001[] = "(//,\002ITERATION \002,i5)";
00678     static char fmt_1005[] = "(/,\002 Singular triangular system detected\
00679 ;\002,/,\002   refresh the lbfgs memory and restart the iteration.\002)";
00680     static char fmt_1006[] = "(/,\002 Nonpositive definiteness in Cholesky f\
00681 actorization in formk;\002,/,\002   refresh the lbfgs memory and restart the\
00682  iteration.\002)";
00683     static char fmt_1008[] = "(/,\002 Bad direction in the line search;\002,\
00684 /,\002   refresh the lbfgs memory and restart the iteration.\002)";
00685     static char fmt_1004[] = "(\002  ys=\002,1p,e10.3,\002  -gs=\002,1p,e10.\
00686 3,\002 BFGS update SKIPPED\002)";
00687     static char fmt_1007[] = "(/,\002 Nonpositive definiteness in Cholesky f\
00688 actorization in formt;\002,/,\002   refresh the lbfgs memory and restart the\
00689  iteration.\002)";*/
00690 
00691     /* System generated locals */
00692     long int ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, 
00693             ss_dim1, ss_offset, yy_dim1, yy_offset, wt_dim1, wt_offset, 
00694             wn_dim1, wn_offset, snd_dim1, snd_offset, i__1;
00695     double d__1, d__2;
00696 //    olist o__1;
00697 
00698     /* Builtin functions */
00699     long int s_cmp(char *, const char *const, long int, long int);
00700     /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
00701     //long int f_open(olist *), s_wsfe(cilist *), do_fio(long int*, char*, long int), e_wsfe();
00702 
00703     /* Local variables */
00704     static long int head;
00705     static double fold;
00706     static long int nact;
00707     static double ddum;
00708     extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
00709     static long int info;
00710     static double time;
00711     static long int nfgv, ifun, iter, nint;
00712     static char word[3];
00713     static double time1, time2;
00714     static long int i__, iback, k;
00715     extern /* Subroutine */ int dscal_(long int *n, double *da, double *dx, long int *incx);
00716     static double gdold;
00717     static long int nfree;
00718     static long int boxed;
00719     static long int itail;
00720     static double theta;
00721     extern /* Subroutine */ int freev_(long int *n, long int *nfree, long int *index, long int *nenter, long int *ileave, long int *indx2, long int *iwhere, 
00722                                  long int *wrk, long int *updatd, long int *cnstnd, long int *iprint, long int *iter),
00723                                  dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
00724     static double dnorm;
00725     extern /* Subroutine */ int timer_(double *ttime), 
00726                           formk_(long int *n, long int *nsub, long int *ind, long int *nenter, long int *ileave, long int *indx2, long int *iupdat, 
00727                                  long int *updatd, double *wn, double *wn1, long int *m, double *ws, double *wy, double *sy, 
00728                                  double *theta, long int *col, long int *head, long int *info);
00729     static long int nskip, iword;
00730     extern /* Subroutine */ int formt_(long int *m, double *wt, double *sy, double *ss, long int *col, double *theta, long int *info), 
00731                            subsm_(long int *n, long int *m, long int *nsub, long int *ind, double *l, double *u, long int *nbd,
00732                                   double *x, double *d__, double *ws, double *wy, double *theta, long int *col,
00733                                   long int *head, long int *iword, double *wv, double *wn, long int *iprint, long int *info);
00734     static double xstep, stpmx;
00735     extern /* Subroutine */ int prn1lb_(long int*, long int*, double*, double*, double*, long int*, long int*, double*),
00736                            prn2lb_(long int *n, double *x, double *f, double *g, long int *iprint, long int *itfile,
00737                          long int *iter, long int *nfgv, long int *nact, double *sbgnrm, long int *nint, char *word,
00738                          long int *iword, long int *iback, double *stp, double *xstep, long int word_len), 
00739                          prn3lb_(long int *n, double *x, double *f, char *task, long int *iprint, long int *info, 
00740                          long int *itfile, long int *iter, long int *nfgv, long int *nintol, long int *nskip, long int *nact,
00741                          double *sbgnrm, double *time, long int *nint, char *word, long int *iback, double *stp,
00742                          double *xstep, long int *k, double *cachyt, double *sbtime, double *lnscht,
00743                          long int task_len, long int word_len);
00744     static double gd, dr, rr;
00745     static long int ileave;
00746     extern /* Subroutine */ int errclb_(long int *n, long int *m, double *factr, double *l, double *u, long int *nbd, char *task,
00747                          long int *info, long int *k, long int task_len);
00748     static long int itfile;
00749     static double cachyt, epsmch;
00750     static long int updatd;
00751     static double sbtime;
00752     extern /* Subroutine */ int active_(long int *n, double *l, double *u, long int *nbd, double *x, long int *iwhere, 
00753                          long int *iprint, long int *prjctd, long int *cnstnd, long int *boxed);
00754     static long int prjctd;
00755     static long int iupdat;
00756     extern double dpmeps_();
00757     static long int cnstnd;
00758     static double sbgnrm;
00759     static long int nenter;
00760     static double lnscht;
00761     extern /* Subroutine */ int cauchy_(long int *n, double *x, double *l, double *u, long int *nbd, double *g, long int *iorder,
00762                          long int *iwhere, double *t, double *d__, double *xcp, long int *m, double *wy, double *ws,
00763                          double *sy, double *wt, double *theta, long int *col, long int *head, double *p, double *c__,
00764                          double *wbp, double *v, long int *nint, double *sg, double *yg, long int *iprint, double *sbgnrm, long int *info),
00765                           cmprlb_(long int *n, long int *m, double *x, double *g, double *ws, double *wy, double *sy,
00766                          double *wt, double *z__, double *r__, double *wa, long int *index, double *theta,
00767                          long int *col, long int *head, long int *nfree, long int *cnstnd, long int *info), 
00768                           lnsrlb_(long int *n, double *l, double *u, long int *nbd, double *x, double *f, double *fold,
00769                          double *gd, double *gdold, double *g, double *d__, double *r__, double *t, 
00770                          double *z__, double *stp, double *dnorm, double *dtd, double *xstep, double *stpmx,
00771                          long int *iter, long int *ifun, long int *iback, long int *nfgv, long int *info, char *task, long int *boxed, 
00772                          long int *cnstnd, char *csave, long int *isave, double *dsave, long int task_len, long int csave_len), 
00773                           matupd_(long int *n, long int *m, double *ws, double *wy, double *sy, double *ss,
00774                          double *d__, double *r__, long int *itail, long int *iupdat, long int *col, long int *head,
00775                          double *theta, double *rr, double *dr, double *stp, double *dtd);
00776     static long int nintol;
00777     extern /* Subroutine */ int projgr_(long int *n, double *l, double *u, long int *nbd, double *x, double *g, double *sbgnrm);
00778     static double dtd;
00779     static long int col;
00780     static double tol;
00781     static long int wrk;
00782     static double stp, cpu1, cpu2;
00783 
00784     /* Fortran I/O blocks */
00785 /*    static cilist io___62 = { 0, 6, 0, fmt_1002, 0 };
00786     static cilist io___63 = { 0, 0, 0, fmt_1003, 0 };
00787     static cilist io___64 = { 0, 6, 0, fmt_1001, 0 };
00788     static cilist io___66 = { 0, 6, 0, fmt_1005, 0 };
00789     static cilist io___68 = { 0, 6, 0, fmt_1006, 0 };
00790     static cilist io___69 = { 0, 6, 0, fmt_1005, 0 };
00791     static cilist io___71 = { 0, 6, 0, fmt_1008, 0 };
00792     static cilist io___75 = { 0, 6, 0, fmt_1004, 0 };
00793     static cilist io___76 = { 0, 6, 0, fmt_1007, 0 };*/
00794 
00795 
00796 /*     ************ */
00797 
00798 /*     Subroutine mainlb */
00799 
00800 /*     This subroutine solves bound constrained optimization problems by */
00801 /*       using the compact formula of the limited memory BFGS updates. */
00802 
00803 /*     n is an long int variable. */
00804 /*       On entry n is the number of variables. */
00805 /*       On exit n is unchanged. */
00806 
00807 /*     m is an long int variable. */
00808 /*       On entry m is the maximum number of variable metric */
00809 /*          corrections allowed in the limited memory matrix. */
00810 /*       On exit m is unchanged. */
00811 
00812 /*     x is a double precision array of dimension n. */
00813 /*       On entry x is an approximation to the solution. */
00814 /*       On exit x is the current approximation. */
00815 
00816 /*     l is a double precision array of dimension n. */
00817 /*       On entry l is the lower bound of x. */
00818 /*       On exit l is unchanged. */
00819 
00820 /*     u is a double precision array of dimension n. */
00821 /*       On entry u is the upper bound of x. */
00822 /*       On exit u is unchanged. */
00823 
00824 /*     nbd is an long int array of dimension n. */
00825 /*       On entry nbd represents the type of bounds imposed on the */
00826 /*         variables, and must be specified as follows: */
00827 /*         nbd(i)=0 if x(i) is unbounded, */
00828 /*                1 if x(i) has only a lower bound, */
00829 /*                2 if x(i) has both lower and upper bounds, */
00830 /*                3 if x(i) has only an upper bound. */
00831 /*       On exit nbd is unchanged. */
00832 
00833 /*     f is a double precision variable. */
00834 /*       On first entry f is unspecified. */
00835 /*       On final exit f is the value of the function at x. */
00836 
00837 /*     g is a double precision array of dimension n. */
00838 /*       On first entry g is unspecified. */
00839 /*       On final exit g is the value of the gradient at x. */
00840 
00841 /*     factr is a double precision variable. */
00842 /*       On entry factr >= 0 is specified by the user.  The iteration */
00843 /*         will stop when */
00844 
00845 /*         (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch */
00846 
00847 /*         where epsmch is the machine precision, which is automatically */
00848 /*         generated by the code. */
00849 /*       On exit factr is unchanged. */
00850 
00851 /*     pgtol is a double precision variable. */
00852 /*       On entry pgtol >= 0 is specified by the user.  The iteration */
00853 /*         will stop when */
00854 
00855 /*                 max{|proj g_i | i = 1, ..., n} <= pgtol */
00856 
00857 /*         where pg_i is the ith component of the projected gradient. */
00858 /*       On exit pgtol is unchanged. */
00859 
00860 /*     ws, wy, sy, and wt are double precision working arrays used to */
00861 /*       store the following information defining the limited memory */
00862 /*          BFGS matrix: */
00863 /*          ws, of dimension n x m, stores S, the matrix of s-vectors; */
00864 /*          wy, of dimension n x m, stores Y, the matrix of y-vectors; */
00865 /*          sy, of dimension m x m, stores S'Y; */
00866 /*          ss, of dimension m x m, stores S'S; */
00867 /*         yy, of dimension m x m, stores Y'Y; */
00868 /*          wt, of dimension m x m, stores the Cholesky factorization */
00869 /*                                  of (theta*S'S+LD^(-1)L'); see eq. */
00870 /*                                  (2.26) in [3]. */
00871 
00872 /*     wn is a double precision working array of dimension 2m x 2m */
00873 /*       used to store the LEL^T factorization of the indefinite matrix */
00874 /*                 K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ] */
00875 /*                     [L_a -R_z           theta*S'AA'S ] */
00876 
00877 /*       where     E = [-I  0] */
00878 /*                     [ 0  I] */
00879 
00880 /*     snd is a double precision working array of dimension 2m x 2m */
00881 /*       used to store the lower triangular part of */
00882 /*                 N = [Y' ZZ'Y   L_a'+R_z'] */
00883 /*                     [L_a +R_z  S'AA'S   ] */
00884 
00885 /*     z(n),r(n),d(n),t(n),wa(8*m) are double precision working arrays. */
00886 /*       z is used at different times to store the Cauchy point and */
00887 /*       the Newton point. */
00888 
00889 /*     sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays. */
00890 
00891 /*     index is an long int working array of dimension n. */
00892 /*       In subroutine freev, index is used to store the free and fixed */
00893 /*          variables at the Generalized Cauchy Point (GCP). */
00894 
00895 /*     iwhere is an long int working array of dimension n used to record */
00896 /*       the status of the vector x for GCP computation. */
00897 /*       iwhere(i)=0 or -3 if x(i) is free and has bounds, */
00898 /*                 1       if x(i) is fixed at l(i), and l(i) .ne. u(i) */
00899 /*                 2       if x(i) is fixed at u(i), and u(i) .ne. l(i) */
00900 /*                 3       if x(i) is always fixed, i.e.,  u(i)=x(i)=l(i) */
00901 /*                -1       if x(i) is always free, i.e., no bounds on it. */
00902 
00903 /*     indx2 is an long int working array of dimension n. */
00904 /*       Within subroutine cauchy, indx2 corresponds to the array iorder. */
00905 /*       In subroutine freev, a list of variables entering and leaving */
00906 /*       the free set is stored in indx2, and it is passed on to */
00907 /*       subroutine formk with this information. */
00908 
00909 /*     task is a working string of characters of length 60 indicating */
00910 /*       the current job when entering and leaving this subroutine. */
00911 
00912 /*     iprint is an long int variable that must be set by the user. */
00913 /*       It controls the frequency and type of output generated: */
00914 /*        iprint<0    no output is generated; */
00915 /*        iprint=0    print only one line at the last iteration; */
00916 /*        0<iprint<99 print also f and |proj g| every iprint iterations; */
00917 /*        iprint=99   print details of every iteration except n-vectors; */
00918 /*        iprint=100  print also the changes of active set and final x; */
00919 /*        iprint>100  print details of every iteration including x and g; */
00920 /*       When iprint > 0, the file iterate.dat will be created to */
00921 /*                        summarize the iteration. */
00922 
00923 /*     csave is a working string of characters of length 60. */
00924 
00925 /*     lsave is a long int working array of dimension 4. */
00926 
00927 /*     isave is an long int working array of dimension 23. */
00928 
00929 /*     dsave is a double precision working array of dimension 29. */
00930 
00931 
00932 /*     Subprograms called */
00933 
00934 /*       L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk, */
00935 
00936 /*        errclb, prn1lb, prn2lb, prn3lb, active, projgr, */
00937 
00938 /*        freev, cmprlb, matupd, formt. */
00939 
00940 /*       Minpack2 Library ... timer, dpmeps. */
00941 
00942 /*       Linpack Library ... dcopy, ddot. */
00943 
00944 
00945 /*     References: */
00946 
00947 /*       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
00948 /*       memory algorithm for bound constrained optimization'', */
00949 /*       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
00950 
00951 /*       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN */
00952 /*       Subroutines for Large Scale Bound Constrained Optimization'' */
00953 /*       Tech. Report, NAM-11, EECS Department, Northwestern University, */
00954 /*       1994. */
00955 
00956 /*       [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of */
00957 /*       Quasi-Newton Matrices and their use in Limited Memory Methods'', */
00958 /*       Mathematical Programming 63 (1994), no. 4, pp. 129-156. */
00959 
00960 /*       (Postscript files of these papers are available via anonymous */
00961 /*        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
00962 
00963 /*                           *  *  * */
00964 
00965 /*     NEOS, November 1994. (Latest revision June 1996.) */
00966 /*     Optimization Technology Center. */
00967 /*     Argonne National Laboratory and Northwestern University. */
00968 /*     Written by */
00969 /*                        Ciyou Zhu */
00970 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
00971 
00972 
00973 /*     ************ */
00974     /* Parameter adjustments */
00975     --indx2;
00976     --iwhere;
00977     --index;
00978     --t;
00979     --d__;
00980     --r__;
00981     --z__;
00982     --g;
00983     --nbd;
00984     --u;
00985     --l;
00986     --x;
00987     --ygo;
00988     --yg;
00989     --sgo;
00990     --sg;
00991     --wa;
00992     snd_dim1 = 2 * *m;
00993     snd_offset = 1 + snd_dim1 * 1;
00994     snd -= snd_offset;
00995     wn_dim1 = 2 * *m;
00996     wn_offset = 1 + wn_dim1 * 1;
00997     wn -= wn_offset;
00998     wt_dim1 = *m;
00999     wt_offset = 1 + wt_dim1 * 1;
01000     wt -= wt_offset;
01001     yy_dim1 = *m;
01002     yy_offset = 1 + yy_dim1 * 1;
01003     yy -= yy_offset;
01004     ss_dim1 = *m;
01005     ss_offset = 1 + ss_dim1 * 1;
01006     ss -= ss_offset;
01007     sy_dim1 = *m;
01008     sy_offset = 1 + sy_dim1 * 1;
01009     sy -= sy_offset;
01010     wy_dim1 = *n;
01011     wy_offset = 1 + wy_dim1 * 1;
01012     wy -= wy_offset;
01013     ws_dim1 = *n;
01014     ws_offset = 1 + ws_dim1 * 1;
01015     ws -= ws_offset;
01016     --lsave;
01017     --isave;
01018     --dsave;
01019 
01020 
01021     /* Function Body */
01022     if (s_cmp(task, "START", (long int)60, (long int)5) == 0) {
01023         timer_(&time1);
01024 /*        Generate the current machine precision. */
01025         epsmch = dpmeps_();
01026 /*        Initialize counters and scalars when task='START'. */
01027 /*           for the limited memory BFGS matrices: */
01028         col = 0;
01029         head = 1;
01030         theta = 1.;
01031         iupdat = 0;
01032         updatd = FALSE_;
01033 /*           for operation counts: */
01034         iter = 0;
01035         nfgv = 0;
01036         nint = 0;
01037         nintol = 0;
01038         nskip = 0;
01039         nfree = *n;
01040 /*           for stopping tolerance: */
01041         tol = *factr * epsmch;
01042 /*           for measuring running time: */
01043         cachyt = 0.;
01044         sbtime = 0.;
01045         lnscht = 0.;
01046 /*           'word' records the status of subspace solutions. */
01047         s_copy(word, "---", (long int)3, (long int)3);
01048 /*           'info' records the termination information. */
01049         info = 0;
01050         if (*iprint >= 1) {
01051 /*                                open a summary file 'iterate.dat' */
01052 /*          o__1.oerr = 0;
01053             o__1.ounit = 8;
01054             o__1.ofnmlen = 11;
01055             o__1.ofnm = "iterate.dat";
01056             o__1.orl = 0;
01057             o__1.osta = "unknown";
01058             o__1.oacc = 0;
01059             o__1.ofm = 0;
01060             o__1.oblnk = 0;
01061             f_open(&o__1); */
01062             itfile = 8;
01063         }
01064 /*        Check the input arguments for errors. */
01065         errclb_(n, m, factr, &l[1], &u[1], &nbd[1], task, &info, &k, (long int)60);
01066         if (s_cmp(task, "ERROR", (long int)5, (long int)5) == 0) {
01067             prn3lb_(n, &x[1], f, task, iprint, &info, &itfile, &iter, &nfgv, &
01068                     nintol, &nskip, &nact, &sbgnrm, &c_b9, &nint, word, &
01069                     iback, &stp, &xstep, &k, &cachyt, &sbtime, &lnscht, (
01070                     long int)60, (long int)3);
01071             return 0;
01072         }
01073         prn1lb_(n, m, &l[1], &u[1], &x[1], iprint, &itfile, &epsmch);
01074 /*        Initialize iwhere & project x onto the feasible set. */
01075         active_(n, &l[1], &u[1], &nbd[1], &x[1], &iwhere[1], iprint, &prjctd, 
01076                 &cnstnd, &boxed);
01077 /*        The end of the initialization. */
01078     } else {
01079 /*          restore local variables. */
01080         prjctd = lsave[1];
01081         cnstnd = lsave[2];
01082         boxed = lsave[3];
01083         updatd = lsave[4];
01084         nintol = isave[1];
01085         itfile = isave[3];
01086         iback = isave[4];
01087         nskip = isave[5];
01088         head = isave[6];
01089         col = isave[7];
01090         itail = isave[8];
01091         iter = isave[9];
01092         iupdat = isave[10];
01093         nint = isave[12];
01094         nfgv = isave[13];
01095         info = isave[14];
01096         ifun = isave[15];
01097         iword = isave[16];
01098         nfree = isave[17];
01099         nact = isave[18];
01100         ileave = isave[19];
01101         nenter = isave[20];
01102         theta = dsave[1];
01103         fold = dsave[2];
01104         tol = dsave[3];
01105         dnorm = dsave[4];
01106         epsmch = dsave[5];
01107         cpu1 = dsave[6];
01108         cachyt = dsave[7];
01109         sbtime = dsave[8];
01110         lnscht = dsave[9];
01111         time1 = dsave[10];
01112         gd = dsave[11];
01113         stpmx = dsave[12];
01114         sbgnrm = dsave[13];
01115         stp = dsave[14];
01116         gdold = dsave[15];
01117         dtd = dsave[16];
01118 /*        After returning from the driver go to the point where execution */
01119 /*        is to resume. */
01120         if (s_cmp(task, "FG_LN", (long int)5, (long int)5) == 0) {
01121             goto L666;
01122         }
01123         if (s_cmp(task, "NEW_X", (long int)5, (long int)5) == 0) {
01124             goto L777;
01125         }
01126         if (s_cmp(task, "FG_ST", (long int)5, (long int)5) == 0) {
01127             goto L111;
01128         }
01129         if (s_cmp(task, "STOP", (long int)4, (long int)4) == 0) {
01130             if (s_cmp(task + 6, "CPU", (long int)3, (long int)3) == 0) {
01131 /*                                          restore the previous iterate. */
01132                 dcopy_(n, &t[1], &c__1, &x[1], &c__1);
01133                 dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
01134                 *f = fold;
01135             }
01136             goto L999;
01137         }
01138     }
01139 /*     Compute f0 and g0. */
01140     s_copy(task, "FG_START", (long int)60, (long int)8);
01141 /*          return to the driver to calculate f and g; reenter at 111. */
01142     goto L1000;
01143 L111:
01144     nfgv = 1;
01145 /*     Compute the infinity norm of the (-) projected gradient. */
01146     projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm);
01147     if (*iprint >= 1) {
01148 /*      s_wsfe(&io___62);
01149         do_fio(&c__1, (char *)&iter, (long int)sizeof(long int));
01150         do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
01151         do_fio(&c__1, (char *)&sbgnrm, (long int)sizeof(double));
01152         e_wsfe();
01153         io___63.ciunit = itfile;
01154         s_wsfe(&io___63);
01155         do_fio(&c__1, (char *)&iter, (long int)sizeof(long int));
01156         do_fio(&c__1, (char *)&nfgv, (long int)sizeof(long int));
01157         do_fio(&c__1, (char *)&sbgnrm, (long int)sizeof(double));
01158         do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
01159         e_wsfe();*/
01160     }
01161     if (sbgnrm <= *pgtol) {
01162 /*                                terminate the algorithm. */
01163         s_copy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL", (
01164                 long int)60, (long int)48);
01165         goto L999;
01166     }
01167 /* ----------------- the beginning of the loop -------------------------- */
01168 L222:
01169     if (*iprint >= 99) {
01170 /*      s_wsfe(&io___64);
01171         i__1 = iter + 1;
01172         do_fio(&c__1, (char *)&i__1, (long int)sizeof(long int));
01173         e_wsfe();*/
01174     }
01175     iword = -1;
01176 
01177     if (! cnstnd && col > 0) {
01178 /*                                            skip the search for GCP. */
01179         dcopy_(n, &x[1], &c__1, &z__[1], &c__1);
01180         wrk = updatd;
01181         nint = 0;
01182         goto L333;
01183     }
01184 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01185 
01186 /*     Compute the Generalized Cauchy Point (GCP). */
01187 
01188 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01189     timer_(&cpu1);
01190     cauchy_(n, &x[1], &l[1], &u[1], &nbd[1], &g[1], &indx2[1], &iwhere[1], &t[
01191             1], &d__[1], &z__[1], m, &wy[wy_offset], &ws[ws_offset], &sy[
01192             sy_offset], &wt[wt_offset], &theta, &col, &head, &wa[1], &wa[(*m 
01193             << 1) + 1], &wa[(*m << 2) + 1], &wa[*m * 6 + 1], &nint, &sg[1], &
01194             yg[1], iprint, &sbgnrm, &info);
01195     if (info != 0) {
01196 /*         singular triangular system detected; refresh the lbfgs memory. */
01197         if (*iprint >= 1) {
01198 /*          s_wsfe(&io___66);
01199             e_wsfe();*/
01200         }
01201         info = 0;
01202         col = 0;
01203         head = 1;
01204         theta = 1.;
01205         iupdat = 0;
01206         updatd = FALSE_;
01207         timer_(&cpu2);
01208         cachyt = cachyt + cpu2 - cpu1;
01209         goto L222;
01210     }
01211     timer_(&cpu2);
01212     cachyt = cachyt + cpu2 - cpu1;
01213     nintol += nint;
01214 /*     Count the entering and leaving variables for iter > 0; */
01215 /*     find the index set of free and active variables at the GCP. */
01216     freev_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iwhere[1], &
01217             wrk, &updatd, &cnstnd, iprint, &iter);
01218     nact = *n - nfree;
01219 L333:
01220 /*     If there are no free variables or B=I, then */
01221 /*                                        skip the subspace minimization. */
01222     if (nfree == 0 || col == 0) {
01223         goto L555;
01224     }
01225 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01226 
01227 /*     Subspace minimization. */
01228 
01229 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01230     timer_(&cpu1);
01231 /*     Form  the LEL^T factorization of the indefinite */
01232 /*       matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ] */
01233 /*                     [L_a -R_z           theta*S'AA'S ] */
01234 /*       where     E = [-I  0] */
01235 /*                     [ 0  I] */
01236     if (wrk) {
01237         formk_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iupdat, &
01238                 updatd, &wn[wn_offset], &snd[snd_offset], m, &ws[ws_offset], &
01239                 wy[wy_offset], &sy[sy_offset], &theta, &col, &head, &info);
01240     }
01241     if (info != 0) {
01242 /*          nonpositive definiteness in Cholesky factorization; */
01243 /*          refresh the lbfgs memory and restart the iteration. */
01244         if (*iprint >= 1) {
01245            /* s_wsfe(&io___68);
01246             e_wsfe();*/
01247         }
01248         info = 0;
01249         col = 0;
01250         head = 1;
01251         theta = 1.;
01252         iupdat = 0;
01253         updatd = FALSE_;
01254         timer_(&cpu2);
01255         sbtime = sbtime + cpu2 - cpu1;
01256         goto L222;
01257     }
01258 /*        compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) */
01259 /*                                                   from 'cauchy'). */
01260     cmprlb_(n, m, &x[1], &g[1], &ws[ws_offset], &wy[wy_offset], &sy[sy_offset]
01261             , &wt[wt_offset], &z__[1], &r__[1], &wa[1], &index[1], &theta, &
01262             col, &head, &nfree, &cnstnd, &info);
01263     if (info != 0) {
01264         goto L444;
01265     }
01266 /*       call the direct method. */
01267     subsm_(n, m, &nfree, &index[1], &l[1], &u[1], &nbd[1], &z__[1], &r__[1], &
01268             ws[ws_offset], &wy[wy_offset], &theta, &col, &head, &iword, &wa[1]
01269             , &wn[wn_offset], iprint, &info);
01270 L444:
01271     if (info != 0) {
01272 /*          singular triangular system detected; */
01273 /*          refresh the lbfgs memory and restart the iteration. */
01274         if (*iprint >= 1) {
01275 /*          s_wsfe(&io___69);
01276             e_wsfe();*/
01277         }
01278         info = 0;
01279         col = 0;
01280         head = 1;
01281         theta = 1.;
01282         iupdat = 0;
01283         updatd = FALSE_;
01284         timer_(&cpu2);
01285         sbtime = sbtime + cpu2 - cpu1;
01286         goto L222;
01287     }
01288     timer_(&cpu2);
01289     sbtime = sbtime + cpu2 - cpu1;
01290 L555:
01291 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01292 
01293 /*     Line search and optimality tests. */
01294 
01295 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01296 /*     Generate the search direction d:=z-x. */
01297     i__1 = *n;
01298     for (i__ = 1; i__ <= i__1; ++i__) {
01299         d__[i__] = z__[i__] - x[i__];
01300 /* L40: */
01301     }
01302     timer_(&cpu1);
01303 L666:
01304     lnsrlb_(n, &l[1], &u[1], &nbd[1], &x[1], f, &fold, &gd, &gdold, &g[1], &
01305             d__[1], &r__[1], &t[1], &z__[1], &stp, &dnorm, &dtd, &xstep, &
01306             stpmx, &iter, &ifun, &iback, &nfgv, &info, task, &boxed, &cnstnd, 
01307             csave, &isave[22], &dsave[17], (long int)60, (long int)60);
01308     if (info != 0 || iback >= 20) {
01309 /*          restore the previous iterate. */
01310         dcopy_(n, &t[1], &c__1, &x[1], &c__1);
01311         dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
01312         *f = fold;
01313         if (col == 0) {
01314 /*             abnormal termination. */
01315             if (info == 0) {
01316                 info = -9;
01317 /*                restore the actual number of f and g evaluations etc. */
01318                 --nfgv;
01319                 --ifun;
01320                 --iback;
01321             }
01322             s_copy(task, "ABNORMAL_TERMINATION_IN_LNSRCH", (long int)60, (
01323                     long int)30);
01324             ++iter;
01325             goto L999;
01326         } else {
01327 /*             refresh the lbfgs memory and restart the iteration. */
01328             if (*iprint >= 1) {
01329                 /*s_wsfe(&io___71);
01330                 e_wsfe();*/
01331             }
01332             if (info == 0) {
01333                 --nfgv;
01334             }
01335             info = 0;
01336             col = 0;
01337             head = 1;
01338             theta = 1.;
01339             iupdat = 0;
01340             updatd = FALSE_;
01341             s_copy(task, "RESTART_FROM_LNSRCH", (long int)60, (long int)19);
01342             timer_(&cpu2);
01343             lnscht = lnscht + cpu2 - cpu1;
01344             goto L222;
01345         }
01346     } else if (s_cmp(task, "FG_LN", (long int)5, (long int)5) == 0) {
01347 /*          return to the driver for calculating f and g; reenter at 666. */
01348         goto L1000;
01349     } else {
01350 /*          calculate and print out the quantities related to the new X. */
01351         timer_(&cpu2);
01352         lnscht = lnscht + cpu2 - cpu1;
01353         ++iter;
01354 /*        Compute the infinity norm of the projected (-)gradient. */
01355         projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm);
01356 /*        Print iteration information. */
01357         prn2lb_(n, &x[1], f, &g[1], iprint, &itfile, &iter, &nfgv, &nact, &
01358                 sbgnrm, &nint, word, &iword, &iback, &stp, &xstep, (long int)3);
01359         goto L1000;
01360     }
01361 L777:
01362 /*     Test for termination. */
01363     if (sbgnrm <= *pgtol) {
01364 /*                                terminate the algorithm. */
01365         s_copy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL", (
01366                 long int)60, (long int)48);
01367         goto L999;
01368     }
01369 /* Computing MAX */
01370     d__1 = abs_(fold), d__2 = abs_(*f), d__1 = max_(d__1,d__2);
01371     ddum = max_(d__1,1.);
01372     if (fold - *f <= tol * ddum) {
01373 /*                                        terminate the algorithm. */
01374         s_copy(task, "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH", (
01375                 long int)60, (long int)47);
01376         if (iback >= 10) {
01377             info = -5;
01378         }
01379 /*           i.e., to issue a warning if iback>10 in the line search. */
01380         goto L999;
01381     }
01382 /*     Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. */
01383     i__1 = *n;
01384     for (i__ = 1; i__ <= i__1; ++i__) {
01385         r__[i__] = g[i__] - r__[i__];
01386 /* L42: */
01387     }
01388     rr = ddot_(n, &r__[1], &c__1, &r__[1], &c__1);
01389     if (stp == 1.) {
01390         dr = gd - gdold;
01391         ddum = -gdold;
01392     } else {
01393         dr = (gd - gdold) * stp;
01394         dscal_(n, &stp, &d__[1], &c__1);
01395         ddum = -gdold * stp;
01396     }
01397     if (dr <= epsmch * ddum) {
01398 /*                            skip the L-BFGS update. */
01399         ++nskip;
01400         updatd = FALSE_;
01401         if (*iprint >= 1) {
01402            /* s_wsfe(&io___75);
01403             do_fio(&c__1, (char *)&dr, (long int)sizeof(double));
01404             do_fio(&c__1, (char *)&ddum, (long int)sizeof(double));
01405             e_wsfe();*/
01406         }
01407         goto L888;
01408     }
01409 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01410 
01411 /*     Update the L-BFGS matrix. */
01412 
01413 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01414     updatd = TRUE_;
01415     ++iupdat;
01416 /*     Update matrices WS and WY and form the middle matrix in B. */
01417     matupd_(n, m, &ws[ws_offset], &wy[wy_offset], &sy[sy_offset], &ss[
01418             ss_offset], &d__[1], &r__[1], &itail, &iupdat, &col, &head, &
01419             theta, &rr, &dr, &stp, &dtd);
01420 /*     Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; */
01421 /*        Store T in the upper triangular of the array wt; */
01422 /*        Cholesky factorize T to J*J' with */
01423 /*           J' stored in the upper triangular of wt. */
01424     formt_(m, &wt[wt_offset], &sy[sy_offset], &ss[ss_offset], &col, &theta, &
01425             info);
01426     if (info != 0) {
01427 /*          nonpositive definiteness in Cholesky factorization; */
01428 /*          refresh the lbfgs memory and restart the iteration. */
01429         if (*iprint >= 1) {
01430             /*s_wsfe(&io___76);
01431             e_wsfe();*/
01432         }
01433         info = 0;
01434         col = 0;
01435         head = 1;
01436         theta = 1.;
01437         iupdat = 0;
01438         updatd = FALSE_;
01439         goto L222;
01440     }
01441 /*     Now the inverse of the middle matrix in B is */
01442 /*       [  D^(1/2)      O ] [ -D^(1/2)  D^(-1/2)*L' ] */
01443 /*       [ -L*D^(-1/2)   J ] [  0        J'          ] */
01444 L888:
01445 /* -------------------- the end of the loop ----------------------------- */
01446     goto L222;
01447 L999:
01448     timer_(&time2);
01449     time = time2 - time1;
01450     prn3lb_(n, &x[1], f, task, iprint, &info, &itfile, &iter, &nfgv, &nintol, 
01451             &nskip, &nact, &sbgnrm, &time, &nint, word, &iback, &stp, &xstep, 
01452             &k, &cachyt, &sbtime, &lnscht, (long int)60, (long int)3);
01453 L1000:
01454 /*     Save local variables. */
01455     lsave[1] = prjctd;
01456     lsave[2] = cnstnd;
01457     lsave[3] = boxed;
01458     lsave[4] = updatd;
01459     isave[1] = nintol;
01460     isave[3] = itfile;
01461     isave[4] = iback;
01462     isave[5] = nskip;
01463     isave[6] = head;
01464     isave[7] = col;
01465     isave[8] = itail;
01466     isave[9] = iter;
01467     isave[10] = iupdat;
01468     isave[12] = nint;
01469     isave[13] = nfgv;
01470     isave[14] = info;
01471     isave[15] = ifun;
01472     isave[16] = iword;
01473     isave[17] = nfree;
01474     isave[18] = nact;
01475     isave[19] = ileave;
01476     isave[20] = nenter;
01477     dsave[1] = theta;
01478     dsave[2] = fold;
01479     dsave[3] = tol;
01480     dsave[4] = dnorm;
01481     dsave[5] = epsmch;
01482     dsave[6] = cpu1;
01483     dsave[7] = cachyt;
01484     dsave[8] = sbtime;
01485     dsave[9] = lnscht;
01486     dsave[10] = time1;
01487     dsave[11] = gd;
01488     dsave[12] = stpmx;
01489     dsave[13] = sbgnrm;
01490     dsave[14] = stp;
01491     dsave[15] = gdold;
01492     dsave[16] = dtd;
01493     return 0;
01494 } /* mainlb_ */
01495 
01496 /* ======================= The end of mainlb ============================= */
01497 /* Subroutine */ int active_(long int *n, double *l, double *u, long int *nbd, double *x, long int *iwhere, 
01498                          long int *iprint, long int *prjctd, long int *cnstnd, long int *boxed)
01499 /*long int *n;
01500 double *l, *u;
01501 long int *nbd;
01502 double *x;
01503 long int *iwhere, *iprint;
01504 long int *prjctd, *cnstnd, *boxed;*/
01505 {
01506     /* Format strings *//*
01507     static char fmt_1001[] = "(/,\002At X0 \002,i9,\002 variables are exactl\
01508 y at the bounds\002)";*/
01509 
01510     /* System generated locals */
01511     long int i__1;
01512 
01513     /* Builtin functions */
01514     //long int s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle(), s_wsfe(cilist *), do_fio(long int*, char*, long int), e_wsfe();
01515 
01516     /* Local variables */
01517     static long int nbdd, i__;
01518 
01519     /* Fortran I/O blocks */
01520 /*    static cilist io___81 = { 0, 6, 0, 0, 0 };
01521     static cilist io___82 = { 0, 6, 0, 0, 0 };
01522     static cilist io___83 = { 0, 6, 0, fmt_1001, 0 };*/
01523 
01524 
01525 /*     ************ */
01526 
01527 /*     Subroutine active */
01528 
01529 /*     This subroutine initializes iwhere and projects the initial x to */
01530 /*       the feasible set if necessary. */
01531 
01532 /*     iwhere is an long int array of dimension n. */
01533 /*       On entry iwhere is unspecified. */
01534 /*       On exit iwhere(i)=-1  if x(i) has no bounds */
01535 /*                         3   if l(i)=u(i) */
01536 /*                         0   otherwise. */
01537 /*       In cauchy, iwhere is given finer gradations. */
01538 
01539 
01540 /*                           *  *  * */
01541 
01542 /*     NEOS, November 1994. (Latest revision June 1996.) */
01543 /*     Optimization Technology Center. */
01544 /*     Argonne National Laboratory and Northwestern University. */
01545 /*     Written by */
01546 /*                        Ciyou Zhu */
01547 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
01548 
01549 
01550 /*     ************ */
01551 /*     Initialize nbdd, prjctd, cnstnd and boxed. */
01552     /* Parameter adjustments */
01553     --iwhere;
01554     --x;
01555     --nbd;
01556     --u;
01557     --l;
01558 
01559     /* Function Body */
01560     nbdd = 0;
01561     *prjctd = FALSE_;
01562     *cnstnd = FALSE_;
01563     *boxed = TRUE_;
01564 /*     Project the initial x to the easible set if necessary. */
01565     i__1 = *n;
01566     for (i__ = 1; i__ <= i__1; ++i__) {
01567         if (nbd[i__] > 0) {
01568             if (nbd[i__] <= 2 && x[i__] <= l[i__]) {
01569                 if (x[i__] < l[i__]) {
01570                     *prjctd = TRUE_;
01571                     x[i__] = l[i__];
01572                 }
01573                 ++nbdd;
01574             } else if (nbd[i__] >= 2 && x[i__] >= u[i__]) {
01575                 if (x[i__] > u[i__]) {
01576                     *prjctd = TRUE_;
01577                     x[i__] = u[i__];
01578                 }
01579                 ++nbdd;
01580             }
01581         }
01582 /* L10: */
01583     }
01584 /*     Initialize iwhere and assign values to cnstnd and boxed. */
01585     i__1 = *n;
01586     for (i__ = 1; i__ <= i__1; ++i__) {
01587         if (nbd[i__] != 2) {
01588             *boxed = FALSE_;
01589         }
01590         if (nbd[i__] == 0) {
01591 /*                                this variable is always free */
01592             iwhere[i__] = -1;
01593 /*           otherwise set x(i)=mid(x(i), u(i), l(i)). */
01594         } else {
01595             *cnstnd = TRUE_;
01596             if (nbd[i__] == 2 && u[i__] - l[i__] <= 0.) {
01597 /*                   this variable is always fixed */
01598                 iwhere[i__] = 3;
01599             } else {
01600                 iwhere[i__] = 0;
01601             }
01602         }
01603 /* L20: */
01604     }
01605     if (*iprint >= 0) {
01606 /*      if (*prjctd) {
01607             s_wsle(&io___81);
01608             do_lio(&c__9, &c__1, "The initial X is infeasible.  Restart with\
01609  its projection.", (long int)58);
01610             e_wsle();
01611         }
01612         if (! (*cnstnd)) {
01613             s_wsle(&io___82);
01614             do_lio(&c__9, &c__1, "This problem is unconstrained.", (long int)30)
01615                     ;
01616             e_wsle();
01617         }*/
01618     }
01619     if (*iprint > 0) {
01620 /*      s_wsfe(&io___83);
01621         do_fio(&c__1, (char *)&nbdd, (long int)sizeof(long int));
01622         e_wsfe();*/
01623     }
01624     return 0;
01625 } /* active_ */
01626 
01627 /* ======================= The end of active ============================= */
01628 /* Subroutine */ int bmv_(long int *m, double *sy, double *wt, long int *col, double *v, double *p, long int *info)
01629 /*long int *m;
01630 double *sy, *wt;
01631 long int *col;
01632 double *v, *p;
01633 long int *info;*/
01634 {
01635     /* System generated locals */
01636     long int sy_dim1, sy_offset, wt_dim1, wt_offset, i__1, i__2;
01637 
01638     /* Builtin functions */
01639     //double sqrt();
01640 
01641     /* Local variables */
01642     static long int i__, k;
01643     extern /* Subroutine */ int dtrsl_(double *t, long int *ldt, long int *n, double *b, long int *job, long int *info);
01644     static long int i2;
01645     static double sum;
01646 
01647 /*     ************ */
01648 
01649 /*     Subroutine bmv */
01650 
01651 /*     This subroutine computes the product of the 2m x 2m middle matrix */
01652 /*      in the compact L-BFGS formula of B and a 2m vector v; */
01653 /*      it returns the product in p. */
01654 
01655 /*     m is an long int variable. */
01656 /*       On entry m is the maximum number of variable metric corrections */
01657 /*         used to define the limited memory matrix. */
01658 /*       On exit m is unchanged. */
01659 
01660 /*     sy is a double precision array of dimension m x m. */
01661 /*       On entry sy specifies the matrix S'Y. */
01662 /*       On exit sy is unchanged. */
01663 
01664 /*     wt is a double precision array of dimension m x m. */
01665 /*       On entry wt specifies the upper triangular matrix J' which is */
01666 /*         the Cholesky factor of (thetaS'S+LD^(-1)L'). */
01667 /*       On exit wt is unchanged. */
01668 
01669 /*     col is an long int variable. */
01670 /*       On entry col specifies the number of s-vectors (or y-vectors) */
01671 /*         stored in the compact L-BFGS formula. */
01672 /*       On exit col is unchanged. */
01673 
01674 /*     v is a double precision array of dimension 2col. */
01675 /*       On entry v specifies vector v. */
01676 /*       On exit v is unchanged. */
01677 
01678 /*     p is a double precision array of dimension 2col. */
01679 /*       On entry p is unspecified. */
01680 /*       On exit p is the product Mv. */
01681 
01682 /*     info is an long int variable. */
01683 /*       On entry info is unspecified. */
01684 /*       On exit info = 0       for normal return, */
01685 /*                    = nonzero for abnormal return when the system */
01686 /*                                to be solved by dtrsl is singular. */
01687 
01688 /*     Subprograms called: */
01689 
01690 /*       Linpack ... dtrsl. */
01691 
01692 
01693 /*                           *  *  * */
01694 
01695 /*     NEOS, November 1994. (Latest revision June 1996.) */
01696 /*     Optimization Technology Center. */
01697 /*     Argonne National Laboratory and Northwestern University. */
01698 /*     Written by */
01699 /*                        Ciyou Zhu */
01700 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
01701 
01702 
01703 /*     ************ */
01704     /* Parameter adjustments */
01705     wt_dim1 = *m;
01706     wt_offset = 1 + wt_dim1 * 1;
01707     wt -= wt_offset;
01708     sy_dim1 = *m;
01709     sy_offset = 1 + sy_dim1 * 1;
01710     sy -= sy_offset;
01711     --p;
01712     --v;
01713 
01714     /* Function Body */
01715     if (*col == 0) {
01716         return 0;
01717     }
01718 /*     PART I: solve [  D^(1/2)      O ] [ p1 ] = [ v1 ] */
01719 /*                   [ -L*D^(-1/2)   J ] [ p2 ]   [ v2 ]. */
01720 /*      solve Jp2=v2+LD^(-1)v1. */
01721     p[*col + 1] = v[*col + 1];
01722     i__1 = *col;
01723     for (i__ = 2; i__ <= i__1; ++i__) {
01724         i2 = *col + i__;
01725         sum = 0.;
01726         i__2 = i__ - 1;
01727         for (k = 1; k <= i__2; ++k) {
01728             sum += sy[i__ + k * sy_dim1] * v[k] / sy[k + k * sy_dim1];
01729 /* L10: */
01730         }
01731         p[i2] = v[i2] + sum;
01732 /* L20: */
01733     }
01734 /*     Solve the triangular system */
01735     dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__11, info);
01736     if (*info != 0) {
01737         return 0;
01738     }
01739 /*      solve D^(1/2)p1=v1. */
01740     i__1 = *col;
01741     for (i__ = 1; i__ <= i__1; ++i__) {
01742         p[i__] = v[i__] / sqrt(sy[i__ + i__ * sy_dim1]);
01743 /* L30: */
01744     }
01745 /*     PART II: solve [ -D^(1/2)   D^(-1/2)*L'  ] [ p1 ] = [ p1 ] */
01746 /*                    [  0         J'           ] [ p2 ]   [ p2 ]. */
01747 /*       solve J^Tp2=p2. */
01748     dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__1, info);
01749     if (*info != 0) {
01750         return 0;
01751     }
01752 /*       compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2) */
01753 /*                 =-D^(-1/2)p1+D^(-1)L'p2. */
01754     i__1 = *col;
01755     for (i__ = 1; i__ <= i__1; ++i__) {
01756         p[i__] = -p[i__] / sqrt(sy[i__ + i__ * sy_dim1]);
01757 /* L40: */
01758     }
01759     i__1 = *col;
01760     for (i__ = 1; i__ <= i__1; ++i__) {
01761         sum = 0.;
01762         i__2 = *col;
01763         for (k = i__ + 1; k <= i__2; ++k) {
01764             sum += sy[k + i__ * sy_dim1] * p[*col + k] / sy[i__ + i__ * 
01765                     sy_dim1];
01766 /* L50: */
01767         }
01768         p[i__] += sum;
01769 /* L60: */
01770     }
01771     return 0;
01772 } /* bmv_ */
01773 
01774 /* ======================== The end of bmv =============================== */
01775 /* Subroutine */ int cauchy_(long int *n, double *x, double *l, double *u, long int *nbd, double *g, long int *iorder,
01776                          long int *iwhere, double *t, double *d__, double *xcp, long int *m, double *wy, double *ws,
01777                          double *sy, double *wt, double *theta, long int *col, long int *head, double *p, double *c__,
01778                          double *wbp, double *v, long int *nint, double *sg, double *yg, long int *iprint, double *sbgnrm, long int *info)
01779 /*long int *n;
01780 double *x, *l, *u;
01781 long int *nbd;
01782 double *g;
01783 long int *iorder, *iwhere;
01784 double *t, *d__, *xcp;
01785 long int *m;
01786 double *wy, *ws, *sy, *wt, *theta;
01787 long int *col, *head;
01788 double *p, *c__, *wbp, *v;
01789 long int *nint;
01790 double *sg, *yg;
01791 long int *iprint;
01792 double *sbgnrm;
01793 long int *info;*/
01794 {
01795     /* Format strings *//*
01796     static char fmt_3010[] = "(/,\002---------------- CAUCHY entered--------\
01797 -----------\002)";
01798     static char fmt_1010[] = "(\002Cauchy X =  \002,/,(4x,1p,6(1x,d11.4)))";
01799     static char fmt_4011[] = "(/,\002Piece    \002,i3,\002 --f1, f2 at start\
01800  point \002,1p,2(1x,d11.4))";
01801     static char fmt_5010[] = "(\002Distance to the next break point =  \002,\
01802 1p,d11.4)";
01803     static char fmt_6010[] = "(\002Distance to the stationary point =  \002,\
01804 1p,d11.4)";
01805     static char fmt_4010[] = "(\002Piece    \002,i3,\002 --f1, f2 at start p\
01806 oint \002,1p,2(1x,d11.4))";
01807     static char fmt_2010[] = "(/,\002---------------- exit CAUCHY-----------\
01808 -----------\002,/)";*/
01809 
01810     /* System generated locals */
01811     long int wy_dim1, wy_offset, ws_dim1, ws_offset, sy_dim1, sy_offset, 
01812             wt_dim1, wt_offset, i__1, i__2;
01813     double d__1;
01814 
01815     /* Builtin functions */
01816 //    long int s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle(), s_wsfe(cilist *), e_wsfe(), do_fio(long int*, char*, long int);
01817 
01818     /* Local variables */
01819     static double dibp;
01820     extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
01821     static long int iter;
01822     static double zibp, tsum, dibp2;
01823     static long int i__, j;
01824     static long int bnded;
01825     extern /* Subroutine */ int dscal_(long int *n, double *da, double *dx, long int *incx);
01826     static double neggi;
01827     static long int nfree;
01828     static double bkmin;
01829     static long int nleft;
01830     extern /* Subroutine */ int dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy),
01831                           daxpy_(long int *n, double *da, double *dx, long int *incx, double *dy, long int *incy);
01832     static double f1, f2, dt, tj, tl;
01833     static long int nbreak, ibkmin;
01834     static double tu;
01835     extern /* Subroutine */ int hpsolb_(long int *n, double *t, long int *iorder, long int *iheap);
01836     static long int pointr;
01837     static double tj0;
01838     static long int xlower, xupper;
01839     static long int ibp;
01840     static double dtm;
01841     extern /* Subroutine */ int bmv_(long int *m, double *sy, double *wt, long int *col, double *v, double *p, long int *info);
01842     static double wmc, wmp, wmw;
01843     static long int col2;
01844 
01845     /* Fortran I/O blocks */
01846 /*    static cilist io___88 = { 0, 6, 0, 0, 0 };
01847     static cilist io___96 = { 0, 6, 0, fmt_3010, 0 };
01848     static cilist io___105 = { 0, 6, 0, fmt_1010, 0 };
01849     static cilist io___109 = { 0, 6, 0, 0, 0 };
01850     static cilist io___116 = { 0, 6, 0, fmt_4011, 0 };
01851     static cilist io___117 = { 0, 6, 0, fmt_5010, 0 };
01852     static cilist io___118 = { 0, 6, 0, fmt_6010, 0 };
01853     static cilist io___121 = { 0, 6, 0, 0, 0 };
01854     static cilist io___126 = { 0, 6, 0, 0, 0 };
01855     static cilist io___127 = { 0, 6, 0, 0, 0 };
01856     static cilist io___128 = { 0, 6, 0, fmt_4010, 0 };
01857     static cilist io___129 = { 0, 6, 0, fmt_6010, 0 };
01858     static cilist io___130 = { 0, 6, 0, fmt_1010, 0 };
01859     static cilist io___131 = { 0, 6, 0, fmt_2010, 0 };*/
01860 
01861 
01862 /*     ************ */
01863 
01864 /*     Subroutine cauchy */
01865 
01866 /*     For given x, l, u, g (with sbgnrm > 0), and a limited memory */
01867 /*       BFGS matrix B defined in terms of matrices WY, WS, WT, and */
01868 /*       scalars head, col, and theta, this subroutine computes the */
01869 /*       generalized Cauchy point (GCP), defined as the first local */
01870 /*       minimizer of the quadratic */
01871 
01872 /*                  Q(x + s) = g's + 1/2 s'Bs */
01873 
01874 /*       along the projected gradient direction P(x-tg,l,u). */
01875 /*       The routine returns the GCP in xcp. */
01876 
01877 /*     n is an long int variable. */
01878 /*       On entry n is the dimension of the problem. */
01879 /*       On exit n is unchanged. */
01880 
01881 /*     x is a double precision array of dimension n. */
01882 /*       On entry x is the starting point for the GCP computation. */
01883 /*       On exit x is unchanged. */
01884 
01885 /*     l is a double precision array of dimension n. */
01886 /*       On entry l is the lower bound of x. */
01887 /*       On exit l is unchanged. */
01888 
01889 /*     u is a double precision array of dimension n. */
01890 /*       On entry u is the upper bound of x. */
01891 /*       On exit u is unchanged. */
01892 
01893 /*     nbd is an long int array of dimension n. */
01894 /*       On entry nbd represents the type of bounds imposed on the */
01895 /*         variables, and must be specified as follows: */
01896 /*         nbd(i)=0 if x(i) is unbounded, */
01897 /*                1 if x(i) has only a lower bound, */
01898 /*                2 if x(i) has both lower and upper bounds, and */
01899 /*                3 if x(i) has only an upper bound. */
01900 /*       On exit nbd is unchanged. */
01901 
01902 /*     g is a double precision array of dimension n. */
01903 /*       On entry g is the gradient of f(x).  g must be a nonzero vector. */
01904 /*       On exit g is unchanged. */
01905 
01906 /*     iorder is an long int working array of dimension n. */
01907 /*       iorder will be used to store the breakpoints in the piecewise */
01908 /*       linear path and free variables encountered. On exit, */
01909 /*         iorder(1),...,iorder(nleft) are indices of breakpoints */
01910 /*                                which have not been encountered; */
01911 /*         iorder(nleft+1),...,iorder(nbreak) are indices of */
01912 /*                                     encountered breakpoints; and */
01913 /*         iorder(nfree),...,iorder(n) are indices of variables which */
01914 /*                 have no bound constraits along the search direction. */
01915 
01916 /*     iwhere is an long int array of dimension n. */
01917 /*       On entry iwhere indicates only the permanently fixed (iwhere=3) */
01918 /*       or free (iwhere= -1) components of x. */
01919 /*       On exit iwhere records the status of the current x variables. */
01920 /*       iwhere(i)=-3  if x(i) is free and has bounds, but is not moved */
01921 /*                 0   if x(i) is free and has bounds, and is moved */
01922 /*                 1   if x(i) is fixed at l(i), and l(i) .ne. u(i) */
01923 /*                 2   if x(i) is fixed at u(i), and u(i) .ne. l(i) */
01924 /*                 3   if x(i) is always fixed, i.e.,  u(i)=x(i)=l(i) */
01925 /*                 -1  if x(i) is always free, i.e., it has no bounds. */
01926 
01927 /*     t is a double precision working array of dimension n. */
01928 /*       t will be used to store the break points. */
01929 
01930 /*     d is a double precision array of dimension n used to store */
01931 /*       the Cauchy direction P(x-tg)-x. */
01932 
01933 /*     xcp is a double precision array of dimension n used to return the */
01934 /*       GCP on exit. */
01935 
01936 /*     m is an long int variable. */
01937 /*       On entry m is the maximum number of variable metric corrections */
01938 /*         used to define the limited memory matrix. */
01939 /*       On exit m is unchanged. */
01940 
01941 /*     ws, wy, sy, and wt are double precision arrays. */
01942 /*       On entry they store information that defines the */
01943 /*                             limited memory BFGS matrix: */
01944 /*         ws(n,m) stores S, a set of s-vectors; */
01945 /*         wy(n,m) stores Y, a set of y-vectors; */
01946 /*         sy(m,m) stores S'Y; */
01947 /*         wt(m,m) stores the */
01948 /*                 Cholesky factorization of (theta*S'S+LD^(-1)L'). */
01949 /*       On exit these arrays are unchanged. */
01950 
01951 /*     theta is a double precision variable. */
01952 /*       On entry theta is the scaling factor specifying B_0 = theta I. */
01953 /*       On exit theta is unchanged. */
01954 
01955 /*     col is an long int variable. */
01956 /*       On entry col is the actual number of variable metric */
01957 /*         corrections stored so far. */
01958 /*       On exit col is unchanged. */
01959 
01960 /*     head is an long int variable. */
01961 /*       On entry head is the location of the first s-vector (or y-vector) */
01962 /*         in S (or Y). */
01963 /*       On exit col is unchanged. */
01964 
01965 /*     p is a double precision working array of dimension 2m. */
01966 /*       p will be used to store the vector p = W^(T)d. */
01967 
01968 /*     c is a double precision working array of dimension 2m. */
01969 /*       c will be used to store the vector c = W^(T)(xcp-x). */
01970 
01971 /*     wbp is a double precision working array of dimension 2m. */
01972 /*       wbp will be used to store the row of W corresponding */
01973 /*         to a breakpoint. */
01974 
01975 /*     v is a double precision working array of dimension 2m. */
01976 
01977 /*     nint is an long int variable. */
01978 /*       On exit nint records the number of quadratic segments explored */
01979 /*         in searching for the GCP. */
01980 
01981 /*     sg and yg are double precision arrays of dimension m. */
01982 /*       On entry sg  and yg store S'g and Y'g correspondingly. */
01983 /*       On exit they are unchanged. */
01984 
01985 /*     iprint is an long int variable that must be set by the user. */
01986 /*       It controls the frequency and type of output generated: */
01987 /*        iprint<0    no output is generated; */
01988 /*        iprint=0    print only one line at the last iteration; */
01989 /*        0<iprint<99 print also f and |proj g| every iprint iterations; */
01990 /*        iprint=99   print details of every iteration except n-vectors; */
01991 /*        iprint=100  print also the changes of active set and final x; */
01992 /*        iprint>100  print details of every iteration including x and g; */
01993 /*       When iprint > 0, the file iterate.dat will be created to */
01994 /*                        summarize the iteration. */
01995 
01996 /*     sbgnrm is a double precision variable. */
01997 /*       On entry sbgnrm is the norm of the projected gradient at x. */
01998 /*       On exit sbgnrm is unchanged. */
01999 
02000 /*     info is an long int variable. */
02001 /*       On entry info is 0. */
02002 /*       On exit info = 0       for normal return, */
02003 /*                    = nonzero for abnormal return when the the system */
02004 /*                              used in routine bmv is singular. */
02005 
02006 /*     Subprograms called: */
02007 
02008 /*       L-BFGS-B Library ... hpsolb, bmv. */
02009 
02010 /*       Linpack ... dscal dcopy, daxpy. */
02011 
02012 
02013 /*     References: */
02014 
02015 /*       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
02016 /*       memory algorithm for bound constrained optimization'', */
02017 /*       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
02018 
02019 /*       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN */
02020 /*       Subroutines for Large Scale Bound Constrained Optimization'' */
02021 /*       Tech. Report, NAM-11, EECS Department, Northwestern University, */
02022 /*       1994. */
02023 
02024 /*       (Postscript files of these papers are available via anonymous */
02025 /*        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
02026 
02027 /*                           *  *  * */
02028 
02029 /*     NEOS, November 1994. (Latest revision June 1996.) */
02030 /*     Optimization Technology Center. */
02031 /*     Argonne National Laboratory and Northwestern University. */
02032 /*     Written by */
02033 /*                        Ciyou Zhu */
02034 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
02035 
02036 
02037 /*     ************ */
02038 /*     Check the status of the variables, reset iwhere(i) if necessary; */
02039 /*       compute the Cauchy direction d and the breakpoints t; initialize */
02040 /*       the derivative f1 and the vector p = W'd (for theta = 1). */
02041     /* Parameter adjustments */
02042     --xcp;
02043     --d__;
02044     --t;
02045     --iwhere;
02046     --iorder;
02047     --g;
02048     --nbd;
02049     --u;
02050     --l;
02051     --x;
02052     --yg;
02053     --sg;
02054     --v;
02055     --wbp;
02056     --c__;
02057     --p;
02058     wt_dim1 = *m;
02059     wt_offset = 1 + wt_dim1 * 1;
02060     wt -= wt_offset;
02061     sy_dim1 = *m;
02062     sy_offset = 1 + sy_dim1 * 1;
02063     sy -= sy_offset;
02064     ws_dim1 = *n;
02065     ws_offset = 1 + ws_dim1 * 1;
02066     ws -= ws_offset;
02067     wy_dim1 = *n;
02068     wy_offset = 1 + wy_dim1 * 1;
02069     wy -= wy_offset;
02070 
02071     /* Function Body */
02072     if (*sbgnrm <= 0.) {
02073         if (*iprint >= 0) {
02074 /*          s_wsle(&io___88);
02075             do_lio(&c__9, &c__1, "Subgnorm = 0.  GCP = X.", (long int)23);
02076             e_wsle();*/
02077         }
02078         dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
02079         return 0;
02080     }
02081     bnded = TRUE_;
02082     nfree = *n + 1;
02083     nbreak = 0;
02084     ibkmin = 0;
02085     bkmin = 0.;
02086     col2 = *col << 1;
02087     f1 = 0.;
02088     if (*iprint >= 99) {
02089 /*      s_wsfe(&io___96);
02090         e_wsfe();*/
02091     }
02092 /*     We set p to zero and build it up as we determine d. */
02093     i__1 = col2;
02094     for (i__ = 1; i__ <= i__1; ++i__) {
02095         p[i__] = 0.;
02096 /* L20: */
02097     }
02098 /*     In the following loop we determine for each variable its bound */
02099 /*        status and its breakpoint, and update p accordingly. */
02100 /*        Smallest breakpoint is identified. */
02101     i__1 = *n;
02102     for (i__ = 1; i__ <= i__1; ++i__) {
02103         neggi = -g[i__];
02104         if (iwhere[i__] != 3 && iwhere[i__] != -1) {
02105 /*             if x(i) is not a constant and has bounds, */
02106 /*             compute the difference between x(i) and its bounds. */
02107             if (nbd[i__] <= 2) {
02108                 tl = x[i__] - l[i__];
02109             }
02110             if (nbd[i__] >= 2) {
02111                 tu = u[i__] - x[i__];
02112             }
02113 /*           If a variable is close enough to a bound */
02114 /*             we treat it as at bound. */
02115             xlower = nbd[i__] <= 2 && tl <= 0.;
02116             xupper = nbd[i__] >= 2 && tu <= 0.;
02117 /*              reset iwhere(i). */
02118             iwhere[i__] = 0;
02119             if (xlower) {
02120                 if (neggi <= 0.) {
02121                     iwhere[i__] = 1;
02122                 }
02123             } else if (xupper) {
02124                 if (neggi >= 0.) {
02125                     iwhere[i__] = 2;
02126                 }
02127             } else {
02128                 if (abs_(neggi) <= 0.) {
02129                     iwhere[i__] = -3;
02130                 }
02131             }
02132         }
02133         pointr = *head;
02134         if (iwhere[i__] != 0 && iwhere[i__] != -1) {
02135             d__[i__] = 0.;
02136         } else {
02137             d__[i__] = neggi;
02138             f1 -= neggi * neggi;
02139 /*             calculate p := p - W'e_i* (g_i). */
02140             i__2 = *col;
02141             for (j = 1; j <= i__2; ++j) {
02142                 p[j] += wy[i__ + pointr * wy_dim1] * neggi;
02143                 p[*col + j] += ws[i__ + pointr * ws_dim1] * neggi;
02144                 pointr = pointr % *m + 1;
02145 /* L40: */
02146             }
02147             if (nbd[i__] <= 2 && nbd[i__] != 0 && neggi < 0.) {
02148 /*                                 x(i) + d(i) is bounded; compute t(i). */
02149                 ++nbreak;
02150                 iorder[nbreak] = i__;
02151                 t[nbreak] = tl / (-neggi);
02152                 if (nbreak == 1 || t[nbreak] < bkmin) {
02153                     bkmin = t[nbreak];
02154                     ibkmin = nbreak;
02155                 }
02156             } else if (nbd[i__] >= 2 && neggi > 0.) {
02157 /*                                 x(i) + d(i) is bounded; compute t(i). */
02158                 ++nbreak;
02159                 iorder[nbreak] = i__;
02160                 t[nbreak] = tu / neggi;
02161                 if (nbreak == 1 || t[nbreak] < bkmin) {
02162                     bkmin = t[nbreak];
02163                     ibkmin = nbreak;
02164                 }
02165             } else {
02166 /*                x(i) + d(i) is not bounded. */
02167                 --nfree;
02168                 iorder[nfree] = i__;
02169                 if (abs_(neggi) > 0.) {
02170                     bnded = FALSE_;
02171                 }
02172             }
02173         }
02174 /* L50: */
02175     }
02176 /*     The indices of the nonzero components of d are now stored */
02177 /*       in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n). */
02178 /*       The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin. */
02179     if (*theta != 1.) {
02180 /*                   complete the initialization of p for theta not= one. */
02181         dscal_(col, theta, &p[*col + 1], &c__1);
02182     }
02183 /*     Initialize GCP xcp = x. */
02184     dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
02185     if (nbreak == 0 && nfree == *n + 1) {
02186 /*                  is a zero vector, return with the initial xcp as GCP. */
02187         if (*iprint > 100) {
02188 /*          s_wsfe(&io___105);
02189             i__1 = *n;
02190             for (i__ = 1; i__ <= i__1; ++i__) {
02191                 do_fio(&c__1, (char *)&xcp[i__], (long int)sizeof(double));
02192             }
02193             e_wsfe();*/
02194         }
02195         return 0;
02196     }
02197 /*     Initialize c = W'(xcp - x) = 0. */
02198     i__1 = col2;
02199     for (j = 1; j <= i__1; ++j) {
02200         c__[j] = 0.;
02201 /* L60: */
02202     }
02203 /*     Initialize derivative f2. */
02204     f2 = -(*theta) * f1;
02205     if (*col > 0) {
02206         bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &p[1], &v[1], info);
02207         if (*info != 0) {
02208             return 0;
02209         }
02210         f2 -= ddot_(&col2, &v[1], &c__1, &p[1], &c__1);
02211     }
02212     dtm = -f1 / f2;
02213     tsum = 0.;
02214     *nint = 1;
02215     if (*iprint >= 99) {
02216 /*      s_wsle(&io___109);
02217         do_lio(&c__9, &c__1, "There are ", (long int)10);
02218         do_lio(&c__3, &c__1, (char *)&nbreak, (long int)sizeof(long int));
02219         do_lio(&c__9, &c__1, "  breakpoints ", (long int)14);
02220         e_wsle();*/
02221     }
02222 /*     If there are no breakpoints, locate the GCP and return. */
02223     if (nbreak == 0) {
02224         goto L888;
02225     }
02226     nleft = nbreak;
02227     iter = 1;
02228     tj = 0.;
02229 /* ------------------- the beginning of the loop ------------------------- */
02230 L777:
02231 /*     Find the next smallest breakpoint; */
02232 /*       compute dt = t(nleft) - t(nleft + 1). */
02233     tj0 = tj;
02234     if (iter == 1) {
02235 /*         Since we already have the smallest breakpoint we need not do */
02236 /*         heapsort yet. Often only one breakpoint is used and the */
02237 /*         cost of heapsort is avoided. */
02238         tj = bkmin;
02239         ibp = iorder[ibkmin];
02240     } else {
02241         if (iter == 2) {
02242 /*             Replace the already used smallest breakpoint with the */
02243 /*             breakpoint numbered nbreak > nlast, before heapsort call. */
02244             if (ibkmin != nbreak) {
02245                 t[ibkmin] = t[nbreak];
02246                 iorder[ibkmin] = iorder[nbreak];
02247             }
02248 /*        Update heap structure of breakpoints */
02249 /*           (if iter=2, initialize heap). */
02250         }
02251         i__1 = iter - 2;
02252         hpsolb_(&nleft, &t[1], &iorder[1], &i__1);
02253         tj = t[nleft];
02254         ibp = iorder[nleft];
02255     }
02256     dt = tj - tj0;
02257     if (dt != 0. && *iprint >= 100) {
02258 /*      s_wsfe(&io___116);
02259         do_fio(&c__1, (char *)&(*nint), (long int)sizeof(long int));
02260         do_fio(&c__1, (char *)&f1, (long int)sizeof(double));
02261         do_fio(&c__1, (char *)&f2, (long int)sizeof(double));
02262         e_wsfe();
02263         s_wsfe(&io___117);
02264         do_fio(&c__1, (char *)&dt, (long int)sizeof(double));
02265         e_wsfe();
02266         s_wsfe(&io___118);
02267         do_fio(&c__1, (char *)&dtm, (long int)sizeof(double));
02268         e_wsfe();*/
02269     }
02270 /*     If a minimizer is within this interval, locate the GCP and return. */
02271     if (dtm < dt) {
02272         goto L888;
02273     }
02274 /*     Otherwise fix one variable and */
02275 /*       reset the corresponding component of d to zero. */
02276     tsum += dt;
02277     --nleft;
02278     ++iter;
02279     dibp = d__[ibp];
02280     d__[ibp] = 0.;
02281     if (dibp > 0.) {
02282         zibp = u[ibp] - x[ibp];
02283         xcp[ibp] = u[ibp];
02284         iwhere[ibp] = 2;
02285     } else {
02286         zibp = l[ibp] - x[ibp];
02287         xcp[ibp] = l[ibp];
02288         iwhere[ibp] = 1;
02289     }
02290     if (*iprint >= 100) {
02291 /*      s_wsle(&io___121);
02292         do_lio(&c__9, &c__1, "Variable  ", (long int)10);
02293         do_lio(&c__3, &c__1, (char *)&ibp, (long int)sizeof(long int));
02294         do_lio(&c__9, &c__1, "  is fixed.", (long int)11);
02295         e_wsle();*/
02296     }
02297     if (nleft == 0 && nbreak == *n) {
02298 /*                                             all n variables are fixed, */
02299 /*                                                return with xcp as GCP. */
02300         dtm = dt;
02301         goto L999;
02302     }
02303 /*     Update the derivative information. */
02304     ++(*nint);
02305 /* Computing 2nd power */
02306     d__1 = dibp;
02307     dibp2 = d__1 * d__1;
02308 /*     Update f1 and f2. */
02309 /*        temporarily set f1 and f2 for col=0. */
02310     f1 = f1 + dt * f2 + dibp2 - *theta * dibp * zibp;
02311     f2 -= *theta * dibp2;
02312     if (*col > 0) {
02313 /*                          update c = c + dt*p. */
02314         daxpy_(&col2, &dt, &p[1], &c__1, &c__[1], &c__1);
02315 /*           choose wbp, */
02316 /*           the row of W corresponding to the breakpoint encountered. */
02317         pointr = *head;
02318         i__1 = *col;
02319         for (j = 1; j <= i__1; ++j) {
02320             wbp[j] = wy[ibp + pointr * wy_dim1];
02321             wbp[*col + j] = *theta * ws[ibp + pointr * ws_dim1];
02322             pointr = pointr % *m + 1;
02323 /* L70: */
02324         }
02325 /*           compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'. */
02326         bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wbp[1], &v[1], info);
02327         if (*info != 0) {
02328             return 0;
02329         }
02330         wmc = ddot_(&col2, &c__[1], &c__1, &v[1], &c__1);
02331         wmp = ddot_(&col2, &p[1], &c__1, &v[1], &c__1);
02332         wmw = ddot_(&col2, &wbp[1], &c__1, &v[1], &c__1);
02333 /*           update p = p - dibp*wbp. */
02334         d__1 = -dibp;
02335         daxpy_(&col2, &d__1, &wbp[1], &c__1, &p[1], &c__1);
02336 /*           complete updating f1 and f2 while col > 0. */
02337         f1 += dibp * wmc;
02338         f2 = f2 + dibp * 2. * wmp - dibp2 * wmw;
02339     }
02340     if (nleft > 0) {
02341         dtm = -f1 / f2;
02342         goto L777;
02343 /*                 to repeat the loop for unsearched intervals. */
02344     } else if (bnded) {
02345         f1 = 0.;
02346         f2 = 0.;
02347         dtm = 0.;
02348     } else {
02349         dtm = -f1 / f2;
02350     }
02351 /* ------------------- the end of the loop ------------------------------- */
02352 L888:
02353     if (*iprint >= 99) {
02354 /*      s_wsle(&io___126);
02355         e_wsle();
02356         s_wsle(&io___127);
02357         do_lio(&c__9, &c__1, "GCP found in this segment", (long int)25);
02358         e_wsle();
02359         s_wsfe(&io___128);
02360         do_fio(&c__1, (char *)&(*nint), (long int)sizeof(long int));
02361         do_fio(&c__1, (char *)&f1, (long int)sizeof(double));
02362         do_fio(&c__1, (char *)&f2, (long int)sizeof(double));
02363         e_wsfe();
02364         s_wsfe(&io___129);
02365         do_fio(&c__1, (char *)&dtm, (long int)sizeof(double));
02366         e_wsfe();*/
02367     }
02368     if (dtm <= 0.) {
02369         dtm = 0.;
02370     }
02371     tsum += dtm;
02372 /*     Move free variables (i.e., the ones w/o breakpoints) and */
02373 /*       the variables whose breakpoints haven't been reached. */
02374     daxpy_(n, &tsum, &d__[1], &c__1, &xcp[1], &c__1);
02375 L999:
02376 /*     Update c = c + dtm*p = W'(x^c - x) */
02377 /*       which will be used in computing r = Z'(B(x^c - x) + g). */
02378     if (*col > 0) {
02379         daxpy_(&col2, &dtm, &p[1], &c__1, &c__[1], &c__1);
02380     }
02381     if (*iprint > 100) {
02382 /*      s_wsfe(&io___130);
02383         i__1 = *n;
02384         for (i__ = 1; i__ <= i__1; ++i__) {
02385             do_fio(&c__1, (char *)&xcp[i__], (long int)sizeof(double));
02386         }
02387         e_wsfe();*/
02388     }
02389     if (*iprint >= 99) {
02390 /*      s_wsfe(&io___131);
02391         e_wsfe();*/
02392     }
02393     return 0;
02394 } /* cauchy_ */
02395 
02396 /* ====================== The end of cauchy ============================== */
02397 /* Subroutine */ int cmprlb_(long int *n, long int *m, double *x, double *g, double *ws, double *wy, double *sy,
02398                          double *wt, double *z__, double *r__, double *wa, long int *index, double *theta,
02399                          long int *col, long int *head, long int *nfree, long int *cnstnd, long int *info)
02400 /*long int *n, *m;
02401 double *x, *g, *ws, *wy, *sy, *wt, *z__, *r__, *wa;
02402 long int *index;
02403 double *theta;
02404 long int *col, *head, *nfree;
02405 long int *cnstnd;
02406 long int *info;*/
02407 {
02408     /* System generated locals */
02409     long int ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, 
02410             wt_dim1, wt_offset, i__1, i__2;
02411 
02412     /* Local variables */
02413     static long int i__, j, k;
02414     static double a1, a2;
02415     static long int pointr;
02416     extern /* Subroutine */ int bmv_(long int *m, double *sy, double *wt, long int *col, double *v, double *p, long int *info);
02417 
02418 /*     ************ */
02419 
02420 /*     Subroutine cmprlb */
02421 
02422 /*       This subroutine computes r=-Z'B(xcp-xk)-Z'g by using */
02423 /*         wa(2m+1)=W'(xcp-x) from subroutine cauchy. */
02424 
02425 /*     Subprograms called: */
02426 
02427 /*       L-BFGS-B Library ... bmv. */
02428 
02429 
02430 /*                           *  *  * */
02431 
02432 /*     NEOS, November 1994. (Latest revision June 1996.) */
02433 /*     Optimization Technology Center. */
02434 /*     Argonne National Laboratory and Northwestern University. */
02435 /*     Written by */
02436 /*                        Ciyou Zhu */
02437 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
02438 
02439 
02440 /*     ************ */
02441     /* Parameter adjustments */
02442     --index;
02443     --r__;
02444     --z__;
02445     --g;
02446     --x;
02447     --wa;
02448     wt_dim1 = *m;
02449     wt_offset = 1 + wt_dim1 * 1;
02450     wt -= wt_offset;
02451     sy_dim1 = *m;
02452     sy_offset = 1 + sy_dim1 * 1;
02453     sy -= sy_offset;
02454     wy_dim1 = *n;
02455     wy_offset = 1 + wy_dim1 * 1;
02456     wy -= wy_offset;
02457     ws_dim1 = *n;
02458     ws_offset = 1 + ws_dim1 * 1;
02459     ws -= ws_offset;
02460 
02461     /* Function Body */
02462     if (! (*cnstnd) && *col > 0) {
02463         i__1 = *n;
02464         for (i__ = 1; i__ <= i__1; ++i__) {
02465             r__[i__] = -g[i__];
02466 /* L26: */
02467         }
02468     } else {
02469         i__1 = *nfree;
02470         for (i__ = 1; i__ <= i__1; ++i__) {
02471             k = index[i__];
02472             r__[i__] = -(*theta) * (z__[k] - x[k]) - g[k];
02473 /* L30: */
02474         }
02475         bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wa[(*m << 1) + 1], &wa[
02476                 1], info);
02477         if (*info != 0) {
02478             *info = -8;
02479             return 0;
02480         }
02481         pointr = *head;
02482         i__1 = *col;
02483         for (j = 1; j <= i__1; ++j) {
02484             a1 = wa[j];
02485             a2 = *theta * wa[*col + j];
02486             i__2 = *nfree;
02487             for (i__ = 1; i__ <= i__2; ++i__) {
02488                 k = index[i__];
02489                 r__[i__] = r__[i__] + wy[k + pointr * wy_dim1] * a1 + ws[k + 
02490                         pointr * ws_dim1] * a2;
02491 /* L32: */
02492             }
02493             pointr = pointr % *m + 1;
02494 /* L34: */
02495         }
02496     }
02497     return 0;
02498 } /* cmprlb_ */
02499 
02500 /* ======================= The end of cmprlb ============================= */
02501 /* Subroutine */ int errclb_(long int *n, long int *m, double *factr, double *l, double *u, long int *nbd, char *task,
02502                          long int *info, long int *k, long int task_len)
02503 /*long int *n, *m;
02504 double *factr, *l, *u;
02505 long int *nbd;
02506 char *task;
02507 long int *info, *k;
02508 long int task_len;*/
02509 {
02510     /* System generated locals */
02511     long int i__1;
02512 
02513     /* Builtin functions */
02514     /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
02515 
02516     /* Local variables */
02517     static long int i__;
02518 
02519 /*     ************ */
02520 
02521 /*     Subroutine errclb */
02522 
02523 /*     This subroutine checks the validity of the input data. */
02524 
02525 
02526 /*                           *  *  * */
02527 
02528 /*     NEOS, November 1994. (Latest revision June 1996.) */
02529 /*     Optimization Technology Center. */
02530 /*     Argonne National Laboratory and Northwestern University. */
02531 /*     Written by */
02532 /*                        Ciyou Zhu */
02533 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
02534 
02535 
02536 /*     ************ */
02537 /*     Check the input arguments for errors. */
02538     /* Parameter adjustments */
02539     --nbd;
02540     --u;
02541     --l;
02542 
02543     /* Function Body */
02544     if (*n <= 0) {
02545         s_copy(task, "ERROR: N .LE. 0", (long int)60, (long int)15);
02546     }
02547     if (*m <= 0) {
02548         s_copy(task, "ERROR: M .LE. 0", (long int)60, (long int)15);
02549     }
02550     if (*factr < 0.) {
02551         s_copy(task, "ERROR: FACTR .LT. 0", (long int)60, (long int)19);
02552     }
02553 /*     Check the validity of the arrays nbd(i), u(i), and l(i). */
02554     i__1 = *n;
02555     for (i__ = 1; i__ <= i__1; ++i__) {
02556         if (nbd[i__] < 0 || nbd[i__] > 3) {
02557 /*                                                   return */
02558             s_copy(task, "ERROR: INVALID NBD", (long int)60, (long int)18);
02559             *info = -6;
02560             *k = i__;
02561         }
02562         if (nbd[i__] == 2) {
02563             if (l[i__] > u[i__]) {
02564 /*                                    return */
02565                 s_copy(task, "ERROR: NO FEASIBLE SOLUTION", (long int)60, (
02566                         long int)27);
02567                 *info = -7;
02568                 *k = i__;
02569             }
02570         }
02571 /* L10: */
02572     }
02573     return 0;
02574 } /* errclb_ */
02575 
02576 /* ======================= The end of errclb ============================= */
02577 /* Subroutine */ int formk_(long int *n, long int *nsub, long int *ind, long int *nenter, long int *ileave, long int *indx2, long int *iupdat, 
02578                         long int *updatd, double *wn, double *wn1, long int *m, double *ws, double *wy, double *sy, 
02579                         double *theta, long int *col, long int *head, long int *info)
02580 /*long int *n, *nsub, *ind, *nenter, *ileave, *indx2, *iupdat;
02581 long int *updatd;
02582 double *wn, *wn1;
02583 long int *m;
02584 double *ws, *wy, *sy, *theta;
02585 long int *col, *head, *info;*/
02586 {
02587     /* System generated locals */
02588     long int wn_dim1, wn_offset, wn1_dim1, wn1_offset, ws_dim1, ws_offset, 
02589             wy_dim1, wy_offset, sy_dim1, sy_offset, i__1, i__2, i__3;
02590 
02591     /* Local variables */
02592     static long int dend, pend;
02593     extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
02594     static long int upcl;
02595     static double temp1, temp2, temp3, temp4;
02596     static long int i__, k;
02597     extern /* Subroutine */ int dpofa_(double *a, long int *lda, long int *n, long int *info), 
02598                         dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy),
02599                         dtrsl_(double *t, long int *ldt, long int *n, double *b, long int *job, long int *info);
02600     static long int ipntr, jpntr, k1, m2, dbegin, is, js, iy, jy, pbegin, is1, 
02601             js1, col2;
02602 
02603 /*     ************ */
02604 
02605 /*     Subroutine formk */
02606 
02607 /*     This subroutine forms  the LEL^T factorization of the indefinite */
02608 
02609 /*       matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ] */
02610 /*                     [L_a -R_z           theta*S'AA'S ] */
02611 /*                                                    where E = [-I  0] */
02612 /*                                                              [ 0  I] */
02613 /*     The matrix K can be shown to be equal to the matrix M^[-1]N */
02614 /*       occurring in section 5.1 of [1], as well as to the matrix */
02615 /*       Mbar^[-1] Nbar in section 5.3. */
02616 
02617 /*     n is an long int variable. */
02618 /*       On entry n is the dimension of the problem. */
02619 /*       On exit n is unchanged. */
02620 
02621 /*     nsub is an long int variable */
02622 /*       On entry nsub is the number of subspace variables in free set. */
02623 /*       On exit nsub is not changed. */
02624 
02625 /*     ind is an long int array of dimension nsub. */
02626 /*       On entry ind specifies the indices of subspace variables. */
02627 /*       On exit ind is unchanged. */
02628 
02629 /*     nenter is an long int variable. */
02630 /*       On entry nenter is the number of variables entering the */
02631 /*         free set. */
02632 /*       On exit nenter is unchanged. */
02633 
02634 /*     ileave is an long int variable. */
02635 /*       On entry indx2(ileave),...,indx2(n) are the variables leaving */
02636 /*         the free set. */
02637 /*       On exit ileave is unchanged. */
02638 
02639 /*     indx2 is an long int array of dimension n. */
02640 /*       On entry indx2(1),...,indx2(nenter) are the variables entering */
02641 /*         the free set, while indx2(ileave),...,indx2(n) are the */
02642 /*         variables leaving the free set. */
02643 /*       On exit indx2 is unchanged. */
02644 
02645 /*     iupdat is an long int variable. */
02646 /*       On entry iupdat is the total number of BFGS updates made so far. */
02647 /*       On exit iupdat is unchanged. */
02648 
02649 /*     updatd is a long int variable. */
02650 /*       On entry 'updatd' is true if the L-BFGS matrix is updatd. */
02651 /*       On exit 'updatd' is unchanged. */
02652 
02653 /*     wn is a double precision array of dimension 2m x 2m. */
02654 /*       On entry wn is unspecified. */
02655 /*       On exit the upper triangle of wn stores the LEL^T factorization */
02656 /*         of the 2*col x 2*col indefinite matrix */
02657 /*                     [-D -Y'ZZ'Y/theta     L_a'-R_z'  ] */
02658 /*                     [L_a -R_z           theta*S'AA'S ] */
02659 
02660 /*     wn1 is a double precision array of dimension 2m x 2m. */
02661 /*       On entry wn1 stores the lower triangular part of */
02662 /*                     [Y' ZZ'Y   L_a'+R_z'] */
02663 /*                     [L_a+R_z   S'AA'S   ] */
02664 /*         in the previous iteration. */
02665 /*       On exit wn1 stores the corresponding updated matrices. */
02666 /*       The purpose of wn1 is just to store these inner products */
02667 /*       so they can be easily updated and inserted into wn. */
02668 
02669 /*     m is an long int variable. */
02670 /*       On entry m is the maximum number of variable metric corrections */
02671 /*         used to define the limited memory matrix. */
02672 /*       On exit m is unchanged. */
02673 
02674 /*     ws, wy, sy, and wtyy are double precision arrays; */
02675 /*     theta is a double precision variable; */
02676 /*     col is an long int variable; */
02677 /*     head is an long int variable. */
02678 /*       On entry they store the information defining the */
02679 /*                                          limited memory BFGS matrix: */
02680 /*         ws(n,m) stores S, a set of s-vectors; */
02681 /*         wy(n,m) stores Y, a set of y-vectors; */
02682 /*         sy(m,m) stores S'Y; */
02683 /*         wtyy(m,m) stores the Cholesky factorization */
02684 /*                                   of (theta*S'S+LD^(-1)L') */
02685 /*         theta is the scaling factor specifying B_0 = theta I; */
02686 /*         col is the number of variable metric corrections stored; */
02687 /*         head is the location of the 1st s- (or y-) vector in S (or Y). */
02688 /*       On exit they are unchanged. */
02689 
02690 /*     info is an long int variable. */
02691 /*       On entry info is unspecified. */
02692 /*       On exit info =  0 for normal return; */
02693 /*                    = -1 when the 1st Cholesky factorization failed; */
02694 /*                    = -2 when the 2st Cholesky factorization failed. */
02695 
02696 /*     Subprograms called: */
02697 
02698 /*       Linpack ... dcopy, dpofa, dtrsl. */
02699 
02700 
02701 /*     References: */
02702 /*       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
02703 /*       memory algorithm for bound constrained optimization'', */
02704 /*       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
02705 
02706 /*       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a */
02707 /*       limited memory FORTRAN code for solving bound constrained */
02708 /*       optimization problems'', Tech. Report, NAM-11, EECS Department, */
02709 /*       Northwestern University, 1994. */
02710 
02711 /*       (Postscript files of these papers are available via anonymous */
02712 /*        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
02713 
02714 /*                           *  *  * */
02715 
02716 /*     NEOS, November 1994. (Latest revision June 1996.) */
02717 /*     Optimization Technology Center. */
02718 /*     Argonne National Laboratory and Northwestern University. */
02719 /*     Written by */
02720 /*                        Ciyou Zhu */
02721 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
02722 
02723 
02724 /*     ************ */
02725 /*     Form the lower triangular part of */
02726 /*               WN1 = [Y' ZZ'Y   L_a'+R_z'] */
02727 /*                     [L_a+R_z   S'AA'S   ] */
02728 /*        where L_a is the strictly lower triangular part of S'AA'Y */
02729 /*              R_z is the upper triangular part of S'ZZ'Y. */
02730     /* Parameter adjustments */
02731     --indx2;
02732     --ind;
02733     sy_dim1 = *m;
02734     sy_offset = 1 + sy_dim1 * 1;
02735     sy -= sy_offset;
02736     wy_dim1 = *n;
02737     wy_offset = 1 + wy_dim1 * 1;
02738     wy -= wy_offset;
02739     ws_dim1 = *n;
02740     ws_offset = 1 + ws_dim1 * 1;
02741     ws -= ws_offset;
02742     wn1_dim1 = 2 * *m;
02743     wn1_offset = 1 + wn1_dim1 * 1;
02744     wn1 -= wn1_offset;
02745     wn_dim1 = 2 * *m;
02746     wn_offset = 1 + wn_dim1 * 1;
02747     wn -= wn_offset;
02748 
02749     /* Function Body */
02750     if (*updatd) {
02751         if (*iupdat > *m) {
02752 /*                                 shift old part of WN1. */
02753             i__1 = *m - 1;
02754             for (jy = 1; jy <= i__1; ++jy) {
02755                 js = *m + jy;
02756                 i__2 = *m - jy;
02757                 dcopy_(&i__2, &wn1[jy + 1 + (jy + 1) * wn1_dim1], &c__1, &wn1[
02758                         jy + jy * wn1_dim1], &c__1);
02759                 i__2 = *m - jy;
02760                 dcopy_(&i__2, &wn1[js + 1 + (js + 1) * wn1_dim1], &c__1, &wn1[
02761                         js + js * wn1_dim1], &c__1);
02762                 i__2 = *m - 1;
02763                 dcopy_(&i__2, &wn1[*m + 2 + (jy + 1) * wn1_dim1], &c__1, &wn1[
02764                         *m + 1 + jy * wn1_dim1], &c__1);
02765 /* L10: */
02766             }
02767         }
02768 /*          put new rows in blocks (1,1), (2,1) and (2,2). */
02769         pbegin = 1;
02770         pend = *nsub;
02771         dbegin = *nsub + 1;
02772         dend = *n;
02773         iy = *col;
02774         is = *m + *col;
02775         ipntr = *head + *col - 1;
02776         if (ipntr > *m) {
02777             ipntr -= *m;
02778         }
02779         jpntr = *head;
02780         i__1 = *col;
02781         for (jy = 1; jy <= i__1; ++jy) {
02782             js = *m + jy;
02783             temp1 = 0.;
02784             temp2 = 0.;
02785             temp3 = 0.;
02786 /*             compute element jy of row 'col' of Y'ZZ'Y */
02787             i__2 = pend;
02788             for (k = pbegin; k <= i__2; ++k) {
02789                 k1 = ind[k];
02790                 temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
02791 /* L15: */
02792             }
02793 /*             compute elements jy of row 'col' of L_a and S'AA'S */
02794             i__2 = dend;
02795             for (k = dbegin; k <= i__2; ++k) {
02796                 k1 = ind[k];
02797                 temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
02798                 temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
02799 /* L16: */
02800             }
02801             wn1[iy + jy * wn1_dim1] = temp1;
02802             wn1[is + js * wn1_dim1] = temp2;
02803             wn1[is + jy * wn1_dim1] = temp3;
02804             jpntr = jpntr % *m + 1;
02805 /* L20: */
02806         }
02807 /*          put new column in block (2,1). */
02808         jy = *col;
02809         jpntr = *head + *col - 1;
02810         if (jpntr > *m) {
02811             jpntr -= *m;
02812         }
02813         ipntr = *head;
02814         i__1 = *col;
02815         for (i__ = 1; i__ <= i__1; ++i__) {
02816             is = *m + i__;
02817             temp3 = 0.;
02818 /*             compute element i of column 'col' of R_z */
02819             i__2 = pend;
02820             for (k = pbegin; k <= i__2; ++k) {
02821                 k1 = ind[k];
02822                 temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
02823 /* L25: */
02824             }
02825             ipntr = ipntr % *m + 1;
02826             wn1[is + jy * wn1_dim1] = temp3;
02827 /* L30: */
02828         }
02829         upcl = *col - 1;
02830     } else {
02831         upcl = *col;
02832     }
02833 /*       modify the old parts in blocks (1,1) and (2,2) due to changes */
02834 /*       in the set of free variables. */
02835     ipntr = *head;
02836     i__1 = upcl;
02837     for (iy = 1; iy <= i__1; ++iy) {
02838         is = *m + iy;
02839         jpntr = *head;
02840         i__2 = iy;
02841         for (jy = 1; jy <= i__2; ++jy) {
02842             js = *m + jy;
02843             temp1 = 0.;
02844             temp2 = 0.;
02845             temp3 = 0.;
02846             temp4 = 0.;
02847             i__3 = *nenter;
02848             for (k = 1; k <= i__3; ++k) {
02849                 k1 = indx2[k];
02850                 temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
02851                 temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
02852 /* L35: */
02853             }
02854             i__3 = *n;
02855             for (k = *ileave; k <= i__3; ++k) {
02856                 k1 = indx2[k];
02857                 temp3 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
02858                 temp4 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
02859 /* L36: */
02860             }
02861             wn1[iy + jy * wn1_dim1] = wn1[iy + jy * wn1_dim1] + temp1 - temp3;
02862             wn1[is + js * wn1_dim1] = wn1[is + js * wn1_dim1] - temp2 + temp4;
02863             jpntr = jpntr % *m + 1;
02864 /* L40: */
02865         }
02866         ipntr = ipntr % *m + 1;
02867 /* L45: */
02868     }
02869 /*       modify the old parts in block (2,1). */
02870     ipntr = *head;
02871     i__1 = *m + upcl;
02872     for (is = *m + 1; is <= i__1; ++is) {
02873         jpntr = *head;
02874         i__2 = upcl;
02875         for (jy = 1; jy <= i__2; ++jy) {
02876             temp1 = 0.;
02877             temp3 = 0.;
02878             i__3 = *nenter;
02879             for (k = 1; k <= i__3; ++k) {
02880                 k1 = indx2[k];
02881                 temp1 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
02882 /* L50: */
02883             }
02884             i__3 = *n;
02885             for (k = *ileave; k <= i__3; ++k) {
02886                 k1 = indx2[k];
02887                 temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
02888 /* L51: */
02889             }
02890             if (is <= jy + *m) {
02891                 wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] + temp1 - 
02892                         temp3;
02893             } else {
02894                 wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] - temp1 + 
02895                         temp3;
02896             }
02897             jpntr = jpntr % *m + 1;
02898 /* L55: */
02899         }
02900         ipntr = ipntr % *m + 1;
02901 /* L60: */
02902     }
02903 /*     Form the upper triangle of WN = [D+Y' ZZ'Y/theta   -L_a'+R_z' ] */
02904 /*                                     [-L_a +R_z        S'AA'S*theta] */
02905     m2 = *m << 1;
02906     i__1 = *col;
02907     for (iy = 1; iy <= i__1; ++iy) {
02908         is = *col + iy;
02909         is1 = *m + iy;
02910         i__2 = iy;
02911         for (jy = 1; jy <= i__2; ++jy) {
02912             js = *col + jy;
02913             js1 = *m + jy;
02914             wn[jy + iy * wn_dim1] = wn1[iy + jy * wn1_dim1] / *theta;
02915             wn[js + is * wn_dim1] = wn1[is1 + js1 * wn1_dim1] * *theta;
02916 /* L65: */
02917         }
02918         i__2 = iy - 1;
02919         for (jy = 1; jy <= i__2; ++jy) {
02920             wn[jy + is * wn_dim1] = -wn1[is1 + jy * wn1_dim1];
02921 /* L66: */
02922         }
02923         i__2 = *col;
02924         for (jy = iy; jy <= i__2; ++jy) {
02925             wn[jy + is * wn_dim1] = wn1[is1 + jy * wn1_dim1];
02926 /* L67: */
02927         }
02928         wn[iy + iy * wn_dim1] += sy[iy + iy * sy_dim1];
02929 /* L70: */
02930     }
02931 /*     Form the upper triangle of WN= [  LL'            L^-1(-L_a'+R_z')] */
02932 /*                                    [(-L_a +R_z)L'^-1   S'AA'S*theta  ] */
02933 /*        first Cholesky factor (1,1) block of wn to get LL' */
02934 /*                          with L' stored in the upper triangle of wn. */
02935     dpofa_(&wn[wn_offset], &m2, col, info);
02936     if (*info != 0) {
02937         *info = -1;
02938         return 0;
02939     }
02940 /*        then form L^-1(-L_a'+R_z') in the (1,2) block. */
02941     col2 = *col << 1;
02942     i__1 = col2;
02943     for (js = *col + 1; js <= i__1; ++js) {
02944         dtrsl_(&wn[wn_offset], &m2, col, &wn[js * wn_dim1 + 1], &c__11, info);
02945 /* L71: */
02946     }
02947 /*     Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the */
02948 /*        upper triangle of (2,2) block of wn. */
02949     i__1 = col2;
02950     for (is = *col + 1; is <= i__1; ++is) {
02951         i__2 = col2;
02952         for (js = is; js <= i__2; ++js) {
02953             wn[is + js * wn_dim1] += ddot_(col, &wn[is * wn_dim1 + 1], &c__1, 
02954                     &wn[js * wn_dim1 + 1], &c__1);
02955 /* L74: */
02956         }
02957 /* L72: */
02958     }
02959 /*     Cholesky factorization of (2,2) block of wn. */
02960     dpofa_(&wn[*col + 1 + (*col + 1) * wn_dim1], &m2, col, info);
02961     if (*info != 0) {
02962         *info = -2;
02963         return 0;
02964     }
02965     return 0;
02966 } /* formk_ */
02967 
02968 /* ======================= The end of formk ============================== */
02969 /* Subroutine */ int formt_(long int *m, double *wt, double *sy, double *ss, long int *col, double *theta, long int *info)
02970 /*long int *m;
02971 double *wt, *sy, *ss;
02972 long int *col;
02973 double *theta;
02974 long int *info;*/
02975 {
02976     /* System generated locals */
02977     long int wt_dim1, wt_offset, sy_dim1, sy_offset, ss_dim1, ss_offset, i__1, 
02978             i__2, i__3;
02979 
02980     /* Local variables */
02981     static double ddum;
02982     static long int i__, j, k;
02983     extern /* Subroutine */ int dpofa_(double *a, long int *lda, long int *n, long int *info);
02984     static long int k1;
02985 
02986 /*     ************ */
02987 
02988 /*     Subroutine formt */
02989 
02990 /*       This subroutine forms the upper half of the pos. def. and symm. */
02991 /*         T = theta*SS + L*D^(-1)*L', stores T in the upper triangle */
02992 /*         of the array wt, and performs the Cholesky factorization of T */
02993 /*         to produce J*J', with J' stored in the upper triangle of wt. */
02994 
02995 /*     Subprograms called: */
02996 
02997 /*       Linpack ... dpofa. */
02998 
02999 
03000 /*                           *  *  * */
03001 
03002 /*     NEOS, November 1994. (Latest revision June 1996.) */
03003 /*     Optimization Technology Center. */
03004 /*     Argonne National Laboratory and Northwestern University. */
03005 /*     Written by */
03006 /*                        Ciyou Zhu */
03007 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03008 
03009 
03010 /*     ************ */
03011 /*     Form the upper half of  T = theta*SS + L*D^(-1)*L', */
03012 /*        store T in the upper triangle of the array wt. */
03013     /* Parameter adjustments */
03014     ss_dim1 = *m;
03015     ss_offset = 1 + ss_dim1 * 1;
03016     ss -= ss_offset;
03017     sy_dim1 = *m;
03018     sy_offset = 1 + sy_dim1 * 1;
03019     sy -= sy_offset;
03020     wt_dim1 = *m;
03021     wt_offset = 1 + wt_dim1 * 1;
03022     wt -= wt_offset;
03023 
03024     /* Function Body */
03025     i__1 = *col;
03026     for (j = 1; j <= i__1; ++j) {
03027         wt[j * wt_dim1 + 1] = *theta * ss[j * ss_dim1 + 1];
03028 /* L52: */
03029     }
03030     i__1 = *col;
03031     for (i__ = 2; i__ <= i__1; ++i__) {
03032         i__2 = *col;
03033         for (j = i__; j <= i__2; ++j) {
03034             k1 = min_(i__,j) - 1;
03035             ddum = 0.;
03036             i__3 = k1;
03037             for (k = 1; k <= i__3; ++k) {
03038                 ddum += sy[i__ + k * sy_dim1] * sy[j + k * sy_dim1] / sy[k + 
03039                         k * sy_dim1];
03040 /* L53: */
03041             }
03042             wt[i__ + j * wt_dim1] = ddum + *theta * ss[i__ + j * ss_dim1];
03043 /* L54: */
03044         }
03045 /* L55: */
03046     }
03047 /*     Cholesky factorize T to J*J' with */
03048 /*        J' stored in the upper triangle of wt. */
03049     dpofa_(&wt[wt_offset], m, col, info);
03050     if (*info != 0) {
03051         *info = -3;
03052     }
03053     return 0;
03054 } /* formt_ */
03055 
03056 /* ======================= The end of formt ============================== */
03057 /* Subroutine */ int freev_(long int *n, long int *nfree, long int *index, long int *nenter, long int *ileave, long int *indx2, long int *iwhere, 
03058                         long int *wrk, long int *updatd, long int *cnstnd, long int *iprint, long int *iter)
03059 /*long int *n, *nfree, *index, *nenter, *ileave, *indx2, *iwhere;
03060 long int *wrk, *updatd, *cnstnd;
03061 long int *iprint, *iter;*/
03062 {
03063     /* System generated locals */
03064     long int i__1;
03065 
03066     /* Builtin functions */
03067 //    long int s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle();
03068 
03069     /* Local variables */
03070     static long int iact, i__, k;
03071 
03072     /* Fortran I/O blocks */
03073 /*    static cilist io___168 = { 0, 6, 0, 0, 0 };
03074     static cilist io___169 = { 0, 6, 0, 0, 0 };
03075     static cilist io___170 = { 0, 6, 0, 0, 0 };
03076     static cilist io___172 = { 0, 6, 0, 0, 0 };*/
03077 
03078 
03079 /*     ************ */
03080 
03081 /*     Subroutine freev */
03082 
03083 /*     This subroutine counts the entering and leaving variables when */
03084 /*       iter > 0, and finds the index set of free and active variables */
03085 /*       at the GCP. */
03086 
03087 /*     cnstnd is a long int variable indicating whether bounds are present */
03088 
03089 /*     index is an long int array of dimension n */
03090 /*       for i=1,...,nfree, index(i) are the indices of free variables */
03091 /*       for i=nfree+1,...,n, index(i) are the indices of bound variables */
03092 /*       On entry after the first iteration, index gives */
03093 /*         the free variables at the previous iteration. */
03094 /*       On exit it gives the free variables based on the determination */
03095 /*         in cauchy using the array iwhere. */
03096 
03097 /*     indx2 is an long int array of dimension n */
03098 /*       On entry indx2 is unspecified. */
03099 /*       On exit with iter>0, indx2 indicates which variables */
03100 /*          have changed status since the previous iteration. */
03101 /*       For i= 1,...,nenter, indx2(i) have changed from bound to free. */
03102 /*       For i= ileave+1,...,n, indx2(i) have changed from free to bound. */
03103 
03104 
03105 /*                           *  *  * */
03106 
03107 /*     NEOS, November 1994. (Latest revision June 1996.) */
03108 /*     Optimization Technology Center. */
03109 /*     Argonne National Laboratory and Northwestern University. */
03110 /*     Written by */
03111 /*                        Ciyou Zhu */
03112 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03113 
03114 
03115 /*     ************ */
03116     /* Parameter adjustments */
03117     --iwhere;
03118     --indx2;
03119     --index;
03120 
03121     /* Function Body */
03122     *nenter = 0;
03123     *ileave = *n + 1;
03124     if (*iter > 0 && *cnstnd) {
03125 /*                           count the entering and leaving variables. */
03126         i__1 = *nfree;
03127         for (i__ = 1; i__ <= i__1; ++i__) {
03128             k = index[i__];
03129             if (iwhere[k] > 0) {
03130                 --(*ileave);
03131                 indx2[*ileave] = k;
03132                 if (*iprint >= 100) {
03133 /*                  s_wsle(&io___168);
03134                     do_lio(&c__9, &c__1, "Variable ", (long int)9);
03135                     do_lio(&c__3, &c__1, (char *)&k, (long int)sizeof(long int));
03136                     do_lio(&c__9, &c__1, " leaves the set of free variables", 
03137                             (long int)33);
03138                     e_wsle();*/
03139                 }
03140             }
03141 /* L20: */
03142         }
03143         i__1 = *n;
03144         for (i__ = *nfree + 1; i__ <= i__1; ++i__) {
03145             k = index[i__];
03146             if (iwhere[k] <= 0) {
03147                 ++(*nenter);
03148                 indx2[*nenter] = k;
03149                 if (*iprint >= 100) {
03150 /*                  s_wsle(&io___169);
03151                     do_lio(&c__9, &c__1, "Variable ", (long int)9);
03152                     do_lio(&c__3, &c__1, (char *)&k, (long int)sizeof(long int));
03153                     do_lio(&c__9, &c__1, " enters the set of free variables", 
03154                             (long int)33);
03155                     e_wsle();*/
03156                 }
03157             }
03158 /* L22: */
03159         }
03160         if (*iprint >= 99) {
03161 /*          s_wsle(&io___170);
03162             i__1 = *n + 1 - *ileave;
03163             do_lio(&c__3, &c__1, (char *)&i__1, (long int)sizeof(long int));
03164             do_lio(&c__9, &c__1, " variables leave; ", (long int)18);
03165             do_lio(&c__3, &c__1, (char *)&(*nenter), (long int)sizeof(long int));
03166             do_lio(&c__9, &c__1, " variables enter", (long int)16);
03167             e_wsle();*/
03168         }
03169     }
03170     *wrk = *ileave < *n + 1 || *nenter > 0 || *updatd;
03171 /*     Find the index set of free and active variables at the GCP. */
03172     *nfree = 0;
03173     iact = *n + 1;
03174     i__1 = *n;
03175     for (i__ = 1; i__ <= i__1; ++i__) {
03176         if (iwhere[i__] <= 0) {
03177             ++(*nfree);
03178             index[*nfree] = i__;
03179         } else {
03180             --iact;
03181             index[iact] = i__;
03182         }
03183 /* L24: */
03184     }
03185     if (*iprint >= 99) {
03186 /*      s_wsle(&io___172);
03187         do_lio(&c__3, &c__1, (char *)&(*nfree), (long int)sizeof(long int));
03188         do_lio(&c__9, &c__1, " variables are free at GCP ", (long int)27);
03189         i__1 = *iter + 1;
03190         do_lio(&c__3, &c__1, (char *)&i__1, (long int)sizeof(long int));
03191         e_wsle();*/
03192     }
03193     return 0;
03194 } /* freev_ */
03195 
03196 /* ======================= The end of freev ============================== */
03197 /* Subroutine */ int hpsolb_(long int *n, double *t, long int *iorder, long int *iheap)
03198 /*long int *n;
03199 double *t;
03200 long int *iorder, *iheap;*/
03201 {
03202     /* System generated locals */
03203     long int i__1;
03204 
03205     /* Local variables */
03206     static double ddum;
03207     static long int i__, j, k, indxin, indxou;
03208     static double out;
03209 
03210 /*     ************ */
03211 
03212 /*     Subroutine hpsolb */
03213 
03214 /*     This subroutine sorts out the least element of t, and puts the */
03215 /*       remaining elements of t in a heap. */
03216 
03217 /*     n is an long int variable. */
03218 /*       On entry n is the dimension of the arrays t and iorder. */
03219 /*       On exit n is unchanged. */
03220 
03221 /*     t is a double precision array of dimension n. */
03222 /*       On entry t stores the elements to be sorted, */
03223 /*       On exit t(n) stores the least elements of t, and t(1) to t(n-1) */
03224 /*         stores the remaining elements in the form of a heap. */
03225 
03226 /*     iorder is an long int array of dimension n. */
03227 /*       On entry iorder(i) is the index of t(i). */
03228 /*       On exit iorder(i) is still the index of t(i), but iorder may be */
03229 /*         permuted in accordance with t. */
03230 
03231 /*     iheap is an long int variable specifying the task. */
03232 /*       On entry iheap should be set as follows: */
03233 /*         iheap .eq. 0 if t(1) to t(n) is not in the form of a heap, */
03234 /*         iheap .ne. 0 if otherwise. */
03235 /*       On exit iheap is unchanged. */
03236 
03237 
03238 /*     References: */
03239 /*       Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT. */
03240 
03241 /*                           *  *  * */
03242 
03243 /*     NEOS, November 1994. (Latest revision June 1996.) */
03244 /*     Optimization Technology Center. */
03245 /*     Argonne National Laboratory and Northwestern University. */
03246 /*     Written by */
03247 /*                        Ciyou Zhu */
03248 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03249 
03250 /*     ************ */
03251     /* Parameter adjustments */
03252     --iorder;
03253     --t;
03254 
03255     /* Function Body */
03256     if (*iheap == 0) {
03257 /*        Rearrange the elements t(1) to t(n) to form a heap. */
03258         i__1 = *n;
03259         for (k = 2; k <= i__1; ++k) {
03260             ddum = t[k];
03261             indxin = iorder[k];
03262 /*           Add ddum to the heap. */
03263             i__ = k;
03264 L10:
03265             if (i__ > 1) {
03266                 j = i__ / 2;
03267                 if (ddum < t[j]) {
03268                     t[i__] = t[j];
03269                     iorder[i__] = iorder[j];
03270                     i__ = j;
03271                     goto L10;
03272                 }
03273             }
03274             t[i__] = ddum;
03275             iorder[i__] = indxin;
03276 /* L20: */
03277         }
03278     }
03279 /*     Assign to 'out' the value of t(1), the least member of the heap, */
03280 /*        and rearrange the remaining members to form a heap as */
03281 /*        elements 1 to n-1 of t. */
03282     if (*n > 1) {
03283         i__ = 1;
03284         out = t[1];
03285         indxou = iorder[1];
03286         ddum = t[*n];
03287         indxin = iorder[*n];
03288 /*        Restore the heap */
03289 L30:
03290         j = i__ + i__;
03291         if (j <= *n - 1) {
03292             if (t[j + 1] < t[j]) {
03293                 ++j;
03294             }
03295             if (t[j] < ddum) {
03296                 t[i__] = t[j];
03297                 iorder[i__] = iorder[j];
03298                 i__ = j;
03299                 goto L30;
03300             }
03301         }
03302         t[i__] = ddum;
03303         iorder[i__] = indxin;
03304 /*     Put the least member in t(n). */
03305         t[*n] = out;
03306         iorder[*n] = indxou;
03307     }
03308     return 0;
03309 } /* hpsolb_ */
03310 
03311 /* ====================== The end of hpsolb ============================== */
03312 /* Subroutine */ int lnsrlb_(long int *n, double *l, double *u, long int *nbd, double *x, double *f, double *fold,
03313                          double *gd, double *gdold, double *g, double *d__, double *r__, double *t, 
03314                          double *z__, double *stp, double *dnorm, double *dtd, double *xstep, double *stpmx,
03315                          long int *iter, long int *ifun, long int *iback, long int *nfgv, long int *info, char *task, long int *boxed, 
03316                          long int *cnstnd, char *csave, long int *isave, double *dsave, long int task_len, long int csave_len)
03317 /*long int *n;
03318 double *l, *u;
03319 long int *nbd;
03320 double *x, *f, *fold, *gd, *gdold, *g, *d__, *r__, *t, *z__, *stp, *dnorm,
03321          *dtd, *xstep, *stpmx;
03322 long int *iter, *ifun, *iback, *nfgv, *info;
03323 char *task;
03324 long int *boxed, *cnstnd;
03325 char *csave;
03326 long int *isave;
03327 double *dsave;
03328 long int task_len;
03329 long int csave_len;*/
03330 {
03331     /* System generated locals */
03332     long int i__1;
03333     double d__1;
03334 
03335     /* Builtin functions */
03336     long int s_cmp(char *, const char *const, long int, long int);
03337     //double sqrt();
03338     /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
03339 
03340     /* Local variables */
03341     extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
03342     static long int i__;
03343     extern /* Subroutine */ int dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
03344     static double a1, a2;
03345     extern /* Subroutine */ int dcsrch_(double *f, double *g, double *stp, double *ftol, double *gtol, double *xtol,
03346                          double *stpmin, double *stpmax, char *task, long int *isave, double *dsave, long int task_len);
03347 
03348 /*     ********** */
03349 
03350 /*     Subroutine lnsrlb */
03351 
03352 /*     This subroutine calls subroutine dcsrch from the Minpack2 library */
03353 /*       to perform the line search.  Subroutine dscrch is safeguarded so */
03354 /*       that all trial points lie within the feasible region. */
03355 
03356 /*     Subprograms called: */
03357 
03358 /*       Minpack2 Library ... dcsrch. */
03359 
03360 /*       Linpack ... dtrsl, ddot. */
03361 
03362 
03363 /*                           *  *  * */
03364 
03365 /*     NEOS, November 1994. (Latest revision June 1996.) */
03366 /*     Optimization Technology Center. */
03367 /*     Argonne National Laboratory and Northwestern University. */
03368 /*     Written by */
03369 /*                        Ciyou Zhu */
03370 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03371 
03372 
03373 /*     ********** */
03374     /* Parameter adjustments */
03375     --z__;
03376     --t;
03377     --r__;
03378     --d__;
03379     --g;
03380     --x;
03381     --nbd;
03382     --u;
03383     --l;
03384     --isave;
03385     --dsave;
03386 
03387     /* Function Body */
03388     if (s_cmp(task, "FG_LN", (long int)5, (long int)5) == 0) {
03389         goto L556;
03390     }
03391     *dtd = ddot_(n, &d__[1], &c__1, &d__[1], &c__1);
03392     *dnorm = sqrt(*dtd);
03393 /*     Determine the maximum step length. */
03394     *stpmx = 1e10;
03395     if (*cnstnd) {
03396         if (*iter == 0) {
03397             *stpmx = 1.;
03398         } else {
03399             i__1 = *n;
03400             for (i__ = 1; i__ <= i__1; ++i__) {
03401                 a1 = d__[i__];
03402                 if (nbd[i__] != 0) {
03403                     if (a1 < 0. && nbd[i__] <= 2) {
03404                         a2 = l[i__] - x[i__];
03405                         if (a2 >= 0.) {
03406                             *stpmx = 0.;
03407                         } else if (a1 * *stpmx < a2) {
03408                             *stpmx = a2 / a1;
03409                         }
03410                     } else if (a1 > 0. && nbd[i__] >= 2) {
03411                         a2 = u[i__] - x[i__];
03412                         if (a2 <= 0.) {
03413                             *stpmx = 0.;
03414                         } else if (a1 * *stpmx > a2) {
03415                             *stpmx = a2 / a1;
03416                         }
03417                     }
03418                 }
03419 /* L43: */
03420             }
03421         }
03422     }
03423     if (*iter == 0 && ! (*boxed)) {
03424 /* Computing MIN */
03425         d__1 = 1. / *dnorm;
03426         *stp = min_(d__1,*stpmx);
03427     } else {
03428         *stp = 1.;
03429     }
03430     dcopy_(n, &x[1], &c__1, &t[1], &c__1);
03431     dcopy_(n, &g[1], &c__1, &r__[1], &c__1);
03432     *fold = *f;
03433     *ifun = 0;
03434     *iback = 0;
03435     s_copy(csave, "START", (long int)60, (long int)5);
03436 L556:
03437     *gd = ddot_(n, &g[1], &c__1, &d__[1], &c__1);
03438     if (*ifun == 0) {
03439         *gdold = *gd;
03440         if (*gd >= 0.) {
03441 /*                               the directional derivative >=0. */
03442 /*                               Line search is impossible. */
03443             *info = -4;
03444             return 0;
03445         }
03446     }
03447     dcsrch_(f, gd, stp, &c_b275, &c_b276, &c_b277, &c_b9, stpmx, csave, &
03448             isave[1], &dsave[1], (long int)60);
03449     *xstep = *stp * *dnorm;
03450     if (s_cmp(csave, "CONV", (long int)4, (long int)4) != 0 && s_cmp(csave, "WARN"
03451             , (long int)4, (long int)4) != 0) {
03452         s_copy(task, "FG_LNSRCH", (long int)60, (long int)9);
03453         ++(*ifun);
03454         ++(*nfgv);
03455         *iback = *ifun - 1;
03456         if (*stp == 1.) {
03457             dcopy_(n, &z__[1], &c__1, &x[1], &c__1);
03458         } else {
03459             i__1 = *n;
03460             for (i__ = 1; i__ <= i__1; ++i__) {
03461                 x[i__] = *stp * d__[i__] + t[i__];
03462 /* L41: */
03463             }
03464         }
03465     } else {
03466         s_copy(task, "NEW_X", (long int)60, (long int)5);
03467     }
03468     return 0;
03469 } /* lnsrlb_ */
03470 
03471 /* ======================= The end of lnsrlb ============================= */
03472 /* Subroutine */ int matupd_(long int *n, long int *m, double *ws, double *wy, double *sy, double *ss,
03473                          double *d__, double *r__, long int *itail, long int *iupdat, long int *col, long int *head,
03474                          double *theta, double *rr, double *dr, double *stp, double *dtd)
03475 /*long int *n, *m;
03476 double *ws, *wy, *sy, *ss, *d__, *r__;
03477 long int *itail, *iupdat, *col, *head;
03478 double *theta, *rr, *dr, *stp, *dtd;*/
03479 {
03480     /* System generated locals */
03481     long int ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, 
03482             ss_dim1, ss_offset, i__1, i__2;
03483 
03484     /* Local variables */
03485     extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
03486     static long int j;
03487     extern /* Subroutine */ int dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
03488     static long int pointr;
03489 
03490 /*     ************ */
03491 
03492 /*     Subroutine matupd */
03493 
03494 /*       This subroutine updates matrices WS and WY, and forms the */
03495 /*         middle matrix in B. */
03496 
03497 /*     Subprograms called: */
03498 
03499 /*       Linpack ... dcopy, ddot. */
03500 
03501 
03502 /*                           *  *  * */
03503 
03504 /*     NEOS, November 1994. (Latest revision June 1996.) */
03505 /*     Optimization Technology Center. */
03506 /*     Argonne National Laboratory and Northwestern University. */
03507 /*     Written by */
03508 /*                        Ciyou Zhu */
03509 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03510 
03511 
03512 /*     ************ */
03513 /*     Set pointers for matrices WS and WY. */
03514     /* Parameter adjustments */
03515     --r__;
03516     --d__;
03517     ss_dim1 = *m;
03518     ss_offset = 1 + ss_dim1 * 1;
03519     ss -= ss_offset;
03520     sy_dim1 = *m;
03521     sy_offset = 1 + sy_dim1 * 1;
03522     sy -= sy_offset;
03523     wy_dim1 = *n;
03524     wy_offset = 1 + wy_dim1 * 1;
03525     wy -= wy_offset;
03526     ws_dim1 = *n;
03527     ws_offset = 1 + ws_dim1 * 1;
03528     ws -= ws_offset;
03529 
03530     /* Function Body */
03531     if (*iupdat <= *m) {
03532         *col = *iupdat;
03533         *itail = (*head + *iupdat - 2) % *m + 1;
03534     } else {
03535         *itail = *itail % *m + 1;
03536         *head = *head % *m + 1;
03537     }
03538 /*     Update matrices WS and WY. */
03539     dcopy_(n, &d__[1], &c__1, &ws[*itail * ws_dim1 + 1], &c__1);
03540     dcopy_(n, &r__[1], &c__1, &wy[*itail * wy_dim1 + 1], &c__1);
03541 /*     Set theta=yy/ys. */
03542     *theta = *rr / *dr;
03543 /*     Form the middle matrix in B. */
03544 /*        update the upper triangle of SS, */
03545 /*                                         and the lower triangle of SY: */
03546     if (*iupdat > *m) {
03547 /*                              move old information */
03548         i__1 = *col - 1;
03549         for (j = 1; j <= i__1; ++j) {
03550             dcopy_(&j, &ss[(j + 1) * ss_dim1 + 2], &c__1, &ss[j * ss_dim1 + 1]
03551                     , &c__1);
03552             i__2 = *col - j;
03553             dcopy_(&i__2, &sy[j + 1 + (j + 1) * sy_dim1], &c__1, &sy[j + j * 
03554                     sy_dim1], &c__1);
03555 /* L50: */
03556         }
03557     }
03558 /*        add new information: the last row of SY */
03559 /*                                             and the last column of SS: */
03560     pointr = *head;
03561     i__1 = *col - 1;
03562     for (j = 1; j <= i__1; ++j) {
03563         sy[*col + j * sy_dim1] = ddot_(n, &d__[1], &c__1, &wy[pointr * 
03564                 wy_dim1 + 1], &c__1);
03565         ss[j + *col * ss_dim1] = ddot_(n, &ws[pointr * ws_dim1 + 1], &c__1, &
03566                 d__[1], &c__1);
03567         pointr = pointr % *m + 1;
03568 /* L51: */
03569     }
03570     if (*stp == 1.) {
03571         ss[*col + *col * ss_dim1] = *dtd;
03572     } else {
03573         ss[*col + *col * ss_dim1] = *stp * *stp * *dtd;
03574     }
03575     sy[*col + *col * sy_dim1] = *dr;
03576     return 0;
03577 } /* matupd_ */
03578 
03579 /* ======================= The end of matupd ============================= */
03580 /* Subroutine */ int prn1lb_(long int *n, long int *m, double *l, double *u, double *x, long int *iprint, long int *itfile, double *epsmch)
03581 /*long int *n, *m;
03582 double *l, *u, *x;
03583 long int *iprint, *itfile;
03584 double *epsmch;*/
03585 {
03586     /* Format strings *//*
03587     static char fmt_7001[] = "(\002RUNNING THE L-BFGS-B CODE\002,/,/,\002   \
03588         * * *\002,/,/,\002Machine precision =\002,1p,d10.3)";
03589     static char fmt_2001[] = "(\002RUNNING THE L-BFGS-B CODE\002,/,/,\002it \
03590    = iteration number\002,/,\002nf    = number of function evaluations\002,/,\
03591 \002nint  = number of segments explored during the Cauchy search\002,/,\002n\
03592 act  = number of active bounds at the generalized Cauchy point\002,/,\002sub\
03593    = manner in which the subspace minimization terminated:\002,/,\002       \
03594  con = converged, bnd = a bound was reached\002,/,\002itls  = number of iter\
03595 ations performed in the line search\002,/,\002stepl = step length used\002,/,\
03596 \002tstep = norm of the displacement (total step)\002,/,\002projg = norm of \
03597 the projected gradient\002,/,\002f     = function value\002,/,/,\002        \
03598    * * *\002,/,/,\002Machine precision =\002,1p,d10.3)";
03599     static char fmt_9001[] = "(/,3x,\002it\002,3x,\002nf\002,2x,\002nint\002\
03600 ,2x,\002nact\002,2x,\002sub\002,2x,\002itls\002,2x,\002stepl\002,4x,\002tstep\
03601 \002,5x,\002projg\002,8x,\002f\002)";
03602     static char fmt_1004[] = "(/,a4,1p,6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))";*/
03603 
03604     /* System generated locals */
03605 //    long int i__1;
03606 
03607     /* Builtin functions */
03608 //    long int s_wsfe(cilist *), do_fio(long int*, char*, long int), e_wsfe(), s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle();
03609 
03610     /* Local variables */
03611 //    static long int i__;
03612 
03613     /* Fortran I/O blocks */
03614 /*    static cilist io___185 = { 0, 6, 0, fmt_7001, 0 };
03615     static cilist io___186 = { 0, 6, 0, 0, 0 };
03616     static cilist io___187 = { 0, 0, 0, fmt_2001, 0 };
03617     static cilist io___188 = { 0, 0, 0, 0, 0 };
03618     static cilist io___189 = { 0, 0, 0, fmt_9001, 0 };
03619     static cilist io___190 = { 0, 6, 0, fmt_1004, 0 };
03620     static cilist io___192 = { 0, 6, 0, fmt_1004, 0 };
03621     static cilist io___193 = { 0, 6, 0, fmt_1004, 0 }; */
03622 
03623 
03624 /*     ************ */
03625 
03626 /*     Subroutine prn1lb */
03627 
03628 /*     This subroutine prints the input data, initial point, upper and */
03629 /*       lower bounds of each variable, machine precision, as well as */
03630 /*       the headings of the output. */
03631 
03632 
03633 /*                           *  *  * */
03634 
03635 /*     NEOS, November 1994. (Latest revision June 1996.) */
03636 /*     Optimization Technology Center. */
03637 /*     Argonne National Laboratory and Northwestern University. */
03638 /*     Written by */
03639 /*                        Ciyou Zhu */
03640 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03641 
03642 
03643 /*     ************ */
03644     /* Parameter adjustments */
03645     --x;
03646     --u;
03647     --l;
03648 
03649     /* Function Body */
03650     if (*iprint >= 0) {
03651 /*      s_wsfe(&io___185);
03652         do_fio(&c__1, (char *)&(*epsmch), (long int)sizeof(double));
03653         e_wsfe();
03654         s_wsle(&io___186);
03655         do_lio(&c__9, &c__1, "N = ", (long int)4);
03656         do_lio(&c__3, &c__1, (char *)&(*n), (long int)sizeof(long int));
03657         do_lio(&c__9, &c__1, "    M = ", (long int)8);
03658         do_lio(&c__3, &c__1, (char *)&(*m), (long int)sizeof(long int));
03659         e_wsle();
03660         if (*iprint >= 1) {
03661             io___187.ciunit = *itfile;
03662             s_wsfe(&io___187);
03663             do_fio(&c__1, (char *)&(*epsmch), (long int)sizeof(double));
03664             e_wsfe();
03665             io___188.ciunit = *itfile;
03666             s_wsle(&io___188);
03667             do_lio(&c__9, &c__1, "N = ", (long int)4);
03668             do_lio(&c__3, &c__1, (char *)&(*n), (long int)sizeof(long int));
03669             do_lio(&c__9, &c__1, "    M = ", (long int)8);
03670             do_lio(&c__3, &c__1, (char *)&(*m), (long int)sizeof(long int));
03671             e_wsle();
03672             io___189.ciunit = *itfile;
03673             s_wsfe(&io___189);
03674             e_wsfe();
03675             if (*iprint > 100) {
03676                 s_wsfe(&io___190);
03677                 do_fio(&c__1, "L =", (long int)3);
03678                 i__1 = *n;
03679                 for (i__ = 1; i__ <= i__1; ++i__) {
03680                     do_fio(&c__1, (char *)&l[i__], (long int)sizeof(double))
03681                             ;
03682                 }
03683                 e_wsfe();
03684                 s_wsfe(&io___192);
03685                 do_fio(&c__1, "X0 =", (long int)4);
03686                 i__1 = *n;
03687                 for (i__ = 1; i__ <= i__1; ++i__) {
03688                     do_fio(&c__1, (char *)&x[i__], (long int)sizeof(double))
03689                             ;
03690                 }
03691                 e_wsfe();
03692                 s_wsfe(&io___193);
03693                 do_fio(&c__1, "U =", (long int)3);
03694                 i__1 = *n;
03695                 for (i__ = 1; i__ <= i__1; ++i__) {
03696                     do_fio(&c__1, (char *)&u[i__], (long int)sizeof(double))
03697                             ;
03698                 }
03699                 e_wsfe();
03700             }
03701         }*/
03702     }
03703     return 0;
03704 } /* prn1lb_ */
03705 
03706 /* ======================= The end of prn1lb ============================= */
03707 /* Subroutine */ int prn2lb_(long int *n, double *x, double *f, double *g, long int *iprint, long int *itfile,
03708                          long int *iter, long int *nfgv, long int *nact, double *sbgnrm, long int *nint, char *word,
03709                          long int *iword, long int *iback, double *stp, double *xstep, long int word_len)
03710 /*long int *n;
03711 double *x, *f, *g;
03712 long int *iprint, *itfile, *iter, *nfgv, *nact;
03713 double *sbgnrm;
03714 long int *nint;
03715 char *word;
03716 long int *iword, *iback;
03717 double *stp, *xstep;
03718 long int word_len;*/
03719 {
03720     /* Format strings *//*
03721     static char fmt_2001[] = "(/,\002At iterate\002,i5,4x,\002f= \002,1p,d12\
03722 .5,4x,\002|proj g|= \002,1p,d12.5)";
03723     static char fmt_1004[] = "(/,a4,1p,6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))";
03724     static char fmt_3001[] = "(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),1\
03725 p,2(1x,d10.3))";*/
03726 
03727     /* System generated locals */
03728 //    long int i__1;
03729 
03730     /* Builtin functions */
03731     /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
03732  //   long int s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle(), s_wsfe(cilist *), do_fio(long int*, char*, long int), e_wsfe();
03733 
03734     /* Local variables */
03735 //    static long int imod, i__;
03736 
03737     /* Fortran I/O blocks */
03738 /*    static cilist io___194 = { 0, 6, 0, 0, 0 };
03739     static cilist io___195 = { 0, 6, 0, fmt_2001, 0 };
03740     static cilist io___196 = { 0, 6, 0, fmt_1004, 0 };
03741     static cilist io___198 = { 0, 6, 0, fmt_1004, 0 };
03742     static cilist io___200 = { 0, 6, 0, fmt_2001, 0 };
03743     static cilist io___201 = { 0, 0, 0, fmt_3001, 0 };*/
03744 
03745 
03746 /*     ************ */
03747 
03748 /*     Subroutine prn2lb */
03749 
03750 /*     This subroutine prints out new information after a successful */
03751 /*       line search. */
03752 
03753 
03754 /*                           *  *  * */
03755 
03756 /*     NEOS, November 1994. (Latest revision June 1996.) */
03757 /*     Optimization Technology Center. */
03758 /*     Argonne National Laboratory and Northwestern University. */
03759 /*     Written by */
03760 /*                        Ciyou Zhu */
03761 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03762 
03763 
03764 /*     ************ */
03765 /*           'word' records the status of subspace solutions. */
03766     /* Parameter adjustments */
03767     --g;
03768     --x;
03769 
03770     /* Function Body */
03771     if (*iword == 0) {
03772 /*                            the subspace minimization converged. */
03773         s_copy(word, "con", (long int)3, (long int)3);
03774     } else if (*iword == 1) {
03775 /*                          the subspace minimization stopped at a bound. */
03776         s_copy(word, "bnd", (long int)3, (long int)3);
03777     } else if (*iword == 5) {
03778 /*                             the truncated Newton step has been used. */
03779         s_copy(word, "TNT", (long int)3, (long int)3);
03780     } else {
03781         s_copy(word, "---", (long int)3, (long int)3);
03782     }
03783     if (*iprint >= 99) {
03784 /*      s_wsle(&io___194);
03785         do_lio(&c__9, &c__1, "LINE SEARCH", (long int)11);
03786         do_lio(&c__3, &c__1, (char *)&(*iback), (long int)sizeof(long int));
03787         do_lio(&c__9, &c__1, " times; norm of step = ", (long int)23);
03788         do_lio(&c__5, &c__1, (char *)&(*xstep), (long int)sizeof(double));
03789         e_wsle();
03790         s_wsfe(&io___195);
03791         do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
03792         do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
03793         do_fio(&c__1, (char *)&(*sbgnrm), (long int)sizeof(double));
03794         e_wsfe();
03795         if (*iprint > 100) {
03796             s_wsfe(&io___196);
03797             do_fio(&c__1, "X =", (long int)3);
03798             i__1 = *n;
03799             for (i__ = 1; i__ <= i__1; ++i__) {
03800                 do_fio(&c__1, (char *)&x[i__], (long int)sizeof(double));
03801             }
03802             e_wsfe();
03803             s_wsfe(&io___198);
03804             do_fio(&c__1, "G =", (long int)3);
03805             i__1 = *n;
03806             for (i__ = 1; i__ <= i__1; ++i__) {
03807                 do_fio(&c__1, (char *)&g[i__], (long int)sizeof(double));
03808             }
03809             e_wsfe();
03810         }*/
03811     } else if (*iprint > 0) {
03812 /*      imod = *iter % *iprint;
03813         if (imod == 0) {
03814             s_wsfe(&io___200);
03815             do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
03816             do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
03817             do_fio(&c__1, (char *)&(*sbgnrm), (long int)sizeof(double));
03818             e_wsfe();
03819         }*/
03820     }
03821     if (*iprint >= 1) {
03822 /*      io___201.ciunit = *itfile;
03823         s_wsfe(&io___201);
03824         do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
03825         do_fio(&c__1, (char *)&(*nfgv), (long int)sizeof(long int));
03826         do_fio(&c__1, (char *)&(*nint), (long int)sizeof(long int));
03827         do_fio(&c__1, (char *)&(*nact), (long int)sizeof(long int));
03828         do_fio(&c__1, word, (long int)3);
03829         do_fio(&c__1, (char *)&(*iback), (long int)sizeof(long int));
03830         do_fio(&c__1, (char *)&(*stp), (long int)sizeof(double));
03831         do_fio(&c__1, (char *)&(*xstep), (long int)sizeof(double));
03832         do_fio(&c__1, (char *)&(*sbgnrm), (long int)sizeof(double));
03833         do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
03834         e_wsfe();*/
03835     }
03836     return 0;
03837 } /* prn2lb_ */
03838 
03839 /* ======================= The end of prn2lb ============================= */
03840 /* Subroutine */ int prn3lb_(long int *n, double *x, double *f, char *task, long int *iprint, long int *info, 
03841                          long int *itfile, long int *iter, long int *nfgv, long int *nintol, long int *nskip, long int *nact,
03842                          double *sbgnrm, double *time, long int *nint, char *word, long int *iback, double *stp,
03843                          double *xstep, long int *k, double *cachyt, double *sbtime, double *lnscht,
03844                          long int task_len, long int word_len)
03845 /*long int *n;
03846 double *x, *f;
03847 char *task;
03848 long int *iprint, *info, *itfile, *iter, *nfgv, *nintol, *nskip, *nact;
03849 double *sbgnrm, *time;
03850 long int *nint;
03851 char *word;
03852 long int *iback;
03853 double *stp, *xstep;
03854 long int *k;
03855 double *cachyt, *sbtime, *lnscht;
03856 long int task_len;
03857 long int word_len;*/
03858 {
03859     /* Format strings *//*
03860     static char fmt_3003[] = "(/,\002           * * *\002,/,/,\002Tit   = to\
03861 tal number of iterations\002,/,\002Tnf   = total number of function evaluati\
03862 ons\002,/,\002Tnint = total number of segments explored during\002,\002 Cauc\
03863 hy searches\002,/,\002Skip  = number of BFGS updates skipped\002,/,\002Nact \
03864  = number of active bounds at final generalized\002,\002 Cauchy point\002,/\
03865 ,\002Projg = norm of the final projected gradient\002,/,\002F     = final fu\
03866 nction value\002,/,/,\002           * * *\002)";
03867     static char fmt_3004[] = "(/,3x,\002N\002,3x,\002Tit\002,2x,\002Tnf\002,\
03868 2x,\002Tnint\002,2x,\002Skip\002,2x,\002Nact\002,5x,\002Projg\002,8x,\002\
03869 F\002)";
03870     static char fmt_3005[] = "(i5,2(1x,i4),(1x,i6),(2x,i4),(1x,i5),1p,2(2x,d\
03871 10.3))";
03872     static char fmt_1004[] = "(/,a4,1p,6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))";
03873     static char fmt_3009[] = "(/,a60)";
03874     static char fmt_9011[] = "(/,\002 Matrix in 1st Cholesky factorization i\
03875 n formk is not Pos. Def.\002)";
03876     static char fmt_9012[] = "(/,\002 Matrix in 2st Cholesky factorization i\
03877 n formk is not Pos. Def.\002)";
03878     static char fmt_9013[] = "(/,\002 Matrix in the Cholesky factorization i\
03879 n formt is not Pos. Def.\002)";
03880     static char fmt_9014[] = "(/,\002 Derivative >= 0, backtracking line sea\
03881 rch impossible.\002,/,\002   Previous x, f and g restored.\002,/,\002 Possib\
03882 le causes: 1 error in function or gradient evaluation;\002,/,\002           \
03883        2 rounding errors dominate computation.\002)";
03884     static char fmt_9015[] = "(/,\002 Warning:  more than 10 function and gr\
03885 adient\002,/,\002   evaluations in the last line search.  Termination\002,/\
03886 ,\002   may possibly be caused by a bad search direction.\002)";
03887     static char fmt_9018[] = "(/,\002 The triangular system is singular.\002)"
03888             ;
03889     static char fmt_9019[] = "(/,\002 Line search cannot locate an adequate \
03890 point after 20 function\002,/,\002  and gradient evaluations.  Previous x, f\
03891  and g restored.\002,/,\002 Possible causes: 1 error in function or gradient\
03892  evaluation;\002,/,\002                  2 rounding error dominate computati\
03893 on.\002)";
03894     static char fmt_3007[] = "(/,\002 Cauchy                time\002,1p,e10.\
03895 3,\002 seconds.\002,/\002 Subspace minimization time\002,1p,e10.3,\002 secon\
03896 ds.\002,/\002 Line search           time\002,1p,e10.3,\002 seconds.\002)";
03897     static char fmt_3008[] = "(/,\002 Total User time\002,1p,e10.3,\002 seco\
03898 nds.\002,/)";
03899     static char fmt_3002[] = "(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),6\
03900 x,\002-\002,10x,\002-\002)";*/
03901 
03902     /* System generated locals */
03903 //    long int i__1;
03904 
03905     /* Builtin functions */
03906     long int s_cmp(char *, const char *const, long int, long int);
03907 //    long int s_wsfe(cilist *), e_wsfe(), do_fio(long int*, char*, long int), s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle();
03908 
03909     /* Local variables */
03910 //    static long int i__;
03911 
03912     /* Fortran I/O blocks */
03913 /*    static cilist io___202 = { 0, 6, 0, fmt_3003, 0 };
03914     static cilist io___203 = { 0, 6, 0, fmt_3004, 0 };
03915     static cilist io___204 = { 0, 6, 0, fmt_3005, 0 };
03916     static cilist io___205 = { 0, 6, 0, fmt_1004, 0 };
03917     static cilist io___207 = { 0, 6, 0, 0, 0 };
03918     static cilist io___208 = { 0, 6, 0, fmt_3009, 0 };
03919     static cilist io___209 = { 0, 6, 0, fmt_9011, 0 };
03920     static cilist io___210 = { 0, 6, 0, fmt_9012, 0 };
03921     static cilist io___211 = { 0, 6, 0, fmt_9013, 0 };
03922     static cilist io___212 = { 0, 6, 0, fmt_9014, 0 };
03923     static cilist io___213 = { 0, 6, 0, fmt_9015, 0 };
03924     static cilist io___214 = { 0, 6, 0, 0, 0 };
03925     static cilist io___215 = { 0, 6, 0, 0, 0 };
03926     static cilist io___216 = { 0, 6, 0, fmt_9018, 0 };
03927     static cilist io___217 = { 0, 6, 0, fmt_9019, 0 };
03928     static cilist io___218 = { 0, 6, 0, fmt_3007, 0 };
03929     static cilist io___219 = { 0, 6, 0, fmt_3008, 0 };
03930     static cilist io___220 = { 0, 0, 0, fmt_3002, 0 };
03931     static cilist io___221 = { 0, 0, 0, fmt_3009, 0 };
03932     static cilist io___222 = { 0, 0, 0, fmt_9011, 0 };
03933     static cilist io___223 = { 0, 0, 0, fmt_9012, 0 };
03934     static cilist io___224 = { 0, 0, 0, fmt_9013, 0 };
03935     static cilist io___225 = { 0, 0, 0, fmt_9014, 0 };
03936     static cilist io___226 = { 0, 0, 0, fmt_9015, 0 };
03937     static cilist io___227 = { 0, 0, 0, fmt_9018, 0 };
03938     static cilist io___228 = { 0, 0, 0, fmt_9019, 0 };
03939     static cilist io___229 = { 0, 0, 0, fmt_3008, 0 };*/
03940 
03941 
03942 /*     ************ */
03943 
03944 /*     Subroutine prn3lb */
03945 
03946 /*     This subroutine prints out information when either a built-in */
03947 /*       convergence test is satisfied or when an error message is */
03948 /*       generated. */
03949 
03950 
03951 /*                           *  *  * */
03952 
03953 /*     NEOS, November 1994. (Latest revision June 1996.) */
03954 /*     Optimization Technology Center. */
03955 /*     Argonne National Laboratory and Northwestern University. */
03956 /*     Written by */
03957 /*                        Ciyou Zhu */
03958 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03959 
03960 
03961 /*     ************ */
03962     /* Parameter adjustments */
03963     --x;
03964 
03965     /* Function Body */
03966     if (s_cmp(task, "ERROR", (long int)5, (long int)5) == 0) {
03967         goto L999;
03968     }
03969     if (*iprint >= 0) {
03970 /*      s_wsfe(&io___202);
03971         e_wsfe();
03972         s_wsfe(&io___203);
03973         e_wsfe();
03974         s_wsfe(&io___204);
03975         do_fio(&c__1, (char *)&(*n), (long int)sizeof(long int));
03976         do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
03977         do_fio(&c__1, (char *)&(*nfgv), (long int)sizeof(long int));
03978         do_fio(&c__1, (char *)&(*nintol), (long int)sizeof(long int));
03979         do_fio(&c__1, (char *)&(*nskip), (long int)sizeof(long int));
03980         do_fio(&c__1, (char *)&(*nact), (long int)sizeof(long int));
03981         do_fio(&c__1, (char *)&(*sbgnrm), (long int)sizeof(double));
03982         do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
03983         e_wsfe();
03984         if (*iprint >= 100) {
03985             s_wsfe(&io___205);
03986             do_fio(&c__1, "X =", (long int)3);
03987             i__1 = *n;
03988             for (i__ = 1; i__ <= i__1; ++i__) {
03989                 do_fio(&c__1, (char *)&x[i__], (long int)sizeof(double));
03990             }
03991             e_wsfe();
03992         }
03993         if (*iprint >= 1) {
03994             s_wsle(&io___207);
03995             do_lio(&c__9, &c__1, " F =", (long int)4);
03996             do_lio(&c__5, &c__1, (char *)&(*f), (long int)sizeof(double));
03997             e_wsle();
03998         }*/
03999     }
04000 L999:
04001     if (*iprint >= 0) {
04002 /*      s_wsfe(&io___208);
04003         do_fio(&c__1, task, (long int)60);
04004         e_wsfe();
04005         if (*info != 0) {
04006             if (*info == -1) {
04007                 s_wsfe(&io___209);
04008                 e_wsfe();
04009             }
04010             if (*info == -2) {
04011                 s_wsfe(&io___210);
04012                 e_wsfe();
04013             }
04014             if (*info == -3) {
04015                 s_wsfe(&io___211);
04016                 e_wsfe();
04017             }
04018             if (*info == -4) {
04019                 s_wsfe(&io___212);
04020                 e_wsfe();
04021             }
04022             if (*info == -5) {
04023                 s_wsfe(&io___213);
04024                 e_wsfe();
04025             }
04026             if (*info == -6) {
04027                 s_wsle(&io___214);
04028                 do_lio(&c__9, &c__1, " Input nbd(", (long int)11);
04029                 do_lio(&c__3, &c__1, (char *)&(*k), (long int)sizeof(long int));
04030                 do_lio(&c__9, &c__1, ") is invalid.", (long int)13);
04031                 e_wsle();
04032             }
04033             if (*info == -7) {
04034                 s_wsle(&io___215);
04035                 do_lio(&c__9, &c__1, " l(", (long int)3);
04036                 do_lio(&c__3, &c__1, (char *)&(*k), (long int)sizeof(long int));
04037                 do_lio(&c__9, &c__1, ") > u(", (long int)6);
04038                 do_lio(&c__3, &c__1, (char *)&(*k), (long int)sizeof(long int));
04039                 do_lio(&c__9, &c__1, ").  No feasible solution.", (long int)25);
04040                 e_wsle();
04041             }
04042             if (*info == -8) {
04043                 s_wsfe(&io___216);
04044                 e_wsfe();
04045             }
04046             if (*info == -9) {
04047                 s_wsfe(&io___217);
04048                 e_wsfe();
04049             }
04050         }
04051         if (*iprint >= 1) {
04052             s_wsfe(&io___218);
04053             do_fio(&c__1, (char *)&(*cachyt), (long int)sizeof(double));
04054             do_fio(&c__1, (char *)&(*sbtime), (long int)sizeof(double));
04055             do_fio(&c__1, (char *)&(*lnscht), (long int)sizeof(double));
04056             e_wsfe();
04057         }
04058         s_wsfe(&io___219);
04059         do_fio(&c__1, (char *)&(*time), (long int)sizeof(double));
04060         e_wsfe();
04061         if (*iprint >= 1) {
04062             if (*info == -4 || *info == -9) {
04063                 io___220.ciunit = *itfile;
04064                 s_wsfe(&io___220);
04065                 do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
04066                 do_fio(&c__1, (char *)&(*nfgv), (long int)sizeof(long int));
04067                 do_fio(&c__1, (char *)&(*nint), (long int)sizeof(long int));
04068                 do_fio(&c__1, (char *)&(*nact), (long int)sizeof(long int));
04069                 do_fio(&c__1, word, (long int)3);
04070                 do_fio(&c__1, (char *)&(*iback), (long int)sizeof(long int));
04071                 do_fio(&c__1, (char *)&(*stp), (long int)sizeof(double));
04072                 do_fio(&c__1, (char *)&(*xstep), (long int)sizeof(double));
04073                 e_wsfe();
04074             }
04075             io___221.ciunit = *itfile;
04076             s_wsfe(&io___221);
04077             do_fio(&c__1, task, (long int)60);
04078             e_wsfe();
04079             if (*info != 0) {
04080                 if (*info == -1) {
04081                     io___222.ciunit = *itfile;
04082                     s_wsfe(&io___222);
04083                     e_wsfe();
04084                 }
04085                 if (*info == -2) {
04086                     io___223.ciunit = *itfile;
04087                     s_wsfe(&io___223);
04088                     e_wsfe();
04089                 }
04090                 if (*info == -3) {
04091                     io___224.ciunit = *itfile;
04092                     s_wsfe(&io___224);
04093                     e_wsfe();
04094                 }
04095                 if (*info == -4) {
04096                     io___225.ciunit = *itfile;
04097                     s_wsfe(&io___225);
04098                     e_wsfe();
04099                 }
04100                 if (*info == -5) {
04101                     io___226.ciunit = *itfile;
04102                     s_wsfe(&io___226);
04103                     e_wsfe();
04104                 }
04105                 if (*info == -8) {
04106                     io___227.ciunit = *itfile;
04107                     s_wsfe(&io___227);
04108                     e_wsfe();
04109                 }
04110                 if (*info == -9) {
04111                     io___228.ciunit = *itfile;
04112                     s_wsfe(&io___228);
04113                     e_wsfe();
04114                 }
04115             }
04116             io___229.ciunit = *itfile;
04117             s_wsfe(&io___229);
04118             do_fio(&c__1, (char *)&(*time), (long int)sizeof(double));
04119             e_wsfe();
04120         }*/
04121     }
04122 /* L3006: */
04123     return 0;
04124 } /* prn3lb_ */
04125 
04126 /* ======================= The end of prn3lb ============================= */
04127 /* Subroutine */ int projgr_(long int *n, double *l, double *u, long int *nbd, double *x, double *g, double *sbgnrm)
04128 /*long int *n;
04129 double *l, *u;
04130 long int *nbd;
04131 double *x, *g, *sbgnrm;*/
04132 {
04133     /* System generated locals */
04134     long int i__1;
04135     double d__1, d__2;
04136 
04137     /* Local variables */
04138     static long int i__;
04139     static double gi;
04140 
04141 /*     ************ */
04142 
04143 /*     Subroutine projgr */
04144 
04145 /*     This subroutine computes the infinity norm of the projected */
04146 /*       gradient. */
04147 
04148 
04149 /*                           *  *  * */
04150 
04151 /*     NEOS, November 1994. (Latest revision June 1996.) */
04152 /*     Optimization Technology Center. */
04153 /*     Argonne National Laboratory and Northwestern University. */
04154 /*     Written by */
04155 /*                        Ciyou Zhu */
04156 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
04157 
04158 
04159 /*     ************ */
04160     /* Parameter adjustments */
04161     --g;
04162     --x;
04163     --nbd;
04164     --u;
04165     --l;
04166 
04167     /* Function Body */
04168     *sbgnrm = 0.;
04169     i__1 = *n;
04170     for (i__ = 1; i__ <= i__1; ++i__) {
04171         gi = g[i__];
04172         if (nbd[i__] != 0) {
04173             if (gi < 0.) {
04174                 if (nbd[i__] >= 2) {
04175 /* Computing MAX */
04176                     d__1 = x[i__] - u[i__];
04177                     gi = max_(d__1,gi);
04178                 }
04179             } else {
04180                 if (nbd[i__] <= 2) {
04181 /* Computing MIN */
04182                     d__1 = x[i__] - l[i__];
04183                     gi = min_(d__1,gi);
04184                 }
04185             }
04186         }
04187 /* Computing MAX */
04188         d__1 = *sbgnrm, d__2 = abs_(gi);
04189         *sbgnrm = max_(d__1,d__2);
04190 /* L15: */
04191     }
04192     return 0;
04193 } /* projgr_ */
04194 
04195 /* ======================= The end of projgr ============================= */
04196 /* Subroutine */ int subsm_(long int *n, long int *m, long int *nsub, long int *ind, double *l, double *u, long int *nbd,
04197                         double *x, double *d__, double *ws, double *wy, double *theta, long int *col,
04198                         long int *head, long int *iword, double *wv, double *wn, long int *iprint, long int *info)
04199 /*long int *n, *m, *nsub, *ind;
04200 double *l, *u;
04201 long int *nbd;
04202 double *x, *d__, *ws, *wy, *theta;
04203 long int *col, *head, *iword;
04204 double *wv, *wn;
04205 long int *iprint, *info;*/
04206 {
04207     /* Format strings *//*
04208     static char fmt_1001[] = "(/,\002----------------SUBSM entered----------\
04209 -------\002,/)";
04210     static char fmt_1002[] = "(\002ALPHA = \002,f7.5,\002 backtrack to the B\
04211 OX\002)";
04212     static char fmt_1003[] = "(\002Subspace solution X =  \002,/,(4x,1p,6(1x\
04213 ,d11.4)))";
04214     static char fmt_1004[] = "(/,\002----------------exit SUBSM ------------\
04215 --------\002,/)";*/
04216 
04217     /* System generated locals */
04218     long int ws_dim1, ws_offset, wy_dim1, wy_offset, wn_dim1, wn_offset, i__1, 
04219             i__2;
04220 
04221     /* Builtin functions */
04222 //    long int s_wsfe(cilist *), e_wsfe(), do_fio(long int*, char*, long int), s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle();
04223 
04224     /* Local variables */
04225     static double temp1, temp2;
04226     static long int i__, j, k;
04227     static double alpha;
04228     //extern /* Subroutine */ int dtrsl_();
04229     extern int dtrsl_(double *t, long int *ldt, long int *n, double *b, long int *job, long int *info);
04230     static long int m2;
04231     static double dk;
04232     static long int js, jy, pointr, ibd, col2;
04233 
04234     /* Fortran I/O blocks */
04235 /*    static cilist io___232 = { 0, 6, 0, fmt_1001, 0 };
04236     static cilist io___246 = { 0, 6, 0, fmt_1002, 0 };
04237     static cilist io___247 = { 0, 6, 0, 0, 0 };
04238     static cilist io___248 = { 0, 6, 0, fmt_1003, 0 };
04239     static cilist io___249 = { 0, 6, 0, fmt_1004, 0 }; */
04240 
04241 
04242 /*     ************ */
04243 
04244 /*     Subroutine subsm */
04245 
04246 /*     Given xcp, l, u, r, an index set that specifies */
04247 /*      the active set at xcp, and an l-BFGS matrix B */
04248 /*      (in terms of WY, WS, SY, WT, head, col, and theta), */
04249 /*      this subroutine computes an approximate solution */
04250 /*      of the subspace problem */
04251 
04252 /*      (P)   min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp) */
04253 
04254 /*             subject to l<=x<=u */
04255 /*                      x_i=xcp_i for all i in A(xcp) */
04256 
04257 /*      along the subspace unconstrained Newton direction */
04258 
04259 /*         d = -(Z'BZ)^(-1) r. */
04260 
04261 /*       The formula for the Newton direction, given the L-BFGS matrix */
04262 /*       and the Sherman-Morrison formula, is */
04263 
04264 /*         d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r. */
04265 
04266 /*       where */
04267 /*                 K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ] */
04268 /*                     [L_a -R_z           theta*S'AA'S ] */
04269 
04270 /*     Note that this procedure for computing d differs */
04271 /*     from that described in [1]. One can show that the matrix K is */
04272 /*     equal to the matrix M^[-1]N in that paper. */
04273 
04274 /*     n is an long int variable. */
04275 /*       On entry n is the dimension of the problem. */
04276 /*       On exit n is unchanged. */
04277 
04278 /*     m is an long int variable. */
04279 /*       On entry m is the maximum number of variable metric corrections */
04280 /*         used to define the limited memory matrix. */
04281 /*       On exit m is unchanged. */
04282 
04283 /*     nsub is an long int variable. */
04284 /*       On entry nsub is the number of free variables. */
04285 /*       On exit nsub is unchanged. */
04286 
04287 /*     ind is an long int array of dimension nsub. */
04288 /*       On entry ind specifies the coordinate indices of free variables. */
04289 /*       On exit ind is unchanged. */
04290 
04291 /*     l is a double precision array of dimension n. */
04292 /*       On entry l is the lower bound of x. */
04293 /*       On exit l is unchanged. */
04294 
04295 /*     u is a double precision array of dimension n. */
04296 /*       On entry u is the upper bound of x. */
04297 /*       On exit u is unchanged. */
04298 
04299 /*     nbd is a long int array of dimension n. */
04300 /*       On entry nbd represents the type of bounds imposed on the */
04301 /*         variables, and must be specified as follows: */
04302 /*         nbd(i)=0 if x(i) is unbounded, */
04303 /*                1 if x(i) has only a lower bound, */
04304 /*                2 if x(i) has both lower and upper bounds, and */
04305 /*                3 if x(i) has only an upper bound. */
04306 /*       On exit nbd is unchanged. */
04307 
04308 /*     x is a double precision array of dimension n. */
04309 /*       On entry x specifies the Cauchy point xcp. */
04310 /*       On exit x(i) is the minimizer of Q over the subspace of */
04311 /*                                                        free variables. */
04312 
04313 /*     d is a double precision array of dimension n. */
04314 /*       On entry d is the reduced gradient of Q at xcp. */
04315 /*       On exit d is the Newton direction of Q. */
04316 
04317 /*     ws and wy are double precision arrays; */
04318 /*     theta is a double precision variable; */
04319 /*     col is an long int variable; */
04320 /*     head is an long int variable. */
04321 /*       On entry they store the information defining the */
04322 /*                                          limited memory BFGS matrix: */
04323 /*         ws(n,m) stores S, a set of s-vectors; */
04324 /*         wy(n,m) stores Y, a set of y-vectors; */
04325 /*         theta is the scaling factor specifying B_0 = theta I; */
04326 /*         col is the number of variable metric corrections stored; */
04327 /*         head is the location of the 1st s- (or y-) vector in S (or Y). */
04328 /*       On exit they are unchanged. */
04329 
04330 /*     iword is an long int variable. */
04331 /*       On entry iword is unspecified. */
04332 /*       On exit iword specifies the status of the subspace solution. */
04333 /*         iword = 0 if the solution is in the box, */
04334 /*                 1 if some bound is encountered. */
04335 
04336 /*     wv is a double precision working array of dimension 2m. */
04337 
04338 /*     wn is a double precision array of dimension 2m x 2m. */
04339 /*       On entry the upper triangle of wn stores the LEL^T factorization */
04340 /*         of the indefinite matrix */
04341 
04342 /*              K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ] */
04343 /*                  [L_a -R_z           theta*S'AA'S ] */
04344 /*                                                    where E = [-I  0] */
04345 /*                                                              [ 0  I] */
04346 /*       On exit wn is unchanged. */
04347 
04348 /*     iprint is an long int variable that must be set by the user. */
04349 /*       It controls the frequency and type of output generated: */
04350 /*        iprint<0    no output is generated; */
04351 /*        iprint=0    print only one line at the last iteration; */
04352 /*        0<iprint<99 print also f and |proj g| every iprint iterations; */
04353 /*        iprint=99   print details of every iteration except n-vectors; */
04354 /*        iprint=100  print also the changes of active set and final x; */
04355 /*        iprint>100  print details of every iteration including x and g; */
04356 /*       When iprint > 0, the file iterate.dat will be created to */
04357 /*                        summarize the iteration. */
04358 
04359 /*     info is an long int variable. */
04360 /*       On entry info is unspecified. */
04361 /*       On exit info = 0       for normal return, */
04362 /*                    = nonzero for abnormal return */
04363 /*                                  when the matrix K is ill-conditioned. */
04364 
04365 /*     Subprograms called: */
04366 
04367 /*       Linpack dtrsl. */
04368 
04369 
04370 /*     References: */
04371 
04372 /*       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
04373 /*       memory algorithm for bound constrained optimization'', */
04374 /*       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
04375 
04376 
04377 
04378 /*                           *  *  * */
04379 
04380 /*     NEOS, November 1994. (Latest revision June 1996.) */
04381 /*     Optimization Technology Center. */
04382 /*     Argonne National Laboratory and Northwestern University. */
04383 /*     Written by */
04384 /*                        Ciyou Zhu */
04385 /*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
04386 
04387 
04388 /*     ************ */
04389     /* Parameter adjustments */
04390     --d__;
04391     --x;
04392     --nbd;
04393     --u;
04394     --l;
04395     wn_dim1 = 2 * *m;
04396     wn_offset = 1 + wn_dim1 * 1;
04397     wn -= wn_offset;
04398     --wv;
04399     wy_dim1 = *n;
04400     wy_offset = 1 + wy_dim1 * 1;
04401     wy -= wy_offset;
04402     ws_dim1 = *n;
04403     ws_offset = 1 + ws_dim1 * 1;
04404     ws -= ws_offset;
04405     --ind;
04406 
04407     /* Function Body */
04408     if (*nsub <= 0) {
04409         return 0;
04410     }
04411     if (*iprint >= 99) {
04412 /*      s_wsfe(&io___232);
04413         e_wsfe();*/
04414     }
04415 /*     Compute wv = W'Zd. */
04416     pointr = *head;
04417     i__1 = *col;
04418     for (i__ = 1; i__ <= i__1; ++i__) {
04419         temp1 = 0.;
04420         temp2 = 0.;
04421         i__2 = *nsub;
04422         for (j = 1; j <= i__2; ++j) {
04423             k = ind[j];
04424             temp1 += wy[k + pointr * wy_dim1] * d__[j];
04425             temp2 += ws[k + pointr * ws_dim1] * d__[j];
04426 /* L10: */
04427         }
04428         wv[i__] = temp1;
04429         wv[*col + i__] = *theta * temp2;
04430         pointr = pointr % *m + 1;
04431 /* L20: */
04432     }
04433 /*     Compute wv:=K^(-1)wv. */
04434     m2 = *m << 1;
04435     col2 = *col << 1;
04436     dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__11, info);
04437     if (*info != 0) {
04438         return 0;
04439     }
04440     i__1 = *col;
04441     for (i__ = 1; i__ <= i__1; ++i__) {
04442         wv[i__] = -wv[i__];
04443 /* L25: */
04444     }
04445     dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__1, info);
04446     if (*info != 0) {
04447         return 0;
04448     }
04449 /*     Compute d = (1/theta)d + (1/theta**2)Z'W wv. */
04450     pointr = *head;
04451     i__1 = *col;
04452     for (jy = 1; jy <= i__1; ++jy) {
04453         js = *col + jy;
04454         i__2 = *nsub;
04455         for (i__ = 1; i__ <= i__2; ++i__) {
04456             k = ind[i__];
04457             d__[i__] = d__[i__] + wy[k + pointr * wy_dim1] * wv[jy] / *theta 
04458                     + ws[k + pointr * ws_dim1] * wv[js];
04459 /* L30: */
04460         }
04461         pointr = pointr % *m + 1;
04462 /* L40: */
04463     }
04464     i__1 = *nsub;
04465     for (i__ = 1; i__ <= i__1; ++i__) {
04466         d__[i__] /= *theta;
04467 /* L50: */
04468     }
04469 /*     Backtrack to the feasible region. */
04470     alpha = 1.;
04471     temp1 = alpha;
04472     i__1 = *nsub;
04473     for (i__ = 1; i__ <= i__1; ++i__) {
04474         k = ind[i__];
04475         dk = d__[i__];
04476         if (nbd[k] != 0) {
04477             if (dk < 0. && nbd[k] <= 2) {
04478                 temp2 = l[k] - x[k];
04479                 if (temp2 >= 0.) {
04480                     temp1 = 0.;
04481                 } else if (dk * alpha < temp2) {
04482                     temp1 = temp2 / dk;
04483                 }
04484             } else if (dk > 0. && nbd[k] >= 2) {
04485                 temp2 = u[k] - x[k];
04486                 if (temp2 <= 0.) {
04487                     temp1 = 0.;
04488                 } else if (dk * alpha > temp2) {
04489                     temp1 = temp2 / dk;
04490                 }
04491             }
04492             if (temp1 < alpha) {
04493                 alpha = temp1;
04494                 ibd = i__;
04495             }
04496         }
04497 /* L60: */
04498     }
04499     if (alpha < 1.) {
04500         dk = d__[ibd];
04501         k = ind[ibd];
04502         if (dk > 0.) {
04503             x[k] = u[k];
04504             d__[ibd] = 0.;
04505         } else if (dk < 0.) {
04506             x[k] = l[k];
04507             d__[ibd] = 0.;
04508         }
04509     }
04510     i__1 = *nsub;
04511     for (i__ = 1; i__ <= i__1; ++i__) {
04512         k = ind[i__];
04513         x[k] += alpha * d__[i__];
04514 /* L70: */
04515     }
04516     if (*iprint >= 99) {
04517 /*      if (alpha < 1.) {
04518             s_wsfe(&io___246);
04519             do_fio(&c__1, (char *)&alpha, (long int)sizeof(double));
04520             e_wsfe();
04521         } else { 
04522             s_wsle(&io___247);
04523             do_lio(&c__9, &c__1, "SM solution inside the box", (long int)26);
04524             e_wsle();
04525         }
04526         if (*iprint > 100) {
04527             s_wsfe(&io___248);
04528             i__1 = *n;
04529             for (i__ = 1; i__ <= i__1; ++i__) {
04530                 do_fio(&c__1, (char *)&x[i__], (long int)sizeof(double));
04531             }
04532             e_wsfe();
04533         }*/
04534     }
04535     if (alpha < 1.) {
04536         *iword = 1;
04537     } else {
04538         *iword = 0;
04539     }
04540     if (*iprint >= 99) {
04541 /*      s_wsfe(&io___249);
04542         e_wsfe();*/
04543     }
04544     return 0;
04545 } /* subsm_ */
04546 
04547 /* ====================== The end of subsm =============================== */
04548 /* Subroutine */ int dcsrch_(double *f, double *g, double *stp, double *ftol, double *gtol, double *xtol,
04549                          double *stpmin, double *stpmax, char *task, long int *isave, double *dsave, long int task_len)
04550 /*double *f, *g, *stp, *ftol, *gtol, *xtol, *stpmin, *stpmax;
04551 char *task;
04552 long int *isave;
04553 double *dsave;
04554 long int task_len;*/
04555 {
04556     /* System generated locals */
04557     double d__1;
04558 
04559     /* Builtin functions */
04560     long int s_cmp(char *, const char *const, long int, long int);
04561     /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
04562 
04563     /* Local variables */
04564     static long int stage;
04565     static double finit, ginit, width, ftest, gtest, stmin, stmax, width1,
04566              fm, gm, fx, fy, gx, gy;
04567     static long int brackt;
04568     //extern /* Subroutine */ int dcstep_();
04569     extern int dcstep_(double *stx, double *fx, double *dx, double *sty, double *fy, double *dy,
04570                      double *stp, double *fp, double *dp, long int *brackt, double *stpmin, double *stpmax);
04571 
04572     static double fxm, fym, gxm, gym, stx, sty;
04573 
04574 /*     ********** */
04575 
04576 /*     Subroutine dcsrch */
04577 
04578 /*     This subroutine finds a step that satisfies a sufficient */
04579 /*     decrease condition and a curvature condition. */
04580 
04581 /*     Each call of the subroutine updates an interval with */
04582 /*     endpoints stx and sty. The interval is initially chosen */
04583 /*     so that it contains a minimizer of the modified function */
04584 
04585 /*           psi(stp) = f(stp) - f(0) - ftol*stp*f'(0). */
04586 
04587 /*     If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the */
04588 /*     interval is chosen so that it contains a minimizer of f. */
04589 
04590 /*     The algorithm is designed to find a step that satisfies */
04591 /*     the sufficient decrease condition */
04592 
04593 /*           f(stp) <= f(0) + ftol*stp*f'(0), */
04594 
04595 /*     and the curvature condition */
04596 
04597 /*           abs_(f'(stp)) <= gtol*abs_(f'(0)). */
04598 
04599 /*     If ftol is less than gtol and if, for example, the function */
04600 /*     is bounded below, then there is always a step which satisfies */
04601 /*     both conditions. */
04602 
04603 /*     If no step can be found that satisfies both conditions, then */
04604 /*     the algorithm stops with a warning. In this case stp only */
04605 /*     satisfies the sufficient decrease condition. */
04606 
04607 /*     A typical invocation of dcsrch has the following outline: */
04608 
04609 /*     task = 'START' */
04610 /*  10 continue */
04611 /*        call dcsrch( ... ) */
04612 /*        if (task .eq. 'FG') then */
04613 /*           Evaluate the function and the gradient at stp */
04614 /*           goto 10 */
04615 /*           end if */
04616 
04617 /*     NOTE: The user must no alter work arrays between calls. */
04618 
04619 /*     The subroutine statement is */
04620 
04621 /*        subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax, */
04622 /*                          task,isave,dsave) */
04623 /*     where */
04624 
04625 /*       f is a double precision variable. */
04626 /*         On initial entry f is the value of the function at 0. */
04627 /*            On subsequent entries f is the value of the */
04628 /*            function at stp. */
04629 /*         On exit f is the value of the function at stp. */
04630 
04631 /*      g is a double precision variable. */
04632 /*         On initial entry g is the derivative of the function at 0. */
04633 /*            On subsequent entries g is the derivative of the */
04634 /*            function at stp. */
04635 /*         On exit g is the derivative of the function at stp. */
04636 
04637 /*      stp is a double precision variable. */
04638 /*         On entry stp is the current estimate of a satisfactory */
04639 /*            step. On initial entry, a positive initial estimate */
04640 /*            must be provided. */
04641 /*         On exit stp is the current estimate of a satisfactory step */
04642 /*            if task = 'FG'. If task = 'CONV' then stp satisfies */
04643 /*            the sufficient decrease and curvature condition. */
04644 
04645 /*       ftol is a double precision variable. */
04646 /*         On entry ftol specifies a nonnegative tolerance for the */
04647 /*            sufficient decrease condition. */
04648 /*         On exit ftol is unchanged. */
04649 
04650 /*       gtol is a double precision variable. */
04651 /*         On entry gtol specifies a nonnegative tolerance for the */
04652 /*            curvature condition. */
04653 /*         On exit gtol is unchanged. */
04654 
04655 /*      xtol is a double precision variable. */
04656 /*         On entry xtol specifies a nonnegative relative tolerance */
04657 /*            for an acceptable step. The subroutine exits with a */
04658 /*            warning if the relative difference between sty and stx */
04659 /*            is less than xtol. */
04660 /*         On exit xtol is unchanged. */
04661 
04662 /*      stpmin is a double precision variable. */
04663 /*         On entry stpmin is a nonnegative lower bound for the step. */
04664 /*         On exit stpmin is unchanged. */
04665 
04666 /*      stpmax is a double precision variable. */
04667 /*         On entry stpmax is a nonnegative upper bound for the step. */
04668 /*         On exit stpmax is unchanged. */
04669 
04670 /*       task is a character variable of length at least 60. */
04671 /*         On initial entry task must be set to 'START'. */
04672 /*         On exit task indicates the required action: */
04673 
04674 /*            If task(1:2) = 'FG' then evaluate the function and */
04675 /*            derivative at stp and call dcsrch again. */
04676 
04677 /*            If task(1:4) = 'CONV' then the search is successful. */
04678 
04679 /*            If task(1:4) = 'WARN' then the subroutine is not able */
04680 /*            to satisfy the convergence conditions. The exit value of */
04681 /*            stp contains the best point found during the search. */
04682 
04683 /*            If task(1:5) = 'ERROR' then there is an error in the */
04684 /*            input arguments. */
04685 
04686 /*         On exit with convergence, a warning or an error, the */
04687 /*            variable task contains additional information. */
04688 
04689 /*       isave is an long int work array of dimension 2. */
04690 
04691 /*       dsave is a double precision work array of dimension 13. */
04692 
04693 /*     Subprograms called */
04694 
04695 /*      MINPACK-2 ... dcstep */
04696 
04697 /*     MINPACK-1 Project. June 1983. */
04698 /*     Argonne National Laboratory. */
04699 /*     Jorge J. More' and David J. Thuente. */
04700 
04701 /*     MINPACK-2 Project. October 1993. */
04702 /*     Argonne National Laboratory and University of Minnesota. */
04703 /*     Brett M. Averick, Richard G. Carter, and Jorge J. More'. */
04704 
04705 /*     ********** */
04706 /*     Initialization block. */
04707     /* Parameter adjustments */
04708     --dsave;
04709     --isave;
04710 
04711     /* Function Body */
04712     if (s_cmp(task, "START", (long int)5, (long int)5) == 0) {
04713 /*        Check the input arguments for errors. */
04714         if (*stp < *stpmin) {
04715             s_copy(task, "ERROR: STP .LT. STPMIN", task_len, (long int)22);
04716         }
04717         if (*stp > *stpmax) {
04718             s_copy(task, "ERROR: STP .GT. STPMAX", task_len, (long int)22);
04719         }
04720         if (*g >= 0.) {
04721             s_copy(task, "ERROR: INITIAL G .GE. ZERO", task_len, (long int)26);
04722         }
04723         if (*ftol < 0.) {
04724             s_copy(task, "ERROR: FTOL .LT. ZERO", task_len, (long int)21);
04725         }
04726         if (*gtol < 0.) {
04727             s_copy(task, "ERROR: GTOL .LT. ZERO", task_len, (long int)21);
04728         }
04729         if (*xtol < 0.) {
04730             s_copy(task, "ERROR: XTOL .LT. ZERO", task_len, (long int)21);
04731         }
04732         if (*stpmin < 0.) {
04733             s_copy(task, "ERROR: STPMIN .LT. ZERO", task_len, (long int)23);
04734         }
04735         if (*stpmax < *stpmin) {
04736             s_copy(task, "ERROR: STPMAX .LT. STPMIN", task_len, (long int)25);
04737         }
04738 /*        Exit if there are errors on input. */
04739         if (s_cmp(task, "ERROR", (long int)5, (long int)5) == 0) {
04740             return 0;
04741         }
04742 /*        Initialize local variables. */
04743         brackt = FALSE_;
04744         stage = 1;
04745         finit = *f;
04746         ginit = *g;
04747         gtest = *ftol * ginit;
04748         width = *stpmax - *stpmin;
04749         width1 = width / .5;
04750 /*        The variables stx, fx, gx contain the values of the step, */
04751 /*        function, and derivative at the best step. */
04752 /*        The variables sty, fy, gy contain the value of the step, */
04753 /*        function, and derivative at sty. */
04754 /*        The variables stp, f, g contain the values of the step, */
04755 /*        function, and derivative at stp. */
04756         stx = 0.;
04757         fx = finit;
04758         gx = ginit;
04759         sty = 0.;
04760         fy = finit;
04761         gy = ginit;
04762         stmin = 0.;
04763         stmax = *stp + *stp * 4.;
04764         s_copy(task, "FG", task_len, (long int)2);
04765         goto L1000;
04766     } else {
04767 /*        Restore local variables. */
04768         if (isave[1] == 1) {
04769             brackt = TRUE_;
04770         } else {
04771             brackt = FALSE_;
04772         }
04773         stage = isave[2];
04774         ginit = dsave[1];
04775         gtest = dsave[2];
04776         gx = dsave[3];
04777         gy = dsave[4];
04778         finit = dsave[5];
04779         fx = dsave[6];
04780         fy = dsave[7];
04781         stx = dsave[8];
04782         sty = dsave[9];
04783         stmin = dsave[10];
04784         stmax = dsave[11];
04785         width = dsave[12];
04786         width1 = dsave[13];
04787     }
04788 /*     If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the */
04789 /*     algorithm enters the second stage. */
04790     ftest = finit + *stp * gtest;
04791     if (stage == 1 && *f <= ftest && *g >= 0.) {
04792         stage = 2;
04793     }
04794 /*     Test for warnings. */
04795     if (brackt && (*stp <= stmin || *stp >= stmax)) {
04796         s_copy(task, "WARNING: ROUNDING ERRORS PREVENT PROGRESS", task_len, (
04797                 long int)41);
04798     }
04799     if (brackt && stmax - stmin <= *xtol * stmax) {
04800         s_copy(task, "WARNING: XTOL TEST SATISFIED", task_len, (long int)28);
04801     }
04802     if (*stp == *stpmax && *f <= ftest && *g <= gtest) {
04803         s_copy(task, "WARNING: STP = STPMAX", task_len, (long int)21);
04804     }
04805     if (*stp == *stpmin && (*f > ftest || *g >= gtest)) {
04806         s_copy(task, "WARNING: STP = STPMIN", task_len, (long int)21);
04807     }
04808 /*     Test for convergence. */
04809     if (*f <= ftest && abs_(*g) <= *gtol * (-ginit)) {
04810         s_copy(task, "CONVERGENCE", task_len, (long int)11);
04811     }
04812 /*     Test for termination. */
04813     if (s_cmp(task, "WARN", (long int)4, (long int)4) == 0 || s_cmp(task, "CONV", 
04814             (long int)4, (long int)4) == 0) {
04815         goto L1000;
04816     }
04817 /*     A modified function is used to predict the step during the */
04818 /*     first stage if a lower function value has been obtained but */
04819 /*     the decrease is not sufficient. */
04820     if (stage == 1 && *f <= fx && *f > ftest) {
04821 /*        Define the modified function and derivative values. */
04822         fm = *f - *stp * gtest;
04823         fxm = fx - stx * gtest;
04824         fym = fy - sty * gtest;
04825         gm = *g - gtest;
04826         gxm = gx - gtest;
04827         gym = gy - gtest;
04828 /*        Call dcstep to update stx, sty, and to compute the new step. */
04829         dcstep_(&stx, &fxm, &gxm, &sty, &fym, &gym, stp, &fm, &gm, &brackt, &
04830                 stmin, &stmax);
04831 /*        Reset the function and derivative values for f. */
04832         fx = fxm + stx * gtest;
04833         fy = fym + sty * gtest;
04834         gx = gxm + gtest;
04835         gy = gym + gtest;
04836     } else {
04837 /*       Call dcstep to update stx, sty, and to compute the new step. */
04838         dcstep_(&stx, &fx, &gx, &sty, &fy, &gy, stp, f, g, &brackt, &stmin, &
04839                 stmax);
04840     }
04841 /*     Decide if a bisection step is needed. */
04842     if (brackt) {
04843         if ((d__1 = sty - stx, abs_(d__1)) >= width1 * .66) {
04844             *stp = stx + (sty - stx) * .5;
04845         }
04846         width1 = width;
04847         width = (d__1 = sty - stx, abs_(d__1));
04848     }
04849 /*     Set the minimum and maximum steps allowed for stp. */
04850     if (brackt) {
04851         stmin = min_(stx,sty);
04852         stmax = max_(stx,sty);
04853     } else {
04854         stmin = *stp + (*stp - stx) * 1.1;
04855         stmax = *stp + (*stp - stx) * 4.;
04856     }
04857 /*     Force the step to be within the bounds stpmax and stpmin. */
04858     *stp = max_(*stp,*stpmin);
04859     *stp = min_(*stp,*stpmax);
04860 /*     If further progress is not possible, let stp be the best */
04861 /*     point obtained during the search. */
04862     if (brackt && (*stp <= stmin || *stp >= stmax) || brackt && stmax - stmin 
04863             <= *xtol * stmax) {
04864         *stp = stx;
04865     }
04866 /*     Obtain another function and derivative. */
04867     s_copy(task, "FG", task_len, (long int)2);
04868 L1000:
04869 /*     Save local variables. */
04870     if (brackt) {
04871         isave[1] = 1;
04872     } else {
04873         isave[1] = 0;
04874     }
04875     isave[2] = stage;
04876     dsave[1] = ginit;
04877     dsave[2] = gtest;
04878     dsave[3] = gx;
04879     dsave[4] = gy;
04880     dsave[5] = finit;
04881     dsave[6] = fx;
04882     dsave[7] = fy;
04883     dsave[8] = stx;
04884     dsave[9] = sty;
04885     dsave[10] = stmin;
04886     dsave[11] = stmax;
04887     dsave[12] = width;
04888     dsave[13] = width1;
04889     return 0;
04890 } /* dcsrch_ */
04891 
04892 /* ====================== The end of dcsrch ============================== */
04893 /* Subroutine */ int dcstep_(double *stx, double *fx, double *dx, double *sty, double *fy, double *dy,
04894                          double *stp, double *fp, double *dp, long int *brackt, double *stpmin, double *stpmax)
04895 /*double *stx, *fx, *dx, *sty, *fy, *dy, *stp, *fp, *dp;
04896 long int *brackt;
04897 double *stpmin, *stpmax;*/
04898 {
04899     /* System generated locals */
04900     double d__1, d__2, d__3;
04901 
04902     /* Builtin functions */
04903     //double sqrt();
04904 
04905     /* Local variables */
04906     static double sgnd, stpc, stpf, stpq, p, q, gamma, r__, s, theta;
04907 
04908 /*     ********** */
04909 
04910 /*     Subroutine dcstep */
04911 
04912 /*     This subroutine computes a safeguarded step for a search */
04913 /*     procedure and updates an interval that contains a step that */
04914 /*     satisfies a sufficient decrease and a curvature condition. */
04915 
04916 /*     The parameter stx contains the step with the least function */
04917 /*     value. If brackt is set to .true. then a minimizer has */
04918 /*     been bracketed in an interval with endpoints stx and sty. */
04919 /*     The parameter stp contains the current step. */
04920 /*     The subroutine assumes that if brackt is set to .true. then */
04921 
04922 /*           min_(stx,sty) < stp < max_(stx,sty), */
04923 
04924 /*     and that the derivative at stx is negative in the direction */
04925 /*     of the step. */
04926 
04927 /*     The subroutine statement is */
04928 
04929 /*       subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, */
04930 /*                         stpmin,stpmax) */
04931 
04932 /*     where */
04933 
04934 /*       stx is a double precision variable. */
04935 /*         On entry stx is the best step obtained so far and is an */
04936 /*            endpoint of the interval that contains the minimizer. */
04937 /*         On exit stx is the updated best step. */
04938 
04939 /*       fx is a double precision variable. */
04940 /*         On entry fx is the function at stx. */
04941 /*         On exit fx is the function at stx. */
04942 
04943 /*       dx is a double precision variable. */
04944 /*         On entry dx is the derivative of the function at */
04945 /*            stx. The derivative must be negative in the direction of */
04946 /*            the step, that is, dx and stp - stx must have opposite */
04947 /*            signs. */
04948 /*         On exit dx is the derivative of the function at stx. */
04949 
04950 /*       sty is a double precision variable. */
04951 /*         On entry sty is the second endpoint of the interval that */
04952 /*            contains the minimizer. */
04953 /*         On exit sty is the updated endpoint of the interval that */
04954 /*            contains the minimizer. */
04955 
04956 /*       fy is a double precision variable. */
04957 /*         On entry fy is the function at sty. */
04958 /*         On exit fy is the function at sty. */
04959 
04960 /*       dy is a double precision variable. */
04961 /*         On entry dy is the derivative of the function at sty. */
04962 /*         On exit dy is the derivative of the function at the exit sty. */
04963 
04964 /*       stp is a double precision variable. */
04965 /*         On entry stp is the current step. If brackt is set to .true. */
04966 /*            then on input stp must be between stx and sty. */
04967 /*         On exit stp is a new trial step. */
04968 
04969 /*       fp is a double precision variable. */
04970 /*         On entry fp is the function at stp */
04971 /*         On exit fp is unchanged. */
04972 
04973 /*       dp is a double precision variable. */
04974 /*         On entry dp is the the derivative of the function at stp. */
04975 /*         On exit dp is unchanged. */
04976 
04977 /*       brackt is an long int variable. */
04978 /*         On entry brackt specifies if a minimizer has been bracketed. */
04979 /*            Initially brackt must be set to .false. */
04980 /*         On exit brackt specifies if a minimizer has been bracketed. */
04981 /*            When a minimizer is bracketed brackt is set to .true. */
04982 
04983 /*       stpmin is a double precision variable. */
04984 /*         On entry stpmin is a lower bound for the step. */
04985 /*         On exit stpmin is unchanged. */
04986 
04987 /*       stpmax is a double precision variable. */
04988 /*         On entry stpmax is an upper bound for the step. */
04989 /*         On exit stpmax is unchanged. */
04990 
04991 /*     MINPACK-1 Project. June 1983 */
04992 /*     Argonne National Laboratory. */
04993 /*     Jorge J. More' and David J. Thuente. */
04994 
04995 /*     MINPACK-2 Project. October 1993. */
04996 /*     Argonne National Laboratory and University of Minnesota. */
04997 /*     Brett M. Averick and Jorge J. More'. */
04998 
04999 /*     ********** */
05000     sgnd = *dp * (*dx / abs_(*dx));
05001 /*     First case: A higher function value. The minimum is bracketed. */
05002 /*     If the cubic step is closer to stx than the quadratic step, the */
05003 /*     cubic step is taken, otherwise the average of the cubic and */
05004 /*     quadratic steps is taken. */
05005     if (*fp > *fx) {
05006         theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp;
05007 /* Computing MAX */
05008         d__1 = abs_(theta), d__2 = abs_(*dx), d__1 = max_(d__1,d__2), d__2 = abs_(
05009                 *dp);
05010         s = max_(d__1,d__2);
05011 /* Computing 2nd power */
05012         d__1 = theta / s;
05013         gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
05014         if (*stp < *stx) {
05015             gamma = -gamma;
05016         }
05017         p = gamma - *dx + theta;
05018         q = gamma - *dx + gamma + *dp;
05019         r__ = p / q;
05020         stpc = *stx + r__ * (*stp - *stx);
05021         stpq = *stx + *dx / ((*fx - *fp) / (*stp - *stx) + *dx) / 2. * (*stp 
05022                 - *stx);
05023         if ((d__1 = stpc - *stx, abs_(d__1)) < (d__2 = stpq - *stx, abs_(d__2)))
05024                  {
05025             stpf = stpc;
05026         } else {
05027             stpf = stpc + (stpq - stpc) / 2.;
05028         }
05029         *brackt = TRUE_;
05030 /*     Second case: A lower function value and derivatives of opposite */
05031 /*     sign. The minimum is bracketed. If the cubic step is farther from */
05032 /*     stp than the secant step, the cubic step is taken, otherwise the */
05033 /*     secant step is taken. */
05034     } else if (sgnd < 0.) {
05035         theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp;
05036 /* Computing MAX */
05037         d__1 = abs_(theta), d__2 = abs_(*dx), d__1 = max_(d__1,d__2), d__2 = abs_(
05038                 *dp);
05039         s = max_(d__1,d__2);
05040 /* Computing 2nd power */
05041         d__1 = theta / s;
05042         gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
05043         if (*stp > *stx) {
05044             gamma = -gamma;
05045         }
05046         p = gamma - *dp + theta;
05047         q = gamma - *dp + gamma + *dx;
05048         r__ = p / q;
05049         stpc = *stp + r__ * (*stx - *stp);
05050         stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
05051         if ((d__1 = stpc - *stp, abs_(d__1)) > (d__2 = stpq - *stp, abs_(d__2)))
05052                  {
05053             stpf = stpc;
05054         } else {
05055             stpf = stpq;
05056         }
05057         *brackt = TRUE_;
05058 /*     Third case: A lower function value, derivatives of the same sign, */
05059 /*     and the magnitude of the derivative decreases. */
05060     } else if (abs_(*dp) < abs_(*dx)) {
05061 /*        The cubic step is computed only if the cubic tends to infinity */
05062 /*        in the direction of the step or if the minimum of the cubic */
05063 /*        is beyond stp. Otherwise the cubic step is defined to be the */
05064 /*        secant step. */
05065         theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp;
05066 /* Computing MAX */
05067         d__1 = abs_(theta), d__2 = abs_(*dx), d__1 = max_(d__1,d__2), d__2 = abs_(
05068                 *dp);
05069         s = max_(d__1,d__2);
05070 /*        The case gamma = 0 only arises if the cubic does not tend */
05071 /*        to infinity in the direction of the step. */
05072 /* Computing MAX */
05073 /* Computing 2nd power */
05074         d__3 = theta / s;
05075         d__1 = 0., d__2 = d__3 * d__3 - *dx / s * (*dp / s);
05076         gamma = s * sqrt((max_(d__1,d__2)));
05077         if (*stp > *stx) {
05078             gamma = -gamma;
05079         }
05080         p = gamma - *dp + theta;
05081         q = gamma + (*dx - *dp) + gamma;
05082         r__ = p / q;
05083         if (r__ < 0. && gamma != 0.) {
05084             stpc = *stp + r__ * (*stx - *stp);
05085         } else if (*stp > *stx) {
05086             stpc = *stpmax;
05087         } else {
05088             stpc = *stpmin;
05089         }
05090         stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
05091         if (*brackt) {
05092 /*           A minimizer has been bracketed. If the cubic step is */
05093 /*           closer to stp than the secant step, the cubic step is */
05094 /*           taken, otherwise the secant step is taken. */
05095             if ((d__1 = stpc - *stp, abs_(d__1)) < (d__2 = stpq - *stp, abs_(
05096                     d__2))) {
05097                 stpf = stpc;
05098             } else {
05099                 stpf = stpq;
05100             }
05101             if (*stp > *stx) {
05102 /* Computing MIN */
05103                 d__1 = *stp + (*sty - *stp) * .66;
05104                 stpf = min_(d__1,stpf);
05105             } else {
05106 /* Computing MAX */
05107                 d__1 = *stp + (*sty - *stp) * .66;
05108                 stpf = max_(d__1,stpf);
05109             }
05110         } else {
05111 /*           A minimizer has not been bracketed. If the cubic step is */
05112 /*           farther from stp than the secant step, the cubic step is */
05113 /*           taken, otherwise the secant step is taken. */
05114             if ((d__1 = stpc - *stp, abs_(d__1)) > (d__2 = stpq - *stp, abs_(
05115                     d__2))) {
05116                 stpf = stpc;
05117             } else {
05118                 stpf = stpq;
05119             }
05120             stpf = min_(*stpmax,stpf);
05121             stpf = max_(*stpmin,stpf);
05122         }
05123 /*     Fourth case: A lower function value, derivatives of the same sign, */
05124 /*     and the magnitude of the derivative does not decrease. If the */
05125 /*     minimum is not bracketed, the step is either stpmin or stpmax, */
05126 /*     otherwise the cubic step is taken. */
05127     } else {
05128         if (*brackt) {
05129             theta = (*fp - *fy) * 3. / (*sty - *stp) + *dy + *dp;
05130 /* Computing MAX */
05131             d__1 = abs_(theta), d__2 = abs_(*dy), d__1 = max_(d__1,d__2), d__2 = 
05132                     abs_(*dp);
05133             s = max_(d__1,d__2);
05134 /* Computing 2nd power */
05135             d__1 = theta / s;
05136             gamma = s * sqrt(d__1 * d__1 - *dy / s * (*dp / s));
05137             if (*stp > *sty) {
05138                 gamma = -gamma;
05139             }
05140             p = gamma - *dp + theta;
05141             q = gamma - *dp + gamma + *dy;
05142             r__ = p / q;
05143             stpc = *stp + r__ * (*sty - *stp);
05144             stpf = stpc;
05145         } else if (*stp > *stx) {
05146             stpf = *stpmax;
05147         } else {
05148             stpf = *stpmin;
05149         }
05150     }
05151 /*     Update the interval which contains a minimizer. */
05152     if (*fp > *fx) {
05153         *sty = *stp;
05154         *fy = *fp;
05155         *dy = *dp;
05156     } else {
05157         if (sgnd < 0.) {
05158             *sty = *stx;
05159             *fy = *fx;
05160             *dy = *dx;
05161         }
05162         *stx = *stp;
05163         *fx = *fp;
05164         *dx = *dp;
05165     }
05166 /*     Compute the new step. */
05167     *stp = stpf;
05168     return 0;
05169 } /* dcstep_ */
05170 
05171 /* ====================== The end of dcstep ============================== */
05172 /* Subroutine */ int timer_(double *ttime)
05173 //double *ttime;
05174 {
05175  //   static float temp;
05176  //    extern double etime_(float *);
05177     static float tarray[2];
05178 
05179 /*     ********* */
05180 
05181 /*     Subroutine timer */
05182 
05183 /*     This subroutine is used to determine user time. In a typical */
05184 /*     application, the user time for a code segment requires calls */
05185 /*     to subroutine timer to determine the initial and final time. */
05186 
05187 /*     The subroutine statement is */
05188 
05189 /*       subroutine timer(ttime) */
05190 
05191 /*     where */
05192 
05193 /*       ttime is an output variable which specifies the user time. */
05194 
05195 /*     Argonne National Laboratory and University of Minnesota. */
05196 /*     MINPACK-2 Project. */
05197 
05198 /*     Modified October 1990 by Brett M. Averick. */
05199 
05200 /*     ********** */
05201 /*     The first element of the array tarray specifies user time */
05202     //temp = etime_(tarray);
05203     tarray[0]=0.0;
05204     tarray[1]=1.0;
05205     *ttime = (double) tarray[0];
05206     return 0;
05207 } /* timer_ */
05208 
05209 /* ====================== The end of timer =============================== */
05210 double dnrm2_(long int *n, double *x, long int *incx)
05211 /*long int *n;
05212 double *x;
05213 long int *incx;*/
05214 {
05215     /* System generated locals */
05216     long int i__1, i__2;
05217     double ret_val, d__1, d__2, d__3;
05218 
05219     /* Builtin functions */
05220     //double sqrt();
05221 
05222     /* Local variables */
05223     static long int i__;
05224     static double scale;
05225 
05226 /*     ********** */
05227 
05228 /*     Function dnrm2 */
05229 
05230 /*     Given a vector x of length n, this function calculates the */
05231 /*     Euclidean norm of x with stride incx. */
05232 
05233 /*     The function statement is */
05234 
05235 /*       double precision function dnrm2(n,x,incx) */
05236 
05237 /*     where */
05238 
05239 /*       n is a positive long int input variable. */
05240 
05241 /*       x is an input array of length n. */
05242 
05243 /*       incx is a positive long int variable that specifies the */
05244 /*         stride of the vector. */
05245 
05246 /*     Subprograms called */
05247 
05248 /*       FORTRAN-supplied ... abs, max, sqrt */
05249 
05250 /*     MINPACK-2 Project. February 1991. */
05251 /*     Argonne National Laboratory. */
05252 /*     Brett M. Averick. */
05253 
05254 /*     ********** */
05255     /* Parameter adjustments */
05256     --x;
05257 
05258     /* Function Body */
05259     ret_val = 0.;
05260     scale = 0.;
05261     i__1 = *n;
05262     i__2 = *incx;
05263     for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
05264 /* Computing MAX */
05265         d__2 = scale, d__3 = (d__1 = x[i__], abs_(d__1));
05266         scale = max_(d__2,d__3);
05267 /* L10: */
05268     }
05269     if (scale == 0.) {
05270         return ret_val;
05271     }
05272     i__2 = *n;
05273     i__1 = *incx;
05274     for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
05275 /* Computing 2nd power */
05276         d__1 = x[i__] / scale;
05277         ret_val += d__1 * d__1;
05278 /* L20: */
05279     }
05280     ret_val = scale * sqrt(ret_val);
05281     return ret_val;
05282 } /* dnrm2_ */
05283 
05284 /* ====================== The end of dnrm2 =============================== */
05285 double dpmeps_()
05286 {
05287     /* Initialized data */
05288 
05289     static double zero = 0.;
05290     static double one = 1.;
05291     static double two = 2.;
05292 
05293     /* System generated locals */
05294     long int i__1;
05295     double ret_val;
05296 
05297     /* Local variables */
05298     static double beta;
05299     static long int irnd;
05300     static volatile double temp, temp1, a, b;
05301     static long int i__;
05302     static double betah;
05303     static long int ibeta, negep;
05304     static double tempa;
05305     static long int itemp, it;
05306     static double betain;
05307 
05308 /*     ********** */
05309 
05310 /*     Subroutine dpeps */
05311 
05312 /*     This subroutine computes the machine precision parameter */
05313 /*     dpmeps as the smallest floating point number such that */
05314 /*     1 + dpmeps differs from 1. */
05315 
05316 /*     This subroutine is based on the subroutine machar described in */
05317 
05318 /*     W. J. Cody, */
05319 /*     MACHAR: A subroutine to dynamically determine machine parameters, */
05320 /*     ACM Transactions on Mathematical Software, 14, 1988, pages 303-311. */
05321 
05322 /*     The subroutine statement is: */
05323 
05324 /*       subroutine dpeps(dpmeps) */
05325 
05326 /*     where */
05327 
05328 /*       dpmeps is a double precision variable. */
05329 /*         On entry dpmeps need not be specified. */
05330 /*         On exit dpmeps is the machine precision. */
05331 
05332 /*     MINPACK-2 Project. February 1991. */
05333 /*     Argonne National Laboratory and University of Minnesota. */
05334 /*     Brett M. Averick. */
05335 
05336 /*     ******* */
05337 /*     determine ibeta, beta ala malcolm. */
05338     a = one;
05339     b = one;
05340 L10:
05341     a += a;
05342     temp = a + one;
05343     temp1 = temp - a;
05344     if (temp1 - one == zero) {
05345         goto L10;
05346     }
05347 L20:
05348     b += b;
05349     temp = a + b;
05350     itemp = (long int) (temp - a);
05351     if (itemp == 0) {
05352         goto L20;
05353     }
05354     ibeta = itemp;
05355     beta = (double) ibeta;
05356 /*     determine it, irnd. */
05357     it = 0;
05358     b = one;
05359 L30:
05360     ++it;
05361     b *= beta;
05362     temp = b + one;
05363     temp1 = temp - b;
05364     if (temp1 - one == zero) {
05365         goto L30;
05366     }
05367     irnd = 0;
05368     betah = beta / two;
05369     temp = a + betah;
05370     if (temp - a != zero) {
05371         irnd = 1;
05372     }
05373     tempa = a + beta;
05374     temp = tempa + betah;
05375     if (irnd == 0 && temp - tempa != zero) {
05376         irnd = 2;
05377     }
05378 /*     determine dpmeps. */
05379     negep = it + 3;
05380     betain = one / beta;
05381     a = one;
05382     i__1 = negep;
05383     for (i__ = 1; i__ <= i__1; ++i__) {
05384         a *= betain;
05385 /* L40: */
05386     }
05387 L50:
05388     temp = one + a;
05389     if (temp - one != zero) {
05390         goto L60;
05391     }
05392     a *= beta;
05393     goto L50;
05394 L60:
05395     ret_val = a;
05396     if (ibeta == 2 || irnd == 0) {
05397         goto L70;
05398     }
05399     a = a * (one + a) / two;
05400     temp = one + a;
05401     if (temp - one != zero) {
05402         ret_val = a;
05403     }
05404 L70:
05405     return ret_val;
05406 } /* dpmeps_ */
05407 
05408 /* ====================== The end of dpmeps ============================== */
05409 /* Subroutine */ int daxpy_(long int *n, double *da, double *dx, long int *incx, double *dy, long int *incy)
05410 /*long int *n;
05411 double *da, *dx;
05412 long int *incx;
05413 double *dy;
05414 long int *incy;*/
05415 {
05416     /* System generated locals */
05417     long int i__1;
05418 
05419     /* Local variables */
05420     static long int i__, m, ix, iy, mp1;
05421 
05422 
05423 /*     constant times a vector plus a vector. */
05424 /*     uses unrolled loops for increments equal to one. */
05425 /*     jack dongarra, linpack, 3/11/78. */
05426 
05427 
05428     /* Parameter adjustments */
05429     --dy;
05430     --dx;
05431 
05432     /* Function Body */
05433     if (*n <= 0) {
05434         return 0;
05435     }
05436     if (*da == 0.) {
05437         return 0;
05438     }
05439     if (*incx == 1 && *incy == 1) {
05440         goto L20;
05441     }
05442 
05443 /*        code for unequal increments or equal increments */
05444 /*          not equal to 1 */
05445 
05446     ix = 1;
05447     iy = 1;
05448     if (*incx < 0) {
05449         ix = (-(*n) + 1) * *incx + 1;
05450     }
05451     if (*incy < 0) {
05452         iy = (-(*n) + 1) * *incy + 1;
05453     }
05454     i__1 = *n;
05455     for (i__ = 1; i__ <= i__1; ++i__) {
05456         dy[iy] += *da * dx[ix];
05457         ix += *incx;
05458         iy += *incy;
05459 /* L10: */
05460     }
05461     return 0;
05462 
05463 /*        code for both increments equal to 1 */
05464 
05465 
05466 /*        clean-up loop */
05467 
05468 L20:
05469     m = *n % 4;
05470     if (m == 0) {
05471         goto L40;
05472     }
05473     i__1 = m;
05474     for (i__ = 1; i__ <= i__1; ++i__) {
05475         dy[i__] += *da * dx[i__];
05476 /* L30: */
05477     }
05478     if (*n < 4) {
05479         return 0;
05480     }
05481 L40:
05482     mp1 = m + 1;
05483     i__1 = *n;
05484     for (i__ = mp1; i__ <= i__1; i__ += 4) {
05485         dy[i__] += *da * dx[i__];
05486         dy[i__ + 1] += *da * dx[i__ + 1];
05487         dy[i__ + 2] += *da * dx[i__ + 2];
05488         dy[i__ + 3] += *da * dx[i__ + 3];
05489 /* L50: */
05490     }
05491     return 0;
05492 } /* daxpy_ */
05493 
05494 /* ====================== The end of daxpy =============================== */
05495 /* Subroutine */ int dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy)
05496 /*long int *n;
05497 double *dx;
05498 long int *incx;
05499 double *dy;
05500 long int *incy;*/
05501 {
05502     /* System generated locals */
05503     long int i__1;
05504 
05505     /* Local variables */
05506     static long int i__, m, ix, iy, mp1;
05507 
05508 
05509 /*     copies a vector, x, to a vector, y. */
05510 /*     uses unrolled loops for increments equal to one. */
05511 /*     jack dongarra, linpack, 3/11/78. */
05512 
05513 
05514     /* Parameter adjustments */
05515     --dy;
05516     --dx;
05517 
05518     /* Function Body */
05519     if (*n <= 0) {
05520         return 0;
05521     }
05522     if (*incx == 1 && *incy == 1) {
05523         goto L20;
05524     }
05525 
05526 /*        code for unequal increments or equal increments */
05527 /*          not equal to 1 */
05528 
05529     ix = 1;
05530     iy = 1;
05531     if (*incx < 0) {
05532         ix = (-(*n) + 1) * *incx + 1;
05533     }
05534     if (*incy < 0) {
05535         iy = (-(*n) + 1) * *incy + 1;
05536     }
05537     i__1 = *n;
05538     for (i__ = 1; i__ <= i__1; ++i__) {
05539         dy[iy] = dx[ix];
05540         ix += *incx;
05541         iy += *incy;
05542 /* L10: */
05543     }
05544     return 0;
05545 
05546 /*        code for both increments equal to 1 */
05547 
05548 
05549 /*        clean-up loop */
05550 
05551 L20:
05552     m = *n % 7;
05553     if (m == 0) {
05554         goto L40;
05555     }
05556     i__1 = m;
05557     for (i__ = 1; i__ <= i__1; ++i__) {
05558         dy[i__] = dx[i__];
05559 /* L30: */
05560     }
05561     if (*n < 7) {
05562         return 0;
05563     }
05564 L40:
05565     mp1 = m + 1;
05566     i__1 = *n;
05567     for (i__ = mp1; i__ <= i__1; i__ += 7) {
05568         dy[i__] = dx[i__];
05569         dy[i__ + 1] = dx[i__ + 1];
05570         dy[i__ + 2] = dx[i__ + 2];
05571         dy[i__ + 3] = dx[i__ + 3];
05572         dy[i__ + 4] = dx[i__ + 4];
05573         dy[i__ + 5] = dx[i__ + 5];
05574         dy[i__ + 6] = dx[i__ + 6];
05575 /* L50: */
05576     }
05577     return 0;
05578 } /* dcopy_ */
05579 
05580 /* ====================== The end of dcopy =============================== */
05581 double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy)
05582 /*long int *n;
05583 double *dx;
05584 long int *incx;
05585 double *dy;
05586 long int *incy;*/
05587 {
05588     /* System generated locals */
05589     long int i__1;
05590     double ret_val;
05591 
05592     /* Local variables */
05593     static long int i__, m;
05594     static double dtemp;
05595     static long int ix, iy, mp1;
05596 
05597 
05598 /*     forms the dot product of two vectors. */
05599 /*     uses unrolled loops for increments equal to one. */
05600 /*     jack dongarra, linpack, 3/11/78. */
05601 
05602 
05603     /* Parameter adjustments */
05604     --dy;
05605     --dx;
05606 
05607     /* Function Body */
05608     ret_val = 0.;
05609     dtemp = 0.;
05610     if (*n <= 0) {
05611         return ret_val;
05612     }
05613     if (*incx == 1 && *incy == 1) {
05614         goto L20;
05615     }
05616 
05617 /*        code for unequal increments or equal increments */
05618 /*          not equal to 1 */
05619 
05620     ix = 1;
05621     iy = 1;
05622     if (*incx < 0) {
05623         ix = (-(*n) + 1) * *incx + 1;
05624     }
05625     if (*incy < 0) {
05626         iy = (-(*n) + 1) * *incy + 1;
05627     }
05628     i__1 = *n;
05629     for (i__ = 1; i__ <= i__1; ++i__) {
05630         dtemp += dx[ix] * dy[iy];
05631         ix += *incx;
05632         iy += *incy;
05633 /* L10: */
05634     }
05635     ret_val = dtemp;
05636     return ret_val;
05637 
05638 /*        code for both increments equal to 1 */
05639 
05640 
05641 /*        clean-up loop */
05642 
05643 L20:
05644     m = *n % 5;
05645     if (m == 0) {
05646         goto L40;
05647     }
05648     i__1 = m;
05649     for (i__ = 1; i__ <= i__1; ++i__) {
05650         dtemp += dx[i__] * dy[i__];
05651 /* L30: */
05652     }
05653     if (*n < 5) {
05654         goto L60;
05655     }
05656 L40:
05657     mp1 = m + 1;
05658     i__1 = *n;
05659     for (i__ = mp1; i__ <= i__1; i__ += 5) {
05660         dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
05661                 i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ + 
05662                 4] * dy[i__ + 4];
05663 /* L50: */
05664     }
05665 L60:
05666     ret_val = dtemp;
05667     return ret_val;
05668 } /* ddot_ */
05669 
05670 /* ====================== The end of ddot ================================ */
05671 /* Subroutine */ int dpofa_(double *a, long int *lda, long int *n, long int *info)
05672 /*double *a;
05673 long int *lda, *n, *info;*/
05674 {
05675     /* System generated locals */
05676     long int a_dim1, a_offset, i__1, i__2, i__3;
05677 
05678     /* Builtin functions */
05679     //double sqrt();
05680 
05681     /* Local variables */
05682     //extern double ddot_();
05683     static long int j, k;
05684     static double s, t;
05685     static long int jm1;
05686 
05687 
05688 /*     dpofa factors a double precision symmetric positive definite */
05689 /*     matrix. */
05690 
05691 /*     dpofa is usually called by dpoco, but it can be called */
05692 /*     directly with a saving in time if  rcond  is not needed. */
05693 /*     (time for dpoco) = (1 + 18/n)*(time for dpofa) . */
05694 
05695 /*     on entry */
05696 
05697 /*        a       double precision(lda, n) */
05698 /*                the symmetric matrix to be factored.  only the */
05699 /*                diagonal and upper triangle are used. */
05700 
05701 /*        lda     long int */
05702 /*                the leading dimension of the array  a . */
05703 
05704 /*        n       long int */
05705 /*                the order of the matrix  a . */
05706 
05707 /*     on return */
05708 
05709 /*        a       an upper triangular matrix  r  so that  a = trans(r)*r */
05710 /*                where  trans(r)  is the transpose. */
05711 /*                the strict lower triangle is unaltered. */
05712 /*                if  info .ne. 0 , the factorization is not complete. */
05713 
05714 /*        info    long int */
05715 /*                = 0  for normal return. */
05716 /*                = k  signals an error condition.  the leading minor */
05717 /*                     of order  k  is not positive definite. */
05718 
05719 /*     linpack.  this version dated 08/14/78 . */
05720 /*     cleve moler, university of new mexico, argonne national lab. */
05721 
05722 /*     subroutines and functions */
05723 
05724 /*     blas ddot */
05725 /*     fortran sqrt */
05726 
05727 /*     internal variables */
05728 
05729 /*     begin block with ...exits to 40 */
05730 
05731 
05732     /* Parameter adjustments */
05733     a_dim1 = *lda;
05734     a_offset = 1 + a_dim1 * 1;
05735     a -= a_offset;
05736 
05737     /* Function Body */
05738     i__1 = *n;
05739     for (j = 1; j <= i__1; ++j) {
05740         *info = j;
05741         s = 0.;
05742         jm1 = j - 1;
05743         if (jm1 < 1) {
05744             goto L20;
05745         }
05746         i__2 = jm1;
05747         for (k = 1; k <= i__2; ++k) {
05748             i__3 = k - 1;
05749             t = a[k + j * a_dim1] - ddot_(&i__3, &a[k * a_dim1 + 1], &c__1, &
05750                     a[j * a_dim1 + 1], &c__1);
05751             t /= a[k + k * a_dim1];
05752             a[k + j * a_dim1] = t;
05753             s += t * t;
05754 /* L10: */
05755         }
05756 L20:
05757         s = a[j + j * a_dim1] - s;
05758 /*     ......exit */
05759         if (s <= 0.) {
05760             goto L40;
05761         }
05762         a[j + j * a_dim1] = sqrt(s);
05763 /* L30: */
05764     }
05765     *info = 0;
05766 L40:
05767     return 0;
05768 } /* dpofa_ */
05769 
05770 /* ====================== The end of dpofa =============================== */
05771 /* Subroutine */ int dscal_(long int *n, double *da, double *dx, long int *incx)
05772 /*long int *n;
05773 double *da, *dx;
05774 long int *incx;*/
05775 {
05776     /* System generated locals */
05777     long int i__1, i__2;
05778 
05779     /* Local variables */
05780     static long int i__, m, nincx, mp1;
05781 
05782 
05783 /*     scales a vector by a constant. */
05784 /*     uses unrolled loops for increment equal to one. */
05785 /*     jack dongarra, linpack, 3/11/78. */
05786 /*     modified 3/93 to return if incx .le. 0. */
05787 
05788 
05789     /* Parameter adjustments */
05790     --dx;
05791 
05792     /* Function Body */
05793     if (*n <= 0 || *incx <= 0) {
05794         return 0;
05795     }
05796     if (*incx == 1) {
05797         goto L20;
05798     }
05799 
05800 /*        code for increment not equal to 1 */
05801 
05802     nincx = *n * *incx;
05803     i__1 = nincx;
05804     i__2 = *incx;
05805     for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
05806         dx[i__] = *da * dx[i__];
05807 /* L10: */
05808     }
05809     return 0;
05810 
05811 /*        code for increment equal to 1 */
05812 
05813 
05814 /*        clean-up loop */
05815 
05816 L20:
05817     m = *n % 5;
05818     if (m == 0) {
05819         goto L40;
05820     }
05821     i__2 = m;
05822     for (i__ = 1; i__ <= i__2; ++i__) {
05823         dx[i__] = *da * dx[i__];
05824 /* L30: */
05825     }
05826     if (*n < 5) {
05827         return 0;
05828     }
05829 L40:
05830     mp1 = m + 1;
05831     i__2 = *n;
05832     for (i__ = mp1; i__ <= i__2; i__ += 5) {
05833         dx[i__] = *da * dx[i__];
05834         dx[i__ + 1] = *da * dx[i__ + 1];
05835         dx[i__ + 2] = *da * dx[i__ + 2];
05836         dx[i__ + 3] = *da * dx[i__ + 3];
05837         dx[i__ + 4] = *da * dx[i__ + 4];
05838 /* L50: */
05839     }
05840     return 0;
05841 } /* dscal_ */
05842 
05843 /* ====================== The end of dscal =============================== */
05844 /* Subroutine */ int dtrsl_(double *t, long int *ldt, long int *n, double *b, long int *job, long int *info)
05845 /*double *t;
05846 long int *ldt, *n;
05847 double *b;
05848 long int *job, *info;*/
05849 {
05850     /* System generated locals */
05851     long int t_dim1, t_offset, i__1, i__2;
05852 
05853     /* Local variables */
05854     static long int case__;
05855     //extern double ddot_(void);
05856     static double temp;
05857     static long int j;
05858     //extern /* Subroutine */ int daxpy_(double *a, long int *n, double *c, long int *d, long int *e);
05859     extern int daxpy_(long int *n, double *da, double *dx, long int *incx, double *dy, long int *incy);
05860     static long int jj;
05861 
05862 
05863 
05864 /*     dtrsl solves systems of the form */
05865 
05866 /*                   t * x = b */
05867 /*     or */
05868 /*                   trans(t) * x = b */
05869 
05870 /*     where t is a triangular matrix of order n. here trans(t) */
05871 /*     denotes the transpose of the matrix t. */
05872 
05873 /*     on entry */
05874 
05875 /*         t         double precision(ldt,n) */
05876 /*                   t contains the matrix of the system. the zero */
05877 /*                   elements of the matrix are not referenced, and */
05878 /*                   the corresponding elements of the array can be */
05879 /*                   used to store other information. */
05880 
05881 /*         ldt       long int */
05882 /*                   ldt is the leading dimension of the array t. */
05883 
05884 /*         n         long int */
05885 /*                   n is the order of the system. */
05886 
05887 /*         b         double precision(n). */
05888 /*                   b contains the right hand side of the system. */
05889 
05890 /*         job       long int */
05891 /*                   job specifies what kind of system is to be solved. */
05892 /*                   if job is */
05893 
05894 /*                        00   solve t*x=b, t lower triangular, */
05895 /*                        01   solve t*x=b, t upper triangular, */
05896 /*                        10   solve trans(t)*x=b, t lower triangular, */
05897 /*                        11   solve trans(t)*x=b, t upper triangular. */
05898 
05899 /*     on return */
05900 
05901 /*         b         b contains the solution, if info .eq. 0. */
05902 /*                   otherwise b is unaltered. */
05903 
05904 /*         info      long int */
05905 /*                   info contains zero if the system is nonsingular. */
05906 /*                   otherwise info contains the index of */
05907 /*                   the first zero diagonal element of t. */
05908 
05909 /*     linpack. this version dated 08/14/78 . */
05910 /*     g. w. stewart, university of maryland, argonne national lab. */
05911 
05912 /*     subroutines and functions */
05913 
05914 /*     blas daxpy,ddot */
05915 /*     fortran mod */
05916 
05917 /*     internal variables */
05918 
05919 
05920 /*     begin block permitting ...exits to 150 */
05921 
05922 /*        check for zero diagonal elements. */
05923 
05924     /* Parameter adjustments */
05925     t_dim1 = *ldt;
05926     t_offset = 1 + t_dim1 * 1;
05927     t -= t_offset;
05928     --b;
05929 
05930     /* Function Body */
05931     i__1 = *n;
05932     for (*info = 1; *info <= i__1; ++(*info)) {
05933 /*     ......exit */
05934         if (t[*info + *info * t_dim1] == 0.) {
05935             goto L150;
05936         }
05937 /* L10: */
05938     }
05939     *info = 0;
05940 
05941 /*        determine the task and go to it. */
05942 
05943     case__ = 1;
05944     if (*job % 10 != 0) {
05945         case__ = 2;
05946     }
05947     if (*job % 100 / 10 != 0) {
05948         case__ += 2;
05949     }
05950     switch ((int)case__) {
05951         case 1:  goto L20;
05952         case 2:  goto L50;
05953         case 3:  goto L80;
05954         case 4:  goto L110;
05955     }
05956 
05957 /*        solve t*x=b for t lower triangular */
05958 
05959 L20:
05960     b[1] /= t[t_dim1 + 1];
05961     if (*n < 2) {
05962         goto L40;
05963     }
05964     i__1 = *n;
05965     for (j = 2; j <= i__1; ++j) {
05966         temp = -b[j - 1];
05967         i__2 = *n - j + 1;
05968         daxpy_(&i__2, &temp, &t[j + (j - 1) * t_dim1], &c__1, &b[j], &c__1);
05969         b[j] /= t[j + j * t_dim1];
05970 /* L30: */
05971     }
05972 L40:
05973     goto L140;
05974 
05975 /*        solve t*x=b for t upper triangular. */
05976 
05977 L50:
05978     b[*n] /= t[*n + *n * t_dim1];
05979     if (*n < 2) {
05980         goto L70;
05981     }
05982     i__1 = *n;
05983     for (jj = 2; jj <= i__1; ++jj) {
05984         j = *n - jj + 1;
05985         temp = -b[j + 1];
05986         daxpy_(&j, &temp, &t[(j + 1) * t_dim1 + 1], &c__1, &b[1], &c__1);
05987         b[j] /= t[j + j * t_dim1];
05988 /* L60: */
05989     }
05990 L70:
05991     goto L140;
05992 
05993 /*        solve trans(t)*x=b for t lower triangular. */
05994 
05995 L80:
05996     b[*n] /= t[*n + *n * t_dim1];
05997     if (*n < 2) {
05998         goto L100;
05999     }
06000     i__1 = *n;
06001     for (jj = 2; jj <= i__1; ++jj) {
06002         j = *n - jj + 1;
06003         i__2 = jj - 1;
06004         b[j] -= ddot_(&i__2, &t[j + 1 + j * t_dim1], &c__1, &b[j + 1], &c__1);
06005         b[j] /= t[j + j * t_dim1];
06006 /* L90: */
06007     }
06008 L100:
06009     goto L140;
06010 
06011 /*        solve trans(t)*x=b for t upper triangular. */
06012 
06013 L110:
06014     b[1] /= t[t_dim1 + 1];
06015     if (*n < 2) {
06016         goto L130;
06017     }
06018     i__1 = *n;
06019     for (j = 2; j <= i__1; ++j) {
06020         i__2 = j - 1;
06021         b[j] -= ddot_(&i__2, &t[j * t_dim1 + 1], &c__1, &b[1], &c__1);
06022         b[j] /= t[j + j * t_dim1];
06023 /* L120: */
06024     }
06025 L130:
06026 L140:
06027 L150:
06028     return 0;
06029 } /* dtrsl_ */
06030 
06031