subroutine daxpy ( n, da, dx, incx, dy, incy ) !*****************************************************************************80 ! !! DAXPY computes constant times a vector plus a vector. ! ! Discussion: ! ! This routine uses double precision real arithmetic. ! ! This routine uses unrolled loops for increments equal to one. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 May 2005 ! ! Author: ! ! Original FORTRAN77 version by Charles Lawson, Richard Hanson, ! David Kincaid, Fred Krogh. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Algorithm 539, ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of elements in DX and DY. ! ! Input, real ( kind = rk ) DA, the multiplier of DX. ! ! Input, real ( kind = rk ) DX(*), the first vector. ! ! Input, integer INCX, the increment between successive ! entries of DX. ! ! Input/output, real ( kind = rk ) DY(*), the second vector. ! On output, DY(*) has been replaced by DY(*) + DA * DX(*). ! ! Input, integer INCY, the increment between successive ! entries of DY. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) da real ( kind = rk ) dx(*) real ( kind = rk ) dy(*) integer i integer incx integer incy integer ix integer iy integer m integer n if ( n <= 0 ) then return end if if ( da == 0.0D+00 ) then return end if ! ! Code for unequal increments or equal increments ! not equal to 1. ! if ( incx /= 1 .or. incy /= 1 ) then if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n dy(iy) = dy(iy) + da * dx(ix) ix = ix + incx iy = iy + incy end do ! ! Code for both increments equal to 1. ! else m = mod ( n, 4 ) dy(1:m) = dy(1:m) + da * dx(1:m) do i = m+1, n, 4 dy(i ) = dy(i ) + da * dx(i ) dy(i+1) = dy(i+1) + da * dx(i+1) dy(i+2) = dy(i+2) + da * dx(i+2) dy(i+3) = dy(i+3) + da * dx(i+3) end do end if return end function ddot ( n, dx, incx, dy, incy ) !*****************************************************************************80 ! !! DDOT forms the dot product of two vectors. ! ! Discussion: ! ! This routine uses double precision real arithmetic. ! ! This routine uses unrolled loops for increments equal to one. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 May 2005 ! ! Author: ! ! Original FORTRAN77 version by Charles Lawson, Richard Hanson, ! David Kincaid, Fred Krogh. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Algorithm 539, ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of entries in the vectors. ! ! Input, real ( kind = rk ) DX(*), the first vector. ! ! Input, integer INCX, the increment between successive ! entries in DX. ! ! Input, real ( kind = rk ) DY(*), the second vector. ! ! Input, integer INCY, the increment between successive ! entries in DY. ! ! Output, real ( kind = rk ) DDOT, the sum of the product of the ! corresponding entries of DX and DY. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) ddot real ( kind = rk ) dtemp real ( kind = rk ) dx(*) real ( kind = rk ) dy(*) integer i integer incx integer incy integer ix integer iy integer m integer n ddot = 0.0D+00 dtemp = 0.0D+00 if ( n <= 0 ) then return end if ! ! Code for unequal increments or equal increments ! not equal to 1. ! if ( incx /= 1 .or. incy /= 1 ) then if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n dtemp = dtemp + dx(ix) * dy(iy) ix = ix + incx iy = iy + incy end do ! ! Code for both increments equal to 1. ! else m = mod ( n, 5 ) do i = 1, m dtemp = dtemp + dx(i) * dy(i) end do do i = m+1, n, 5 dtemp = dtemp + dx(i ) * dy(i ) & + dx(i+1) * dy(i+1) & + dx(i+2) * dy(i+2) & + dx(i+3) * dy(i+3) & + dx(i+4) * dy(i+4) end do end if ddot = dtemp return end function dnrm2 ( n, x, incx ) !*****************************************************************************80 ! !! DNRM2 returns the euclidean norm of a vector. ! ! Discussion: ! ! This routine uses double precision real arithmetic. ! ! DNRM2 ( X ) = sqrt ( X' * X ) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 May 2005 ! ! Author: ! ! Original FORTRAN77 version by Charles Lawson, Richard Hanson, ! David Kincaid, Fred Krogh. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Algorithm 539, ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real ( kind = rk ) X(*), the vector whose norm is to be computed. ! ! Input, integer INCX, the increment between successive ! entries of X. ! ! Output, real ( kind = rk ) DNRM2, the Euclidean norm of X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) absxi real ( kind = rk ) dnrm2 integer incx integer ix integer n real ( kind = rk ) norm real ( kind = rk ) scale real ( kind = rk ) ssq real ( kind = rk ) x(*) if ( n < 1 .or. incx < 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, 1 + ( n - 1 )*incx, incx if ( x(ix) /= 0.0D+00 ) then absxi = abs ( x(ix) ) if ( scale < 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 dnrm2 = norm return end subroutine drot ( n, x, incx, y, incy, c, s ) !*****************************************************************************80 ! !! DROT applies a plane rotation. ! ! Discussion: ! ! This routine uses double precision real arithmetic. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 April 1999 ! ! Author: ! ! Original FORTRAN77 version by Charles Lawson, Richard Hanson, ! David Kincaid, Fred Krogh. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Algorithm 539, ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of entries in the vectors. ! ! Input/output, real ( kind = rk ) X(*), one of the vectors to be rotated. ! ! Input, integer INCX, the increment between successive ! entries of X. ! ! Input/output, real ( kind = rk ) Y(*), one of the vectors to be rotated. ! ! Input, integer INCY, the increment between successive ! elements of Y. ! ! Input, real ( kind = rk ) C, S, parameters (presumably the cosine and ! sine of some angle) that define a plane rotation. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) c integer i integer incx integer incy integer ix integer iy integer n real ( kind = rk ) s real ( kind = rk ) stemp real ( kind = rk ) x(*) real ( kind = rk ) y(*) if ( n <= 0 ) then else if ( incx == 1 .and. incy == 1 ) then do i = 1, n stemp = c * x(i) + s * y(i) y(i) = c * y(i) - s * x(i) x(i) = stemp end do else if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n stemp = c * x(ix) + s * y(iy) y(iy) = c * y(iy) - s * x(ix) x(ix) = stemp ix = ix + incx iy = iy + incy end do end if return end subroutine drotg ( sa, sb, c, s ) !*****************************************************************************80 ! !! DROTG constructs a Givens plane rotation. ! ! Discussion: ! ! This routine uses double precision real arithmetic. ! ! Given values A and B, this routine computes ! ! SIGMA = sign ( A ) if abs ( A ) > abs ( B ) ! = sign ( B ) if abs ( A ) <= abs ( B ); ! ! R = SIGMA * ( A * A + B * B ); ! ! C = A / R if R is not 0 ! = 1 if R is 0; ! ! S = B / R if R is not 0, ! 0 if R is 0. ! ! The computed numbers then satisfy the equation ! ! ( C S ) ( A ) = ( R ) ! ( -S C ) ( B ) = ( 0 ) ! ! The routine also computes ! ! Z = S if abs ( A ) > abs ( B ), ! = 1 / C if abs ( A ) <= abs ( B ) and C is not 0, ! = 1 if C is 0. ! ! The single value Z encodes C and S, and hence the rotation: ! ! If Z = 1, set C = 0 and S = 1; ! If abs ( Z ) < 1, set C = sqrt ( 1 - Z * Z ) and S = Z; ! if abs ( Z ) > 1, set C = 1/ Z and S = sqrt ( 1 - C * C ); ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 May 2006 ! ! Author: ! ! Original FORTRAN77 version by Charles Lawson, Richard Hanson, ! David Kincaid, Fred Krogh. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Algorithm 539, ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input/output, real ( kind = rk ) SA, SB. On input, SA and SB are the values ! A and B. On output, SA is overwritten with R, and SB is ! overwritten with Z. ! ! Output, real ( kind = rk ) C, S, the cosine and sine of the ! Givens rotation. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) c real ( kind = rk ) r real ( kind = rk ) roe real ( kind = rk ) s real ( kind = rk ) sa real ( kind = rk ) sb real ( kind = rk ) scale real ( kind = rk ) z if ( abs ( sb ) < abs ( sa ) ) then roe = sa else roe = sb end if scale = abs ( sa ) + abs ( sb ) if ( scale == 0.0D+00 ) then c = 1.0D+00 s = 0.0D+00 r = 0.0D+00 else r = scale * sqrt ( ( sa / scale )**2 + ( sb / scale )**2 ) r = sign ( 1.0D+00, roe ) * r c = sa / r s = sb / r end if if ( 0.0D+00 < abs ( c ) .and. abs ( c ) <= s ) then z = 1.0D+00 / c else z = s end if sa = r sb = z return end subroutine dscal ( n, sa, x, incx ) !*****************************************************************************80 ! !! DSCAL scales a vector by a constant. ! ! Discussion: ! ! This routine uses double precision real arithmetic. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 April 1999 ! ! Author: ! ! Original FORTRAN77 version by Charles Lawson, Richard Hanson, ! David Kincaid, Fred Krogh. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Algorithm 539, ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real ( kind = rk ) SA, the multiplier. ! ! Input/output, real ( kind = rk ) X(*), the vector to be scaled. ! ! Input, integer INCX, the increment between successive ! entries of X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i integer incx integer ix integer m integer n real ( kind = rk ) sa real ( kind = rk ) x(*) if ( n <= 0 ) then else if ( incx == 1 ) then m = mod ( n, 5 ) x(1:m) = sa * x(1:m) do i = m+1, n, 5 x(i) = sa * x(i) x(i+1) = sa * x(i+1) x(i+2) = sa * x(i+2) x(i+3) = sa * x(i+3) x(i+4) = sa * x(i+4) end do else if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if do i = 1, n x(ix) = sa * x(ix) ix = ix + incx end do end if return end subroutine dsvdc ( a, lda, m, n, s, e, u, ldu, v, ldv, work, job, info ) !*****************************************************************************80 ! !! DSVDC computes the singular value decomposition of a real rectangular matrix. ! ! Discussion: ! ! This routine reduces an M by N matrix A to diagonal form by orthogonal ! transformations U and V. The diagonal elements S(I) are the singular ! values of A. The columns of U are the corresponding left singular ! vectors, and the columns of V the right singular vectors. ! ! The form of the singular value decomposition is then ! ! A(MxN) = U(MxM) * S(MxN) * V(NxN)' ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 September 2006 ! ! Author: ! ! Original FORTRAN77 version by Jack Dongarra, Jim Bunch, Cleve Moler, ! Pete Stewart. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Parameters: ! ! Input/output, real ( kind = rk ) A(LDA,N). On input, the M by N ! matrix whose singular value decomposition is to be computed. ! On output, the matrix has been destroyed. Depending on the user's ! requests, the matrix may contain other useful information. ! ! Input, integer LDA, the leading dimension of the array A. ! LDA must be at least N. ! ! Input, integer M, the number of rows of the matrix. ! ! Input, integer N, the number of columns of the matrix A. ! ! Output, real ( kind = rk ) S(MM), where MM = max(M+1,N). The first ! min(M,N) entries of S contain the singular values of A arranged in ! descending order of magnitude. ! ! Output, real ( kind = rk ) E(MM), where MM = max(M+1,N). Ordinarily ! contains zeros. However see the discussion of INFO for exceptions. ! ! Output, real ( kind = rk ) U(LDU,K). If JOBA = 1 then K = M; ! if 2 <= JOBA, then K = min(M,N). U contains the M by M matrix of ! left singular vectors. U is not referenced if JOBA = 0. If M <= N ! or if JOBA = 2, then U may be identified with A in the subroutine call. ! ! Input, integer LDU, the leading dimension of the array U. ! LDU must be at least M. ! ! Output, real ( kind = rk ) V(LDV,N), the N by N matrix of right singular ! vectors. V is not referenced if JOB is 0. If N <= M, then V may be ! identified with A in the subroutine call. ! ! Input, integer LDV, the leading dimension of the array V. ! LDV must be at least N. ! ! Workspace, real ( kind = rk ) WORK(M). ! ! Input, integer JOB, controls the computation of the singular ! vectors. It has the decimal expansion AB with the following meaning: ! A = 0, do not compute the left singular vectors. ! A = 1, return the M left singular vectors in U. ! A >= 2, return the first min(M,N) singular vectors in U. ! B = 0, do not compute the right singular vectors. ! B = 1, return the right singular vectors in V. ! ! Output, integer INFO, status indicator. ! The singular values (and their corresponding singular vectors) ! S(INFO+1), S(INFO+2),...,S(MN) are correct. Here MN = min ( M, N ). ! Thus if INFO is 0, all the singular values and their vectors are ! correct. In any event, the matrix B = U' * A * V is the bidiagonal ! matrix with the elements of S on its diagonal and the elements of E on ! its superdiagonal. Thus the singular values of A and B are the same. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer lda integer ldu integer ldv integer m integer n real ( kind = rk ) a(lda,n) real ( kind = rk ) b real ( kind = rk ) c real ( kind = rk ) cs real ( kind = rk ) e(*) real ( kind = rk ) el real ( kind = rk ) emm1 real ( kind = rk ) f real ( kind = rk ) g integer info integer iter integer j integer job integer jobu integer k integer kase integer kk integer l integer ll integer lls integer ls integer lu integer, parameter :: maxit = 30 integer mm integer mm1 integer mn integer nct integer nctp1 integer ncu integer nrt integer nrtp1 real ( kind = rk ) s(*) real ( kind = rk ) scale real ( kind = rk ) ddot real ( kind = rk ) shift real ( kind = rk ) sl real ( kind = rk ) sm real ( kind = rk ) smm1 real ( kind = rk ) sn real ( kind = rk ) dnrm2 real ( kind = rk ) t real ( kind = rk ) t1 real ( kind = rk ) test real ( kind = rk ) u(ldu,m) real ( kind = rk ) v(ldv,n) logical wantu logical wantv real ( kind = rk ) work(m) real ( kind = rk ) ztest ! ! Determine what is to be computed. ! wantu = .false. wantv = .false. jobu = mod ( job, 100 ) / 10 if ( 1 < jobu ) then ncu = min ( m, n ) else ncu = m end if if ( jobu /= 0 ) then wantu = .true. end if if ( mod ( job, 10 ) /= 0 ) then wantv = .true. end if ! ! Reduce A to bidiagonal form, storing the diagonal elements ! in S and the super-diagonal elements in E. ! info = 0 nct = min ( m-1, n ) nrt = max ( 0, min ( m, n-2 ) ) lu = max ( nct, nrt ) do l = 1, lu ! ! Compute the transformation for the L-th column and ! place the L-th diagonal in S(L). ! if ( l <= nct ) then s(l) = dnrm2 ( m-l+1, a(l,l), 1 ) if ( s(l) /= 0.0D+00 ) then if ( a(l,l) /= 0.0D+00 ) then s(l) = sign ( s(l), a(l,l) ) end if call dscal ( m-l+1, 1.0D+00 / s(l), a(l,l), 1 ) a(l,l) = 1.0D+00 + a(l,l) end if s(l) = -s(l) end if do j = l+1, n ! ! Apply the transformation. ! if ( l <= nct .and. s(l) /= 0.0D+00 ) then t = -ddot ( m-l+1, a(l,l), 1, a(l,j), 1 ) / a(l,l) call daxpy ( m-l+1, t, a(l,l), 1, a(l,j), 1 ) end if ! ! Place the L-th row of A into E for the ! subsequent calculation of the row transformation. ! e(j) = a(l,j) end do ! ! Place the transformation in U for subsequent back multiplication. ! if ( wantu .and. l <= nct ) then u(l:m,l) = a(l:m,l) end if ! ! Compute the L-th row transformation and place the ! L-th superdiagonal in E(L). ! if ( l <= nrt ) then e(l) = dnrm2 ( n-l, e(l+1), 1 ) if ( e(l) /= 0.0D+00 ) then if ( e(l+1) /= 0.0D+00 ) then e(l) = sign ( e(l), e(l+1) ) end if call dscal ( n-l, 1.0D+00 / e(l), e(l+1), 1 ) e(l+1) = 1.0D+00 + e(l+1) end if e(l) = -e(l) ! ! Apply the transformation. ! if ( l + 1 <= m .and. e(l) /= 0.0D+00 ) then work(l+1:m) = 0.0D+00 do j = l+1, n call daxpy ( m-l, e(j), a(l+1,j), 1, work(l+1), 1 ) end do do j = l+1, n call daxpy ( m-l, -e(j)/e(l+1), work(l+1), 1, a(l+1,j), 1 ) end do end if ! ! Place the transformation in V for subsequent back multiplication. ! if ( wantv ) then v(l+1:n,l) = e(l+1:n) end if end if end do ! ! Set up the final bidiagonal matrix of order MN. ! mn = min ( m + 1, n ) nctp1 = nct + 1 nrtp1 = nrt + 1 if ( nct < n ) then s(nctp1) = a(nctp1,nctp1) end if if ( m < mn ) then s(mn) = 0.0D+00 end if if ( nrtp1 < mn ) then e(nrtp1) = a(nrtp1,mn) end if e(mn) = 0.0D+00 ! ! If required, generate U. ! if ( wantu ) then u(1:m,nctp1:ncu) = 0.0D+00 do j = nctp1, ncu u(j,j) = 1.0D+00 end do do ll = 1, nct l = nct - ll + 1 if ( s(l) /= 0.0D+00 ) then do j = l+1, ncu t = -ddot ( m-l+1, u(l,l), 1, u(l,j), 1 ) / u(l,l) call daxpy ( m-l+1, t, u(l,l), 1, u(l,j), 1 ) end do u(l:m,l) = -u(l:m,l) u(l,l) = 1.0D+00 + u(l,l) u(1:l-1,l) = 0.0D+00 else u(1:m,l) = 0.0D+00 u(l,l) = 1.0D+00 end if end do end if ! ! If it is required, generate V. ! if ( wantv ) then do ll = 1, n l = n - ll + 1 if ( l <= nrt .and. e(l) /= 0.0D+00 ) then do j = l + 1, n t = -ddot ( n-l, v(l+1,l), 1, v(l+1,j), 1 ) / v(l+1,l) call daxpy ( n-l, t, v(l+1,l), 1, v(l+1,j), 1 ) end do end if v(1:n,l) = 0.0D+00 v(l,l) = 1.0D+00 end do end if ! ! Main iteration loop for the singular values. ! mm = mn iter = 0 do while ( 0 < mn ) ! ! If too many iterations have been performed, set flag and return. ! if ( maxit <= iter ) then info = mn return end if ! ! This section of the program inspects for ! negligible elements in the S and E arrays. ! ! On completion the variables KASE and L are set as follows: ! ! KASE = 1 if S(MN) and E(L-1) are negligible and L < MN ! KASE = 2 if S(L) is negligible and L < MN ! KASE = 3 if E(L-1) is negligible, L < MN, and ! S(L), ..., S(MN) are not negligible (QR step). ! KASE = 4 if E(MN-1) is negligible (convergence). ! do ll = 1, mn l = mn - ll if ( l == 0 ) then exit end if test = abs ( s(l) ) + abs ( s(l+1) ) ztest = test + abs ( e(l) ) if ( ztest == test ) then e(l) = 0.0D+00 exit end if end do if ( l == mn - 1 ) then kase = 4 else do lls = l + 1, mn + 1 ls = mn - lls + l + 1 if ( ls == l ) then exit end if test = 0.0D+00 if ( ls /= mn ) then test = test + abs ( e(ls) ) end if if ( ls /= l + 1 ) then test = test + abs ( e(ls-1) ) end if ztest = test + abs ( s(ls) ) if ( ztest == test ) then s(ls) = 0.0D+00 exit end if end do if ( ls == l ) then kase = 3 else if ( ls == mn ) then kase = 1 else kase = 2 l = ls end if end if l = l + 1 ! ! Deflate negligible S(MN). ! if ( kase == 1 ) then mm1 = mn - 1 f = e(mn-1) e(mn-1) = 0.0D+00 do kk = l, mm1 k = mm1 - kk + l t1 = s(k) call drotg ( t1, f, cs, sn ) s(k) = t1 if ( k /= l ) then f = -sn * e(k-1) e(k-1) = cs * e(k-1) end if if ( wantv ) then call drot ( n, v(1,k), 1, v(1,mn), 1, cs, sn ) end if end do ! ! Split at negligible S(L). ! else if ( kase == 2 ) then f = e(l-1) e(l-1) = 0.0D+00 do k = l, mn t1 = s(k) call drotg ( t1, f, cs, sn ) s(k) = t1 f = -sn * e(k) e(k) = cs * e(k) if ( wantu ) then call drot ( m, u(1,k), 1, u(1,l-1), 1, cs, sn ) end if end do ! ! Perform one QR step. ! else if ( kase == 3 ) then ! ! Calculate the shift. ! scale = max ( abs ( s(mn) ), abs ( s(mn-1) ), abs ( e(mn-1) ), & abs ( s(l) ), abs ( e(l) ) ) sm = s(mn) / scale smm1 = s(mn-1) / scale emm1 = e(mn-1) / scale sl = s(l) / scale el = e(l) / scale b = ( ( smm1 + sm ) * ( smm1 - sm ) + emm1 * emm1 ) / 2.0D+00 c = sm * sm * emm1 * emm1 shift = 0.0D+00 if ( b /= 0.0D+00 .or. c /= 0.0D+00 ) then shift = sqrt ( b * b + c ) if ( b < 0.0D+00 ) then shift = -shift end if shift = c / ( b + shift ) end if f = ( sl + sm ) * ( sl - sm ) + shift g = sl * el ! ! Chase zeros. ! mm1 = mn - 1 do k = l, mm1 call drotg ( f, g, cs, sn ) if ( k /= l ) then e(k-1) = f end if f = cs * s(k) + sn * e(k) e(k) = cs * e(k) - sn * s(k) g = sn * s(k+1) s(k+1) = cs * s(k+1) if ( wantv ) then call drot ( n, v(1,k), 1, v(1,k+1), 1, cs, sn ) end if call drotg ( f, g, cs, sn ) s(k) = f f = cs * e(k) + sn * s(k+1) s(k+1) = -sn * e(k) + cs * s(k+1) g = sn * e(k+1) e(k+1) = cs * e(k+1) if ( wantu .and. k < m ) then call drot ( m, u(1,k), 1, u(1,k+1), 1, cs, sn ) end if end do e(mn-1) = f iter = iter + 1 ! ! Convergence. ! else if ( kase == 4 ) then ! ! Make the singular value nonnegative. ! if ( s(l) < 0.0D+00 ) then s(l) = -s(l) if ( wantv ) then v(1:n,l) = -v(1:n,l) end if end if ! ! Order the singular value. ! do if ( l == mm ) then exit end if if ( s(l+1) <= s(l) ) then exit end if t = s(l) s(l) = s(l+1) s(l+1) = t if ( wantv .and. l < n ) then call dswap ( n, v(1,l), 1, v(1,l+1), 1 ) end if if ( wantu .and. l < m ) then call dswap ( m, u(1,l), 1, u(1,l+1), 1 ) end if l = l + 1 end do iter = 0 mn = mn - 1 end if end do return end subroutine dswap ( n, x, incx, y, incy ) !*****************************************************************************80 ! !! DSWAP interchanges two vectors. ! ! Discussion: ! ! This routine uses double precision real arithmetic. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 April 1999 ! ! Author: ! ! Original FORTRAN77 version by Charles Lawson, Richard Hanson, ! David Kincaid, Fred Krogh. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Algorithm 539, ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of entries in the vectors. ! ! Input/output, real ( kind = rk ) X(*), one of the vectors to swap. ! ! Input, integer INCX, the increment between successive ! entries of X. ! ! Input/output, real ( kind = rk ) Y(*), one of the vectors to swap. ! ! Input, integer INCY, the increment between successive ! elements of Y. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i integer incx integer incy integer ix integer iy integer m integer n real ( kind = rk ) temp real ( kind = rk ) x(*) real ( kind = rk ) y(*) if ( n <= 0 ) then else if ( incx == 1 .and. incy == 1 ) then m = mod ( n, 3 ) do i = 1, m temp = x(i) x(i) = y(i) y(i) = temp end do do i = m+1, n, 3 temp = x(i) x(i) = y(i) y(i) = temp temp = x(i+1) x(i+1) = y(i+1) y(i+1) = temp temp = x(i+2) x(i+2) = y(i+2) y(i+2) = temp end do else if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n temp = x(ix) x(ix) = y(iy) y(iy) = temp ix = ix + incx iy = iy + incy end do end if return end subroutine phi1 ( n, r, r0, v ) !*****************************************************************************80 ! !! PHI1 evaluates the multiquadric radial basis function. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 June 2012 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! William Press, Brian Flannery, Saul Teukolsky, William Vetterling, ! Numerical Recipes in FORTRAN: The Art of Scientific Computing, ! Third Edition, ! Cambridge University Press, 2007, ! ISBN13: 978-0-521-88068-8, ! LC: QA297.N866. ! ! Parameters: ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) R(N), the radial separation. ! 0 < R. ! ! Input, real ( kind = rk ) R0, a scale factor. ! ! Output, real ( kind = rk ) V(N), the value of the radial basis function. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) r(n) real ( kind = rk ) r0 real ( kind = rk ) v(n) v(1:n) = sqrt ( r(1:n) ** 2 + r0 ** 2 ) return end subroutine phi2 ( n, r, r0, v ) !*****************************************************************************80 ! !! PHI2 evaluates the inverse multiquadric radial basis function. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 June 2012 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! William Press, Brian Flannery, Saul Teukolsky, William Vetterling, ! Numerical Recipes in FORTRAN: The Art of Scientific Computing, ! Third Edition, ! Cambridge University Press, 2007, ! ISBN13: 978-0-521-88068-8, ! LC: QA297.N866. ! ! Parameters: ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) R(N), the radial separation. ! 0 < R. ! ! Input, real ( kind = rk ) R0, a scale factor. ! ! Output, real ( kind = rk ) V(N), the value of the radial basis function. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) r(n) real ( kind = rk ) r0 real ( kind = rk ) v(n) v = 1.0D+00 / sqrt ( r ** 2 + r0 ** 2 ) return end subroutine phi3 ( n, r, r0, v ) !*****************************************************************************80 ! !! PHI3 evaluates the thin-plate spline radial basis function. ! ! Discussion: ! ! Note that PHI3(R,R0) is negative if R < R0. Thus, for this basis function, ! it may be desirable to choose a value of R0 smaller than any possible R. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 June 2012 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! William Press, Brian Flannery, Saul Teukolsky, William Vetterling, ! Numerical Recipes in FORTRAN: The Art of Scientific Computing, ! Third Edition, ! Cambridge University Press, 2007, ! ISBN13: 978-0-521-88068-8, ! LC: QA297.N866. ! ! Parameters: ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) R(N), the radial separation. ! 0 < R. ! ! Input, real ( kind = rk ) R0, a scale factor. ! ! Output, real ( kind = rk ) V(N), the value of the radial basis function. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer i real ( kind = rk ) r(n) real ( kind = rk ) r0 real ( kind = rk ) v(n) do i = 1, n if ( r(i) <= 0.0D+00 ) then v(i) = 0.0D+00 else v(i) = r(i) ** 2 * log ( r(i) / r0 ) end if end do return end subroutine phi4 ( n, r, r0, v ) !*****************************************************************************80 ! !! PHI4 evaluates the gaussian radial basis function. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 June 2012 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! William Press, Brian Flannery, Saul Teukolsky, William Vetterling, ! Numerical Recipes in FORTRAN: The Art of Scientific Computing, ! Third Edition, ! Cambridge University Press, 2007, ! ISBN13: 978-0-521-88068-8, ! LC: QA297.N866. ! ! Parameters: ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) R(N), the radial separation. ! 0 < R. ! ! Input, real ( kind = rk ) R0, a scale factor. ! ! Output, real ( kind = rk ) V(N), the value of the radial basis function. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) r(n) real ( kind = rk ) r0 real ( kind = rk ) v(n) v(1:n) = exp ( - 0.5D+00 * r(1:n) ** 2 / r0 ** 2 ) return end subroutine phi5 ( n, r, r0, v ) !*****************************************************************************80 ! !! PHI5 evaluates a Wendland radial basis function. ! ! Discussion: ! ! Since each basis function is zero beyond a radius of R0, care must ! be taken so that every point in the region is no further than R0 from ! a basis point! ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 July 2018 ! ! Author: ! ! John Burkardt ! ! Reference ! ! Holger Wendland, ! Piecewise polynomial, positive definite and compactly supported radial ! functions of minimal degree, ! Advances in Computational Mathematics, ! Volume 4, Number 1, December 1995, pages 389–396. ! ! Parameters: ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) R(N), the radial separation. ! 0 < R. ! ! Input, real ( kind = rk ) R0, a scale factor. ! ! Output, real ( kind = rk ) V(N), the value of the radial basis function. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) r(n) real ( kind = rk ) r0 real ( kind = rk ) v(n) v(1:n) = max ( & 1.0D+00 - ( r(1:n) / r0 ) ** 2, & 0.0D+00 ) return end subroutine phi6 ( n, r, r0, v ) !*****************************************************************************80 ! !! PHI6 evaluates a Wendland radial basis function. ! ! Discussion: ! ! Since each basis function is zero beyond a radius of R0, care must ! be taken so that every point in the region is no further than R0 from ! a basis point! ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 July 2018 ! ! Author: ! ! John Burkardt ! ! Reference ! ! Holger Wendland, ! Piecewise polynomial, positive definite and compactly supported radial ! functions of minimal degree, ! Advances in Computational Mathematics, ! Volume 4, Number 1, December 1995, pages 389–396. ! ! Parameters: ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) R(N), the radial separation. ! 0 < R. ! ! Input, real ( kind = rk ) R0, a scale factor. ! ! Output, real ( kind = rk ) V(N), the value of the radial basis function. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) r(n) real ( kind = rk ) r0 real ( kind = rk ) v(n) v(1:n) = max ( & 1.0D+00 - ( r(1:n) / r0 ) ** 3, & 0.0D+00 ) return end subroutine phi7 ( n, r, r0, v ) !*****************************************************************************80 ! !! PHI7 evaluates a Wendland radial basis function. ! ! Discussion: ! ! Since each basis function is zero beyond a radius of R0, care must ! be taken so that every point in the region is no further than R0 from ! a basis point! ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 03 July 2018 ! ! Author: ! ! John Burkardt ! ! Reference ! ! Holger Wendland, ! Piecewise polynomial, positive definite and compactly supported radial ! functions of minimal degree, ! Advances in Computational Mathematics, ! Volume 4, Number 1, December 1995, pages 389–396. ! ! Parameters: ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) R(N), the radial separation. ! 0 < R. ! ! Input, real ( kind = rk ) R0, a scale factor. ! ! Output, real ( kind = rk ) V(N), the value of the radial basis function. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) r(n) real ( kind = rk ) r0 real ( kind = rk ) v(n) v(1:n) = max ( & 1.0D+00 - ( r(1:n) / r0 ) ** 4 * ( 4.0D+00 * r(1:n) / r0 + 1.0D+00 ), & 0.0D+00 ) return end subroutine r8mat_solve_svd ( m, n, a, b, x ) !*****************************************************************************80 ! !! R8MAT_SOLVE_SVD solves a linear system A*x=b using the SVD. ! ! Discussion: ! ! When the system is determined, the solution is the solution in the ! ordinary sense, and A*x = b. ! ! When the system is overdetermined, the solution minimizes the ! L2 norm of the residual ||A*x-b||. ! ! When the system is underdetermined, ||A*x-b|| should be zero, and ! the solution is the solution of minimum L2 norm, ||x||. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 June 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns ! in the matrix A. ! ! Input, real ( kind = rk ) A(M,N), the matrix. ! ! Input, real ( kind = rk ) B(M), the right hand side. ! ! Output, real ( kind = rk ) X(N), the solution. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) a_copy(m,n) real ( kind = rk ) a_pseudo(n,m) real ( kind = rk ) b(m) real ( kind = rk ) e(max(m+1,n)) integer i integer info integer lda integer ldu integer ldv integer job real ( kind = rk ) s(m,n) real ( kind = rk ) sp(n,m) real ( kind = rk ) sdiag(max(m+1,n)) real ( kind = rk ) u(m,m) real ( kind = rk ) v(n,n) real ( kind = rk ) work(m) real ( kind = rk ) x(n) ! ! Compute the SVD decomposition. ! job = 11 lda = m ldu = m ldv = n a_copy(1:m,1:n) = a(1:m,1:n) call dsvdc ( a_copy, lda, m, n, sdiag, e, u, ldu, v, ldv, work, job, info ) if ( info /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_SOLVE_SVD - Failure!' write ( *, '(a)' ) ' The SVD could not be calculated.' write ( *, '(a)' ) ' LINPACK routine DSVDC returned a nonzero' write ( *, '(a,i8)' ) ' value of the error flag, INFO = ', info return end if s(1:m,1:n) = 0.0D+00 do i = 1, min ( m, n ) s(i,i) = sdiag(i) end do ! ! Compute the pseudo inverse. ! sp(1:n,1:m) = 0.0D+00 do i = 1, min ( m, n ) if ( s(i,i) /= 0.0D+00 ) then sp(i,i) = 1.0D+00 / s(i,i) end if end do a_pseudo(1:n,1:m) = matmul ( v(1:n,1:n), & matmul ( sp(1:n,1:m), transpose ( u(1:m,1:m) ) ) ) ! ! Compute x = A_pseudo * b ! x(1:n) = matmul ( a_pseudo(1:n,1:m), b(1:m) ) return end subroutine rbf_interp ( m, nd, xd, r0, phi, w, ni, xi, fi ) !*****************************************************************************80 ! !! RBF_INTERP evaluates a radial basis function interpolant. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 June 2012 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! William Press, Brian Flannery, Saul Teukolsky, William Vetterling, ! Numerical Recipes in FORTRAN: The Art of Scientific Computing, ! Third Edition, ! Cambridge University Press, 2007, ! ISBN13: 978-0-521-88068-8, ! LC: QA297.N866. ! ! Parameters: ! ! Input, integer M, the spatial dimension. ! ! Input, integer ND, the number of data points. ! ! Input, real ( kind = rk ) XD(M,ND), the data points. ! ! Input, real ( kind = rk ) R0, a scale factor. R0 should be larger than ! the typical separation between points, but smaller than the maximum ! separation. The value of R0 has a significant effect on the resulting ! interpolant. ! ! Input, subroutine PHI ( N, R, R0, V ), a subroutine to evaluate the radial ! basis functions. ! ! Input, real ( kind = rk ) W(ND), the weights, as computed by RBF_WEIGHTS. ! ! Input, integer NI, the number of interpolation points. ! ! Input, real ( kind = rk ) XI(M,NI), the interpolation points. ! ! Output, real ( kind = rk ) FI(NI), the interpolated values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer nd integer ni real ( kind = rk ) fi(ni) integer i integer j external phi real ( kind = rk ) r(nd) real ( kind = rk ) r0 real ( kind = rk ) v(nd) real ( kind = rk ) w(nd) real ( kind = rk ) xd(m,nd) real ( kind = rk ) xi(m,ni) do i = 1, ni do j = 1, nd r(j) = sqrt ( sum ( ( xi(1:m,i) - xd(1:m,j) ) ** 2 ) ) end do call phi ( nd, r, r0, v ) fi(i) = dot_product ( v, w ) end do return end subroutine rbf_weight ( m, nd, xd, r0, phi, fd, w ) !*****************************************************************************80 ! !! RBF_WEIGHT computes weights for radial basis function interpolation. ! ! Discussion: ! ! We assume that there are N (nonsingular) equations in N unknowns. ! ! However, it should be clear that, if we are willing to do some kind ! of least squares calculation, we could allow for singularity, ! inconsistency, or underdetermine systems. This could be associated ! with data points that are very close or repeated, a smaller number ! of data points than function values, or some other ill-conditioning ! of the system arising from a peculiarity in the point spacing. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 June 2012 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! William Press, Brian Flannery, Saul Teukolsky, William Vetterling, ! Numerical Recipes in FORTRAN: The Art of Scientific Computing, ! Third Edition, ! Cambridge University Press, 2007, ! ISBN13: 978-0-521-88068-8, ! LC: QA297.N866. ! ! Parameters: ! ! Input, integer M, the spatial dimension. ! ! Input, integer ND, the number of data points. ! ! Input, real ( kind = rk ) XD(M,ND), the data points. ! ! Input, real ( kind = rk ) R0, a scale factor. R0 should be larger than ! the typical separation between points, but smaller than the maximum ! separation. The value of R0 has a significant effect on the resulting ! interpolant. ! ! Input, subroutine PHI ( N, R, R0, V ), a subroutine to evaluate the radial ! basis functions. ! ! Input, real ( kind = rk ) FD(ND), the function values at the data points. ! ! Output, real ( kind = rk ) W(ND), the weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer nd real ( kind = rk ) a(nd,nd) real ( kind = rk ) fd(nd) integer i integer j external phi real ( kind = rk ) r(nd) real ( kind = rk ) r0 real ( kind = rk ) v(nd) real ( kind = rk ) w(nd) real ( kind = rk ) xd(m,nd) do i = 1, nd do j = 1, nd r(j) = sqrt ( sum ( ( xd(1:m,i) - xd(1:m,j) ) ** 2 ) ) end do call phi ( nd, r, r0, v ) a(i,1:nd) = v(1:nd) end do ! ! Solve for the weights. ! call r8mat_solve_svd ( nd, nd, a, fd, w ) return end