Procedures

This subroutine sets the bounds on the coefficient changes in slnpro.

Given m n-dimensional vectors p(j) as the first m columns of the matrix pmat1 and an n-vector r, this subroutine returns in ptnr the nearest point to r in the cone of points summation( coef(j)*p(j)), where coef(j) >= 0.0 for j=1,...,m (unless jflag > 0, which indicates failure). the number of vectors p(j) in the final corral is returned in ncor with their indices in icor, the distance is returned in dist, the number of major cycles (i.e. adding a vector) is returned in nmaj, and the number of minor cycles (i.e. removing a vector) is returned in nmin. if the user sets istrt1=0 the subroutne starts from scratch, but the user can set istrt1=1 and initially specify ncor, icor, and coef (noting that ncor must be <= n, and if j does not occur in icor, then coef(j) should be set to 0.0.)

CONMAX consists of two programs for solving the problem

This subroutine determines whether parprj violates any type -2 or type -1 (i.e. standard) constraints by more than tolcon, and if so it attempts to correct back to the feasible region. if it is successful it sets icorct=0 and replaces parprj by the corrected vector. if it is not successful it sets icorct=1 and leaves parprj unchanged. if no correction was needed it sets icorct=-1 and leaves parprj unchanged.

This subroutine uses fnset to compute confun(i,1) and the partial derivatives of the function whose value is in confun(i,1) for certain value(s) of i. note that we do not want the icntyp computed by fnset to override the icntyp (or jcntyp) carried into this subroutine in icntyp, so we use kcntyp when we call fnset. (the icntyp computed by fnset was stored earlier through a call to ercmp1 from conmax.)

This subprogram computes the dot product of vectors vec1 and vec2 of length lgth. vec1 and vec2 do not appear in function iloc since they are used only as input names for this subprogram, and so they don't need to have space reserved for them in the array work.

Given ncor n dimensional vectors as columns of the n by ncor matrix picor and an n dimensional vector r, this subroutine uses householder transformations to find the best least squares solution vec to the linear system of equations picor*vec = r, where vec is an ncor dimensional vector. if the rank of picor is (computationally) 0, the subroutine will return with the failure warning ihouse=1, otherwise it will return with ihouse=0. if the rank is > 0 but < ncor, then (ncor - rank) of the components of vec will be set to 0.0. the arays picor and r will not be changed by this subroutine. the subroutine will attempt up to numref iterative refinements of the solution, where the user can set numref as any nonnegative integer, but to get the most out of the iterative refinement process, the computation of the residual summ near the end of this subroutine should be done in higher precision than the other computations in the subroutine.

This function subprogram returns the subscript of the first element of array iarr relative to iwork (if the array is integer, i.e. 13 <= iarr <= 23) or relative to work (if the array is floating point, i.e. 1 <= iarr <= 12 or 24 <= iarr <= 48).

In this subroutine we are given a base vector zwork, a direction vector dvec, a scalar procor with emin = f(procor) = (the maximum type -2 and -1 error with parameters zwork + procor*dvec) < -tolcon, and a scalar p1 with p1 < procor and f1 = f(p1) > tolcon. we do a revised mullers method approach (with a solution contained in a shrinking interval) to attempt to adjust procor so that -tolcon <= f(procor) <= tolcon, but if we are not successful we return with the leftmost procor found satisfying emin = f(procor) < -tolcon on the theory that overcorrection is better than no correction. note that when corrct calls this subroutine it will have lumped the type -1 constraints in with the type -2 constraints using jcntyp, which is carried through this subroutine into subroutine ercmp1 in iwork.

This subroutine sets up the (nparm+1) by mactrk matrix pmat. for 1 <= j <= mactrk, the top nparm elements of column j of pmat will contain the negative of the gradient of active constraint j (if constraint j is of type 2, i.e. of the form abs(f(x) - f(parwrk,x)) <= w, the left side will be treated as f(x) - f(parwrk,x) if this quantity is nonnegative and will be treated as f(parwrk,x) - f(x) otherwise). the (nparm+1)st row of pmat will contain actdif, the right side of the inequalities gradient.vector >= actdif.

This subroutine increases rchdwn or rchin if it appears some constraints which should have been declared active were not so declared.

