 # document.write(pagetitle) Home Download Tutorial 1 Tutorial 2 Tutorial 3 Tutorial 4 Tutorial 5
```;+
; NAME:              BVLS_LH
;
; PURPOSE:           Bounded value least-square minimization (BVLS)
;
; USAGE:
;          bvls_lh, a, b, bnd, x, STATUS=status, /DOUBLE, /QUIET
;                   RNORM=rnorm,
;                   NSETP=nsetp,
;                   W=w,
;                   INDEX=index,
;                   ITMAX=itmax,
;                   EPS=eps
;
; ARGUMENTS:
;       A
;          On entry A contains the M by N matrix, A.
;          On return A contains the product matrix, Q*A, where Q is an M by M
;          orthogonal matrix generated by this subroutine.
;       B
;          On entry B contains the M-vector, B.
;          On return, B contains Q*B.  The same Q multiplies A.
;       BND
;          Input, unchanged on return
;          BND[0,J] is the lower bound for X[J].
;          BND[1,J] is the upper bound for X[J].
;          Require:  BND[0,J]  <=  BND[1,J].
;       X
;          On return, X will contain the solution N-vector.
;
; KEYWORDS:   [Optional input]
;       ITMAX=itmax,
;          Set to 3*N.  Maximum number of iterations permitted. This is usually
;          larger than required.
;       EPS=eps,
;          Determines the relative linear dependence of a column vector for a
;          variable moved from its initial value.  The default value is
;          EPS=(machar()).eps (single precision).  Other values, larger or
;          smaller may be needed for some problems.
;       /DOUBLE,
;          Setting this keyword forces some computations to be performed in
;          double precision. This may be required if: (i) N is large (>10^6)
;          or (ii) if the A or B values span a large amplitude (e.g. 10^20).
;          This slows the computation, and is not required in most cases.
;       /QUIET
;          Verbosity control. By default the procedure prints the messages.
;
; KEYWORDS:   [Optional output]
;       RNORM=rnorm,
;         Euclidean norm of the residual vector, b - A*X.
;       NSETP=nsetp,
;         Indicates the number of components of the solution vector, X, that
;         are not at their constraint values.
;       W=w,
;         An N-array.  On return, W() will contain the dual solution vector.
;         Using the following definitions:
;         W[J] = 0 for all j in Set P (i.e. X not at a limit),
;         W[J]  <=  0 for all j in Set Z, such that X[J] is at its lower bound,
;         W[J]  >=  0 for all j in Set Z, such that X[J] is at its upper bound.
;         If BND[0,J] = BND[1,J], so the variable X[J] is fixed, then W[J] will
;         have an arbitrary value.
;       INDEX=index,
;         An INTEGER N-array.  On exit the contents of this array define
;         the sets P, Z, and F as follows:
;            INDEX   through INDEX[NSETP-1] =  Set P.
;            INDEX[NSETP] through INDEX[IZ2]   = Set Z.
;            INDEX[IZ2+1] through INDEX[N-1]   = Set F.
;         Any of these sets may be empty.  Set F is those components that are
;         constrained to a unique value by the given constraints.
;         Sets P and Z are those that are allowed a non-zero range of values.
;         Of these, set Z are those whose final value is a constraint value,
;         while set P are those whose final value is not a constraint.
;         The value of IZ2 is computable as the number of bounds constraining
;         a component of X uniquely.
;       STATUS=status,
;         Indicates status on return.
;            0   Solution completed.
;            1   M  <=  0 or N  <=  0
;            2   B, X, BND, W, or INDEX size or shape violation.
;            3   Input bounds are inconsistent.
;            4   Exceed maximum number of iterations.
;
; DESCRIPTION:
;   Given an M by N matrix, A, and an M-vector, B,  compute an N-vector, X,
;   that solves the least-squares problem A *X = B
;   subject to X[J] satisfying  BND[0,J]  <=  X[J]  <=  BND[1,J]
;
;   Values BND[0,J] <= -(machar()).xmax and BND[1,J] >= (machar()).xmax
;   designate that there is no constraint in that direction.
;   Note that this flags are the single precision largest values, therefore
;   double precision bounds values exceeding this range will be taken
;   as no bound.
;   Non finite values in BND are also interpreted as absence of
;   constraint in the corresponding direction.
;
; AUTHOR: R. J. Hanson (1995 April)
;
; HISTORY:
;   This algorithm is a generalization of  NNLS, that solves the least-squares
;   problem,  A * X = B,  subject to all X[J]  >=  0.
;   The subroutine NNLS appeared in 'SOLVING LEAST SQUARES PROBLEMS', by
;   Lawson and Hanson, Prentice-Hall, 1974.  Work on BVLS was started by
;   C. L. Lawson and R. J. Hanson at Jet Propulsion Laboratory, 1973 June 12.
;   Many modifications were subsequently made.
;   The Fortran 90 code was completed in April, 1995 by R. J. Hanson.
;   The BVLS package is an additional item for the reprinting of the book by
;   SIAM Publications and the distribution of the code package using netlib
;   and Internet or network facilities.
;
;   The GDL/IDL transcoding was completed by Philippe Prugniel on
;   2010 November 10th. A similar transcoding by M. Cappellari was used
;   in earlier releases of the ULySS package. There are some differences in the
;   arguments of both versions, and the present one was found slightly faster.
;
;   See also the implementation of BVSL by R. Parker and P. Stark. It allows
;   a "warm" start and may be more efficient than the present one.
;-
; CATEGORY:    ULY_UTIL
;------------------------------------------------------------------------------

function lh_nrm2, X, DOUBLE=double
compile_opt idl2, hidden

if keyword_set(double) then return, sqrt(total(double(X)^2)) \$
else return, sqrt(total(X^2))

end

pro lh_termination
compile_opt idl2, hidden

common bvls_common, \$
FIND, HITBND, FREE1, FREE2, FREE, \$
IERR, M, N, \$
I, IBOUND, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, \$
J, JJ, JZ, L, LBOUND, NPP1, NSETP, \$
INDEX, \$
ABV, BBV, S, X, W, Z, BNDBV, \$
ALPHA, ASAVE, CC, EPS, RANGE, RNORM, \$
NORM, SM, SS, T, UNORM, UP, ZTEST

IF IERR le 0 then begin

;  Compute the norm of the residual vector.
SM = 0.
IF NPP1 le M then  SM = lh_NRM2(BBV[NPP1-1:*]) \$
else  W[0:N-1] = 0.
RNORM = SM

ENDIF                         ; ( IERR...)
END                             ; ( TERMINATION )

pro lh_move_coef_i_from_set_p_to_set_z
compile_opt idl2, hidden

common bvls_common

X[I-1] = BNDBV[IBOUND-1,I-1]
IF  abs(X[I-1]) gt 0 and JJ gt 0 then BBV[0:JJ-1] -= ABV[0:JJ-1,I-1] * X[I-1]

;   The following loop can be null.
for J = JJ, NSETP-1 do begin
II = INDEX[J]
INDEX[J-1] = II

;  The plane rotation is applied to two rows of ABV and the right-hand
;  side.  One row is moved to the scratch array S and then the updates
;  are computed.  The intent is for array operations to be performed
;  and minimal extra data movement.  One extra rotation is applied
;  to column II in this approach.
SCALE = ABS(ABV[J-1,II]) + ABS(ABV[J,II])
IF SCALE gt 0. THEN begin
ROE = ABV[J,II]
IF ABS(ABV[J-1,II]) GT ABS(ABV[J,II]) then ROE = ABV[J-1,II]
R = SCALE*SQRT((ABV[J-1,II]/SCALE)^2 + (ABV[J,II]/SCALE)^2)
IF ROE lt 0. then R = -R
CC = ABV[J-1,II]/R
SS = ABV[J,II]/R
ABV[J-1,II] = R
SM = ABV[J-1,II]
S = ABV[J-1,0:N-1]
ABV[J-1,0:*] = CC*S + SS*ABV[J,0:*]
ABV[J,0:*] = CC*ABV[J,0:*] - SS*S
ABV[J-1,II] = SM

SM = BBV[J-1]
BBV[J-1] = CC*SM + SS*BBV[J]
BBV[J] = CC*BBV[J] - SS*SM
endif

ABV[J,II] = 0.

Endfor

NPP1 = NSETP
NSETP --
IZ1 = IZ1-1
INDEX[IZ1] = I - 1
END; ( MOVE COEF I FROM SET P TO SET Z )

pro lh_see_if_all_constrained_coeffs_are_feasible
compile_opt idl2, hidden

common bvls_common

;  See if each coefficient in set P is strictly interior to its constraint region.
;  If so, set HITBND = false.
;  If not, set HITBND = true, and also set ALPHA, JJ, and IBOUND.
;  Then ALPHA will satisfy  0.  < ALPHA  <=  1.

ALPHA = 2.
for IP=0,NSETP-1 do begin
L = INDEX[IP]
IF Z[IP] le BNDBV[0,L] then \$ ;   Z(IP) HITS LOWER BOUND
LBOUND = 1   \$
ELSE IF  Z[IP] ge BNDBV[1,L] then \$ ; Z(IP) HITS UPPER BOUND
LBOUND = 2   \$
else  LBOUND = 0

IF LBOUND ne 0 then begin
T = (BNDBV[LBOUND-1,L]-X[L]) / (Z[IP] -X[L])
IF  ALPHA gt T then begin
ALPHA = T
JJ = IP+1
IBOUND = LBOUND
ENDIF
ENDIF
ENDfor
HITBND = abs(ALPHA   -  2) gt 0.
END ; SEE IF ALL CONSTRAINED COEFFS ARE FEASIBLE

pro lh_test_set_p_against_constraints
compile_opt idl2, hidden

common bvls_common

while 1 do begin
;  The solution obtained by solving the current set P is in the array Z().

ITER ++
if  ITER gt ITMAX then begin
IERR = 4
if not keyword_set(quiet) then \$
message, /INFO, 'Maximum number of iterations reached', itmax
break
endif

lh_SEE_IF_ALL_CONSTRAINED_COEFFS_ARE_FEASIBLE

;  The above call sets HITBND.  If HITBND = true then it also sets
;  ALPHA, JJ, and IBOUND.
IF  not HITBND then break

;  Here ALPHA will be between 0 and 1 for interpolation
;  between the old X() and the new Z().
if nsetp gt 0 then begin
l = index[0:nsetp-1]
x[l] += alpha * (z[0:nsetp-1] - x[l])
endif

I = INDEX[JJ-1] + 1
;  Note:  The exit test is done at the end of the loop, so the loop
;  will always be executed at least once.
while 1 DO begin

;  Modify A(*,*), B(*) and the index arrays to move coefficient I
;  from set P to set Z.

lh_MOVE_COEF_I_FROM_SET_P_TO_SET_Z

IF  NSETP le 0 then break

;  See if the remaining coefficients in set P are feasible.  They should
;  be because of the way ALPHA was determined.  If any are infeasible
;  it is due to round-off error.  Any that are infeasible or on a boundary
;  will be set to the boundary value and moved from set P to set Z.

IBOUND = 0
for JJ=0L,NSETP-1 do begin
I = INDEX[JJ]
IF  X[I] le  BNDBV[0,I] then begin
IBOUND = 1
break
endif ELSE IF X[I] ge BNDBV[1,I] then begin
IBOUND = 2
break
ENDIF
ENDfor
IF  IBOUND le  0 then break
ENDwhile

;  Copy BBV( ) into Z( ).  Then solve again and loop back.
Z = BBV

for I=NSETP-1,0, -1 do begin
IF  I ne NSETP-1 then Z[0:I] -= ABV[0:I,II]*Z[I+1]
II = INDEX[I]
Z[I] /= ABV[I,II]
Endfor

ENDwhile

if nsetp ge 1 then begin
i = index[0:nsetp-1]
x[i] = z[0:nsetp-1]
endif

;
END; ( TEST SET P AGAINST CONSTRAINTS)

pro lh_move_j_from_set_z_to_set_p
compile_opt idl2, hidden

common bvls_common

;  The index  J=index(IZ)  has been selected to be moved from
;  set Z to set P.  Z() contains the old BBV() adjusted as though X(J) = 0.
;  A(*,J) contains the new Householder transformation vector.
BBV = Z

INDEX[IZ] = INDEX[IZ1]
INDEX[IZ1] = J
IZ1 = IZ1+1
NSETP = NPP1
NPP1 = NPP1 + 1
;  The following loop can be null or not required.
NORM = ABV[NSETP-1,J]
ABV[NSETP-1,J] = UP
IF ABS(NORM) gt 0 THEN begin
for JZ=IZ1,IZ2 do begin
JJ = INDEX[JZ]
SM = total(ABV[NSETP-1:*,J]/NORM * ABV[NSETP-1:*,JJ],/DOUBLE)/UP
ABV[NSETP-1:*,JJ] += SM*ABV[NSETP-1:*,J]
endfor
ABV[NSETP-1,J] = NORM
ENDIF

if npp1 le m then ABV[npp1-1:*, j] = 0

W[J] = 0

;  Solve the triangular system.  Store this solution temporarily in Z().
for I = NSETP-1, 0, -1 do begin
IF I ne  NSETP-1 then  Z[0:I] -= ABV[0:I,II] * Z[I+1]
II = INDEX[I]
Z[I] /= ABV[I,II]
endfor

END; ( MOVE J FROM SET Z TO SET P )

pro lh_test_coef_j_for_diag_elt_and_direction_of_change
compile_opt idl2, hidden

common bvls_common

;  The sign of W(J) is OK for J to be moved to set P.
;  Begin the transformation and check new diagonal element to avoid
;  near linear dependence.

ASAVE = ABV[NPP1-1,J]

;  Construct a Householder transormation.(routine HTC in the F90 version)
VNORM = lh_NRM2(abv[NPP1-1:*,j])
IF abv[NPP1-1,j] gt 0 then VNORM *= -1.
UP = abv[NPP1-1,j] - VNORM
abv[NPP1-1,j] = VNORM

IF NSETP LT 1 then UNORM = 0.0 ELSE UNORM = lh_NRM2(ABV[0:NSETP-1,J])
IF  abs(ABV[NPP1-1,J]) gt EPS * UNORM then begin

;  Column J is sufficiently independent.  Copy b into Z, update Z.
Z = BBV
; Compute product of transormation and updated right-hand side.
NORM = ABV[NPP1-1,J]
ABV[NPP1-1,J] = UP
IF ABS(NORM) gt 0 THEN begin
sm = total(ABV[NPP1-1:*,J]/NORM * Z[NPP1-1:*]) / UP
Z[NPP1-1:*] += SM*ABV[NPP1-1:*,J]
ABV[NPP1-1,J] = NORM
ENDIF

IF abs(X[J]) gt 0 then Z[0:NPP1-1] = Z[0:NPP1-1] + ABV[0:NPP1-1,J] * X[J]
IF FREE then  FIND = 1 \$
else begin

;  Solve for ZTEST ( proposed new value for X(J) ).
;  Then set FIND to indicate if ZTEST has moved away from X(J) in
;  the expected direction indicated by the sign of W(J).
ZTEST = Z[NPP1-1] / ABV[NPP1-1,J]
FIND = ( W[J] lt 0 and  ZTEST lt  X[j] ) or \$
( W[j] gt  0 and ZTEST gt  X[J] )
endelse
ENDIF

;  If J was not accepted to be moved from set Z to set P,
;  restore A(NNP1,J).  Failing these tests may mean the computed
;  sign of W(J) is suspect, so here we set W(J) = 0.  This will
;  not affect subsequent computation, but cleans up the W() array.
IF  not FIND then begin
ABV[NPP1-1,J] = ASAVE
W[J] = 0
ENDIF                      ; ( .not. FIND )
END; TEST_COEF_J_FOR_DIAG_ELT_AND_DIRECTION_OF_CHANGE

pro lh_select_another_coeff_to_solve_for
compile_opt idl2, hidden

common bvls_common

;  1. Search through set z for a new coefficient to solve for.
;     First select a candidate that is either an unconstrained
;     coefficient or else a constrained coefficient that has room
;     to move in the direction consistent with the sign of its dual
;     vector component.  Components of the dual (negative gradient)
;     vector will be computed as needed.
;  2. For each candidate start the transformation to bring this
;     candidate into the triangle, and then do two tests:  Test size
;     of new diagonal value to avoid extreme ill-conditioning, and
;     the value of this new coefficient to be sure it moved in the
;     expected direction.
;  3. If some coefficient passes all these conditions, set FIND = true,
;     The index of the selected coefficient is J = INDEX(IZ).
;  4. If no coefficient is selected, set FIND = false.

FIND = 0

for iz=iz1,iz2 do begin
J = INDEX[IZ]

;  Set FREE1 = true if X(J) is not at the left end-point of its
;  constraint region.
FREE1 = X[J] gt  BNDBV[0,J]
;  Set FREE2 = true if X(J) is not at the right end-point of its
;  constraint region.
FREE2 = X[J] lt  BNDBV[1,J]
;  Set FREE = true if X(J) is not at either end-point of its
;  constraint region.
FREE = FREE1 and FREE2

IF  FREE  then \$
lh_TEST_COEF_J_FOR_DIAG_ELT_AND_DIRECTION_OF_CHANGE \$
else begin
;  Compute dual coefficient W(J).
W[J] = total(ABV[NPP1-1:*,J] * BBV[NPP1-1:*])

;  Can X(J) move in the direction indicated by the sign of W(J)?
if  W[J] lt 0 and  FREE1 then \$
lh_TEST_COEF_J_FOR_DIAG_ELT_AND_DIRECTION_OF_CHANGE \$
else if W[J] gt 0 and FREE2 then \$
lh_TEST_COEF_J_FOR_DIAG_ELT_AND_DIRECTION_OF_CHANGE
endelse

IF  FIND then RETURN
Endfor;  iz
END    ; ( SELECT ANOTHER COEF TO SOLVE FOR )
;

pro lh_initialize, QUIET=quiet, ITMAX=itmaxl, EPS=epsl
compile_opt idl2, hidden

common bvls_common

HITBND=0
i=0 & IBOUND=0 & ii=0 & ip=0
jj=0 & jz=0 & l=0 & lbound=0
alpha=0 & asave=0 & cc=0 & rnorm=0
norm=0 & ss=0 & t=0 & unorm=0 & up=0 & ztest=0

if n_elements(abv) eq  0 then begin
ierr = 1
if not keyword_set(quiet) then \$
message, /INOF, 'The A matrix has no element'
return
endif

sz = size(abv)
m = sz & n = sz

; Check array sizes for consistency and with M and N.
if (size(BBV)) ne 1 then begin
ierr=2
if not keyword_set(quiet) then \$
message, /INFO, 'The B array must be a vector '
return
endif

IF (SIZE(BBV)) ne M THEN begin
IERR=2
if not keyword_set(quiet) then \$
message, 'The B vector must have '+strtrim(m,2)+' elements ... it has '+strtrim((SIZE(BBV)),2)
return
endif

X = make_array(n, TYPE=sz) & W = X & S = X
index = lindgen(n)              ;   Initialize the array index().
Z = make_array(m, TYPE=sz)

IERR = 0

huge = (machar()).xmax
if n_elements(epsl) eq 0 then epsl = (machar()).eps
eps = epsl                      ; load in the common

if n_elements(itmaxl) eq 0 then itmaxl = 3*N
itmax = itmaxl                  ; load in the common
iter=0

IZ1 = 0
IZ2 = N-1
NSETP = 0
NPP1 = 1

nl = where(finite(BNDBV[0,*]) eq 0, cnt)
if cnt gt 0 then BNDBV[0,nl] = -huge
nl = where(finite(BNDBV[1,*]) eq 0, cnt)
if cnt gt 0 then BNDBV[1,nl] = huge

;   Begin:  Loop on IZ to initialize  X().
IZ=IZ1
while 1 do begin
IF  IZ  gt  IZ2  then break
J = INDEX[IZ]
IF  BNDBV[0,J] le -huge then begin
IF  BNDBV[1,J] ge huge then X[J] = 0. \$
else X[j] = 0. < BNDBV[1,j]
endif ELSE IF BNDBV[1,J] ge huge then  X[j] = 0 > BNDBV[0,J] \$
else begin
RANGE = BNDBV[1,J] - BNDBV[0,J]
IF  RANGE le 0. then begin
;   Here X(J) is constrained to a single value.
INDEX[IZ] = INDEX[IZ2]
INDEX[IZ2] = J + 1
IZ --
IZ2 --
X[J] = BNDBV[0,J]
W[J] = 0.
endif else IF  RANGE  gt  0 then  \$
;
;   The following statement sets X(J) to 0 if the constraint interval
;   includes 0, and otherwise sets X(J) to the endpoint of the
;   constraint interval that is closest to 0.

X[J] = BNDBV[0,J] > (BNDBV[1,J] < 0.) \$
else begin
ierr = 3
if not keyword_set(quiet) then \$
message, /INFO, 'Input bounds are inconsistent'
return
endelse;  RANGE
endelse

IF  abs(X[J]) gt 0 then begin
;   Change B() to reflect a nonzero starting value for X(J).
BBV -= ABV[0:*,J] * X[J]
ENDIF
IZ ++
endwhile; ( IZ   <=   IZ2 )

end ; INITIALIZE

pro bvls_lh, A, B, BND, X1, \$
RNORM=RNORM1,   \$
NSETP=NSETP1,   \$
W=W1,           \$
INDEX=INDEX1,   \$
ITMAX=itmaxl,   \$
EPS=epsl,       \$
STATUS=status,  \$
QUIET=quiet

; Another variable that we may consider to define as an argument:
;   ITER   [integer] Iteration counter.

;   Z()     [Scratch]  An M-array of working space.
;   S()     [Scratch]  An N-array of working space.

COMMON bvls_common

abv = temporary(a)  ; copy the argument A into the common
bbv = temporary(b)  ; copy the argument B into the common
bndbv = bnd         ; copy the argument BND into the common

lh_initialize, QUIET=quiet, ITMAX=itmaxl, EPS=epsl  ;   This will set IERR.

while 1 do begin  ; LOOPA
;   Quit on error flag, or if all coefficients are already in the
;   solution, .or. if M columns of A have been triangularized.
IF  IERR  ne  0  or  IZ1  gt IZ2 or NSETP ge M then break

lh_select_another_coeff_to_solve_for

;  See if no index was found to be moved from set Z to set P.
;  Then go to termination.
IF  not FIND then break

lh_move_j_from_set_z_to_set_p

lh_test_set_p_against_constraints;  This call may set IERR.

;!   All coefficients in set P are strictly feasible.  Loop back.
endwhile

lh_termination

a = temporary(abv)  ; copy to the output argument
b = temporary(bbv)  ; copy to the output argument

x1 = temporary(x)  ; copy to the output argument
rnorm1 = temporary(rnorm)  ; copy to the output argument
nsetp1 = temporary(nsetp)  ; copy to the output argument
w1 = temporary(w)  ; copy to the output argument
index1 = temporary(index)  ; copy to the output argument

status = ierr

RETURN
END; BVLS

pro prog7

; DEMONSTRATE THE USE OF THE SUBROUTINE BVLS  FOR LEAST
; SQUARES SOLVING WITH BOUNDS ON THE VARIABLES.
;
; The original version of this code was developed by
; Charles L. Lawson and Richard J. Hanson and published in the book
; "SOLVING LEAST SQUARES PROBLEMS." REVISED APRIL, 1995 to accompany
; reprinting of the book by SIAM.

;    Test driver for BVLS.  Bounded Variables Least Squares.
;    C.L.Lawson, & R.J.Hanson, Jet Propulsion Laboratory, 1973 June 12
;    Changes made in September 1982, June 1986, October 1987.
;    Conversion made to Fortran 90 by R. J. Hanson, April 1995.
;    ------------------------------------------------------------------
;
mm = 10
nn = 10
mxcase = 6
jstep = 5

MTAB = [2, 5000000, 4,  5, 10, 6]
NTAB = [2, 4, 2, 10,  5, 4]
UNBTAB = [1.0E6, 1.0E6, 1.0E6, 1.0E6, 1.0E6, 999.0E0]

bndtab = fltarr (2,NN,MXCASE)
BNDTAB[0:1,0:1,0] = [[1.,2.],[3.,4.]]
BNDTAB[0:1,0:3,1] = [[0,10], [0,10], [0,10], [0,10]]
BNDTAB[0:1,0:1,2] = [ 0,100,-100,100]
BNDTAB[0:1,0:9,3] = [[0,0],[-.3994E0,-.3994E0],[-1,1],[-.3E0,-.2E0],[21,22], \$
[-4,-3], [45,46], [100,101], [1.E6,1.E6], [-1,1]]
BNDTAB[0:1,0:4,4] = [[0,1], [-1,0], [0,1], [.3E0,.4E0], [.048E0,.049E0]]
BNDTAB[0:1,0:3,5] = [[-100.,100.],  [999.,999.],   [999.,999.],  [999.,999.]]

;!     ------------------------------------------------------------------

print, \$
' TESTBVLS.. Test driver for BVLS,';,\$
;     ' Bounded Variables Least Squares.',\$
;     '          If the algorithm succeeds the solution vector, X(),',\$
;     '           and the dual vector, W(),',\$
;     '           should be related as follows:',\$
;     '        X(i) not at a bound       =>  W(i) = 0',\$
;     '        X(i) at its lower bound  =>  W(i) .le. 0',\$
;     '        X(i) at its upper bound  =>  W(i) .ge. 0',\$
;     ' except that if an upper bound and lower bound are equal then',\$
;     ' the corresponding X(i) must take that value and W(i) may have',\$
;     ' any value.'

for ICASE = 0,MXCASE-1 do begin
M = MTAB[ICASE]
N = NTAB[ICASE]
UNBND = UNBTAB[ICASE]
BND = BNDTAB[0:1,0:n-1,ICASE]

huge = (machar()).xmax

for nn=0,n-1 do begin
if bnd[0,nn] eq UNBND then bnd[0,nn] = -huge
if bnd[1,nn] eq UNBND then bnd[1,nn] = huge
endfor

print, 'Case No.',ICASE, '   M =',M,',   N =',N,',   UNBND =',UNBND
;     print, ' Bounds ='
;        print, BND[0,0:*]
;        print, BND[1,0:*]

b = randomu(1, m)
A = randomu(2, m, n)

;     print, 'A=', A[0:M-1,0:N-1]
;     print, 'B=', b[0:m-1]

B2 = B[0:M-1]
A2 = A[0:M-1,0:N-1]
bvls_lh, A2[0:M-1,0:N-1], B2[0:m-1],BND,X,RNORM=rnorm,NSETP=nsetp,W=w,INDEX=index, STATUS=ierr
print, 'X=', x
print, 'No. of components not at constraints =',NSETP
R = B[0:M-1] - A[0:M-1,0:N-1] # X[0:N-1]

RNORM2 = sqrt(total(R[0:M-1] * R[0:M-1]))

D = R[0:M-1] # A[0:M-1,0:N-1]
;     print, ' R = B - A*X Computed by the driver:', r[0:m-1]
print, 'RNORM2 computed by the driver =',RNORM2
print, 'RNORM computed by BVLS       = ',RNORM

;    B2 = B[0:M-1]
;    A2 = A[0:M-1,0:N-1]
;    t  = systime(1)
;    x = mregress(a2, b2, MEAS=b2*0+1)
;    print, 'mregress X= ', reform(x,n_elements(x))
;    print<, 'mregress time= ', systime(1)-t
;    print, ' '

Endfor                        ; ICASE
end                             ; program PROG7
```
Part of ULySS package