| 1 | c======================================================================= |
|---|
| 2 | subroutine bvls(key, m, n, a, b, bl, bu, maxIter, x, w, act, zz, |
|---|
| 3 | + istate, istop, msg, loopA) |
|---|
| 4 | c======================================================================= |
|---|
| 5 | c |
|---|
| 6 | c$$$$ calls qr |
|---|
| 7 | c--------------------Bounded Variable Least Squares--------------------- |
|---|
| 8 | c |
|---|
| 9 | c Robert L. Parker and Philip B. Stark Version 3/19/90 |
|---|
| 10 | c |
|---|
| 11 | c Robert L. Parker Philip B. Stark |
|---|
| 12 | c Scripps Institution of Oceanography Department of Statistics |
|---|
| 13 | c University of California, San Diego University of California |
|---|
| 14 | c La Jolla CA 92093 Berkeley CA 94720-3860 |
|---|
| 15 | c rlparker@ucsd.edu stark@stat.berkeley.edu |
|---|
| 16 | c |
|---|
| 17 | c Copyright of this software is reserved by the authors; however, this |
|---|
| 18 | c algorithm and subroutine may be used freely for non-commercial |
|---|
| 19 | c purposes, and may be distributed freely for non-commercial purposes. |
|---|
| 20 | c |
|---|
| 21 | c The authors do not warrant this software in any way: use it at your |
|---|
| 22 | c own risk. |
|---|
| 23 | c |
|---|
| 24 | c |
|---|
| 25 | c See the article ``Bounded Variable Least Squares: An Algorithm and |
|---|
| 26 | c Applications'' by P.B. Stark and R.L. Parker, in the journal |
|---|
| 27 | c Computational Statistics, in press (1995) for further description |
|---|
| 28 | c and applications to minimum l-1, l-2 and l-infinity fitting problems, |
|---|
| 29 | c as well as finding bounds on linear functionals subject to bounds on |
|---|
| 30 | c variables and fitting linear data within l-1, l-2 or l-infinity |
|---|
| 31 | c measures of misfit. |
|---|
| 32 | c |
|---|
| 33 | c BVLS solves the problem: |
|---|
| 34 | c |
|---|
| 35 | c min || a.x - b || such that bl <= x <= bu |
|---|
| 36 | c 2 |
|---|
| 37 | c where |
|---|
| 38 | c x is an unknown n-vector |
|---|
| 39 | c a is a given m by n matrix |
|---|
| 40 | c b is a given m-vector |
|---|
| 41 | c bl is a given n-vector of lower bounds on the |
|---|
| 42 | c components of x. |
|---|
| 43 | c bu is a given n-vector of upper bounds on the |
|---|
| 44 | c components of x. |
|---|
| 45 | c |
|---|
| 46 | c |
|---|
| 47 | c----------------------------------------------------------------------- |
|---|
| 48 | c Input parameters: |
|---|
| 49 | c |
|---|
| 50 | c m, n, a, b, bl, bu see above. Let mm=min(m,n). |
|---|
| 51 | c maxIter - max number of iterations. Default value by routine authors |
|---|
| 52 | c was 3*n (changed by Dmitrey) |
|---|
| 53 | c |
|---|
| 54 | c If key = 0, the subroutine solves the problem from scratch. |
|---|
| 55 | c |
|---|
| 56 | c If key > 0 the routine initializes using the user's guess about |
|---|
| 57 | c which components of x are `active', i.e. are stricly within their |
|---|
| 58 | c bounds, which are at their lower bounds, and which are at their |
|---|
| 59 | c upper bounds. This information is supplied through the array |
|---|
| 60 | c istate. istate(n+1) should contain the total number of components |
|---|
| 61 | c at their bounds (the `bound variables'). The absolute values of the |
|---|
| 62 | c first nbound=istate(n+1) entries of istate are the indices |
|---|
| 63 | c of these `bound' components of x. The sign of istate(j), j=1,..., |
|---|
| 64 | c nbound, indicates whether x(|istate(j)|) is at its upper or lower |
|---|
| 65 | c bound. istate(j) is positive if the component is at its upper |
|---|
| 66 | c bound, negative if the component is at its lower bound. |
|---|
| 67 | c istate(j), j=nbound+1,...,n contain the indices of the components |
|---|
| 68 | c of x that are active (i.e. are expected to lie strictly within |
|---|
| 69 | c their bounds). When key > 0, the routine initially sets the active |
|---|
| 70 | c components to the averages of their upper and lower bounds: |
|---|
| 71 | c x(j)=(bl(j)+bu(j))/2, for j in the active set. |
|---|
| 72 | c |
|---|
| 73 | c----------------------------------------------------------------------- |
|---|
| 74 | c Output parameters: |
|---|
| 75 | c |
|---|
| 76 | c x the solution vector. |
|---|
| 77 | c |
|---|
| 78 | c w(1) the minimum 2-norm || a.x-b ||. |
|---|
| 79 | c |
|---|
| 80 | c istate vector indicating which components of x are active and |
|---|
| 81 | c which are at their bounds (see the previous paragraph). |
|---|
| 82 | c istate can be supplied to the routine to give it a good |
|---|
| 83 | c starting guess for the solution. |
|---|
| 84 | c |
|---|
| 85 | c loopA number of iterations taken in the main loop, Loop A. |
|---|
| 86 | c istop stop condition case, should be positive if converged |
|---|
| 87 | c msg stop condition case text description |
|---|
| 88 | c |
|---|
| 89 | c----------------------------------------------------------------------- |
|---|
| 90 | c Working arrays: |
|---|
| 91 | c |
|---|
| 92 | c w dimension n. act dimension m*(mm+2). |
|---|
| 93 | c zz dimension m. istate dimension n+1. |
|---|
| 94 | c |
|---|
| 95 | c----------------------------------------------------------------------- |
|---|
| 96 | c Method: active variable method along the general plan of NNLS by |
|---|
| 97 | c Lawson & Hanson, "Solving Least Squares Problems," 1974. See |
|---|
| 98 | c Algorithm 23.10. Step numbers in comment statements refer to their |
|---|
| 99 | c scheme. |
|---|
| 100 | c For more details and further uses, see the article |
|---|
| 101 | c "Bounded Variable Least Squares: An Algorithm and Applications" |
|---|
| 102 | c by Stark and Parker in 1995 Computational Statistics. |
|---|
| 103 | c |
|---|
| 104 | c----------------------------------------------------------------------- |
|---|
| 105 | c A number of measures are taken to enhance numerical reliability: |
|---|
| 106 | c |
|---|
| 107 | c 1. As noted by Lawson and Hanson, roundoff errors in the computation |
|---|
| 108 | c of the gradient of the misfit may cause a component on the bounds |
|---|
| 109 | c to appear to want to become active, yet when the component is added |
|---|
| 110 | c to the active set, it moves away from the feasible region. In this |
|---|
| 111 | c case the component is not made active, the gradient of the misfit |
|---|
| 112 | c with respect to a change in that component is set to zero, and the |
|---|
| 113 | c program returns to the Kuhn-Tucker test. Flag ifrom5 is used in |
|---|
| 114 | c this test, which occurs at the end of Step 6. |
|---|
| 115 | c |
|---|
| 116 | c |
|---|
| 117 | c 2. When the least-squares minimizer after Step 6 is infeasible, it |
|---|
| 118 | c is used in a convex interpolation with the previous solution to |
|---|
| 119 | c obtain a feasible vector. The constant in this interpolation is |
|---|
| 120 | c supposed to put at least one component of x on a bound. There can |
|---|
| 121 | c be difficulties: |
|---|
| 122 | c |
|---|
| 123 | c 2a. Sometimes, due to roundoff, no interpolated component ends up on |
|---|
| 124 | c a bound. The code in Step 11 uses the flag jj, computed in Step 8, |
|---|
| 125 | c to ensure that at least the component that determined the |
|---|
| 126 | c interpolation constant alpha is moved to the appropriate bound. |
|---|
| 127 | c This guarantees that what Lawson and Hanson call `Loop B' is finite. |
|---|
| 128 | c |
|---|
| 129 | c 2b. The code in Step 11 also incorporates Lawson and Hanson's feature |
|---|
| 130 | c that any components remaining infeasible at this stage (which must |
|---|
| 131 | c be due to roundoff) are moved to their nearer bound. |
|---|
| 132 | c |
|---|
| 133 | c |
|---|
| 134 | c 3. If the columns of a passed to qr are linearly dependent, the new |
|---|
| 135 | c potentially active component is not introduced: the gradient of the |
|---|
| 136 | c misfit with respect to that component is set to zero, and control |
|---|
| 137 | c returns to the Kuhn-Tucker test. |
|---|
| 138 | c |
|---|
| 139 | c |
|---|
| 140 | c 4. When some of the columns of a are approximately linearly |
|---|
| 141 | c dependent, we have observed cycling of active components: a |
|---|
| 142 | c component just moved to a bound desires immediately to become |
|---|
| 143 | c active again; qr allows it to become active and a different |
|---|
| 144 | c component is moved to its bound. This component immediately wants |
|---|
| 145 | c to become active, which qr allows, and the original component is |
|---|
| 146 | c moved back to its bound. We have taken two steps to avoid this |
|---|
| 147 | c problem: |
|---|
| 148 | c |
|---|
| 149 | c 4a. First, the column of the matrix a corresponding to the new |
|---|
| 150 | c potentially active component is passed to qr as the last column of |
|---|
| 151 | c its matrix. This ordering tends to make a component recently moved |
|---|
| 152 | c to a bound fail the test mentioned in (1), above. |
|---|
| 153 | c |
|---|
| 154 | c 4b. Second, we have incorporated a test that prohibits short cycles. |
|---|
| 155 | c If the most recent successful change to the active set was to move |
|---|
| 156 | c the component x(jj) to a bound, x(jj) is not permitted to reenter |
|---|
| 157 | c the solution at this stage. This test occurs just after checking |
|---|
| 158 | c the Kuhn-Tucker conditions, and uses the flag jj, set in Step 8. |
|---|
| 159 | c The flag jj is reset after Step 6 if Step 6 was entered from |
|---|
| 160 | c Step 5 indicating that a new component has successfully entered the |
|---|
| 161 | c active set. The test for resetting jj uses the flag ifrom5, |
|---|
| 162 | c which will not equal zero in case Step 6 was entered from Step 5. |
|---|
| 163 | c |
|---|
| 164 | c |
|---|
| 165 | c----------------------------------------------------------------------- |
|---|
| 166 | |
|---|
| 167 | cf2py intent(out) x, w, istop, msg, loopa |
|---|
| 168 | c-f-2-p-y intent(-in) key, a, b, bl, bu, maxIter, istate |
|---|
| 169 | |
|---|
| 170 | c Initialize flags, etc. |
|---|
| 171 | implicit double precision (a-h, o-z) |
|---|
| 172 | character *(64) msg |
|---|
| 173 | dimension a(m,n), b(m), x(n), bl(n), bu(n) |
|---|
| 174 | dimension w(n), act(m,m+2), zz(m), istate(n+1) |
|---|
| 175 | |
|---|
| 176 | data eps/1.0d-11/ |
|---|
| 177 | c |
|---|
| 178 | c----------------------First Executable Statement----------------------- |
|---|
| 179 | c |
|---|
| 180 | c Step 1. Initialize everything--active and bound sets, initial |
|---|
| 181 | c values, etc. |
|---|
| 182 | c |
|---|
| 183 | |
|---|
| 184 | istop = 1000 |
|---|
| 185 | msg = 'Converged' |
|---|
| 186 | |
|---|
| 187 | mm=min(m,n) |
|---|
| 188 | mm1 = mm + 1 |
|---|
| 189 | jj = 0 |
|---|
| 190 | ifrom5 = 0 |
|---|
| 191 | c Check consistency of given bounds bl, bu. |
|---|
| 192 | bdiff = 0.0 |
|---|
| 193 | do 1005 j=1, n |
|---|
| 194 | bdiff=max(bdiff, bu(j)-bl(j)) |
|---|
| 195 | if (bl(j) .gt. bu(j)) then |
|---|
| 196 | istop = -101 |
|---|
| 197 | msg = 'Inconsistent bounds in BVLS.' |
|---|
| 198 | return |
|---|
| 199 | c stop |
|---|
| 200 | endif |
|---|
| 201 | 1005 continue |
|---|
| 202 | if (bdiff .eq. 0.0) then |
|---|
| 203 | istop = -102 |
|---|
| 204 | msg = 'No free variables in BVLS--check input bounds' |
|---|
| 205 | return |
|---|
| 206 | c stop |
|---|
| 207 | endif |
|---|
| 208 | c |
|---|
| 209 | c In a fresh initialization (key = 0) bind all variables at their lower |
|---|
| 210 | c bounds. If (key != 0), use the supplied istate vector to |
|---|
| 211 | c initialize the variables. istate(n+1) contains the number of |
|---|
| 212 | c bound variables. The absolute values of the first |
|---|
| 213 | c nbound=istate(n+1) entries of istate are the indices of the bound |
|---|
| 214 | c variables. The sign of each entry determines whether the indicated |
|---|
| 215 | c variable is at its upper (positive) or lower (negative) bound. |
|---|
| 216 | if (key .eq. 0) then |
|---|
| 217 | nbound=n |
|---|
| 218 | nact=0 |
|---|
| 219 | do 1010 j=1, nbound |
|---|
| 220 | istate(j)=-j |
|---|
| 221 | 1010 continue |
|---|
| 222 | else |
|---|
| 223 | nbound=istate(n+1) |
|---|
| 224 | endif |
|---|
| 225 | nact=n - nbound |
|---|
| 226 | if ( nact .gt. mm ) then |
|---|
| 227 | istop = -103 |
|---|
| 228 | msg = 'Too many active variables in BVLS starting solution!' |
|---|
| 229 | return |
|---|
| 230 | c stop |
|---|
| 231 | endif |
|---|
| 232 | do 1100 k=1, nbound |
|---|
| 233 | j=abs(istate(k)) |
|---|
| 234 | if (istate(k) .lt. 0) x(j)=bl(j) |
|---|
| 235 | if (istate(k) .gt. 0) x(j)=bu(j) |
|---|
| 236 | 1100 continue |
|---|
| 237 | c |
|---|
| 238 | c In a warm start (key != 0) initialize the active variables to |
|---|
| 239 | c (bl+bu)/2. This is needed in case the initial qr results in |
|---|
| 240 | c active variables out-of-bounds and Steps 8-11 get executed the |
|---|
| 241 | c first time through. |
|---|
| 242 | do 1150 k=nbound+1,n |
|---|
| 243 | kk=istate(k) |
|---|
| 244 | x(kk)=(bu(kk)+bl(kk))/2 |
|---|
| 245 | 1150 continue |
|---|
| 246 | c |
|---|
| 247 | c Compute bnorm, the norm of the data vector b, for reference. |
|---|
| 248 | bsq=0.0 |
|---|
| 249 | do 1200 i=1, m |
|---|
| 250 | bsq=bsq + b(i)**2 |
|---|
| 251 | 1200 continue |
|---|
| 252 | bnorm=sqrt(bsq) |
|---|
| 253 | c |
|---|
| 254 | c-----------------------------Main Loop--------------------------------- |
|---|
| 255 | c |
|---|
| 256 | c Initialization complete. Begin major loop (Loop A). |
|---|
| 257 | do 15000 loopA=1, maxIter |
|---|
| 258 | c |
|---|
| 259 | c Step 2. |
|---|
| 260 | c Initialize the negative gradient vector w(*). |
|---|
| 261 | 2000 obj=0.0 |
|---|
| 262 | do 2050 j=1, n |
|---|
| 263 | w(j)=0.0 |
|---|
| 264 | 2050 continue |
|---|
| 265 | c |
|---|
| 266 | c Compute the residual vector b-a.x , the negative gradient vector |
|---|
| 267 | c w(*), and the current objective value obj = || a.x - b ||. |
|---|
| 268 | c The residual vector is stored in the mm+1'st column of act(*,*). |
|---|
| 269 | do 2300 i=1, m |
|---|
| 270 | ri=b(i) |
|---|
| 271 | do 2100 j=1, n |
|---|
| 272 | ri=ri - a(i,j)*x(j) |
|---|
| 273 | 2100 continue |
|---|
| 274 | obj=obj + ri**2 |
|---|
| 275 | do 2200 j=1, n |
|---|
| 276 | w(j)=w(j) + a(i,j)*ri |
|---|
| 277 | 2200 continue |
|---|
| 278 | act(i,mm1)=ri |
|---|
| 279 | 2300 continue |
|---|
| 280 | c |
|---|
| 281 | c Converged? Stop if the misfit << || b ||, or if all components are |
|---|
| 282 | c active (unless this is the first iteration from a warm start). |
|---|
| 283 | if (sqrt(obj) .le. bnorm*eps .or. |
|---|
| 284 | + (loopA .gt. 1 .and. nbound .eq. 0)) then |
|---|
| 285 | istate(n+1)=nbound |
|---|
| 286 | w(1)=sqrt(obj) |
|---|
| 287 | return |
|---|
| 288 | endif |
|---|
| 289 | c |
|---|
| 290 | c Add the contribution of the active components back into the residual. |
|---|
| 291 | do 2500 k=nbound+1, n |
|---|
| 292 | j=istate(k) |
|---|
| 293 | do 2400 i=1, m |
|---|
| 294 | act(i,mm1)=act(i,mm1) + a(i,j)*x(j) |
|---|
| 295 | 2400 continue |
|---|
| 296 | 2500 continue |
|---|
| 297 | c |
|---|
| 298 | c The first iteration in a warm start requires immediate qr. |
|---|
| 299 | if (loopA .eq. 1 .and. key .ne. 0) goto 6000 |
|---|
| 300 | c |
|---|
| 301 | c Steps 3, 4. |
|---|
| 302 | c Find the bound element that most wants to be active. |
|---|
| 303 | 3000 worst=0.0 |
|---|
| 304 | it=1 |
|---|
| 305 | do 3100 j=1, nbound |
|---|
| 306 | ks=abs(istate(j)) |
|---|
| 307 | bad=w(ks)*sign(1, istate(j)) |
|---|
| 308 | if (bad .lt. worst) then |
|---|
| 309 | it=j |
|---|
| 310 | worst=bad |
|---|
| 311 | iact=ks |
|---|
| 312 | endif |
|---|
| 313 | 3100 continue |
|---|
| 314 | c |
|---|
| 315 | c Test whether the Kuhn-Tucker condition is met. |
|---|
| 316 | if (worst .ge. 0.0 ) then |
|---|
| 317 | istate(n+1)=nbound |
|---|
| 318 | w(1)=sqrt(obj) |
|---|
| 319 | return |
|---|
| 320 | endif |
|---|
| 321 | c |
|---|
| 322 | c The component x(iact) is the one that most wants to become active. |
|---|
| 323 | c If the last successful change in the active set was to move x(iact) |
|---|
| 324 | c to a bound, don't let x(iact) in now: set the derivative of the |
|---|
| 325 | c misfit with respect to x(iact) to zero and return to the Kuhn-Tucker |
|---|
| 326 | c test. |
|---|
| 327 | if ( iact .eq. jj ) then |
|---|
| 328 | w(jj)=0.0 |
|---|
| 329 | goto 3000 |
|---|
| 330 | endif |
|---|
| 331 | c |
|---|
| 332 | c Step 5. |
|---|
| 333 | c Undo the effect of the new (potentially) active variable on the |
|---|
| 334 | c residual vector. |
|---|
| 335 | if (istate(it) .gt. 0) bound=bu(iact) |
|---|
| 336 | if (istate(it) .lt. 0) bound=bl(iact) |
|---|
| 337 | do 5100 i=1, m |
|---|
| 338 | act(i,mm1)=act(i,mm1) + bound*a(i,iact) |
|---|
| 339 | 5100 continue |
|---|
| 340 | c |
|---|
| 341 | c Set flag ifrom5, indicating that Step 6 was entered from Step 5. |
|---|
| 342 | c This forms the basis of a test for instability: the gradient |
|---|
| 343 | c calculation shows that x(iact) wants to join the active set; if |
|---|
| 344 | c qr puts x(iact) beyond the bound from which it came, the gradient |
|---|
| 345 | c calculation was in error and the variable should not have been |
|---|
| 346 | c introduced. |
|---|
| 347 | ifrom5=istate(it) |
|---|
| 348 | c |
|---|
| 349 | c Swap the indices (in istate) of the new active variable and the |
|---|
| 350 | c rightmost bound variable; `unbind' that location by decrementing |
|---|
| 351 | c nbound. |
|---|
| 352 | istate(it)=istate(nbound) |
|---|
| 353 | nbound=nbound - 1 |
|---|
| 354 | nact=nact + 1 |
|---|
| 355 | istate(nbound+1)=iact |
|---|
| 356 | c |
|---|
| 357 | if (mm .lt. nact) then |
|---|
| 358 | istop = -104 |
|---|
| 359 | msg = 'Too many free variables in BVLS!' |
|---|
| 360 | return |
|---|
| 361 | c stop |
|---|
| 362 | endif |
|---|
| 363 | c |
|---|
| 364 | c Step 6. |
|---|
| 365 | c Load array act with the appropriate columns of a for qr. For |
|---|
| 366 | c added stability, reverse the column ordering so that the most |
|---|
| 367 | c recent addition to the active set is in the last column. Also |
|---|
| 368 | c copy the residual vector from act(., mm1) into act(., mm1+1). |
|---|
| 369 | 6000 do 6200 i=1, m |
|---|
| 370 | act(i,mm1+1)=act(i,mm1) |
|---|
| 371 | do 6100 k=nbound+1, n |
|---|
| 372 | j=istate(k) |
|---|
| 373 | act(i,nact+1-k+nbound)=a(i,j) |
|---|
| 374 | 6100 continue |
|---|
| 375 | 6200 continue |
|---|
| 376 | c |
|---|
| 377 | call qr(m, nact, act, act(1,mm1+1), zz, resq) |
|---|
| 378 | c |
|---|
| 379 | c Test for linear dependence in qr, and for an instability that moves |
|---|
| 380 | c the variable just introduced away from the feasible region |
|---|
| 381 | c (rather than into the region or all the way through it). |
|---|
| 382 | c In either case, remove the latest vector introduced from the |
|---|
| 383 | c active set and adjust the residual vector accordingly. |
|---|
| 384 | c Set the gradient component (w(iact)) to zero and return to |
|---|
| 385 | c the Kuhn-Tucker test. |
|---|
| 386 | if (resq .lt. 0.0 |
|---|
| 387 | + .or. (ifrom5 .gt. 0 .and. zz(nact) .gt. bu(iact)) |
|---|
| 388 | + .or. (ifrom5 .lt. 0 .and. zz(nact) .lt. bl(iact))) then |
|---|
| 389 | nbound=nbound + 1 |
|---|
| 390 | istate(nbound)=istate(nbound)*sign(1.0d0, x(iact)-bu(iact)) |
|---|
| 391 | nact=nact - 1 |
|---|
| 392 | do 6500 i=1, m |
|---|
| 393 | act(i,mm1)=act(i,mm1) - x(iact)*a(i,iact) |
|---|
| 394 | 6500 continue |
|---|
| 395 | ifrom5=0 |
|---|
| 396 | w(iact)=0.0 |
|---|
| 397 | goto 3000 |
|---|
| 398 | endif |
|---|
| 399 | c |
|---|
| 400 | c If Step 6 was entered from Step 5 and we are here, a new variable |
|---|
| 401 | c has been successfully introduced into the active set; the last |
|---|
| 402 | c variable that was fixed at a bound is again permitted to become |
|---|
| 403 | c active. |
|---|
| 404 | if ( ifrom5 .ne. 0 ) jj=0 |
|---|
| 405 | ifrom5=0 |
|---|
| 406 | c |
|---|
| 407 | c Step 7. Check for strict feasibility of the new qr solution. |
|---|
| 408 | do 7100 k=1, nact |
|---|
| 409 | k1=k |
|---|
| 410 | j=istate(k+nbound) |
|---|
| 411 | if (zz(nact+1-k).lt.bl(j) .or. zz(nact+1-k).gt.bu(j)) goto 8000 |
|---|
| 412 | 7100 continue |
|---|
| 413 | do 7200 k=1, nact |
|---|
| 414 | j=istate(k+nbound) |
|---|
| 415 | x(j)=zz(nact+1-k) |
|---|
| 416 | 7200 continue |
|---|
| 417 | c New iterate is feasible; back to the top. |
|---|
| 418 | goto 15000 |
|---|
| 419 | c |
|---|
| 420 | c Steps 8, 9. |
|---|
| 421 | 8000 alpha=2.0 |
|---|
| 422 | alf=alpha |
|---|
| 423 | do 8200 k=k1, nact |
|---|
| 424 | j=istate(k+nbound) |
|---|
| 425 | if (zz(nact+1-k) .gt. bu(j)) |
|---|
| 426 | + alf=(bu(j)-x(j))/(zz(nact+1-k)-x(j)) |
|---|
| 427 | if (zz(nact+1-k) .lt. bl(j)) |
|---|
| 428 | + alf=(bl(j)-x(j))/(zz(nact+1-k)-x(j)) |
|---|
| 429 | if (alf .lt. alpha) then |
|---|
| 430 | alpha=alf |
|---|
| 431 | jj=j |
|---|
| 432 | sj=sign(1.0, zz(nact+1-k)-bl(j)) |
|---|
| 433 | endif |
|---|
| 434 | 8200 continue |
|---|
| 435 | c |
|---|
| 436 | c Step 10 |
|---|
| 437 | do 10000 k=1, nact |
|---|
| 438 | j=istate(k+nbound) |
|---|
| 439 | x(j)=x(j) + alpha*(zz(nact+1-k)-x(j)) |
|---|
| 440 | 10000 continue |
|---|
| 441 | c |
|---|
| 442 | c Step 11. |
|---|
| 443 | c Move the variable that determined alpha to the appropriate bound. |
|---|
| 444 | c (jj is its index; sj is + if zz(jj)> bu(jj), - if zz(jj)<bl(jj) ). |
|---|
| 445 | c If any other component of x is infeasible at this stage, it must |
|---|
| 446 | c be due to roundoff. Bind every infeasible component and every |
|---|
| 447 | c component at a bound to the appropriate bound. Correct the |
|---|
| 448 | c residual vector for any variables moved to bounds. Since at least |
|---|
| 449 | c one variable is removed from the active set in this step, Loop B |
|---|
| 450 | c (Steps 6-11) terminates after at most nact steps. |
|---|
| 451 | noldb=nbound |
|---|
| 452 | do 11200 k=1, nact |
|---|
| 453 | j=istate(k+noldb) |
|---|
| 454 | if (((bu(j)-x(j)) .le. 0.0) .or. |
|---|
| 455 | + (j .eq. jj .and. sj .gt. 0.0)) then |
|---|
| 456 | c Move x(j) to its upper bound. |
|---|
| 457 | x(j)=bu(j) |
|---|
| 458 | istate(k+noldb)=istate(nbound+1) |
|---|
| 459 | istate(nbound+1)=j |
|---|
| 460 | nbound=nbound+1 |
|---|
| 461 | do 11100 i=1, m |
|---|
| 462 | act(i,mm1)=act(i,mm1) - bu(j)*a(i,j) |
|---|
| 463 | 11100 continue |
|---|
| 464 | else if (((x(j)-bl(j)) .le. 0.0) .or. |
|---|
| 465 | + (j .eq. jj .and. sj .lt. 0.0)) then |
|---|
| 466 | c Move x(j) to its lower bound. |
|---|
| 467 | x(j)=bl(j) |
|---|
| 468 | istate(k+noldb)=istate(nbound+1) |
|---|
| 469 | istate(nbound+1)=-j |
|---|
| 470 | nbound=nbound+1 |
|---|
| 471 | do 11150 i=1, m |
|---|
| 472 | act(i,mm1)=act(i,mm1) - bl(j)*a(i,j) |
|---|
| 473 | 11150 continue |
|---|
| 474 | endif |
|---|
| 475 | 11200 continue |
|---|
| 476 | nact=n - nbound |
|---|
| 477 | c |
|---|
| 478 | c If there are still active variables left repeat the qr; if not, |
|---|
| 479 | c go back to the top. |
|---|
| 480 | if (nact .gt. 0 ) goto 6000 |
|---|
| 481 | |
|---|
| 482 | c |
|---|
| 483 | 15000 continue |
|---|
| 484 | c |
|---|
| 485 | msg = 'BVLS fails to converge' |
|---|
| 486 | istop = -105 |
|---|
| 487 | return |
|---|
| 488 | c stop |
|---|
| 489 | end |
|---|
| 490 | c====================================================================== |
|---|
| 491 | subroutine qr(m, n, a, b, x, resq) |
|---|
| 492 | c====================================================================== |
|---|
| 493 | implicit double precision (a-h, o-z) *doub |
|---|
| 494 | c$$$$ calls no other routines |
|---|
| 495 | c Relies on FORTRAN77 do-loop conventions! |
|---|
| 496 | c Solves over-determined least-squares problem ax ~ b |
|---|
| 497 | c where a is an m by n matrix, b is an m-vector . |
|---|
| 498 | c resq is the sum of squared residuals of optimal solution. Also used |
|---|
| 499 | c to signal error conditions - if -2 , system is underdetermined, if |
|---|
| 500 | c -1, system is singular. |
|---|
| 501 | c Method - successive Householder rotations. See Lawson & Hanson - |
|---|
| 502 | c Solving Least Squares Problems (1974). |
|---|
| 503 | c Routine will also work when m=n. |
|---|
| 504 | c***** CAUTION - a and b are overwritten by this routine. |
|---|
| 505 | dimension a(m,n),b(m),x(n) |
|---|
| 506 | double precision sum, dot |
|---|
| 507 | c |
|---|
| 508 | resq=-2.0 |
|---|
| 509 | if (m .lt. n) return |
|---|
| 510 | resq=-1.0 |
|---|
| 511 | c Loop ending on 1800 rotates a into upper triangular form. |
|---|
| 512 | do 1800 j=1, n |
|---|
| 513 | c Find constants for rotation and diagonal entry. |
|---|
| 514 | sq=0.0 |
|---|
| 515 | do 1100 i=j, m |
|---|
| 516 | sq=a(i,j)**2 + sq |
|---|
| 517 | 1100 continue |
|---|
| 518 | if (sq .eq. 0.0) return |
|---|
| 519 | qv1=-sign(sqrt(sq), a(j,j)) |
|---|
| 520 | u1=a(j,j) - qv1 |
|---|
| 521 | a(j,j)=qv1 |
|---|
| 522 | j1=j + 1 |
|---|
| 523 | c Rotate remaining columns of sub-matrix. |
|---|
| 524 | do 1400 jj=j1, n |
|---|
| 525 | dot=u1*a(j,jj) |
|---|
| 526 | do 1200 i=j1, m |
|---|
| 527 | dot=a(i,jj)*a(i,j) + dot |
|---|
| 528 | 1200 continue |
|---|
| 529 | const=dot/abs(qv1*u1) |
|---|
| 530 | do 1300 i=j1, m |
|---|
| 531 | a(i,jj)=a(i,jj) - const*a(i,j) |
|---|
| 532 | 1300 continue |
|---|
| 533 | a(j,jj)=a(j,jj) - const*u1 |
|---|
| 534 | 1400 continue |
|---|
| 535 | c Rotate b vector. |
|---|
| 536 | dot=u1*b(j) |
|---|
| 537 | do 1600 i=j1, m |
|---|
| 538 | dot=b(i)*a(i,j) + dot |
|---|
| 539 | 1600 continue |
|---|
| 540 | const=dot/abs(qv1*u1) |
|---|
| 541 | b(j)=b(j) - const*u1 |
|---|
| 542 | do 1700 i=j1, m |
|---|
| 543 | b(i)=b(i) - const*a(i,j) |
|---|
| 544 | 1700 continue |
|---|
| 545 | 1800 continue |
|---|
| 546 | c Solve triangular system by back-substitution. |
|---|
| 547 | do 2200 ii=1, n |
|---|
| 548 | i=n-ii+1 |
|---|
| 549 | sum=b(i) |
|---|
| 550 | do 2100 j=i+1, n |
|---|
| 551 | sum=sum - a(i,j)*x(j) |
|---|
| 552 | 2100 continue |
|---|
| 553 | if (a(i,i).eq. 0.0) return |
|---|
| 554 | x(i)=sum/a(i,i) |
|---|
| 555 | 2200 continue |
|---|
| 556 | c Find residual in overdetermined case. |
|---|
| 557 | resq=0.0 |
|---|
| 558 | do 2300 i=n+1, m |
|---|
| 559 | resq=b(i)**2 + resq |
|---|
| 560 | 2300 continue |
|---|
| 561 | return |
|---|
| 562 | end qr |
|---|
| 563 | c______________________________________________________________________ |
|---|
| 564 | |
|---|
