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.
To run the program, write a driver program which dimensions the arrays in the calling sequence for wolfe and sets the input variables as specified in the list below, then call subroutine wolfe.
This program was developed by ed kaufman, david leeming, and jerry taylor. the method used is an enhanced version of the method described in (wolfe, philip, finding the nearest point in a polytope, mathematical programming 11 (1976), 128-149).
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | Ndm |
The number of variables. it must be less than or equal to |
||
| integer, | intent(in) | :: | m |
The number of inequalities defining the polytope. it
must be less than or equal to |
||
| real(kind=wp), | intent(in) | :: | Pmat(Nparm+1,Numgr) |
This is an array whose kth column contains the vector (a(k),b(k)) for k=1,...,m, where the m inequalities a(k).x + b(k) <= 0.0 define the polytope whose nearest point to the origin we seek. the first dimension of pmat in the driver program must be exactly nparm+1, while the second dimension of pmat in the driver program must be at least numgr. if we actually want the nearest point in the polytope to some point y other than the origin, we translate y to the origin before calling wolfe, that is, call wolfe to find the nearest point z to the origin in the polytope defined by a(k).z + (b(k) + a(k).y) <= 0.0, then compute x = y + z. |
||
| integer, | intent(in) | :: | Istrt |
Set this equal to zero unless a hot start is desired-- see next paragraph of comments for more details. if istrt is set equal to 1, then s, wcoef, ncor, and icor must also be assigned values initially. |
||
| real(kind=wp), | intent(out) | :: | s |
You may ignore this scale factor unless you want to use the hot start option. |
||
| integer, | intent(out) | :: | Ncor |
This is the number of vectors (i.e. colunns of pmat) in the final corral. |
||
| integer, | intent(out) | :: | Icor(Nparm+1) |
This array contains the ncor indices of the vectors in the final corral. its dimension in the driver program must be at least nparm+1. |
||
| integer | :: | Iwork(Liwrk) |
Work array |
|||
| integer, | intent(in) | :: | Liwrk |
This is the dimension of |
||
| real(kind=wp) | :: | Work(Lwrk) |
Work array |
|||
| integer, | intent(in) | :: | Lwrk |
This is the dimension of |
||
| real(kind=wp) | :: | r(Nparm+1) |
Work array |
|||
| real(kind=wp) | :: | Coef(Numgr) |
Work array |
|||
| real(kind=wp) | :: | Ptnr(Nparm+1) |
Work array |
|||
| real(kind=wp) | :: | Pmat1(Nparm+1,Numgr) |
Work array |
|||
| integer, | intent(in) | :: | Nparm |
This is basically a dimension parameter here. it must
be greater than or equal to |
||
| integer, | intent(in) | :: | Numgr |
This is basically a dimension parameter here. it must
be greater than or equal to |
||
| real(kind=wp), | intent(out) | :: | Wcoef(Numgr) |
This will give the coefficients of the vectors a(k) needed to form a linear combination equal to the solution in wpt. its dimension in the driver program must be at least numgr. wcoef may not be accurate if it was necessary to call refwl to refine wpt, which rarely happens. |
||
| real(kind=wp), | intent(out) | :: | Wpt(Nparm) |
This will give the coordinates of the point we are seeking, namely the nearest point in the polytope to the origin. its dimension in the driver program must be at least nparm. |
||
| real(kind=wp), | intent(out) | :: | Wdist |
This will be the (minimized) euclidean distance of |
||
| integer, | intent(out) | :: | Nmaj |
This will be the number of major cycles used in |
||
| integer, | intent(out) | :: | Nmin |
This will be the number of major cycles used in |
||
| integer, | intent(out) | :: | Jflag |
This is a flag variable which is 0 in case of a normal solution and is positive otherwise (in which case the returned solution may be no good). |
subroutine wolfe(Ndm, m, Pmat, Istrt, s, Ncor, Icor, Iwork, Liwrk, Work, & Lwrk, r, Coef, Ptnr, Pmat1, Nparm, Numgr, Wcoef, Wpt, & Wdist, Nmaj, Nmin, Jflag) implicit none integer, intent(in) :: Ndm !! The number of variables. it must be less than or equal to `nparm`. integer, intent(in) :: m !! The number of inequalities defining the polytope. it !! must be less than or equal to `numgr`. integer, intent(in) :: Istrt !! Set this equal to zero unless a hot start is desired-- !! see next paragraph of comments for more details. if istrt is set !! equal to 1, then s, wcoef, ncor, and icor must also be assigned !! values initially. integer, intent(out) :: Ncor !! This is the number of vectors (i.e. colunns of pmat) in !! the final corral. integer, intent(in) :: Liwrk !! This is the dimension of `iwork`. It must be at least 7*nparm + 7*numgr + 3. integer, intent(in) :: Lwrk !! This is the dimension of `work`. it must be at least !! 2*nparm**2 + 4*numgr*nparm + 11*numgr + 27*nparm + 13. !! note that some storage could be saved by rewriting function !! subprogram iloc to take out all but the arrays needed (namely 1, 3, !! 4, 9, 28, 32, 34, 39 for work, 18, 23 for iwork) and scrunching !! work and iwork in iloc so the remaining arrays follow one after !! another. integer, intent(in) :: Nparm !! This is basically a dimension parameter here. it must !! be greater than or equal to `ndm`. integer, intent(in) :: Numgr !! This is basically a dimension parameter here. it must !! be greater than or equal to `m`. integer, intent(out) :: Nmaj !! This will be the number of major cycles used in `wolfe`. integer, intent(out) :: Nmin !! This will be the number of major cycles used in `wolfe`. integer, intent(out) :: Jflag !! This is a flag variable which is 0 in case of a normal !! solution and is positive otherwise (in which case the returned !! solution may be no good). real(wp), intent(out) :: s !! You may ignore this scale factor unless you want to use the hot start option. real(wp), intent(out) :: Wdist !! This will be the (minimized) euclidean distance of `wpt` from the origin. integer, intent(out) :: Icor(Nparm + 1) !! This array contains the ncor indices of the vectors in !! the final corral. its dimension in the driver program must be at !! least nparm+1. real(wp), intent(in) :: Pmat(Nparm + 1, Numgr) !! This is an array whose kth column contains the vector !! (a(k),b(k)) for k=1,...,m, where the m inequalities a(k).x + b(k) !! <= 0.0 define the polytope whose nearest point to the origin we !! seek. the first dimension of pmat in the driver program must be !! exactly nparm+1, while the second dimension of pmat in the driver !! program must be at least numgr. !! if we actually want the nearest point in the polytope to some point !! y other than the origin, we translate y to the origin before calling !! wolfe, that is, call wolfe to find the nearest point z to the origin !! in the polytope defined by a(k).z + (b(k) + a(k).y) <= 0.0, then !! compute x = y + z. integer :: Iwork(Liwrk) !! Work array real(wp) :: Work(Lwrk) !! Work array real(wp) :: r(Nparm + 1) !! Work array real(wp) :: Coef(Numgr) !! Work array real(wp) :: Ptnr(Nparm + 1) !! Work array real(wp) :: Pmat1(Nparm + 1, Numgr) !! Work array real(wp), intent(out) :: Wcoef(Numgr) !! This will give the coefficients of the vectors a(k) !! needed to form a linear combination equal to the solution in wpt. !! its dimension in the driver program must be at least numgr. !! wcoef may not be accurate if it was necessary to call refwl to !! refine wpt, which rarely happens. real(wp), intent(out) :: Wpt(Nparm) !! This will give the coordinates of the point we are seeking, !! namely the nearest point in the polytope to the origin. its dimension !! in the driver program must be at least nparm. real(wp) :: ab, bk, dist, fackp, facsc, & fact, quot, s1, s1hi, s1low, & scl, scl1, scl1a, tol, tol1, & tols, v1, violm, vmax integer :: i, ilc18, ilc28, ilc32, ilc34, ilc39, & ind, iref, istrt1, itcon, iup, j, & jmax, k, l, lmcon, n ! set machine and precision dependent constants for wolfe. tol = ten*ten*spcmn tol1 = (ten**4)*spcmn tols = sqrt(spcmn) iref = 0 violm = one/two lmcon = 3 itcon = 0 iup = 0 s1low = ten*ten*ten*spcmn s1hi = one - s1low ! make sure s1low <= one third and s1hi >= two thirds to avoid ! squeezing the allowable region for s1 too tightly (or even making it ! empty). if (s1low > one/three) s1low = one/three if (s1hi < two/three) s1hi = two/three facsc = ten*ten*ten*ten fackp = facsc ilc18 = iloc(18, nparm, numgr) ilc28 = iloc(28, nparm, numgr) ilc32 = iloc(32, nparm, numgr) ilc34 = iloc(34, nparm, numgr) ilc39 = iloc(39, nparm, numgr) n = ndm + 1 istrt1 = istrt do i = 1, ndm r(i) = zero end do r(n) = one ! now compute the scale factor scl, whose main purpose is to avoid ! having all vectors in pmat with positive last component form an angle ! close to 90 degrees with r = (0...0 1), which can cause numerical ! problems. we will compute scl = min(max(abs(a(i,k)): 1 <= i <= ! ndm)/b(k), b(k) >= tols, 1 <= k <= m) unless no b(k) is >= ! tols, in which case we set scl=1.0, or some b(k) is >= tols but ! scl would be < tol, in which case we set scl = tol. 100 scl = one ind = 0 do k = 1, m bk = pmat(n, k) if (bk >= tols) then quot = zero do i = 1, ndm ab = abs(pmat(i, k)) if (ab > quot) quot = ab end do quot = quot/bk if (ind > 0) then if (quot >= scl) cycle end if ind = 1 scl = quot end if end do 300 if (scl < tol) scl = tol ! put scaled pmat into pmat1 for use in conenr. pmat itself will remain ! unchanged. 400 do j = 1, m do i = 1, ndm pmat1(i, j) = pmat(i, j)/scl end do pmat1(n, j) = pmat(n, j) end do ! now do a normal scaling on each column of pmat1 which has an element ! with absolute value >= tol1. do j = 1, m scl1 = zero do i = 1, n ab = abs(pmat1(i, j)) if (ab > scl1) scl1 = ab end do if (scl1 >= tol1) then do i = 1, n pmat1(i, j) = pmat1(i, j)/scl1 end do if (istrt1 > 0) coef(j) = wcoef(j)*scl1 ! also put a scaled version of wcoef into coef if istrt1=1. else if (istrt1 > 0) then coef(j) = wcoef(j) end if end do ! if istrt1=1, for use in conenr set coef = (-s1*scl**2)*coef, where ! s1 = s/(s + (1.0-s)*scl**2) is the s value in the scaled situation. ! note that a partly scaled version of wcoef (see loop ending with the ! statement numbered 190 above) is already in coef if istrt1=1. if (istrt1 > 0) then ! if we had ncor > n, reset ncor to n. if (ncor > n) ncor = n fact = -(s/(s + (one - s)*scl**2))*scl**2 do j = 1, m coef(j) = fact*coef(j) end do end if ! call conenr to compute the nearest point to r in the cone of ! nonnegative linear combinations of columns of pmat1. call conenr(n, m, pmat1, r, istrt1, ncor, icor, tol, iwork, liwrk, work, & lwrk, work(ilc39), work(ilc32), work(ilc28), nparm, numgr, & coef, ptnr, dist, nmaj, nmin, jflag) ! if jflag=3 then conenr has failed, possibly because scl was too large. if (jflag /= 3) then ! here jflag /= 3 and we compute s1 = 1.0 - ptnr(n). s1 = one - ptnr(n) if (s1 >= s1low) then ! here jflag /= 3 and s1 >= s1low, so if also s1 <= s1hi we accept ! the result from conenr and move on. if (s1 <= s1hi) goto 700 ! here jflag /= 3 and s1 > s1hi, so if itcon < lmcon we try ! again with larger scl. ! if here jflag=0 and ncor=0 the nearest point to the origin in the ! polytope appears to be the origin so we forego adjusting scl. if (jflag /= 0) goto 600 if (ncor > 0) goto 600 goto 700 end if end if ! here jflag=3 or s1 < s1low, so if itcon < lmcon we try again with ! smaller scl. if (itcon < lmcon) then ! here we increment itcon and if scl was not already very small we ! decrease it and try conenr again. itcon = itcon + 1 if (iup < 0) then else if (iup == 0) then ! here iup=0 so either we are just starting (in which case we set iup=-1 ! to indicate we are in a phase of decreasing scl) or we are oscillating. if (itcon <= 1) then iup = -1 else facsc = sqrt(facsc) end if else ! here iup=1 and we have oscillation in the search for a usable scl so ! we replace the correction factor by its square root and reset iup to ! 0 to indicate oscillation. iup = 0 facsc = sqrt(facsc) end if ! here we decrease scl if it was not already very small. if (scl >= (one + one/ten)*tol) then scl = scl/facsc goto 300 end if end if ! here we were unable to get an acceptable s1 from conenr so we set ! jflag=4 as a warning and return. first try again from scratch if this ! has not been done. 500 if (istrt1 <= 0) then jflag = 4 return else istrt1 = 0 itcon = 0 iref = 0 iup = 0 facsc = fackp goto 100 end if 600 if (itcon >= lmcon) goto 500 itcon = itcon + 1 if (iup < 0) then ! here iup=-1 and we have oscillation in the search for a usable scl so ! we replace the correction factor by its square root and set iup=0 ! to indicate oscillation. iup = 0 facsc = sqrt(facsc) scl = scl*facsc else if (iup == 0) then ! here iup=0 so either we are just starting (in which case we set iup=1 ! to indicate we are in a phase of increasing scl) or we are oscillating. if (itcon <= 1) then iup = 1 scl = scl*facsc else facsc = sqrt(facsc) scl = scl*facsc end if else scl = scl*facsc end if goto 400 ! here conenr may have succeeded and we compute the nearest point ! (wpt,s1)=r-ptnr to r from the dual of the cone described earlier. ! this new cone is the set of (x,t) such that (a(k)/scl,b(k)).(x,t) <= ! 0.0 for k=1,...,m. 700 do i = 1, ndm wpt(i) = -ptnr(i) end do ! divide wpt by s1*scl. do i = 1, ndm wpt(i) = wpt(i)/(s1*scl) end do ! compute the maximum wolfe constraint violation as a check. 800 do j = 1, m v1 = pmat(n, j) do i = 1, ndm v1 = v1 + pmat(i, j)*wpt(i) end do if (j > 1) then if (v1 <= vmax) goto 900 end if jmax = j vmax = v1 900 end do ! if vmax <= violm we reset jflag to 0 and accept the result. if (vmax <= violm) then jflag = 0 ! divide the coefficients by -s1*scl**2. do j = 1, m wcoef(j) = -coef(j)/(s1*scl**2) end do ! we now reconstruct the normal scaling factors computed in the loop ! ending with the statement labelled 190 in this subroutine. in a later ! version of this subroutine an array may be created to store these in ! that loop, but for now we avoid the extra storage and programming work ! of fiddling with the variable dimensioning. to recreate the factor ! scl1 corresponding to column j, we compute the maximum absolute value ! of the first ndm elements of pmat in this column, divide it by scl, take ! the maximum of this and abs(pmat(ndm+1,j)), and take scl1 to be this ! value unless it is less than tol1, in which we (in effect) take scl1=1.0. ! finally, since wcoef(j) was computed with the jth column of pmat divided ! by scl1 it contains a hidden factor of scl1, which we divide out. do j = 1, m scl1a = zero do i = 1, ndm ab = abs(pmat(i, j)) if (ab > scl1a) scl1a = ab end do scl1 = scl1a/scl ab = abs(pmat(ndm + 1, j)) if (ab > scl1) scl1 = ab if (scl1 >= tol1) wcoef(j) = wcoef(j)/scl1 end do ! compute the s value for the unscaled situation. s = s1/(s1 + (one - s1)/scl**2) ! copy wpt into ptnr to get the right dimension for dotprd and compute ! the distance. do i = 1, ndm ptnr(i) = wpt(i) end do ptnr(n) = zero wdist = sqrt(dotprd(ndm, ptnr, ptnr, nparm)) return else ! here vmax is too large. if (iref > 0) then ! here we have unsuccessfully tried to refine wpt with refwl at least ! once. if ncor < ndm and the worst violation occurred outside ! icor we will put it in icor and try refwl again, otherwise we will ! set jflag=7 and return (first trying from scratch if this has not ! been done). if (ncor >= ndm) goto 1000 if (ncor > 0) then do l = 1, ncor if (jmax == icor(l)) goto 1000 end do end if ncor = ncor + 1 icor(ncor) = jmax end if ! increment iref and call refwl to attempt to refine wpt, then go back ! and recheck the maximum constraint violation. iref = iref + 1 call refwl(ndm, ncor, icor, pmat, pmat1, nparm, numgr, iwork(ilc18), & work(ilc34), wpt) goto 800 end if 1000 if (istrt1 > 0) then istrt1 = 0 itcon = 0 iref = 0 iup = 0 facsc = fackp goto 100 end if jflag = 7 end subroutine wolfe