|
|
|
|
|
;+ ; 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[0] 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] ; Adjust Z() as though X(J) had been reset to zero. 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[1] & n = sz[2] ; Check array sizes for consistency and with M and N. if (size(BBV))[0] 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))[1] 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))[1],2) return endif X = make_array(n, TYPE=sz[3]) & W = X & S = X index = lindgen(n) ; Initialize the array index(). Z = make_array(m, TYPE=sz[3]) 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