subroutine bvls ( m, n, a, b, bnd, x, rnorm, nsetp, w, index, & ierr ) c*********************************************************************72 c cc bvls() solves a least squares problem with bounds on the variables. c c Discussion: c c Given an M by N matrix, A, and an M-vector, B, compute an c N-vector, X, that solves the least-squares problem c A * X = B c subject to the lower and upper bounds: c BND(1,1:N) .le. X(J) .le. BND(2,1:N) c c In cases where no bound is intended, simply set the corresponding entry c of BND to a very negative or very positive value. c c This algorithm is a generalization of NNLS, which solves c the least-squares problem, c A * X = B c subject to the constraints c 0 .le. X(1:N). c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 22 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer M, N, the number of rows and columns in A. c c Input/output, double precision A(M,N). c On input, A() contains the M by N matrix, A. c On output, A() contains the product matrix, Q*A, where c Q is an M by M orthogonal matrix generated by this c subroutine. The dimensions are M=size(A,1) and N=size(A,2). c c Input/output, double precision B(M). c On input, B() contains the M-vector, B. c On output, B() contains Q*B. The same Q multiplies A. c c Input, double precision BND(2,N), lower and upper bounds for X. c BND(1,J) .le. X(J) .le. BND(2,J). c c Output, double precision X(N), the solution vector. c c Output, double precision RNORM, the Euclidean norm of the residual c vector, b - A*X. c c Output, integer NSETP, the number of components of the c solution vector, X(), that are not at their constraint values. c c Output, double precision W(N), the dual solution vector. c Using Set definitions below: c * W(J) = 0 for all j in Set P, c * W(J) .le. 0 for all j in Set Z, such that X(J) is at its c lower bound, and c * W(J) >= 0 for all j in Set Z, such that X(J) is at its c upper bound. c If BND(1,J) = BND(2,J), so the variable X(J) is fixed, c then W(J) will have an arbitrary value. c c Output, integer INDEX(N), defines the sets P, Z, and F: c * INDEX(1) through INDEX(NSETP) = Set P. c * INDEX(IZ1) through INDEX(IZ2) = Set Z. c * INDEX(IZ2+1) through INDEX(N) = Set F. c * IZ1 = NSETP + 1 = NPP1 c Any of these sets may be empty. Set F is those components c that are constrained to a unique value by the given c constraints. Sets P and Z are those that are allowed a non-zero c range of values. Of these, set Z are those whose final c value is a constraint value, while set P are those whose c final value is not a constraint. The value of IZ2 is not returned. c It is computable as the number of bounds constraining a component c of X uniquely. c c Output, integer IERR, status flag. c * 0, Solution completed. c * 1, M .le. 0 or N .le. 0 c * 2, B(:), X(:), BND(:,:), W(:), or INDEX(:) size or shape violation. c * 3, Input bounds are inconsistent. c * 4, Exceed maximum number of iterations. c c Local Parameters: c c ITMAX [integer] c Set to 3*N. Maximum number of iterations permitted. c This is usually larger than required. Library software will c likely make this an optional argument. c c ITER [integer] c Iteration counter. c implicit none integer m integer n double precision a(m,n) double precision b(m) double precision bnd(2,n) logical find logical free integer ierr integer index(n) integer iter integer itmax integer iz integer iz1 integer iz2 integer j integer npp1 integer nsetp double precision rnorm double precision up double precision w(n) double precision x(n) double precision z(m) call initialize ( m, n, a, b, bnd, ierr, index, iter, & itmax, iz, iz1, iz2, j, npp1, nsetp, w, x ) 10 continue c c Quit (failure) on error flag. c if ( ierr .ne. 0 ) then go to 20 end if c c Quit (success) if all coefficients are already in the solution c if ( iz2 .lt. iz1 ) then go to 20 end if c c Quit (success) if M columns of A have been triangularized. c if ( m .le. nsetp ) then go to 20 end if call select_another_coeff ( m, n, iz, iz1, iz2, j, npp1, & nsetp, free, up, b, bnd, index, x, a, w, find, z ) c c Quit if no index was found to move from set Z to set P. c if ( .not. find ) then go to 20 end if call move_index_j ( m, n, iz, iz2, j, up, iz1, npp1, a, index, & w, z, nsetp, b ) call test_set_p_against_constraints ( m, n, a, b, bnd, index, & iter, itmax, iz1, npp1, nsetp, x, z, ierr ) go to 10 20 continue call termination ( ierr, m, n, npp1, b, w, rnorm ) return end subroutine bvls_report ( m, n, a, b, bnd, x, rnorm, nsetp, w, & index, ierr ) c*********************************************************************72 c cc BVLS_REPORT reports on the results of a successful call to BVLS. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 25 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer M, N, the number of rows and columns in A. c c Input, double precision A(M,N), the matrix. c c Input, double precision B(M), the right-hand-side vector, B. c c Input, double precision BND(2,N), lower and upper bounds on X. c c Input, double precision X(N), the solution vector. c c Input, double precision RNORM, the Euclidean norm of the residual c vector, b - A*X. c c Input, integer NSETP, the number of components of the c solution vector, X(), that are not at their constraint values. c c Input, double precision W(N), the dual solution vector. c c Input, integer INDEX(N), defines the sets P, Z, and F: c * INDEX(1) through INDEX(NSETP) = Set P. c * INDEX(IZ1) through INDEX(IZ2) = Set Z. c * INDEX(IZ2+1) through INDEX(N) = Set F. c * IZ1 = NSETP + 1 = NPP1 c Any of these sets may be empty. Set F is those components c that are constrained to a unique value by the given c constraints. Sets P and Z are those that are allowed a non- c zero range of values. Of these, set Z are those whose final c value is a constraint value, while set P are those whose c final value is not a constraint. The value of IZ2 is not returned. c It is computable as the number of bounds constraining a component c of X uniquely. c c Input, integer IERR, status flag. c * 0, Solution completed. c * 1, M <= 0 or N <= 0 c * 2, B(:), X(:), BND(:,:), W(:), or INDEX(:) size or shape violation. c * 3, Input bounds are inconsistent. c * 4, Exceed maximum number of iterations. c implicit none integer m integer n double precision a(m,n) double precision b(m) double precision bnd(2,n) integer i integer ierr integer index(n) double precision nrm2 external nrm2 integer nsetp double precision r(m) double precision rnorm double precision rnorm2 double precision w(n) double precision w2(n) double precision x(n) write ( *, '(a)' ) '' write ( *, '(a)' ) 'BVLS_REPORT:' if ( ierr .ne. 0 ) then write ( *, '(a,i8)' ) ' Abnormal error flag, IERR = ', ierr return end if write ( *, '(a,i8)' ) & ' Number of components not at constraints =', nsetp c c X c write ( *, '(a)' ) '' write ( *, '(a)' ) ' Solution vector, X:' write ( *, '(a)' ) '' write ( *, '(2x,5g14.6)' ) x(1:n) c c INDEX c write ( *, '(a)' ) '' write ( *, '(a)' ) ' Variable index INDEX:' write ( *, '(a)' ) '' write ( *, '(2x,5i14)' ) index(1:n) c c R c call r8mat_mv ( m, n, a, x, r ) do i = 1, m r(i) = b(i) - r(i) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Residual R = B - A*X:' write ( *, '(a)' ) ' ' write ( *, '(2x,5g14.6)' ) r(1:m) c c ||R|| c rnorm2 = nrm2 ( m, r ) write ( *, '(a)' ) ' ' write ( *, '(a,g17.5)') ' Residual norm = ', rnorm2 write ( *, '(a,g17.5)') ' Residual norm from BVLS = ', rnorm c c W2 = A'*R c call r8mat_mtv ( m, n, a, r, w2 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Dual vector: W = (A'')*R:' write ( *, '(a)' ) ' ' write ( *, '(2x,5g14.6)' ) w2(1:n) c c W: Dual vector. c write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Dual vector from BVLS: W' write ( *, '(a)' ) ' ' write ( *, '(2x,5g14.6)' ) w(1:n) return end subroutine constrained_feasible ( nsetp, index, bnd, x, z, & hitbnd, alpha, jj, ibound ) c*********************************************************************72 c cc CONSTRAINED_FEASIBLE determines if each coefficient is interior. c c Discussion: c c See if each coefficient in set P is strictly interior to its constraint c region. c c If so, set HITBND = false. c c If not, set HITBND = true, and also set ALPHA, JJ, and IBOUND. c Then ALPHA will satisfy 0 .lt. ALPHA .le. 1. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 11 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer NSETP. c c Input, integer INDEX(N), defines the sets P, Z, and F. c c Input, double precision BND(2,N). c c Input, double precision X(N). c c Input, double precision Z(M). c c Output, logical HITBND. c c Output, double precision ALPHA. c c Output, integer JJ. c c Output, integer IBOUND. c implicit none integer nsetp double precision alpha double precision bnd(2,*) logical hitbnd integer ibound integer index(*) integer ip integer jj integer l integer lbound double precision t double precision x(*) double precision z(*) alpha = 2.0D+00 do ip = 1, nsetp l = index(ip) c c Z(IP) hits lower bound c if ( z(ip) .le. bnd(1,l) ) then lbound = 1 c c Z(IP) hits upper bound c else if ( bnd(2,l) .le. z(ip) ) then lbound = 2 c c Z(IP) is in the interior: c else lbound = 0 end if if ( lbound .ne. 0 ) then t = ( bnd(lbound,l) - x(l) ) / ( z(ip) - x(l) ) if ( t .lt. alpha ) then alpha = t jj = ip ibound = lbound end if end if end do hitbnd = ( 0.0D+00 .lt. abs ( alpha - 2.0D+00 ) ) return end subroutine htc ( p, m, u, up ) c*********************************************************************72 c cc HTC constructs a Householder transformation. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer P, an index into the U vector. c c Input, integer M, the length of U. c c Input/output, double precision U(M). c c Output, double precision UP. c implicit none integer m double precision nrm2 integer p double precision u(m) double precision up double precision vnorm vnorm = nrm2 ( m + 1 - p, u(p:m) ) if ( 0.0D+00 .lt. u(p) ) then vnorm = - vnorm end if up = u(p) - vnorm u(p) = vnorm return end subroutine initialize ( m, n, a, b, bnd, ierr, index, iter, & itmax, iz, iz1, iz2, j, npp1, nsetp, w, x ) c*********************************************************************72 c cc INITIALIZE initializes internal data. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 22 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer M, N, the dimensions of the matrix. c c Input, double precision A(M,N). c c Input, double precision B(M). c c Input, double precision BND(2,N). c c Output, integer IERR. c c Output, integer INDEX(N), defines the sets P, Z, and F. c c Output, integer ITER. c c Output, integer ITMAX. c c Output, integer IZ. c c Output, integer IZ1. c c Output, integer IZ2. c c Output, integer J. c c Output, integer NPP1. c c Output, integer NSETP. c c Output, double precision W(N). c c Output, double precision X(N). c implicit none integer m integer n double precision a(m,n) double precision b(m) double precision bnd(2,n) integer i integer ierr integer index(n) integer iter integer itmax integer iz integer iz1 integer iz2 integer j integer npp1 integer nsetp double precision r double precision r8_huge parameter ( r8_huge = 1.0D+30 ) double precision w(n) double precision x(n) ierr = 0 itmax = 3 * n iter = 0 c c Initialize the array index(). c do i = 1, n index(i) = i end do iz2 = n iz1 = 1 nsetp = 0 npp1 = 1 c c Loop on IZ to initialize X(). c iz = iz1 10 continue if ( iz2 .lt. iz ) then go to 20 end if j = index(iz) if ( bnd(1,j) .le. - r8_huge ) then if ( r8_huge .le. bnd(2,j) ) then x(j) = 0.0D+00 else x(j) = min ( 0.0D+00, bnd(2,j) ) end if else if ( r8_huge .le. bnd(2,j) ) then x(j) = max ( 0.0D+00, bnd(1,j) ) else r = bnd(2,j) - bnd(1,j) c c Here X(J) is constrained to a single value. c if ( r .le. 0.0D+00 ) then index(iz) = index(iz2) index(iz2) = j iz = iz - 1 iz2 = iz2 - 1 x(j) = bnd(1,j) w(j) = 0.0D+00 c c Sets X(J) to 0 if the constraint interval includes 0; c otherwise set X(J) to the endpoint of the constraint interval closest to 0. c else if ( 0.0D+00 .lt. r ) then x(j) = max ( bnd(1,j), min ( bnd(2,j), 0.0D+00 ) ) else ierr = 3 return end if end if c c Change B() to reflect a nonzero starting value for X(J). c if ( 0.0D+00 .lt. abs ( x(j) ) ) then do i = 1, m b(i) = b(i) - a(i,j) * x(j) end do end if iz = iz + 1 go to 10 20 continue return end subroutine move_coef_i ( m, n, i, ibound, bnd, jj, b, a, x, & nsetp, index, npp1, iz1 ) c*********************************************************************72 c cc MOVE_COEF_I moves coefficient I from set P to set Z. c c Discussion: c c Modify A(*,*), B(*) and the index arrays to move coefficient I c from set P to set Z. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer M, the number of rows of A. c c Input, integer N, the number of columns of A. c c Input, integer I, the index of the coefficient to be moved. c c Input, integer IBOUND. c c Input, double precision BND(2,N). c c Input, integer JJ. c c Input/output, double precision B(M). c c Input/output, double precision A(M,N). c c Input/output, double precision X(N). c c Input/output, integer NSETP. c c Input/output, integer INDEX(N), defines the sets P, Z, and F. c c Input/output, integer NPP1. c c Input/output, integer IZ1. c implicit none integer m integer n double precision a(m,n) double precision b(m) double precision bnd(2,n) double precision cc integer i integer i2 integer ibound integer ii integer index(n) integer iz1 integer j integer j2 integer jj integer npp1 integer nsetp double precision s(n) double precision sm double precision ss double precision x(n) x(i) = bnd(ibound,i) if ( 0.0D+00 .lt. abs ( x(i) ) .and. 0 .lt. jj ) then do i2 = 1, jj b(i2) = b(i2) - a(i2,i) * x(i) end do end if do j = jj + 1, nsetp ii = index(j) index(j-1) = ii call rotg ( a(j-1,ii), a(j,ii), cc, ss ) sm = a(j-1,ii) c c The plane rotation is applied to two rows of A and the right-hand c side. One row is moved to the scratch array S and then the updates c are computed. The intent is for array operations to be performed c and minimal extra data movement. One extra rotation is applied c to column II in this approach. c do j2 = 1, n s(j2) = a(j-1,j2) end do do j2 = 1, n a(j-1,j2) = cc * s(j2) + ss * a(j,j2) end do do j2 = 1, n a(j,j2) = cc * a(j,j2) - ss * s(j2) end do a(j-1,ii) = sm a(j,ii) = 0.0D+00 sm = b(j-1) b(j-1) = cc * sm + ss * b(j) b(j) = cc * b(j) - ss * sm end do npp1 = nsetp nsetp = nsetp - 1 iz1 = iz1 - 1 index(iz1) = i return end subroutine move_index_j ( m, n, iz, iz2, j, up, iz1, npp1, a, & index, w, z, nsetp, b ) c*********************************************************************72 c cc MOVE_INDEX_J moves the index J from set Z to set P. c c Discussion: c c The index J=index(IZ) has been selected to be moved from c set Z to set P. Z() contains the old B() adjusted as though X(J) = 0. c c A(*,J) contains the new Householder transformation vector. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 11 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer M. c c Input, integer N. c c Input, integer IZ. c c Input, integer IZ2. c c Input, integer J, the index to be moved. c c Input, double precision UP. c c Input/output, integer IZ1. c c Input/output, integer NPP1. c c Input/output, double precision A(M,N). c c Input/output, integer INDEX(N), defines the sets P, Z, and F. c c Input/output, double precision W(N). c c Input/output, double precision Z(M). c c Output, integer NSETP. c c Output, double precision B(M). c implicit none integer m integer n double precision a(m,n) double precision b(m) integer i integer i2 integer ii integer index(n) integer iz integer iz1 integer iz2 integer j integer jj integer jz integer l double precision norm integer npp1 integer nsetp double precision sm double precision up double precision w(n) double precision z(m) do i2 = 1, m b(i2) = z(i2) end do index(iz) = index(iz1) index(iz1) = j iz1 = iz1 + 1 nsetp = npp1 npp1 = npp1 + 1 norm = a(nsetp,j) a(nsetp,j) = up if ( 0.0D+00 .lt. abs ( norm ) ) then do jz = iz1, iz2 jj = index(jz) sm = 0.0D+00 do i2 = nsetp, m sm = sm + a(i2,j) * a(i2,jj) end do sm = sm / norm / up do i2 = nsetp, m a(i2,jj) = a(i2,jj) + sm * a(i2,j) end do end do a(nsetp,j) = norm end if do l = npp1, m a(l,j) = 0.0D+00 end do w(j) = 0.0D+00 c c Solve the triangular system. c Store the solution temporarily in Z(). c do i = nsetp, 1, -1 if ( i .ne. nsetp ) then do i2 = 1, i z(i2) = z(i2) - a(i2,ii) * z(i+1) end do end if ii = index(i) z(i) = z(i) / a(i,ii) end do return end function nrm2 ( n, x ) c*********************************************************************72 c cc NRM2 returns the Euclidean norm of a vector. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 11 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer N, the length of the vector. c c Input, double precision X(N), the vector. c c Output, double precision NRM2, the norm of the vector. c implicit none integer n double precision absxi integer ix double precision norm double precision nrm2 double precision scale double precision ssq double precision x(n) if ( n .lt. 1 ) then norm = 0.0D+00 else if ( n == 1 ) then norm = abs ( x(1) ) else scale = 0.0D+00 ssq = 1.0D+00 do ix = 1, n absxi = abs ( x(ix) ) if ( 0.0D+00 .lt. absxi ) then if ( scale .lt. absxi ) then ssq = 1.0D+00 + ssq * ( scale / absxi )**2 scale = absxi else ssq = ssq + ( absxi / scale )**2 end if end if end do norm = scale * sqrt ( ssq ) end if nrm2 = norm return end subroutine r8mat_copy ( m, n, a1, a2 ) c*********************************************************************72 c cc R8MAT_COPY copies an R8MAT. c c Discussion: c c An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 26 July 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, N, the order of the matrix. c c Input, double precision A1(M,N), the matrix to be copied. c c Output, double precision A2(M,N), a copy of the matrix. c implicit none integer m integer n double precision a1(m,n) double precision a2(m,n) integer i integer j do j = 1, n do i = 1, m a2(i,j) = a1(i,j) end do end do return end subroutine r8mat_mv ( m, n, a, x, y ) c*********************************************************************72 c cc R8MAT_MV multiplies a matrix times a vector. c c Discussion: c c An R8MAT is an array of R8's. c c In FORTRAN90, this operation can be more efficiently carried c out by the command c c Y(1:M) = MATMUL ( A(1:M,1:N), X(1:N) ) c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 December 2004 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, N, the number of rows and columns of the matrix. c c Input, double precision A(M,N), the M by N matrix. c c Input, double precision X(N), the vector to be multiplied by A. c c Output, double precision Y(M), the product A*X. c implicit none integer m integer n double precision a(m,n) integer i integer j double precision x(n) double precision y(m) double precision y1(m) do i = 1, m y1(i) = 0.0D+00 do j = 1, n y1(i) = y1(i) + a(i,j) * x(j) end do end do do i = 1, m y(i) = y1(i) end do return end subroutine r8mat_mtv ( m, n, a, x, y ) c*****************************************************************************80 c cc R8MAT_MTV multiplies a transposed matrix times a vector c c Discussion: c c An R8MAT is an array of R8 values. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 08 August 2013 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, N, the number of rows and columns of c the matrix. c c Input, double precision A(M,N), the M by N matrix. c c Input, double precision X(M), the vector to be multiplied by A. c c Output, double precision Y(N), the product A'*X. c implicit none integer m integer n double precision a(m,n) integer i integer j double precision x(m) double precision y(n) double precision y1(n) do i = 1, n y1(i) = 0.0D+00 do j = 1, m y1(i) = y1(i) + a(j,i) * x(j) end do end do do i = 1, n y(i) = y1(i) end do return end subroutine r8mat_uniform_01 ( m, n, seed, r ) c*********************************************************************72 c cc R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT. c c Discussion: c c An R8MAT is an array of R8's. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 11 August 2004 c c Author: c c John Burkardt c c Reference: c c Paul Bratley, Bennett Fox, Linus Schrage, c A Guide to Simulation, c Springer Verlag, pages 201-202, 1983. c c Bennett Fox, c Algorithm 647: c Implementation and Relative Efficiency of Quasirandom c Sequence Generators, c ACM Transactions on Mathematical Software, c Volume 12, Number 4, pages 362-376, 1986. c c Peter Lewis, Allen Goodman, James Miller, c A Pseudo-Random Number Generator for the System/360, c IBM Systems Journal, c Volume 8, pages 136-143, 1969. c c Parameters: c c Input, integer M, N, the number of rows and columns in the array. c c Input/output, integer SEED, the "seed" value, which should NOT be 0. c On output, SEED has been updated. c c Output, double precision R(M,N), the array of pseudorandom values. c implicit none integer m integer n integer i integer j integer k integer seed double precision r(m,n) if ( seed .eq. 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do j = 1, n do i = 1, m k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed .lt. 0 ) then seed = seed + 2147483647 end if r(i,j) = dble ( seed ) * 4.656612875D-10 end do end do return end subroutine r8vec_copy ( n, a1, a2 ) c*********************************************************************72 c cc R8VEC_COPY copies an R8VEC. c c Discussion: c c An R8VEC is a vector of R8 values. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 13 August 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the length of the vectors. c c Input, double precision A1(N), the vector to be copied. c c Output, double precision A2(N), a copy of A1. c implicit none integer n double precision a1(n) double precision a2(n) integer i do i = 1, n a2(i) = a1(i) end do return end subroutine r8vec_uniform_01 ( n, seed, r ) c*********************************************************************72 c cc R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC. c c Discussion: c c An R8VEC is a vector of R8's. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 17 July 2006 c c Author: c c John Burkardt c c Reference: c c Paul Bratley, Bennett Fox, Linus Schrage, c A Guide to Simulation, c Springer Verlag, pages 201-202, 1983. c c Bennett Fox, c Algorithm 647: c Implementation and Relative Efficiency of Quasirandom c Sequence Generators, c ACM Transactions on Mathematical Software, c Volume 12, Number 4, pages 362-376, 1986. c c Peter Lewis, Allen Goodman, James Miller, c A Pseudo-Random Number Generator for the System/360, c IBM Systems Journal, c Volume 8, pages 136-143, 1969. c c Parameters: c c Input, integer N, the number of entries in the vector. c c Input/output, integer SEED, the "seed" value, which should NOT be 0. c On output, SEED has been updated. c c Output, double precision R(N), the vector of pseudorandom values. c implicit none integer n integer i integer k integer seed double precision r(n) do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed .lt. 0 ) then seed = seed + 2147483647 end if r(i) = dble ( seed ) * 4.656612875D-10 end do return end subroutine rotg ( sa, sb, c, s ) c*********************************************************************72 c cc ROTG defines a Givens rotation. c c Discussion: c c Given values A and B, this routine computes c c SIGMA = sign ( A ) if abs ( A ) > abs ( B ) c = sign ( B ) if abs ( A ) .le. abs ( B ); c c R = SIGMA * ( A * A + B * B ); c c C = A / R if R is not 0 c = 1 if R is 0; c c S = B / R if R is not 0, c 0 if R is 0. c c The computed numbers then satisfy the equation c c ( C S ) ( A ) = ( R ) c ( -S C ) ( B ) = ( 0 ) c c The routine also computes c c Z = S if abs ( A ) > abs ( B ), c = 1 / C if abs ( A ) .le. abs ( B ) and C is not 0, c = 1 if C is 0. c c The single value Z encodes C and S, and hence the rotation: c c If Z = 1, set C = 0 and S = 1; c If abs ( Z ) .lt. 1, set C = sqrt ( 1 - Z * Z ) and S = Z; c if abs ( Z ) > 1, set C = 1/ Z and S = sqrt ( 1 - C * C ); c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 15 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input/output, double precision SA, SB. On input, SA and SB are the values c A and B. On output, SA is overwritten with R, and SB is c overwritten with Z. c c Output, double precision C, S, the cosine and sine of the c Givens rotation. c implicit none double precision c double precision r double precision roe double precision s double precision sa double precision sb double precision scale roe = sb if ( abs ( sb ) .lt. abs ( sa ) ) then roe = sa end if scale = abs ( sa ) + abs ( sb ) if ( scale .le. 0.0D+00 ) then c = 1.0D+00 s = 0.0D+00 return end if r = scale * sqrt ( ( sa / scale ) ** 2 + ( sb / scale ) ** 2 ) if ( roe .lt. 0.0D+00 ) then r = - r end if c = sa / r s = sb / r sa = r return end subroutine select_another_coeff ( m, n, iz, iz1, iz2, j, npp1, & nsetp, free, up, b, bnd, index, x, a, w, find, z ) c*********************************************************************72 c cc SELECT_ANOTHER_COEFF searches set Z for a new coefficient to solve for. c c Discussion: c c 1. Search through set z for a new coefficient to solve for. c First select a candidate that is either an unconstrained c coefficient or else a constrained coefficient that has room c to move in the direction consistent with the sign of its dual c vector component. Components of the dual (negative gradient) c vector will be computed as needed. c c 2. For each candidate start the transformation to bring this c candidate into the triangle, and then do two tests: Test size c of new diagonal value to avoid extreme ill-conditioning, and c the value of this new coefficient to be sure it moved in the c expected direction. c c 3. if some coefficient passes all these conditions, set FIND = true, c The index of the selected coefficient is J = INDEX(IZ). c c 4. if no coefficient is selected, set FIND = false. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 18 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer M, the number of rows of A. c c Input, integer N, the number of columns of A. c c Output, integer IZ. c c Input, integer IZ1, IZ2. c c Input, integer J. c c Input, integer NPP1. c c Input, integer NSETP. c c Input, logical FREE. c c Input, double precision UP. c c Input, double precision B(M). c c Input, double precision BND(2,N). c c Input, integer INDEX(N), defines the sets P, Z, and F. c c Input, double precision X(N). c c Input/output, double precision A(M,N). c c Input/output, double precision W(N). c c Output, logical FIND. c c Output, double precision Z(M). c implicit none integer m integer n double precision a(m,n) double precision b(m) double precision bnd(2,n) logical find logical free logical free1 logical free2 integer i2 integer index(n) integer iz integer iz_index integer iz1 integer iz2 integer j integer npp1 integer nsetp double precision t double precision up double precision w(n) double precision x(n) double precision z(m) find = .false. do iz_index = iz1, iz2 iz = iz_index j = index(iz) c c FREE1 = true if X(J) is not at the left end-point of its constraint region. c FREE2 = true if X(J) is not at the right end-point of its constraint region. c FREE = true if X(J) is not at either end-point of its constraint region. c free1 = ( bnd(1,j) .lt. x(j) ) free2 = ( x(j) .lt. bnd(2,j) ) free = ( free1 .and. free2 ) if ( free ) then call test_coef_j ( m, n, j, npp1, nsetp, free, up, b, & x, a, w, find, z ) c c Compute dual coefficient W(J). c else t = 0.0D+00 do i2 = npp1, m t = t + a(i2,j) * b(i2) end do w(j) = t c c Can X(J) move in the direction indicated by the sign of W(J)? c if ( w(j) .lt. 0.0D+00 ) then if ( free1 ) then call test_coef_j ( m, n, j, npp1, nsetp, free, up, & b, x, a, w, find, z ) end if else if ( 0.0D+00 .lt. w(j) ) then if ( free2 ) then call test_coef_j ( m, n, j, npp1, nsetp, free, up, & b, x, a, w, find, z ) end if end if end if if ( find ) then return end if end do return end subroutine termination ( ierr, m, n, npp1, b, w, rnorm ) c*********************************************************************72 c cc TERMINATION handles the termination of the BVLS algorithm. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 11 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer IERR, error flag. c c Input, integer M. c c Input, integer N. c c Input, integer NPP1. c c Input, double precision B(M). c c Output, double precision W(N). c c Output, double precision RNORM. c implicit none integer m integer n double precision b(m) integer i2 integer ierr integer npp1 double precision nrm2 double precision rnorm double precision w(n) c c Compute the norm of the residual vector. c rnorm = 0.0D+00 if ( ierr .le. 0 ) then if ( npp1 .le. m ) then rnorm = nrm2 ( m + 1 - npp1, b(npp1:m) ) else do i2 = 1, n w(i2) = 0.0D+00 end do end if end if return end subroutine test_coef_j ( m, n, j, npp1, nsetp, free, up, b, & x, a, w, find, z ) c*********************************************************************72 c cc TEST_COEF_J c c Discussion: c c The sign of W(J) is OK for J to be moved to set P. c c Begin the transformation and check new diagonal element to avoid c near linear dependence. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 11 June 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer M, the number of rows of A. c c Input, integer N, the number of columns of A. c c Input, integer J. c c Input, integer NPP1. c c Input, integer NSETP. c c Input, logical FREE. c c Input, double precision UP. c c Input, double precision B(M). c c Input, double precision X(N). c c Input/output, double precision A(M,N). c c Input/output, double precision W(N). c c Output, logical FIND. c c Output, double precision Z(M). c implicit none integer m integer n double precision a(m,n) double precision asave double precision b(m) double precision eps parameter ( eps = 2.220446049250313D-016 ) logical find logical free integer i2 integer j double precision norm integer npp1 double precision nrm2 integer nsetp double precision sm double precision unorm double precision up double precision w(n) double precision x(n) double precision z(m) double precision ztest asave = a(npp1,j) call htc ( npp1, m, a(1:m,j), up ) unorm = nrm2 ( nsetp, a(1:nsetp,j) ) if ( eps * unorm .lt. abs ( a(npp1,j) ) ) then c c Column J is sufficiently independent. Copy B into Z, update Z. c do i2 = 1, m z(i2) = b(i2) end do c c Compute product of transformation and updated right-hand side. c norm = a(npp1,j) a(npp1,j) = up if ( 0.0D+00 .lt. abs ( norm ) ) then sm = 0.0D+00 do i2 = npp1, m sm = sm + a(i2,j) * z(i2) end do sm = sm / norm / up do i2 = npp1, m z(i2) = z(i2) + sm * a(i2,j) end do a(npp1,j) = norm end if if ( 0.0D+00 .lt. abs ( x(j) ) ) then do i2 = 1, npp1 z(i2) = z(i2) + a(i2,j) * x(j) end do end if c c Adjust Z as though X(J) had been reset to zero. c if ( free ) then find = .true. else c c Solve for ZTEST, the proposed new value for X(J). c FIND indicates if ZTEST has moved away from X(J) in c the expected direction indicated by the sign of W(J). c ztest = z(npp1) / a(npp1,j) find = ( w(j) .lt. 0.0D+00 .and. ztest .lt. x(j) ) .or. & ( 0.0D+00 .lt. w(j) .and. x(j) .lt. ztest ) end if end if c c If J was not accepted to be moved from set Z to set P, c restore A(NNP1,J). Failing these tests may mean the computed c sign of W(J) is suspect, so here we set W(J) = 0. This will c not affect subsequent computation, but cleans up the W() array. c if ( .not. find ) then a(npp1,j) = asave w(j) = 0.0D+00 end if return end subroutine test_set_p_against_constraints ( m, n, a, b, bnd, & index, iter, itmax, iz1, npp1, nsetp, x, z, ierr ) c*********************************************************************72 c cc TEST_SET_P_AGAINST_CONSTRAINTS c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 13 August 2014 c c Author: c c Original FORTRAN90 version by Charles Lawson, Richard Hanson. c This version by John Burkardt. c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c c Parameters: c c Input, integer M, N. c c Input, double precision A(M,N), the matrix. c c Input, double precision B(M). c c Input, double precision BND(2,N). c c Input/output, integer INDEX(N), defines the sets P, Z, and F. c c Input/output, integer ITER. c c Input, integer ITMAX. c c Input/output, integer IZ1. c c Input/output, integer NPP1. c c Input/output, integer NSETP. c c Input/output, double precision X(N). c c Input/output, double precision Z(M). c c Output, integer IERR. c c Local Parameters: c c ALPHA c HITBND c I c IBOUND c II c IP c JJ c L c implicit none integer m integer n double precision a(m,n) double precision alpha double precision b(m) double precision bnd(2,n) logical hitbnd integer i integer i1 integer i2 integer ibound integer ierr integer ii integer index(n) integer ip integer iter integer itmax integer iz1 integer jj integer jjj integer l integer npp1 integer nsetp double precision x(n) double precision z(m) 10 continue c c The solution obtained by solving the current set P is in the array Z(). c iter = iter + 1 if ( itmax .lt. iter ) then ierr = 4 go to 50 end if c c Set HITBND. c If HITBND is true, also set ALPHA, JJ, and IBOUND. c call constrained_feasible ( nsetp, index, bnd, x, z, hitbnd, & alpha, jj, ibound ) if ( .not. hitbnd ) then go to 50 end if c c ALPHA will be between 0 and 1 for interpolation c between the old X and the new Z. c do ip = 1, nsetp l = index(ip) x(l) = x(l) + alpha * ( z(ip) - x(l) ) end do i = index(jj) c c The exit test is done at the end of the loop, so the loop c will always be executed at least once. c 20 continue c c Modify A, B, and the index arrays to move coefficient I c from set P to set Z. c call move_coef_i ( m, n, i, ibound, bnd, jj, b, a, x, nsetp, & index, npp1, iz1 ) if ( nsetp .le. 0 ) then go to 50 end if c c See if the remaining coefficients in set P are feasible. They should c be because of the way ALPHA was determined. If any are infeasible c it is due to round-off error. Any that are infeasible or on a boundary c will be set to the boundary value and moved from set P to set Z. c ibound = 0 jj = nsetp do jjj = 1, nsetp i = index(jjj) if ( x(i) .le. bnd(1,i) ) then ibound = 1 jj = jjj go to 30 else if ( bnd(2,i) .le. x(i) ) then ibound = 2 jj = jjj go to 30 end if end do 30 continue if ( ibound .le. 0 ) then go to 40 end if go to 20 40 continue c c Copy B into Z. c do i2 = 1, m z(i2) = b(i2) end do c c Solve and loop back again. c do i = nsetp, 1, -1 if ( i .ne. nsetp ) then do i1 = 1, i z(i1) = z(i1) - a(i1,ii) * z(i+1) end do end if ii = index(i) z(i) = z(i) / a(i,ii) end do go to 10 50 continue do ip = 1, nsetp i = index(ip) x(i) = z(ip) end do return end subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, month(m), y, h, ':', n, ':', s, '.', mm, ampm return end