!     DGESVD Example Program Text
!     NAG Copyright 2005.
!     .. Parameters ..
  Integer nin, nout
  Parameter (nin=5, nout=6)
  Integer mmax, nb, nmax
  Parameter (mmax=10, nb=64, nmax=8)
  Integer lda, ldvt, lwork
  Parameter (lda=mmax, ldvt=nmax, lwork=mmax+4*nmax+nb*(mmax+nmax))
!     .. Local Scalars ..
  Double Precision eps, serrbd
  Integer i, ifail, info, j, lwkopt, m, n
!     .. Local Arrays ..
  Double Precision a(lda, nmax), dummy(1, 1), rcondu(nmax), rcondv(nmax), &
    s(nmax), uerrbd(nmax), verrbd(nmax), vt(ldvt, nmax), work(lwork)
!     .. External Functions ..
  Double Precision dlamch
  External dlamch
!     .. External Subroutines ..
  External ddisna, dgesvd, x04caf
!     .. Executable Statements ..
  Write (nout, *) 'DGESVD Example Program Results'
  Write (nout, *)
!     Skip heading in data file
  Read (nin, *)
  Read (nin, *) m, n
  If (m<=mmax .And. n<=nmax) Then
!
!        Read the m by n matrix A from data file
!
    Read (nin, *)((a(i,j),j=1,n), i=1, m)
!
!        Compute the singular values and left and right singular vectors
!        of A (A = U*S*(V**T), m.ge.n)
!
    Call dgesvd('Overwrite A by U', 'Singular vectors (V)', m, n, a, lda, s, &
      dummy, 1, vt, ldvt, work, lwork, info)
    lwkopt = work(1)
!
    If (info==0) Then
!
!           Print solution
!
      Write (nout, *) 'Singular values'
      Write (nout, 100)(s(j), j=1, n)
!
      ifail = 0
!      Call x04caf('General', ' ', m, n, a, lda, &
!        'Left singular vectors (first n columns of U)', ifail)
      print *,a
      Write (nout, *)
!      Call x04caf('General', ' ', n, n, vt, ldvt, &
!        'Right singular vectors by row (V**T)', ifail)
      print *,vt
!
!           Get the machine precision, EPS and compute the approximate
!           error bound for the computed singular values.  Note that for
!           the 2-norm, S(1) = norm(A)
!
      eps = dlamch('Eps')
      serrbd = eps*s(1)
!
!           Call DDISNA (F08FLF) to estimate reciprocal condition
!           numbers for the singular vectors
!
      Call ddisna('Left', m, n, s, rcondu, info)
      Call ddisna('Right', m, n, s, rcondv, info)
!
!           Compute the error estimates for the singular vectors
!
      Do i = 1, n
        uerrbd(i) = serrbd/rcondu(i)
        verrbd(i) = serrbd/rcondv(i)
      End Do
!
!           Print the approximate error bounds for the singular values
!           and vectors
!
      Write (nout, *)
      Write (nout, *) 'Error estimate for the singular values'
      Write (nout, 110) serrbd
      Write (nout, *)
      Write (nout, *) 'Error estimates for the left singular vectors'
      Write (nout, 110)(uerrbd(i), i=1, n)
      Write (nout, *)
      Write (nout, *) 'Error estimates for the right singular vectors'
      Write (nout, 110)(verrbd(i), i=1, n)
    Else
      Write (nout, 120) 'Failure in DGESVD. INFO =', info
    End If
!
!        Print workspace information
!
    If (lwork<lwkopt) Then
      Write (nout, *)
      Write (nout, 130) 'Optimum workspace required = ', lwkopt, &
        'Workspace provided         = ', lwork
    End If
  Else
    Write (nout, *) 'MMAX and/or NMAX too small'
  End If
  Stop
!
100 Format (3X, (8F8.4))
110 Format (4X, 1P, 6E11.1)
120 Format (1X, A, I4)
130 Format (1X, A, I5, /1X, A, I5)
End Program
