subroutine bndacc ( g, mdg, nb, ip, ir, mt, jt ) ! !******************************************************************************* ! !! BNDACC accumulates information for a banded least squares problem. ! ! ! Discussion: ! ! BNDACC is used for the accumulation phase. Then BNDSOL should be called ! to solve linear systems. ! ! The user must set IR = 1 and IP = 1 before the first call ! to BNDACC for a new case. ! ! BNDACC is to be called once for each block of data [ C(I), B(I) ] ! to be introduced into the problem. For each block of data, the ! user must assign values to MT and JT, and copy the MT by (NB+1) ! array of data [ C(I), B(I) ] into rows IR through IR+MT-1 of ! the working array G. ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Input/output, double precision G(MDG,NB+1), the banded matrix being ! accumulated. ! ! Input, integer MDG, the leading dimension of G. ! MDG must be at least equal to the number of rows. ! ! Input, integer NB, the bandwidth of G, not counting the diagonal. ! ! Input/output, integer IP, must be set to 1 by the user before ! the first call. Thereafter, its value is controlled by the routine. ! ! Input/output, integer IR, indicates the index of the first row ! of G into which new data is to be placed. The user sets IR to ! 1 before the first call, but does not alter it thereafter. ! Instead, the program keeps this quantity up to date as more ! data is read in. ! ! Input, integer MT, is set by the user to indicate the number of ! new rows of data being introduced by the current call. ! MT must be at least 0. ! ! Input, integer JT, set by the user to indicate the column of the ! submatrix A(I) that is identified with the first column of C(I). ! This means that JT must be at least 1. ! implicit none ! integer mdg integer nb ! double precision g(mdg,nb+1) integer i integer ie integer ig integer ig1 integer ig2 integer ip integer ir integer j integer jg integer jt integer k integer kh integer l integer mh integer mode integer mt integer mu double precision rho ! if ( mt <= 0 ) then return end if if ( jt /= ip ) then if ( jt > ir ) then do i = 1, mt ig1 = jt+mt-i ig2 = ir+mt-i do j = 1, nb+1 g(ig1,j) = g(ig2,j) end do end do do i = 1, jt - ir ig = ir+i-1 g(ig,1:nb+1) = 0.0D+00 end do ir = jt end if mu = min ( nb-1, ir-ip-1 ) do l = 1, mu k = min ( l, jt-ip ) ig = ip+l do i = l+1, nb g(ig,i-k) = g(ig,i) end do g(ig,nb+1-k:nb) = 0.0D+00 end do ip = jt end if mh = ir + mt - ip kh = min ( nb+1, mh ) mode = 1 do i = 1, kh if ( i+1 <= nb ) then call h12 ( mode, i, max ( i+1, ir-ip+1 ), mh, g(ip,i), 1, rho, & g(ip,i+1), 1, mdg, nb+1-i) else call h12 ( mode, i, max ( i+1, ir-ip+1 ), mh, g(ip,i), 1, rho, & g(ip,1), 1, mdg, 0 ) end if end do ir = ip + kh if ( kh >= nb+1 ) then g(ip+kh-1,1:nb) = 0.0D+00 end if return end subroutine bndsol ( mode, g, mdg, nb, ip, ir, x, n, rnorm ) ! !******************************************************************************* ! !! BNDSOL solves a banded least squares problem accumulated by BNDACC. ! ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Input, integer MODE, specifies what system is to be solved. ! 1, solve R * X = Y, where R and Y are in stored in G. ! 2, solve R' * X = Y, where R is in G, Y in X. ! 3, solve R * X = Y, where R is in G, Y is in X. ! ! Input, double precision G(MDG,NB+1), the banded matrix. ! ! Input, integer MDG, the leading dimension of G. ! MDG must be at least equal to the number of rows. ! ! Input, integer NB, the bandwidth of G, not counting the diagonal. ! ! Input, integer IP, IR, must have the values returned by BNDACC ! from its last call. ! ! Output, double precision X(N), the computed solution. ! ! Input, integer N, the dimension of the solution. BNDSOL will ! use just the N by N logical submatrix of G. ! ! Output, double precision RNORM, is set if MODE was 1 on input. ! Its value is the L2 norm of G(N+1:IR-1,NB+1). ! implicit none ! integer mdg integer n integer nb ! double precision g(mdg,nb+1) integer i integer i1 integer i2 integer ie integer ii integer ip integer ir integer ix integer j integer jg integer l integer mode double precision rnorm double precision rsq double precision s double precision x(n) ! rnorm = 0.0D+00 if ( mode == 1 ) then x(1:n) = g(1:n,nb+1) rsq = sqrt ( sum ( g(n+1:ir-1,nb+1)**2 ) ) end if if ( mode == 1 .or. mode == 3 ) then do ii = 1, n i = n+1-ii s = 0.0D+00 l = max ( 0, i-ip ) ie = min ( n+1-i, nb ) do j = 2, ie jg = j+l ix = i-1+j s = s + g(i,jg) * x(ix) end do if ( g(i,l+1) == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BNDSOL - Fatal error!' write ( *, '(a)' ) ' Zero diagonal term.' write ( *, '(a,i6)' ) ' MODE = ', mode write ( *, '(a,i6)' ) ' I = ', i write ( *, '(a,i6)' ) ' J = ', j write ( *, '(a,i6)' ) ' L = ', l stop end if x(i) = ( x(i) - s ) / g(i,l+1) end do else if ( mode == 2 ) then do j = 1, n s = 0.0D+00 if ( j /= 1 ) then i1 = max ( 1, j-nb+1 ) do i = i1, j-1 l = j - i + 1 + max ( 0, i-ip ) s = s + x(i) * g(i,l) end do end if l = max ( 0, j-ip ) if ( g(j,l+1) == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BNDSOL - Fatal error!' write ( *, '(a)' ) ' Zero diagonal term.' write ( *, '(a,i6)' ) ' I = ', i write ( *, '(a,i6)' ) ' J = ', j write ( *, '(a,i6)' ) ' L = ', l stop end if x(j) = ( x(j) - s ) / g(j,l+1) end do end if return end subroutine d_swap ( x, y ) ! !******************************************************************************* ! !! D_SWAP swaps two double precision values. ! ! ! Modified: ! ! 01 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, double precision X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none ! double precision x double precision y double precision z ! z = x x = y y = z return end function diff ( x, y ) ! !******************************************************************************* ! !! DIFF is used in tests that depend on machine precision. ! ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Input, double precision X, Y, the numbers whose difference is desired. ! ! Output, double precision DIFF, the difference of X and Y. ! implicit none ! double precision diff double precision x double precision y ! diff = x - y return end subroutine g1 ( a, b, c, s, sig ) ! !******************************************************************************* ! !! G1 computes an orthogonal rotation matrix. ! ! ! Discussion: ! ! The routine is given a vector ! ! ( A ) ! ( B ) ! ! and computes quantities C and S so that ! ! ( C S ) ( A ) ( SQRT ( A**2 + B**2 ) ) ! ( -S C ) * ( B ) = ( 0 ) ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Modified: ! ! 30 November 2000 ! ! Parameters: ! ! Input, double precision A, B, the entries of the vector. ! ! Output, double precision C, S, the entries of the rotation matrix. ! ! Output, double precision SIG, the value of SQRT ( A**2 + B**2 ). ! implicit none ! double precision a double precision b double precision c double precision s double precision sig double precision xr double precision yr ! if ( abs ( a ) > abs ( b ) ) then xr = b / a yr = sqrt ( 1.0D+00 + xr**2 ) c = sign ( 1.0D+00 / yr, a ) s = c * xr sig = abs ( a ) * yr else if ( b /= 0.0D+00 ) then xr = a / b yr = sqrt ( 1.0D+00 + xr**2 ) s = sign ( 1.0D+00 / yr, b ) c = s * xr sig = abs ( b ) * yr else c = 0.0D+00 s = 1.0D+00 sig = 0.0D+00 end if return end subroutine g2 ( c, s, x, y ) ! !******************************************************************************* ! !! G2 applies a rotation matrix to a vector (X,Y). ! ! ! Discussion: ! ! This routine is no longer used by other routines in the package, ! but is included for completeness. ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Input, double precision C, S, the cosine and sine of the rotation. ! ! Input/output, double precision X, Y, the components of the vector ! to be rotated. ! implicit none ! double precision c double precision s double precision x double precision x_new double precision y double precision y_new ! x_new = c * x + s * y y_new = - s * x + c * y x = x_new y = y_new return end function gen ( anoise ) ! !******************************************************************************* ! !! GEN generates numbers for construction of test cases. ! ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Input, double precision ANOISE, ... ! ! Output, double precision GEN, ... ! implicit none ! double precision, save :: aj = 0.0D+00 double precision anoise double precision gen integer, save :: i = 5 integer, save :: j = 7 integer, parameter :: mi = 891 integer, parameter :: mj = 457 ! if ( anoise < 0.0D+00 ) then aj = 0.0D+00 gen = 0.0D+00 return end if ! ! The sequence of values of J is bounded between 1 and 996. ! If initial j = 1,2,3,4,5,6,7,8, or 9, the period is 332 ! if ( anoise > 0.0D+00 ) then j = j * mj j = j - 997 * ( j / 997 ) aj = dble ( j - 498 ) end if ! ! The sequence of values of I is bounded between 1 and 999. ! If initial i = 1,2,3,6,7, or 9, the period will be 50. ! If initial i = 4 or 8 the period will be 25. ! If initial i = 5 the period will be 10. ! i = i * mi i = i - 1000 * ( i / 1000 ) gen = dble ( i - 500 ) + aj * anoise return end subroutine h12 ( mode, lpivot, l1, m, u, iue, up, c, ice, icv, ncv ) ! !******************************************************************************* ! !! H12 constructs or applies a Householder transformation. ! ! ! Discussion: ! ! The symmetric orthogonal Householder transformation matrix ! has the form ! ! Q = I + u * u' / b ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Input, integer MODE, selects the algorithm. ! 1, algorithm H1 to construct and apply a Householder transformation, ! 2, algorithm H2 to apply a previously constructed transformation. ! ! Input, integer LPIVOT, is the index of the pivot element, which ! should satisfy 1 <= LPIVOT < L1. ! ! l1,m if l1 <= m the transformation will be constructed to ! zero elements indexed from l1 through m. if l1 gt. m ! the subroutine does an identity transformation. ! ! u(),iue,up on entry with mode = 1, U contains the pivot ! vector. iue is the storage increment between elements. ! on exit when mode = 1, U and up contain quantities ! defining the vector u of the householder transformation. ! on entry with mode = 2, U and up should contain ! quantities previously computed with mode = 1. these will ! not be modified during the entry with mode = 2. ! ! Input/output, double precision C(*). On input, C contains a matrix ! which will be regarded as a set of vectors to which the ! Householder transformation is to be applied. ! On output, C contains the set of transformed vectors. ! ! Input, integer ICE, the storage increment between successive elements ! of a single vector in C. ! ! Input, integer ICV, the storage increment between the first element ! of succesive vectors in C. ! ! Input, integer NCV, the number of vectors in C to be transformed. ! If NCV <= 0, no operations are done on C. ! implicit none ! integer iue integer m ! double precision b double precision c(*) double precision cl double precision clinv integer i integer i2 integer i3 integer i4 integer ice integer icv integer incr integer j integer l1 integer lpivot integer mode integer ncv double precision sm double precision u(iue,m) double precision up ! if ( lpivot <= 0 ) then return else if ( lpivot >= l1 ) then return else if ( l1 > m ) then return end if cl = abs ( u(1,lpivot) ) ! ! Construct the transformation. ! if ( mode == 1 ) then do j = l1, m cl = max ( cl, abs ( u(1,j) ) ) end do if ( cl <= 0.0D+00 ) then return end if clinv = 1.0D+00 / cl sm = ( u(1,lpivot) * clinv )**2 do j = l1, m sm = sm + ( u(1,j) * clinv )**2 end do cl = cl * sqrt ( sm ) if ( u(1,lpivot) > 0.0D+00 ) then cl = - cl end if up = u(1,lpivot) - cl u(1,lpivot) = cl else if ( mode == 2 ) then if ( cl <= 0.0D+00 ) then return end if end if ! ! Apply the transformation i+u*(u**t)/b to the NCV vectors in C. ! if ( ncv <= 0 ) then return end if b = up * u(1,lpivot) ! ! B must be nonpositive. ! if ( b >= 0.0D+00 ) then return end if b = 1.0D+00 / b i2 = 1 - icv + ice * ( lpivot - 1 ) incr = ice * ( l1 - lpivot ) do j = 1, ncv i2 = i2 + icv i3 = i2 + incr i4 = i2 + incr sm = c(i2) * up do i = l1, m sm = sm + c(i3) * u(1,i) i3 = i3 + ice end do if ( sm /= 0.0D+00 ) then sm = sm * b c(i2) = c(i2) + sm * up do i = l1, m c(i4) = c(i4) + sm * u(1,i) i4 = i4 + ice end do end if end do return end subroutine hfti ( a, mda, m, n, b, mdb, nb, tau, krank, rnorm, h, g, ip ) ! !******************************************************************************* ! !! HFTI: Householder forward triangulation with column interchanges. ! ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Modified: ! ! 15 May 2002 ! ! Parameters: ! ! Input/output, double precision A(MDA,N), the M by N matrix. ! ! Input, integer MDA, the leading dimension of A. MDA must be at least M. ! ! Input, integer M, N, the number of rows and columns in A. ! ! Input/output, double precision B(MDB,NB), on input, the NB right hand ! side vectors of length M, and on output, the NB least squares solution ! vectors of length N. ! ! Input, integer MDB, the leading dimension of B, which must be at least M. ! ! Input, integer NB, the number of columns in B. ! ! Input, double precision TAU, a tolerance used in the pseudorank test. ! ! Output, integer KRANK, the numerically determined rank of the ! matrix A. ! ! Output, double precision RNORM(NB), the Euclidean norms of the residual ! vectors for the least squares problems. ! ! Workspace, double precision H(N). ! ! Workspace, double precision G(N). ! ! Workspace, integer IP(N). ! implicit none ! integer mda integer mdb integer n integer nb ! double precision a(mda,n) double precision b(mdb,nb) double precision diff double precision, parameter :: factor = 0.001D+00 double precision g(n) double precision h(n) double precision hmax integer i integer ii integer ip(n) integer j integer jb integer jj integer k integer krank integer l integer ldiag integer lmax integer m integer mode double precision rnorm(nb) double precision sm double precision tau double precision tmp ! k = 0 ldiag = min ( m, n ) if ( ldiag <= 0 ) then krank = 0 return end if do j = 1, ldiag ! ! Update the squared column lengths and find LMAX. ! In the original code, LMAX was not set if J <= 1, causing errors. ! if ( j <= 1 ) then lmax = 1 else lmax = j do l = j, n h(l) = h(l) - a(j-1,l)**2 if ( h(l) > h(lmax) ) then lmax = l end if end do end if ! ! Compute the squared column lengths and find LMAX. ! if ( diff ( hmax + factor * h(lmax), hmax ) <= 0.0D+00 ) then lmax = j do l = j, n h(l) = 0.0D+00 do i = j, m h(l) = h(l) + a(i,l)**2 end do if ( h(l) > h(lmax) ) then lmax = l end if end do hmax = h(lmax) end if ! ! Do column interchanges if needed. ! ip(j) = lmax if ( ip(j) /= j ) then do i = 1, m call d_swap ( a(i,j), a(i,lmax) ) end do h(lmax) = h(j) end if ! ! Compute the J-th transformation and apply it to A and B. ! mode = 1 if ( j < n ) then call h12 ( mode, j, j+1, m, a(1,j), 1, h(j), a(1,j+1), 1, mda, n-j ) else call h12 ( mode, j, j+1, m, a(1,j), 1, h(j), a(1,1), 1, mda, n-j ) end if mode = 2 call h12 ( mode, j, j+1, m, a(1,j), 1, h(j), b, 1, mdb, nb ) end do ! ! Determine the pseudorank K using the tolerance TAU. ! k = ldiag do j = 1, ldiag if ( abs ( a(j,j) ) <= tau ) then k = j - 1 exit end if end do ! ! Compute the norms of the residual vectors. ! do jb = 1, nb rnorm(jb) = sqrt ( sum ( b(k+1:m,jb)**2 ) ) end do ! ! Special for pseudorank 0. ! if ( k <= 0 ) then b(1:n,1:nb) = 0.0D+00 krank = k return end if ! ! If the pseudorank is less than N, compute the Householder ! decomposition of the first K rows. ! if ( k /= n ) then mode = 1 do ii = 1, k i = k+1-ii call h12 ( mode, i, k+1, n, a(i,1), mda, g(i), a, mda, 1, i-1 ) end do end if do jb = 1, nb ! ! Solve the K by K triangular system. ! do l = 1, k sm = 0.0D+00 i = k+1-l do j = i+1, k sm = sm + a(i,j) * b(j,jb) end do b(i,jb) = ( b(i,jb) - sm ) / a(i,i) end do ! ! Complete computation of solution vector. ! if ( k < n ) then b(k+1:n,jb) = 0.0D+00 mode = 2 do i = 1, k call h12 ( mode, i, k+1, n, a(i,1), mda, g(i), b(1,jb), 1, mdb, 1 ) end do end if ! ! Reorder the solution vector to compensate for the column interchanges. ! do jj = 1, ldiag j = ldiag+1-jj if ( ip(j) /= j) then l = ip(j) call d_swap ( b(l,jb), b(j,jb) ) end if end do end do ! ! The solution vectors X are now in the first N rows of B. ! krank = k return end subroutine ldp ( g, m, n, h, x, xnorm, mode ) ! !******************************************************************************* ! !! LDP implements least distance programming algorithm. ! ! ! Discussion: ! ! The problem is to determine an N vector X with minimum norm, ! subject to ! ! G * X >= H ! ! where G is an M by N matrix, and H is an M vector. ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Modified: ! ! 03 May 2002 ! ! Parameters: ! ! Input, double precision G(M,N), the M by N matrix. ! There is no restriction on the rank of G. ! ! Input, integer M, N, the number of rows and columns in G. ! M may be greater or less than N. ! ! Input, double precision H(M), the right hand side of the constraint. ! ! Output, double precision X(N), the computed solution. ! ! Output, double precision XNORM, the value of the norm of X. ! ! Output, integer MODE, error flag. ! 1, no error. ! 2, bad dimension. ! 3, maximum number of iterations taken without convergence. ! 4, the inequality constraints G * X >= H are incompatible. ! ! Local Parameters: ! ! double precision W((n+1)*(m+1)+m). ! ! first (n+1)*m locs of w() = matrix e for problem nnls. ! next n+1 locs of w() = vector f for problem nnls. ! next m locs of w() = vector y for problem nnls. ! implicit none ! integer m integer n ! double precision diff double precision fac double precision g(m,n) double precision h(m) integer i integer iw integer iy integer j integer jf integer mode double precision rnorm double precision w(m*n+2*m+n+1) double precision x(n) double precision xnorm ! if ( n <= 0 ) then mode = 2 return end if x(1:n) = 0.0D+00 xnorm = 0.0D+00 if ( m <= 0 ) then mode = 1 return end if iw = 0 do j = 1, m do i = 1, n iw = iw + 1 w(iw) = g(j,i) end do iw = iw + 1 w(iw) = h(j) end do jf = iw + 1 ! ! Store N zeros followed by a 1 into F. ! do i = 1, n iw = iw + 1 w(iw) = 0.0D+00 end do w(iw+1) = 1.0D+00 iy = iw + 2 call nnls ( w, n+1, m, w(jf), w(iy), rnorm, mode ) if ( mode /= 1 ) then return end if if ( rnorm <= 0.0D+00 ) then mode = 4 return end if fac = 1.0D+00 iw = iy - 1 do i = 1, m iw = iw + 1 fac = fac - h(i) * w(iw) end do if ( diff ( 1.0D+00 + fac, 1.0D+00 ) == 0.0D+00 ) then mode = 4 return end if fac = 1.0D+00 / fac do j = 1, n iw = iy - 1 do i = 1, m iw = iw + 1 x(j) = x(j) + g(i,j) * w(iw) end do end do x(1:n) = x(1:n) * fac xnorm = sqrt ( sum ( x(1:n)**2 ) ) mode = 1 return end subroutine mfeout ( a, mda, m, n, names, mode, unit, width ) ! !******************************************************************************* ! !! MFEOUT labeled matrix output for use with singular value analysis. ! ! ! Discussion: ! ! this 1995 version has additional arguments, unit and width, ! to support user options regarding the output unit and the width of ! print lines. also allows user to choose length of names in names(). ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Input, double precision A(MDA,N), the M by N matrix to be output. ! ! Input, integer MDA, the leading dimension of A. MDA must be at least M. ! ! Input, integer M, N, the number of rows and columns in A. ! ! Input, names() [character array] array of names. ! if names(1) contains only blanks, the rest of the names() ! array will be ignored. ! ! Input, mode ! = 1 write header for v matrix and use an f format. ! = 2 write header for candidate solutions and use p format. ! ! Input, unit [integer] selects output unit. if unit .ge. 0 then unit ! is the output unit number. if unit = -1, output to ! the '*' unit. ! ! Input, width [integer] selects width of output lines. ! each output line from this subroutine will have at most ! max ( 26, min ( 124, width ) ) characters plus one additional ! leading character for fortran "carriage control". the ! carriage control character will always be a blank. ! implicit none ! integer mda integer n ! double precision a(mda,n) logical blknam character ( len = 26 ) fmt1(2) character ( len = 26 ) fmt2(2) character ( len = 4 ) head(2) integer i integer j integer j1 integer j2 integer kblock integer l integer lennam integer m integer maxcol integer mode character ( len = * ) names(m) integer namsiz integer nblock logical star integer unit integer width ! data head(1)/' col'/ data head(2)/'soln'/ data fmt1 / '(/7x,00x,8(5x,a4,i4,1x)/)', & '(/7x,00x,8(2x,a4,i4,4x)/)'/ data fmt2 / '(1x,i4,1x,a00,1x,4p8f14.0)', & '(1x,i4,1x,a00,1x,8g14.6 )'/ ! if ( m <= 0 ) then return end if if ( n <= 0 ) then return end if star = unit < 0 ! ! The len function returns the char length of a single element of ! the names() array. ! lennam = len ( names(1) ) blknam = names(1) == ' ' namsiz = 1 if ( .not. blknam) then do i = 1, m namsiz = max ( namsiz, len_trim ( names(i) ) ) end do end if write(fmt1(mode)(6:7),'(i2.2)') namsiz write(fmt2(mode)(12:13),'(i2.2)') namsiz 70 format(/' v-matrix of the singular value decomposition of a*d.'/ & ' (elements of v scaled up by a factor of 10**4)') 80 format(/' sequence of candidate solutions, x') if ( star ) then if (mode == 1) then write (*,70) else write (*,80) end if else if (mode == 1) then write (unit,70) else write (unit,80) end if end if ! ! With namsiz characters allowed for the "name" and maxcol ! columns of numbers, the total line width, exclusive of a ! carriage control character, will be 6 + lennam + 14 * maxcol. ! maxcol = max ( 1, min ( 8, ( width - 6 - namsiz ) / 14 ) ) nblock = (n + maxcol -1) / maxcol j2 = 0 do kblock = 1, nblock j1 = j2 + 1 j2 = min ( n, j2 + maxcol ) if ( star ) then write (*,fmt1(mode)) (head(mode),j,j = j1,j2) else write (unit,fmt1(mode)) (head(mode),j,j = j1,j2) end if do i = 1,m if ( star ) then if ( blknam ) then write (*,fmt2(mode)) i,' ',(a(i,j),j = j1,j2) else write (*,fmt2(mode)) i,names(i),(a(i,j),j = j1,j2) end if else if ( blknam ) then write (unit,fmt2(mode)) i,' ',(a(i,j),j = j1,j2) else write (unit,fmt2(mode)) i,names(i),(a(i,j),j = j1,j2) end if end if end do end do return end subroutine nnls ( a, m, n, b, x, rnorm, mode ) ! !******************************************************************************* ! !! NNLS implements the nonnegative least squares algorithm. ! ! ! Discussion: ! ! Given an M by N matrix A and an M vector B, compute an ! N vector X that solves the least squares problem (that is, ! minimizes the L2 norm of the residual) ! ! A * X = B ! ! subject to ! ! X >= 0. ! ! Householder transformations are applied, which essentially ! multiply the system by an orthogonal matrix Q: ! ! Q * A * X = Q * B ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Modified: ! ! 03 May 2002 ! ! Parameters: ! ! Input/output, double precision A(M,N). On input, the M by N matrix. ! On output, the product matrix Q * A, where Q is an M by M orthogonal ! matrix generated implicitly by this routine. ! ! Input, integer M, N, the number of rows and columns in A. ! ! Input/output, double precision B(M). On input, contains the right ! hand side vector. On output, contains Q * B. ! ! Output, double precision X(N), the computed solution. ! ! Output, double precision RNORM, the L2 or Euclidean norm of the residual. ! ! Output, integer MODE, error flag. ! 1, the solution has been computed successfully. ! 2, the dimensions of the problem are bad, M <= 0 or N <= 0. ! 3, iteration count exceeded. More than 3*N iterations. ! implicit none ! integer m integer n ! double precision a(m,n) double precision alpha double precision asave double precision b(m) double precision cc double precision diff double precision dummy double precision, parameter :: factor = 0.01D+00 integer i integer ii integer index(n) integer ip integer iter integer itmax integer iz integer iz1 integer iz2 integer izmax integer j integer jj integer jz integer l integer mode integer npp1 integer nsetp double precision rnorm integer rtnkey double precision sm double precision ss double precision t double precision temp double precision unorm double precision up double precision w(n) double precision wmax double precision x(n) double precision ztest double precision zz(m) ! mode = 1 if ( m <= 0 .or. n <= 0 ) then mode = 2 return end if iter = 0 itmax = 3 * n ! ! Initialize INDEX and X. ! x(1:n) = 0.0D+00 do i = 1, n index(i) = i end do iz2 = n iz1 = 1 nsetp = 0 npp1 = 1 ! ! The main loop begins here. ! 30 continue ! ! Quit if all coefficients are already in the solution. ! or if M columns of A have been triangularized. ! if ( iz1 > iz2 .or. nsetp >= m ) then go to 350 end if ! ! Compute components of the dual (negative gradient) vector W. ! do iz = iz1, iz2 j = index(iz) w(j) = dot_product ( b(npp1:m), a(npp1:m,j) ) end do ! ! Find the largest positive W. ! 60 continue wmax = 0.0D+00 do iz = iz1, iz2 j = index(iz) if ( w(j) > wmax) then wmax = w(j) izmax = iz end if end do ! ! If WMAX <= 0, go to termination. ! This indicates satisfaction of the Kuhn-Tucker conditions. ! if ( wmax <= 0.0D+00 ) then go to 350 end if iz = izmax j = index(iz) ! ! 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 = a(npp1,j) call h12 ( 1, npp1, npp1+1, m, a(1,j), 1, up, dummy, 1, 1, 0 ) unorm = sqrt ( sum ( a(1:nsetp,j)**2 ) ) if ( diff ( unorm+abs(a(npp1,j))*factor, unorm ) > 0.0D+00 ) then ! ! Column J is sufficiently independent. ! ! Copy B into ZZ, update ZZ ! and solve for ztest ( = proposed new value for x(j) ). ! zz(1:m) = b(1:m) call h12 ( 2, npp1, npp1+1, m, a(1,j), 1, up, zz, 1, 1, 1 ) ztest = zz(npp1) / a(npp1,j) ! ! See if ZTEST is positive. ! if ( ztest > 0.0D+00 ) then go to 140 end if end if ! ! Reject j as a candidate to be moved from set z to set p. ! restore a(npp1,j), set w(j) = 0., and loop back to test dual ! coefficients again. ! a(npp1,j) = asave w(j) = 0.0D+00 go to 60 ! ! The index J = index(iz) has been selected to be moved from ! set z to set p. ! ! Update b, update indices, Apply Householder ! transformations to columns in new set z, zero subdiagonal elements in ! column j, set w(j) = 0. ! 140 continue b(1:m) = zz(1:m) index(iz) = index(iz1) index(iz1) = j iz1 = iz1+1 nsetp = npp1 npp1 = npp1+1 if ( iz1 <= iz2) then do jz = iz1,iz2 jj = index(jz) call h12 ( 2, nsetp, npp1, m, a(1,j), 1, up, a(1,jj), 1, m, 1 ) end do end if if ( nsetp /= m) then a(npp1:m,j) = 0.0D+00 end if w(j) = 0.0D+00 ! ! Solve the triangular system. ! Store the solution temporarily in ZZ. ! rtnkey = 1 go to 400 200 continue ! ! Secondary loop begins here. ! ! Iteration counter. ! 210 continue iter = iter + 1 if ( iter > itmax) then mode = 3 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NNLS - Warning!' write ( *, '(a)' ) ' Quitting on iteration count.' go to 350 end if ! ! See if all new constrained coefficients are feasible. ! If not, compute alpha. ! alpha = 2.0D+00 do ip = 1, nsetp l = index(ip) if ( zz(ip) <= 0.0D+00 ) then t = -x(l)/(zz(ip)-x(l)) if ( alpha > t) then alpha = t jj = ip end if end if end do ! ! If all new constrained coefficients are feasible then ALPHA will ! still be 2. If so exit from secondary loop to main loop. ! if ( alpha == 2.0D+00 ) go to 330 ! ! Otherwise use ALPHA which will be between 0.0D+00 and 1.0 to ! interpolate between the old X and the new ZZ. ! do ip = 1, nsetp l = index(ip) x(l) = x(l) + alpha * ( zz(ip) - x(l) ) end do ! ! Modify A and b and the index arrays to move coefficient i ! from set p to set Z. ! i = index(jj) 260 continue x(i) = 0.0D+00 if ( jj /= nsetp ) then jj = jj+1 do j = jj, nsetp ii = index(j) index(j-1) = ii call g1 ( a(j-1,ii), a(j,ii), cc, ss, a(j-1,ii) ) ! ! Apply procedure g2 (cc,ss,a(j-1,l),a(j,l)) ! a(j,ii) = 0.0D+00 do l = 1, n if ( l /= ii ) then temp = a(j-1,l) a(j-1,l) = cc*temp + ss*a(j,l) a(j,l) = -ss*temp + cc*a(j,l) end if end do ! ! Apply procedure g2 (cc,ss,b(j-1),b(j)) ! temp = b(j-1) b(j-1) = cc*temp + ss*b(j) b(j) = -ss*temp + cc*b(j) end do end if npp1 = nsetp nsetp = nsetp-1 iz1 = iz1-1 index(iz1) = i ! ! 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 nonpositive will be set to zero ! and moved from set p to set z. ! do jj = 1, nsetp i = index(jj) if ( x(i) <= 0.0D+00 ) then go to 260 end if end do ! ! Copy B into ZZ, then solve again and loop back. ! zz(1:m) = b(1:m) rtnkey = 2 go to 400 ! ! End of secondary loop. ! 330 continue do ip = 1, nsetp i = index(ip) x(i) = zz(ip) end do ! ! All new coefficients are positive. Loop back to beginning. ! go to 30 ! ! End of the main loop ! ! Come here for termination. ! Compute the norm of the final residual vector. ! 350 continue sm = sum ( b(npp1:m)**2 ) if ( npp1 > m ) then w(1:n) = 0.0D+00 end if rnorm = sqrt ( sm ) return ! ! The following block of code is used as an internal subroutine ! to solve the triangular system, putting the solution in ZZ. ! 400 continue do ip = nsetp, 1, -1 if ( ip /= nsetp ) then do ii = 1, ip zz(ii) = zz(ii) - a(ii,jj) * zz(ip+1) end do end if jj = index(ip) zz(ip) = zz(ip) / a(ip,jj) end do ! ! Something wrong here...Maybe 210 should be 330??? ! if ( rtnkey == 1 ) then go to 200 else go to 210 end if end subroutine qrbd ( ipass, q, e, nn, v, mdv, nrv, c, mdc, ncc ) ! !******************************************************************************* ! !! QRBD uses the QR algorithm for the singular values of a bidiagonal matrix. ! ! ! Discussion: ! ! The bidiagonal matrix D, with diagonal entries Q(1:N) and ! superdiagonal entries E(2:N) is pre- and post-multiplied by ! elementary rotation matrices R and P so that ! ! Rk*...*R2*R1 * D * P1'*P2'*...*Pk' = diag(S1,...,Sn) ! ! Reference: ! ! Reinsch and Golub, ! Singular Value Decomposition and Least Squares Solutions, ! Numerische Mathematik, ! Volume 14, 1970. ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Output, integer IPASS, error flag. ! 1, no error was detected. ! 2, an error occurred. ! ! Input/output, double precision Q(NN). On input, the diagonal ! entries of the matrix. On output, the singular values, which ! are nonnegative, listed in nonincreasing order. ! ! Input/output, double precision E(NN). On input, E(2:N) contains ! the superdiagonal entries of the matrix. ! ! Input, integer NN, the order of the matrix. ! ! Input/output, double precision V(MDV,NN), an NRV by NN matrix ! whose input value is to be postmultiplied by the transformations ! P1' * P2' * ... * PK. ! ! Input, integer MDV, the leading dimension of V. ! ! Input, integer NRV, the number of rows in V. ! ! Input/output, double precision C(MDC,NCC), an NN by NCC matrix ! whose input value is to be premultiplied by the transformations ! Rm * ... * R2 * R1. ! ! Input, integer MDC, the leading dimension of C, which must be ! at least NN. ! ! Input, integer NCC, the number of columns in C. ! implicit none ! integer mdc integer mdv integer ncc integer nn ! double precision c(mdc,ncc) double precision cs double precision diff double precision dnorm double precision e(nn) double precision f logical fail double precision g double precision h logical havers integer i integer ii integer ipass integer j integer k integer kk integer l integer ll integer n integer n10 integer nqrs integer nrv double precision q(nn) double precision small double precision sn double precision t double precision temp double precision v(mdv,nn) logical wntv double precision x double precision y double precision z ! n = nn ipass = 1 if ( n <= 0 ) then return end if n10 = 10 * n wntv = nrv > 0 havers = ncc > 0 fail = .false. nqrs = 0 e(1) = 0.0D+00 dnorm = 0.0D+00 do j = 1, n dnorm = max ( dnorm, abs ( q(j) ) + abs ( e(j) ) ) end do do k = n, 1, -1 ! ! Test for splitting or rank deficiencies. ! First make test for last diagonal term, Q(K), being small. ! 20 continue if ( k == 1 ) then go to 50 end if if ( diff ( dnorm+q(k), dnorm ) /= 0.0D+00 ) then go to 50 end if ! ! Since Q(K) is small we will make a special pass to ! transform E(k) to zero. ! ! Transformation constructed to zero out position (I,K). ! cs = 0.0D+00 sn = -1.0D+00 do i = k-1, 1, -1 f = -sn * e(i+1) e(i+1) = cs * e(i+1) call g1 ( q(i), f, cs, sn, q(i) ) ! ! Accumulate right hand transformations in V. ! if ( wntv ) then do j = 1, nrv temp = v(j,i) v(j,i) = cs * temp + sn * v(j,k) v(j,k) = -sn * temp + cs * v(j,k) end do end if end do ! ! The matrix is now bidiagonal, and of lower order ! since E(K) == zero. ! 50 continue do l = k, 1, -1 if ( diff ( dnorm+e(l), dnorm ) == 0.0D+00 ) then go to 100 end if if ( diff ( dnorm+q(l-1), dnorm ) == 0.0D+00 ) then go to 70 end if end do ! ! This loop can't complete since E(1) = zero. ! go to 100 ! ! Cancellation of E(L), L > 1. ! 70 continue cs = 0.0D+00 sn = -1.0D+00 do i = l, k f = -sn * e(i) e(i) = cs * e(i) if ( diff ( dnorm+f, dnorm ) == 0.0D+00 ) then go to 100 end if call g1 ( q(i), f, cs, sn, q(i) ) ! ! Accumulate left hand transformations in C. ! if ( havers ) then do j = 1, ncc temp = c(i,j) c(i,j) = cs * temp + sn * c(l-1,j) c(l-1,j) = -sn * temp + cs * c(l-1,j) end do end if end do ! ! Test for convergence. ! 100 continue z = q(k) if ( l == k ) then go to 170 end if ! ! Shift from bottom 2 by 2 minor of B'*B. ! x = q(l) y = q(k-1) g = e(k-1) h = e(k) f = ( (y-z) * (y+z) + (g-h) * (g+h) ) / ( 2.0D+00 * h * y ) g = sqrt ( 1.0D+00 + f**2 ) if ( f >= 0.0D+00 ) then t = f + g else t = f - g end if f = ( (x-z) * (x+z) + h * ( y / t - h ) ) / x ! ! Next QR sweep. ! cs = 1.0D+00 sn = 1.0D+00 do i = l+1, k g = e(i) y = q(i) h = sn * g g = cs * g call g1 ( f, h, cs, sn, e(i-1) ) f = x * cs + g * sn g = -x * sn + g * cs h = y * sn y = y * cs ! ! Accumulate rotations from the right in V. ! if ( wntv ) then do j = 1, nrv temp = v(j,i-1) v(j,i-1) = cs * temp + sn * v(j,i) v(j,i) = -sn * temp + cs * v(j,i) end do end if call g1 ( f, h, cs, sn, q(i-1) ) f = cs * g + sn * y x = -sn * g + cs * y ! ! Accumulate rotations from the left in C. ! if ( havers ) then do j = 1, ncc temp = c(i-1,j) c(i-1,j) = cs * temp + sn * c(i,j) c(i,j) = -sn * temp + cs * c(i,j) end do end if end do e(l) = 0.0D+00 e(k) = f q(k) = x nqrs = nqrs + 1 if ( nqrs <= n10 ) then go to 20 end if ! ! Return to 'test for splitting'. ! small = abs(e(k)) i = k ! ! If failure to converge set smallest magnitude ! term in off-diagonal to zero. Continue on. ! do j = l, k temp = abs ( e(j) ) if ( temp /= 0.0D+00 ) then if ( temp < small ) then small = temp i = j end if end if end do e(i) = 0.0D+00 nqrs = 0 fail = .true. go to 20 ! ! Cutoff for convergence failure. NQRS will be 2*N usually. ! 170 continue if ( z < 0.0D+00 ) then q(k) = -z if ( wntv ) then v(1:nrv,k) = - v(1:nrv,k) end if end if 190 continue end do ! ! Convergence. Q(K) is made nonnegative. ! if ( n == 1 ) then return end if do i = 2,n if ( q(i) > q(i-1) ) then go to 220 end if end do if ( fail ) then ipass = 2 end if return ! ! Every singular value is in order. ! 220 continue do i = 2, n t = q(i-1) k = i - 1 do j = i, n if ( t < q(j) ) then t = q(j) k = j end if end do if ( k == i - 1 ) then cycle end if q(k) = q(i-1) q(i-1) = t if ( havers ) then do j = 1, ncc call d_swap ( c(i-1,j), c(k,j) ) end do end if if ( wntv ) then do j = 1, nrv call d_swap ( v(j,i-1), v(j,k) ) end do end if end do ! ! End of ordering algorithm. ! if ( fail ) then ipass = 2 end if return end subroutine sva ( a, mda, m, n, mdata, b, sing, kpvec, names, iscale, d, work ) ! !******************************************************************************* ! !! SVA carries out a singular value analysis. ! ! ! Discussion: ! ! This routine computes the singular value decomposition of the matrix ! of a least squares problem, and produces a printed report. ! ! This 1995 version differs from the original 1973 version by the ! addition of the arguments kpvec() and work(), and by allowing user to ! choose the length of names in names(). ! ! kpvec() allows the user to exercise options regarding printing. ! work() provides 2*n locations of work space. originally sing() was ! required to have 3*n elements, of which the last 2*n were used for ! work space. now sing() only needs n elements. ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Input/output, double precision A(MDA,N). On input, the M by N ! matrix of the least squares problem to be analyzed. This could be ! a matrix obtained by preliminary orthogonal transformations ! applied to the actual problem matrix which may have had more rows. ! See MDATA below. ! ! Input, integer MDA, the first dimensioning parameter for A. ! It must be the case that MDA >= max ( M, N ). ! ! Input, integer M, N, the number of rows and columns in A. ! m > n or m <= n is permitted. ! require m > 0 and n > 0. ! ! Input, integer mdata, number of rows in actual least squares problem. ! generally mdata .ge. m. mdata is used only in computing ! statistics for the report and is not used as a loop ! count or array dimension. ! ! Input/output, b(). on entry, contains the right-side vector, b, of the ! least squares problem. this vector is of length, m. ! on return, contains the vector, g = (u**t)*b, where u ! comes from the singular value decomposition of a. the ! vector , g, is also of length m. ! ! sing() [out] on return, contains the singular values of A, in ! descending order, in locations indexed 1 thru min ( m, n ). ! if m < n, locations indexed from m+1 through n will be ! set to zero. ! ! kpvec() [integer array, in] array of integers to select print ! options. kpvec(1) determines whether the rest of ! the array is to be used or ignored. ! if kpvec(1) = 1, the contents of (kpvec(i), i=2,4) ! will be used to set internal variables as follows: ! prblk = kpvec(2) ! unit = kpvec(3) ! width = kpvec(4) ! if kpvec(1) = 0 default settings will be used. the user ! need not dimension kpvec() greater than 1. the subr will ! set prblk = 111111, unit = -1, and width = 69. ! ! the internal variables prblk, unit, and width are ! interpreted as follows: ! ! prblk the decimal representation of prblk must be ! representable as at most 6 digits, each being 0 or 1. ! the decimal digits will be interpreted as independant ! on/off flags for the 6 possible blocks of printed output. ! examples: 111111 selects all blocks, 0 suppresses all ! printing, 101010 selects the 1st, 3rd, and 5th blocks, ! etc. ! the six blocks are: ! 1. header, with size and scaling option parameters. ! 2. v-matrix. amount of output depends on m and n. ! 3. singular values and related quantities. amount of ! output depends on n. ! 4. listing of ynorm and rnorm and their logarithms. ! amount of output depends on n. ! 5. levenberg-marquart analysis. ! 6. candidate solutions. amount of output depends on ! m and n. ! ! unit selects the output unit. if unit .ge. 0, ! unit will be used as the output unit number. ! if unit = -1, output will be written to the "*" output ! unit, i.e., the standard system output unit. ! the calling program unit is responsible for opening ! and/or closing the selected output unit if the host ! system requires these actions. ! ! width default value is 79. determines the width of ! blocks 2, 3, and 6 of the output report. ! block 3 will use 95(+1) columns if width .ge. 95, and otherwise ! 69(+1) columns. ! blocks 2 and 6 are printed by subroutine mfeout. these blocks ! generally use at most width(+1) columns, but will use more if ! the names are so long that more space is needed to print one ! name and one numeric column. the (+1)'s above are reminders ! that in all cases there is one extra initial column for fortran ! "carriage control". the carriage control character will always ! be a blank. ! blocks 1, 4, and 5 have fixed widths of 63(+1), 66(+1) and ! 66(+1), respectively. ! ! names() [in] names(j), for j = 1, ..., n, may contain a ! name for the jth component of the solution ! vector. the declared length of the elements of the ! names() array is not specifically limited, but a ! greater length reduces the space available for columns ! of the matris to be printed. ! if names(1) contains only blank characters, ! it will be assumed that no names have been provided, ! and this subr will not access the names() array beyond ! the first element. ! ! iscale [in] set by the user to 1, 2, or 3 to select the column ! scaling option. ! 1 subr will use identity scaling and ignore the d() ! array. ! 2 subr will scale nonzero rolumns to have unit euclid- ! ean length and will store reciprocal lengths of ! original nonzero columns in d(). ! 3 user supplies col scale factors in d(). subr ! will mult col j by d(j) and remove the scaling ! from the soln at the end. ! ! d() [ignored or out or in] usage of d() depends on iscale as ! described above. when used, its length must be ! at least n. ! ! work() [scratch] work space of length at least 2*n. used ! directly in this subroutine. ! implicit none ! integer m integer mda integer n ! double precision a(mda,n) double precision a1 double precision a2 double precision a3 double precision a4 double precision alamb double precision aln10 double precision b(m) logical blk(6) double precision d(n) double precision del double precision el double precision el2 integer i integer ie integer ipass integer iscale integer j integer k integer kpvec(4) integer mdata integer minmn integer minmn1 integer mpass character ( len = * ) names(n) logical narrow integer nsol double precision pcoef integer prblk double precision rl double precision rnorm double precision rs double precision sb double precision sing(n) double precision sl logical star integer unit integer width double precision work(2*n) double precision yl double precision ynorm double precision ys double precision ysq ! 220 format (1x/' index sing. val. p coef ', & ' reciprocal g coef g**2 ', & ' cumulative scaled sqrt'/ & 31x,' sing. val.',26x, & ' sum of sqrs of cum.s.s.') 221 format (1x/' index sing. val. p coef ', & ' reciprocal g coef scaled sqrt'/ & 31x,' sing. val.',13x,' of cum.s.s.') 222 format (1x/' index sing. val. g coef g**2 ', & ' cumulative scaled sqrt'/ & 44x,' sum of sqrs of cum.s.s.') 230 format (' ',4x,'0',64x,2g13.4) 231 format (' ',4x,'0',51x, g13.4) 232 format (' ',4x,'0',38x,2g13.4) 240 format (' ',i5,g12.4,6g13.4) 260 format (1x,' m = ',i6,', n =',i4,', mdata =',i8) 270 format (1x/' singular value analysis of the least squares', & ' problem, a*x = b,'/' scaled as (a*d)*y = b.') 280 format (1x/' scaling option number',i2,'. d is a diagonal', & ' matrix with the following diagonal elements..'/(5x,10e12.4)) 290 format (1x/' scaling option number 1. d is the identity matrix.'/1x) 300 format (1x/' index',12x,'ynorm rnorm',11x, & ' log10 log10'/ & 45x,' ynorm rnorm'/1x) 310 format (' ',i5,6x,2e11.3,11x,2f11.3) 320 format (1x/ & ' norms of solution and residual vectors for a range of values'/ & ' of the levenberg-marquardt parameter, lambda.'// & ' lambda ynorm rnorm', & ' log10 log10 log10'/ & 34x,' lambda ynorm rnorm') 330 format (1x, 3e11.3, 3f11.3) ! if ( m <= 0 .or. n <= 0 ) then return end if minmn = min ( m, n ) minmn1 = minmn + 1 if ( kpvec(1) == 0 ) then prblk = 111111 unit = -1 width = 79 else prblk = kpvec(2) unit = kpvec(3) width = kpvec(4) end if star = unit < 0 ! ! Build logical array blk() by testing decimal digits of PRBLK. ! do i = 6, 1, -1 j = prblk/10 blk(i) = (prblk - 10*j) > 0 prblk = j end do ! ! Optionally print the header and m, n, mdata ! if ( blk(1) ) then if ( star ) then write (*,270) write (*,260) m,n,mdata else write (unit,270) write (unit,260) m,n,mdata end if end if ! ! Handle scaling as selected by ISCALE. ! if ( iscale == 1 ) then if ( blk(1) ) then if ( star ) then write (*,290) else write (unit,290) end if end if else ! ! Apply column scaling to A. ! do j = 1, n a1 = d(j) if ( iscale <= 2 ) then sb = 0.0D+00 do i = 1, m sb = sb + a(i,j)**2 end do a1 = sqrt ( sb ) if ( a1 == 0.0D+00 ) a1 = 1.0D+00 a1 = 1.0D+00 / a1 d(j) = a1 end if do i = 1, m a(i,j) = a(i,j)*a1 end do end do if ( blk(1) ) then if ( star ) then write (*,280) iscale,(d(j),j = 1, n) else write (unit,280) iscale,(d(j),j = 1, n) end if end if end if ! ! Compute the singular value decomposition of the scaled matrix. ! call svdrs ( a, mda, m, n, b, m, 1, sing ) ! ! Determine nsol. ! nsol = minmn do j = 1, minmn if ( sing(j) == 0.0D+00 ) then nsol = j-1 exit end if end do ! ! The array b() contains the vector g. ! Compute cumulative sums of squares of components of ! g and store them in work(i), i = 1,...,minmn+1 ! sb = sum ( b(minmn1:m)**2 ) work(minmn+1) = sb do j = minmn, 1, -1 sb = sb + b(j)**2 work(j) = sb end do ! ! Print the V matrix. ! if ( blk(2) ) then call mfeout ( a, mda, n, n, names, 1, unit, width ) end if ! ! Replace v by d*v in the array a() ! if ( iscale > 1 ) then do i = 1, n a(i,1:n) = d(i) * a(i,1:n) end do end if if ( blk(3) ) then ! ! Print singular values and other summary results. ! ! Output will be done using one of two layouts. the narrow ! layout uses 69 columns + 1 for carriage control, and makes two passes ! through the computation. ! ! The wide layout uses 95 columns + 1 for carriage control, and makes ! only one pass through the computation. ! ! G is now in b() array. v now in a(,) array. ! narrow = width < 95 if ( narrow ) then mpass = 2 else mpass = 1 end if do 170 ipass = 1, mpass if ( star ) then if ( narrow ) then if ( ipass == 1 ) then write(*,221) else write(*,222) end if else write (*,220) end if else if ( narrow ) then if ( ipass == 1 ) then write(unit,221) else write(unit,222) end if else write (unit,220) end if end if ! ! The following statement converts from integer to floating-point. ! a3 = work(1) a4 = sqrt ( a3 / max ( 1, mdata ) ) if ( star ) then if ( narrow ) then if ( ipass == 1 ) then write(*,231) a4 else write(*,232) a3, a4 end if else write (*,230) a3,a4 end if else if ( narrow ) then if ( ipass == 1 ) then write(unit,231) a4 else write(unit,232) a3, a4 end if else write (unit,230) a3,a4 end if end if do k = 1, minmn if ( sing(k) == 0.0D+00 ) then pcoef = 0.0D+00 if ( star ) then write (*,240) k,sing(k) else write (unit,240) k,sing(k) end if else pcoef = b(k) / sing(k) a1 = 1.0D+00 / sing(k) a2 = b(k)**2 a3 = work(k+1) a4 = sqrt ( a3 / max ( 1, mdata-k ) ) if ( star ) then if ( narrow ) then if ( ipass == 1 ) then write(*,240) k,sing(k),pcoef,a1,b(k), a4 else write(*,240) k,sing(k), b(k),a2,a3,a4 end if else write (*,240) k,sing(k),pcoef,a1,b(k),a2,a3,a4 end if else if ( narrow ) then if ( ipass == 1 ) then write(unit,240) k,sing(k),pcoef,a1,b(k), a4 else write(unit,240) k,sing(k), b(k),a2,a3,a4 end if else write (unit,240) k,sing(k),pcoef,a1,b(k),a2,a3,a4 end if end if end if end do 170 continue end if if ( blk(4) ) then ! ! Compute and print values of ynorm, rnorm and their logarithms. ! if ( star ) then write (*,300) else write (unit,300) end if ysq = 0.0D+00 do j = 0, nsol if ( j /= 0 ) ysq = ysq + (b(j) / sing(j))**2 ynorm = sqrt(ysq) rnorm = sqrt(work(j+1)) yl = -1000.0D+00 if ( ynorm > 0.0D+00 ) yl = log10(ynorm) rl = -1000.0D+00 if ( rnorm > 0.0D+00 ) rl = log10(rnorm) if ( star ) then write (*,310) j,ynorm,rnorm,yl,rl else write (unit,310) j,ynorm,rnorm,yl,rl end if end do end if if ( blk(5) .and. sing(1) /= 0.0D+00 ) then ! ! Compute values of xnorm and rnorm for a sequence of values of ! the levenberg-marquardt parameter. ! el = log10 ( sing(1) ) + 1.0D+00 el2 = log10 ( sing(nsol) ) - 1.0D+00 del = ( el2 - el ) / 20.0D+00 aln10 = log ( 10.0D+00 ) if ( star ) then write (*,320) else write (unit,320) end if do ie = 1, 21 ! ! Compute alamb = 10.0D+00 **el ! alamb = exp(aln10*el) ys = 0.0D+00 rs = work(nsol+1) do i = 1, minmn sl = sing(i)**2 + alamb**2 ys = ys + (b(i)*sing(i)/sl)**2 rs = rs + (b(i)*(alamb**2)/sl)**2 end do ynorm = sqrt(ys) rnorm = sqrt(rs) rl = -1000.0D+00 if ( rnorm > 0.0D+00 ) rl = log10(rnorm) yl = -1000.0D+00 if ( ynorm > 0.0D+00 ) yl = log10(ynorm) if ( star ) then write (*,330) alamb,ynorm,rnorm,el,yl,rl else write (unit,330) alamb,ynorm,rnorm,el,yl,rl end if el = el + del end do end if ! ! Compute and optionally print candidate solutions. ! do k = 1, nsol pcoef = b(k) / sing(k) do i = 1, n a(i,k) = a(i,k) * pcoef if ( k > 1) a(i,k) = a(i,k) + a(i,k-1) end do end do if ( blk(6) .and. nsol >= 1 ) then call mfeout ( a, mda, n, nsol, names, 2, unit, width ) end if return end subroutine svdrs ( a, mda, m, n, b, mdb, nb, s ) ! !******************************************************************************* ! !! SVDRS: singular value decomposition with a right side vector. ! ! ! Discussion: ! ! This subroutine computes the singular value decomposition of the ! M by N matrix A and optionally applies the transformations ! from the left to the NB column vectors of the M by NB matrix B. ! ! Either M >= N or M < N is permitted. ! ! The singular value decomposition of A is of the form ! ! A = U * S * V' ! ! where U is M by M orthogonal, S is M by N diagonal with the ! diagonal terms nonnegative and ordered from large to small, and ! V is N by N orthogonal. Note that these matrices also satisfy ! ! S = U' * A * V ! ! The matrix V is returned in the leading N rows and ! columns of the array A. ! ! The singular values, i.e. the diagonal terms of the matrix S, ! are returned in the array S. If M < N, positions M+1 ! through N of S will be set to zero. ! ! The product matrix G = U' * B overwrites the input matrix B. ! ! If the user wishes to obtain a minimum length least squares ! solution of the linear system ! ! A * x ~ = ~ b ! ! the solution X can be constructed, following use of this subroutine, ! by computing the sum for I = 1, ..., R of the outer products ! ! (column I of V) * 1/S(I) * (row I of G) ! ! Here R denotes the pseudorank of A which the user may choose ! in the range 0 through min ( M, N ) based on the sizes of the ! singular values. ! ! This code gives special treatment to rows and columns that are ! entirely zero. This causes certain zero singular values to appear as ! exact zeros rather than as about macheps times the largest singular value. ! It similarly cleans up the associated columns of U and V. ! ! Algorithm: ! ! 1. Exchange columns of A to pack nonzero columns to the left. ! Set N_COPY = number of nonzero columns. ! Use locations A(1,N), A(1,N-1),..., A(1,N_COPY+1) to record the ! column permutations. ! ! 2. Exchange rows of A to pack nonzero rows to the top. ! Quit packing if find N_COPY nonzero rows. Make same row exchanges ! in B. Set M_COPY so that all nonzero rows of the permuted A ! are in first M_COPY rows. If M_COPY <= N_COPY then all M_COPY rows are ! nonzero. If M_COPY > N_COPY then the first N_COPY rows are known ! to be nonzero, and rows NCOPY+1 through M_COPY may be zero or nonzero. ! ! 3. Apply the original algorithm to the M_COPY by N_COPY problem. ! ! 4. Move permutation record from A(,) to S(I), I = N_COPY+1,...,N. ! ! 5. Build V up from N_COPY by N_COPY to N by N by placing ones on ! the diagonal and zeros elsewhere. This is only partly done ! explicitly. It is completed during step 6. ! ! 6. Exchange rows of V to compensate for column exchanges of step 2. ! ! 7. Place zeros in S(N_COPY+1:N) to represent zero singular values. ! ! Reference: ! ! Charles Lawson and Richard Hanson, ! Solving Least Squares Problems, ! Prentice-Hall, 1974, ! Revised edition, SIAM, 1995. ! QA275.L38 ! ! Parameters: ! ! Input/output, double precision A(MDA,N). ! On input, the M by N matrix A. ! On output contains the N by N matrix V. ! ! Input, integer MDA, the leading dimension of A. ! MDA must be at least the maximum of M and N. ! ! Input, integer M, the number of rows of A and B. ! M must be at least 1. ! ! Input, integer N, the number of columns of the matrix A, ! the number of rows and columns of the matrix V. ! N must be greater than 0. ! ! Input/output, double precision B(MDB,NB). If NB > 0, this array ! must contain an M by NB matrix on input and will contain the ! M by NB product matrix, G = U' * B on output. ! ! Input, integer LDB, the leading dimension of B, which must be ! at least M. ! ! Input, integer NB, the number of right hand sides; that is, the ! number of columns of data in B. If no right hand sides are being ! supplied, set NB = 0. ! ! Output, real S(N), the singular values of A, with the ordering ! S(1) >= S(2) >= ... >= S(N) >= 0. ! If M < N the singular values indexed from M+1 through N will be zero. ! ! Local parameters: ! ! double precision WORK(N,2). ! locations 1 thru N will hold the off-diagonal terms of ! the bidiagonal matrix for subroutine QRBD. Locations N+1 ! thru 2*N will save info from one call to the next of H12. ! implicit none ! integer mda integer mdb integer n integer nb ! double precision a(mda,n) double precision b(mdb,nb) integer i integer ipass integer j integer k integer l integer m integer m_copy integer mode integer n_copy integer ns integer nsp1 double precision s(n) double precision t double precision work(n,2) ! n_copy = n if ( m <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SVDRS - Fatal error!' write ( *, '(a,i6)' ) ' Illegal value of M = ', m return end if if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SVDRS - Fatal error!' write ( *, '(a,i6)' ) ' Illegal value of N = ', n return end if if ( mda < m .or. mda < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SVDRS - Fatal error!' write ( *, '(a,i6)' ) ' Illegal value of MDA = ', mda return end if if ( mdb < m ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SVDRS - Fatal error!' write ( *, '(a,i6)' ) ' Illegal value of MDB = ', mdb return end if ! ! If column J is entirely zero, exchange it with column N. ! j = n_copy do while ( j >= 1 ) if ( all ( a(1:m,j) == 0.0D+00 ) ) then if ( j /= n_copy ) then a(1:m,j) = a(1:m,n_copy) end if a(1,n_copy) = j n_copy = n_copy - 1 end if j = j - 1 end do ns = 0 ! ! If N_COPY = 0 then A is entirely zero and the SVD computation can be skipped. ! if ( n_copy /= 0 ) then ! ! Pack nonzero rows to the top. ! Quit packing if find N_COPY nonzero rows. ! i = 1 m_copy = m do while ( i <= n_copy .and. i < m_copy ) if ( a(i,i) == 0.0D+00 ) then if ( all ( a(i,1:n_copy) == 0.0D+00 ) ) then do j = 1, nb call d_swap ( b(i,j), b(m_copy,j) ) end do a(i,1:n_copy) = a(m_copy,1:n_copy) if ( m_copy <= n_copy ) then a(m_copy,1:n_copy) = 0.0D+00 end if m_copy = m_copy - 1 else i = i + 1 end if else i = i + 1 end if end do ! ! Begin the SVD algorithm: ! ! 1: Reduce the matrix to upper bidiagonal form with Householder ! transformations ! H(N_COPY)...H(1) * A * Q(1)...Q(N_COPY-2) = (D',0)' ! where D is upper bidiagonal. ! ! 2: Apply H(N_COPY)...H(1) to B. Here H(N_COPY)...H(1)*B replaces ! B in storage. ! ! 3: The matrix product W = Q(1)...Q(N_COPY-2) overwrites the first ! N_COPY rows of A in storage. ! ! 4: An SVD for D is computed. Here K rotations Ri and Pi are ! computed so that ! Rk...R1*D*P1'...Pk' = diag(S1,...,Sm_copy) ! to working accuracy. The S(I) are nonnegative and nonincreasing. ! Here Rk...R1*B overwrites B in storage while ! A*P1'...Pk' overwrites A in storage. ! ! 5: It follows that, with the proper definitions, U'*B overwrites B, ! while V overwrites the first N_COPY row and columns of A. ! l = min ( m_copy, n_copy ) ! ! The following loop reduces A to upper bidiagonal and ! also applies the premultiplying transformations to B. ! do j = 1, l ! ! Note: In some cases, the memory references may be illegal. But ! this should be OK. For instance, if J = N, then A(1,J+1) is ! an illegal memory reference, but in fact N_COPY - J will be 0, ! indicating that there will be no operations on that memory. ! Just a bit of careless programming from the old days. ! if ( j < m_copy ) then mode = 1 call h12 ( mode, j, j+1, m_copy, a(1,j), 1, t, a(1,j+1), 1, mda, & n_copy-j ) mode = 2 call h12 ( mode, j, j+1, m_copy, a(1,j), 1, t, b, 1, mdb, nb ) end if if ( j < n_copy-1 ) then mode = 1 call h12 ( mode, j+1, j+2, n_copy, a(j,1), mda, work(j,2), a(j+1,1), & mda, 1, m_copy-j ) end if end do ! ! Copy the bidiagonal matrix into S and WORK for QRBD. ! do j = 2, l s(j) = a(j,j) work(j,1) = a(j-1,j) end do s(1) = a(1,1) ns = n_copy if ( m_copy < n_copy ) then ns = m_copy + 1 s(ns) = 0.0D+00 work(ns,1) = a(m_copy,m_copy+1) end if ! ! Construct the explicit N_COPY by N_COPY product matrix, ! W = Q1*Q2*...*QL*I in A. ! do i = n_copy, 1, -1 if ( i <= min ( m_copy, n_copy-2 ) ) then mode = 2 call h12 ( mode, i+1, i+2, n_copy, a(i,1), mda, work(i,2), a(1,i+1), & 1, mda, n_copy-i ) end if a(i,1:n_copy) = 0.0D+00 a(i,i) = 1.0D+00 end do ! ! Compute the singular value decomposition of the bidiagonal matrix. ! call qrbd ( ipass, s(1), work(1,1), ns, a, mda, n_copy, b, mdb, nb ) if ( ipass == 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SVDRS - Warning:' write ( *, '(a)' ) ' Full accuracy was not attained in the singular' write ( *, '(a)' ) ' value decomposition of the bidiagonal matrix.' end if end if s(ns+1:n_copy) = 0.0D+00 ! ! Extract the record of permutations and store zeros. ! s(n_copy+1:n) = a(1,n_copy+1:n) a(1:n_copy,n_copy+1:n) = 0.0D+00 ! ! Permute rows and set zero singular values. ! do k = n_copy+1, n i = int ( s(k) ) s(k) = 0.0D+00 a(k,1:n) = a(i,1:n) a(i,1:n) = 0.0D+00 a(i,k) = 1.0D+00 end do return end subroutine timestamp ( ) ! !******************************************************************************* ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! character ( len = 8 ) ampm integer d character ( len = 8 ) date integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone ! call date_and_time ( date, time, zone, 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 ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end