This subroutine attempts to refine the ndm dimensional vector wpt produced by wolfe by directly solving the system summation(pmat(i,j)*wpt(i), i=1,...,ndm) = -pmat(ndm+1,j) for j = icor(l), l=1,...,ncor. nresl resolvents are chosen by total pivoting. if nresl < ndm then the remaining ndm-nresl elements of wpt are kept form the old wpt. itrlm steps of iterative refinement are attempted at the end.

This subroutine computes a parameter vector parprj using fourth order runge kutta with h=-projct. h is negative since we want to approximate the parameters resulting from decreasing w by abs(h). if we do nstep steps then h=-projct/nstep.

This subroutine puts the (signed) indices of the mactrk active constraints in iact. it also sets the right side vector actdif for the wolfe subproblem.

This subroutine uses a modified quadratic fitting process to search for a projection factor procor for which the maximum of the left sides of the type -2 and -1 constraints evaluated at zwork + procor*dvec is <= tolcon. note that when corrct calls this subroutine it will have lumped the type -1 constraints in with the type -2 constraints using jcntyp, which is carried through this subroutine into subroutine ercmp1 in iwork. if searcr is able to force this maximum <= tolcon it will return with isrcr=0, with the minimum value found for the maximum in emin, with the corresponding projection factor in procor, with the number of times the maximum was computed in nsrch, and with the closest point found to the left with the maximum > tolcon in (p1,f1). the subroutine will begin by computing the maxima for procor = 1.0, 0.5, and 2.0, and if none of these maxima is <= tolcon and it is not the case that the maximum at 1.0 is <= the other two maxima the subroutine will return with the warning isrcr=1. the subroutine will also return with isrcr=1 if it would need to compute f more than limscr times, or the search interval length drops below tol1, or the quadratic fit becomes too flat. even in the event of a return with isrcr=1, emin, procor, and nsrch will be as above.

This subroutine uses a modified quadratic fitting process to search for the minimum of a function f. it requres an initial guess in projct, a tolerance tol1 on the search interval length, an upper bound prjlim on the minimizing point (which should be set very large if no upper bound is desired), and a way to compute f(x) for a given x. the subroutine will print a warning and return if it would need to compute f more than initlm times in the initialization or more than nadd additional times in the main part of the program. when the subroutine returns, it will have put the minimum location in projct, the minimum f value in emin, the f value for the initial projct in emin1, and the number of times it computed f in nsrch.

This subroutine sets up v for slnpro to solve a modified linearized (about the old parameters in param) version of our problem.

This subroutine performs a modified jordan elimination on the l-ll+1 by k matrix consisting of rows ll through l of v and columns 1 through k of v. The resolvent is v(ir,is).

This subroutine solves the linear programming problem

given m inequalities of the form a(k).x + b(k) <= 0.0 for k=1, ...,m, where a(k) and x are ndm dimensional vectors and b(k) are numbers, this subroutine returns the nearest point to the origin in the polytope defined by these inequalities (unless jflag > 0, which indicates failure). the user should put the mdm+1 dimensional vectors (a(k),b(k)) in the columns of pmat. the solution point will be returned in wpt, and will also be a linear combination of the a(k) vectors with (nonpositive) coefficients in the m dimensional vector wcoef. wcoef may not be accurate if refwl was used to refine wpt, which rarely happens. the number of vectors in the final corral will be returned in ncor with their indices in icor, and all entries of wcoef not corresponding to indices in icor will be zero. the distance will be returned in wdist, and the numbers of major and minor cycles in the cone subproblem will be returned in nmaj and nmin respectively. if the user sets istrt=0 the program will start from scratch, but the user can set istrt=1 (hot start) and specify ncor, icor, wcoef, and the factor s. (see later comments; set s=1.0 if no better value is available. set wcoef(j)=0.0 if icor(i) /= j for i=1,...,ncor.) (if inaccurate wcoef or s is used in a hot start attempt little will be lost, since ncor and icor are more important for a successful hot start than wcoef and s.) we must always have ncor <= ndm+1 in theory since the ncor ndm+1 dimensional vectors in a corral should be linearly independent, and in practice we will always require ncor <= ndm+1. if the user sets istrt=1 but the program fails, it will automatically try from scratch before giving up.