program main !*****************************************************************************80 ! !! MAIN is the main program for BVLS_PRB. ! ! Discussion: ! ! This program demonstrates the use of BVLS for solving least squares ! problems with include bounds on the variables. ! ! Modified: ! ! 19 October 2008 ! ! Author: ! ! Charles Lawson, Richard Hanson ! ! Reference: ! ! Charles Lawson, Richard Hanson, ! Solving Least Squares Problems, ! SIAM, 1995, ! ISBN: 0898713560, ! LC: QA275.L38. ! implicit none interface subroutine bvls ( a, b, bnd, x, rnorm, nsetp, w, index, ierr ) real ( kind ( 1e0 ) ) a(:,:) real ( kind ( 1e0 ) ) b(:) real ( kind ( 1e0 ) ) bnd(:,:) real ( kind ( 1e0 ) ) x(:) real ( kind ( 1e0 ) ) rnorm integer nsetp real ( kind ( 1e0 ) ) w(:) integer index(:) integer ierr end subroutine end interface integer, parameter :: mm = 10 integer, parameter :: nn = 10 integer, parameter :: mxcase = 6 integer, parameter :: jstep = 5 integer i integer icase integer ierr integer index(nn) integer j integer j1 integer j2 integer m integer, dimension(mxcase) :: mtab = (/ & 2, 2, 4, 5, 10, 6 /) integer n integer nsetp integer, dimension(mxcase) :: ntab = (/ & 2, 4, 2, 10, 5, 4 /) real(kind(1e0)) unbtab(mxcase), bndtab(2,nn,mxcase) real(kind(1e0)) a(mm,nn),b(mm),x(nn),w(nn) real(kind(1e0)) a2(mm,nn),b2(mm), r(mm), d(nn),bnd(2,nn) real(kind(1e0)) rnorm, rnorm2, unbnd data unbtab / 5 * 1.0e6, 999.0e0 / data ((bndtab(i,j,1),i=1,2),j=1,2)/ 1.,2., 3.,4. / data ((bndtab(i,j,2),i=1,2),j=1,4)/ 0,10, 0,10, 0,10, 0,10/ data ((bndtab(i,j,3),i=1,2),j=1,2)/ 0,100, -100,100/ data ((bndtab(i,j,4),i=1,2),j=1,10)/& 0,0, -.3994e0,-.3994e0, -1,1, -.3e0,-.2e0, 21,22,& -4,-3, 45,46, 100,101, 1.e6,1.e6, -1,1/ data ((bndtab(i,j,5),i=1,2),j=1,5)/& 0,1, -1,0, 0,1, .3e0,.4e0, .048e0,.049e0/ data ((bndtab(i,j,6),i=1,2),j=1,4)/& -100.,100., 999.,999., 999.,999., 999.,999. / write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TESTBVLS' write ( *, '(a)' ) ' Test driver for BVLS' write ( *, '(a)' ) ' Bounded Variables Least Squares.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' If the algorithm succeeds the solution vector, X(),' write ( *, '(a)' ) ' and the dual vector, W(), should be related as follows:' write ( *, '(a)' ) ' X(i) not at a bound => W(i) = 0' write ( *, '(a)' ) ' X(i) at its lower bound => W(i) <= 0' write ( *, '(a)' ) ' X(i) at its upper bound => 0 <= W(i)' write ( *, '(a)' ) ' except that if an upper bound and lower bound are equal, then' write ( *, '(a)' ) ' the corresponding X(i) must take that value and W(i) may have' write ( *, '(a)' ) ' any value.' write ( *, '(a)' ) ' ' do icase = 1, mxcase m = mtab(icase) n = ntab(icase) unbnd = unbtab(icase) do j = 1, n bnd(1,j) = bndtab(1,j,icase) bnd(2,j) = bndtab(2,j,icase) end do where ( bnd(1,1:n) == unbnd ) bnd(1,1:n)=-huge(1e0) where ( bnd(2,1:n) == unbnd ) bnd(2,1:n)= huge(1e0) write ( *, '(a)' ) ' ' write(*,'(////1x,a,i3/1x)') 'Case No.',ICASE write(*,'(1x,a,i5,a,i5,a,g17.5)')& 'M =',M,', N =',N,', UNBND =',UNBND write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Bounds:' write ( *, '(a)' ) ' ' do j1 = 1, n, jstep j2 = min ( j1 - 1 + jstep, n ) write(*,*) ' ' write(*,1001) (bnd(1,j),j=j1,j2) write(*,1001) (bnd(2,j),j=j1,j2) end do call random_number ( harvest = b(1:m) ) call random_number ( harvest = a(1:m,1:n) ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Matrix A:' write ( *, '(a)' ) ' ' do j1 = 1, n, jstep j2 = min(j1 - 1 + jstep, n) write(*,*) ' ' do i = 1,m write(*,1001) (a(i,j),j=j1,j2) end do end do b2(1:m) = b(1:m) a2(1:m,1:n) = a(1:m,1:n) write(*,'(1x/'' B() =''/1X)') write(*,1001) (b(i),i=1,m) call bvls ( a2(1:m,1:n), b2, bnd, x, rnorm, nsetp, w, index, ierr ) if ( ierr > 0) then write(*,'(1x, "Abnormal error flag, IERR = ")') ierr stop end if write(*,'(1x/1x,a,i4)') & 'After BVLS: No. of components not at constraints =',nsetp write(*,'(1x/'' Solution vector, X() =''/1x)') write(*,1001) (x(j),j=1,n) r(1:m) = b(1:m)-matmul(a(1:m,1:n),x(1:n)) rnorm2 = sqrt(dot_product(r(1:m),r(1:m))) d(1:n) = matmul(r(1:m),a(1:m,1:n)) write(*,'(1X/'' R = B - A*X Computed by the driver:''/1X)') write(*,1001) (R(I),I=1,M) write(*,'(1x,a,g17.5)') 'RNORM2 computed by the driver =',RNORM2 write(*,'(1x,a,g17.5)') 'RNORM computed by BVLS = ',RNORM write(*,'(1X/'' W = (A**T)*R Computed by the driver:''/1X)') write(*,1001) (D(J),J=1,N) write(*,'(1X/'' Dual vector from BVLS, W() =''/1X)') write(*,1001) (W(J),J=1,N) end do 1001 format(1x,5g14.5) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BVLS_PRB:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end program main subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 August 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer ( kind = 4 ) d integer ( kind = 4 ) h integer ( kind = 4 ) m integer ( kind = 4 ) mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer ( kind = 4 ) n integer ( kind = 4 ) s integer ( kind = 4 ) values(8) integer ( kind = 4 ) y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 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, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end