hybsvd.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	hybsvd (NA, NU, NV, NZ, NB, M, N, A, W, MATU, U, MATV, V, Z, B, IRHS, IERR, RV1)
subroutine	mgnsvd (NU, NV, NZ, NB, M, N, W, MATU, U, MATV, V, Z,
subroutine	grsvd (NU, NV, NB, M, N, W, MATU, U, MATV, V, B, IRHS,
real *8 function	srelpr (DUMMY)

Function/Subroutine Documentation

grsvd()

subroutine grsvd	(	integer	NU,
		integer	NV,
		integer	NB,
		integer	M,
		integer	N,
		real*8, dimension(1)	W,
		logical	MATU,
		real*8, dimension(nu,1)	U,
		logical	MATV,
		real*8, dimension(nv,1)	V,
		real*8, dimension(nb,irhs)	B,
		integer	IRHS
	)

     * ierr, rv1)
c
      INTEGER i, j, k, l, m, n, ii, i1, kk, k1, ll, l1, mn, nu, nv, nb,
     * its, ierr, irhs
      REAL*8 w(1), u(nu,1), v(nv,1), b(nb,irhs), rv1(1)
      REAL*8 c, f, g, h, s, x, y, z, eps, scale, srelpr, dummy
      REAL*8 dsqrt, dmax1, dabs, dsign
      LOGICAL matu, matv
c
c     this SUBROUTINE is a translation of the algol procedure svd,
c     num. math. 14, 403-420(1970) by golub and reinsch.
c     handbook for auto. comp., vol ii-linear algebra, 134-151(1971).
c
c     this SUBROUTINE determines the singular value decomposition
c          t
c     a=usv  of a REAL m by n rectangular matrix.  householder
c     bidiagonalization and a variant of the qr algorithm are used.
c     grsvd assumes that a copy of the matrix a is in the array u. it
c     also assumes m .GE. n.  IF m .LT. n, THEN compute the singular
c                             t       t    t             t
c     VALUE decomposition of a .  IF a =uwv  , THEN a=vwu  .
c
c     grsvd can also be used to compute the minimal length least squares
c     solution to the overdetermined linear system a*x=b.
c
c     on input-
c
c        nu must be set to the row dimension of two-dimensional
c          array parameters u as declared in the calling PROGRAM
c          dimension statement.  note that nu must be at least
c          as large as m,
c
c        nv must be set to the row dimension of the two-dimensional
c          array PARAMETER v as declared in the calling PROGRAM
c          dimension statement.  nv must be at least as large as n,
c
c        nb must be set to the row dimension of two-dimensional
c          array parameters b as declared in the calling PROGRAM
c          dimension statement.  note that nb must be at least
c          as large as m,
c
c        m is the number of rows of a(and u),
c
c        n is the number of columns of a(and u) and the order of v,
c
c        a CONTAINS the rectangular input matrix to be decomposed,
c
c        b CONTAINS the irhs right-hand-sides of the overdetermined
c          linear system a*x=b.  IF irhs .GT. 0,  THEN on output,
c          the first n components of these irhs columns of b
c                        t
c          will contain u b.  thus, to compute the minimal length least
c                                                +
c          squares solution, one must compute v*w  times the columns of
c                    +                        +
c          b, WHERE w  is a diagonal matrix, w(i)=0 IF w(i) is
c          negligible, otherwise is 1/w(i).  IF irhs=0, b may coincide
c          with a or u and will not be referenced,
c
c        irhs is the number of right-hand-sides of the overdetermined
c          system a*x=b.  irhs should be set to zero IF only the singula
c          VALUE decomposition of a is desired,
c
c        matu should be set to .true. IF the u matrix in the
c          decomposition is desired, and to .false. otherwise,
c
c        matv should be set to .true. IF the v matrix in the
c          decomposition is desired, and to .false. otherwise.
c
c     on output-
c
c        w CONTAINS the n(non-negative) singular values of a(the
c          diagonal elements of s).  they are unordered.  IF an
c          error EXIT is made, the singular values should be correct
c          for indices ierr+1,ierr+2,...,n,
c
c        u CONTAINS the matrix u(orthogonal column vectors) of the
c          decomposition IF matu has been set to .true.  otherwise
c          u is used as a temporary array.
c          IF an error EXIT is made, the columns of u corresponding
c          to indices of correct singular values should be correct,
c
c        v CONTAINS the matrix v(orthogonal) of the decomposition IF
c          matv has been set to .true.  otherwise v is not referenced.
c          IF an error EXIT is made, the columns of v corresponding to
c          indices of correct singular values should be correct,
c
c        ierr is set to
c          zero       for normal RETURN,
c          k          IF the k-th singular VALUE has not been
c                     determined after 30 iterations,
c          -1         IF irhs .LT. 0 ,
c          -2         IF m .LT. n ,
c          -3         IF nu .LT. m .OR. nb .LT. m,
c          -4         IF nv .LT. n .
c
c        rv1 is a temporary storage array.
c
c        this SUBROUTINE has been checked by the pfort verifier
c(ryder, b.g. 'THE PFORT VERIFIER', software - practice and
c        experience, vol.4, 359-377, 1974) for adherence to a large,
c        carefully defined, portable subset of american national standar
c        fortran called pfort.
c
c        original version of this code is SUBROUTINE svd in release 2 of
c        eispack.
c
c        modified by tony f. chan,
c                    comp. sci. dept, yale univ.,
c                    box 2158, yale station,
c                    ct 06520
c        last modified : january, 1982.
c
c     ------------------------------------------------------------------
c
c     ********** srelpr is a machine-dependent FUNCTION specifying
c                the relative precision of floating point arithmetic.
c
c                **********
c
      ierr = 0
      IF (irhs.GE.0) GO TO 10
      ierr = -1
      RETURN
   10 IF (m.GE.n) GO TO 20
      ierr = -2
      RETURN
   20 IF (nu.GE.m .AND. nb.GE.m) GO TO 30
      ierr = -3
      RETURN
   30 IF (nv.GE.n) GO TO 40
      ierr = -4
      RETURN
   40 CONTINUE
c
c     ********** householder reduction to bidiagonal form **********
      g = 0.d0
      scale = 0.d0
      x = 0.d0
c
      DO 260 i=1,n
        l = i + 1
        rv1(i) = scale*g
        g = 0.d0
        s = 0.d0
        scale = 0.d0
c
c     compute left transformations that zero the subdiagonal elements
c     of the i-th column.
c
        DO 50 k=i,m
          scale = scale + dabs(u(k,i))
   50   CONTINUE
c
        IF (scale.EQ.0.d0) GO TO 160
c
        DO 60 k=i,m
          u(k,i) = u(k,i)/scale
          s = s + u(k,i)**2
   60   CONTINUE
c
        f = u(i,i)
        g = -dsign(dsqrt(s),f)
        h = f*g - s
        u(i,i) = f - g
        IF (i.EQ.n) GO TO 100
c
c     apply left transformations to remaining columns of a.
c
        DO 90 j=l,n
          s = 0.d0
c
          DO 70 k=i,m
            s = s + u(k,i)*u(k,j)
   70     CONTINUE
c
          f = s/h
c
          DO 80 k=i,m
            u(k,j) = u(k,j) + f*u(k,i)
   80     CONTINUE
   90   CONTINUE
c
c      apply left transformations to the columns of b IF irhs .GT. 0
c
  100   IF (irhs.EQ.0) GO TO 140
        DO 130 j=1,irhs
          s = 0.d0
          DO 110 k=i,m
            s = s + u(k,i)*b(k,j)
  110     CONTINUE
          f = s/h
          DO 120 k=i,m
            b(k,j) = b(k,j) + f*u(k,i)
  120     CONTINUE
  130   CONTINUE
c
c     compute right transformations.
c
  140   DO 150 k=i,m
          u(k,i) = scale*u(k,i)
  150   CONTINUE
c
  160   w(i) = scale*g
        g = 0.d0
        s = 0.d0
        scale = 0.d0
        IF (i.GT.m .OR. i.EQ.n) GO TO 250
c
        DO 170 k=l,n
          scale = scale + dabs(u(i,k))
  170   CONTINUE
c
        IF (scale.EQ.0.d0) GO TO 250
c
        DO 180 k=l,n
          u(i,k) = u(i,k)/scale
          s = s + u(i,k)**2
  180   CONTINUE
c
        f = u(i,l)
        g = -dsign(dsqrt(s),f)
        h = f*g - s
        u(i,l) = f - g
c
        DO 190 k=l,n
          rv1(k) = u(i,k)/h
  190   CONTINUE
c
        IF (i.EQ.m) GO TO 230
c
        DO 220 j=l,m
          s = 0.d0
c
          DO 200 k=l,n
            s = s + u(j,k)*u(i,k)
  200     CONTINUE
c
          DO 210 k=l,n
            u(j,k) = u(j,k) + s*rv1(k)
  210     CONTINUE
  220   CONTINUE
c
  230   DO 240 k=l,n
          u(i,k) = scale*u(i,k)
  240   CONTINUE
c
  250   x = dmax1(x,dabs(w(i))+dabs(rv1(i)))
  260 CONTINUE
c     ********** accumulation of right-hand transformations **********
      IF (.NOT.matv) GO TO 350
c     ********** for i=n step -1 until 1 DO -- **********
      DO 340 ii=1,n
        i = n + 1 - ii
        IF (i.EQ.n) GO TO 330
        IF (g.EQ.0.d0) GO TO 310
c
        DO 270 j=l,n
c     ********** double division avoids possible underflow **********
          v(j,i) = (u(i,j)/u(i,l))/g
  270   CONTINUE
c
        DO 300 j=l,n
          s = 0.d0
c
          DO 280 k=l,n
            s = s + u(i,k)*v(k,j)
  280     CONTINUE
c
          DO 290 k=l,n
            v(k,j) = v(k,j) + s*v(k,i)
  290     CONTINUE
  300   CONTINUE
c
  310   DO 320 j=l,n
          v(i,j) = 0.d0
          v(j,i) = 0.d0
  320   CONTINUE
c
  330   v(i,i) = 1.d0
        g = rv1(i)
        l = i
  340 CONTINUE
c     ********** accumulation of left-hand transformations **********
  350 IF (.NOT.matu) GO TO 470
c     **********for i=min(m,n) step -1 until 1 DO -- **********
      mn = n
      IF (m.LT.n) mn = m
c
      DO 460 ii=1,mn
        i = mn + 1 - ii
        l = i + 1
        g = w(i)
        IF (i.EQ.n) GO TO 370
c
        DO 360 j=l,n
          u(i,j) = 0.d0
  360   CONTINUE
c
  370   IF (g.EQ.0.d0) GO TO 430
        IF (i.EQ.mn) GO TO 410
c
        DO 400 j=l,n
          s = 0.d0
c
          DO 380 k=l,m
            s = s + u(k,i)*u(k,j)
  380     CONTINUE
c     ********** double division avoids possible underflow **********
          f = (s/u(i,i))/g
c
          DO 390 k=i,m
            u(k,j) = u(k,j) + f*u(k,i)
  390     CONTINUE
  400   CONTINUE
c
  410   DO 420 j=i,m
          u(j,i) = u(j,i)/g
  420   CONTINUE
c
        GO TO 450
c
  430   DO 440 j=i,m
          u(j,i) = 0.d0
  440   CONTINUE
c
  450   u(i,i) = u(i,i) + 1.d0
  460 CONTINUE
c     ********** diagonalization of the bidiagonal form **********
  470 eps = srelpr(dummy)*x
c     ********** for k=n step -1 until 1 DO -- **********
      DO 650 kk=1,n
        k1 = n - kk
        k = k1 + 1
        its = 0
c     ********** test for splitting.
c                for l=k step -1 until 1 DO -- **********
  480   DO 490 ll=1,k
          l1 = k - ll
          l = l1 + 1
          IF (dabs(rv1(l)).LE.eps) GO TO 550
c     ********** rv1(1) is always zero, so there is no EXIT
c                through the bottom of the loop **********
          IF (dabs(w(l1)).LE.eps) GO TO 500
  490   CONTINUE
c     ********** cancellation of rv1(l) IF l greater than 1 **********
  500   c = 0.d0
        s = 1.d0
c
        DO 540 i=l,k
          f = s*rv1(i)
          rv1(i) = c*rv1(i)
          IF (dabs(f).LE.eps) GO TO 550
          g = w(i)
          h = dsqrt(f*f+g*g)
          w(i) = h
          c = g/h
          s = -f/h
c
c     apply left transformations to b IF irhs .GT. 0
c
          IF (irhs.EQ.0) GO TO 520
          DO 510 j=1,irhs
            y = b(l1,j)
            z = b(i,j)
            b(l1,j) = y*c + z*s
            b(i,j) = -y*s + z*c
  510     CONTINUE
  520     CONTINUE
c
          IF (.NOT.matu) GO TO 540
c
          DO 530 j=1,m
            y = u(j,l1)
            z = u(j,i)
            u(j,l1) = y*c + z*s
            u(j,i) = -y*s + z*c
  530     CONTINUE
c
  540   CONTINUE
c     ********** test for convergence **********
  550   z = w(k)
        IF (l.EQ.k) GO TO 630
c     ********** shift from bottom 2 by 2 minor **********
        IF (its.EQ.30) GO TO 660
        its = its + 1
        x = w(l)
        y = w(k1)
        g = rv1(k1)
        h = rv1(k)
        f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.d0*h*y)
        g = dsqrt(f*f+1.0)
        f = ((x-z)*(x+z)+h*(y/(f+dsign(g,f))-h))/x
c     ********** next qr transformation **********
        c = 1.0
        s = 1.0
c
        DO 620 i1=l,k1
          i = i1 + 1
          g = rv1(i)
          y = w(i)
          h = s*g
          g = c*g
          z = dsqrt(f*f+h*h)
          rv1(i1) = z
          c = f/z
          s = h/z
          f = x*c + g*s
          g = -x*s + g*c
          h = y*s
          y = y*c
          IF (.NOT.matv) GO TO 570
c
          DO 560 j=1,n
            x = v(j,i1)
            z = v(j,i)
            v(j,i1) = x*c + z*s
            v(j,i) = -x*s + z*c
  560     CONTINUE
c
  570     z = dsqrt(f*f+h*h)
          w(i1) = z
c     ********** rotation can be arbitrary IF z is zero **********
          IF (z.EQ.0.d0) GO TO 580
          c = f/z
          s = h/z
  580     f = c*g + s*y
          x = -s*g + c*y
c
c     apply left transformations to b IF irhs .GT. 0
c
          IF (irhs.EQ.0) GO TO 600
          DO 590 j=1,irhs
            y = b(i1,j)
            z = b(i,j)
            b(i1,j) = y*c + z*s
            b(i,j) = -y*s + z*c
  590     CONTINUE
  600     CONTINUE
c
          IF (.NOT.matu) GO TO 620
c
          DO 610 j=1,m
            y = u(j,i1)
            z = u(j,i)
            u(j,i1) = y*c + z*s
            u(j,i) = -y*s + z*c
  610     CONTINUE
c
  620   CONTINUE
c
        rv1(l) = 0.d0
        rv1(k) = f
        w(k) = x
        GO TO 480
c     ********** convergence **********
  630   IF (z.GE.0.d0) GO TO 650
c     ********** w(k) is made non-negative **********
        w(k) = -z
        IF (.NOT.matv) GO TO 650
c
        DO 640 j=1,n
          v(j,k) = -v(j,k)
  640   CONTINUE
c
  650 CONTINUE
c
      GO TO 670
c     ********** set error -- no convergence to a
c                singular VALUE after 30 iterations **********
  660 ierr = k
  670 RETURN
c     ********** last card of grsvd **********

hybsvd()

subroutine hybsvd	(	integer	NA,
		integer	NU,
		integer	NV,
		integer	NZ,
			NB,
		integer	M,
		integer	N,
		real*8, dimension(na,1)	A,
		real*8, dimension(1)	W,
		logical	MATU,
		real*8, dimension(nu,1)	U,
		logical	MATV,
		real*8, dimension(nv,1)	V,
		real*8, dimension(nz,1)	Z,
		real*8, dimension(nb,irhs)	B,
		integer	IRHS,
		integer	IERR,
		real*8, dimension(1)	RV1
	)

      INTEGER na, nu, nv, nz, m, n, irhs, ierr, min0
      REAL*8 a(na,1), w(1), u(nu,1), v(nv,1), z(nz,1), b(nb,irhs)
      REAL*8 rv1(1)
      LOGICAL matu, matv
C
C     THIS ROUTINE IS A MODIFICATION OF THE GOLUB-REINSCH PROCEDURE (1)
C                                                           T
C     FOR COMPUTING THE SINGULAR VALUE DECOMPOSITION A = UWV  OF A
C     REAL M BY N RECTANGULAR MATRIX. U IS M BY MIN(M,N) CONTAINING
C     THE LEFT SINGULAR VECTORS, W IS A MIN(M,N) BY MIN(M,N) DIAGONAL
C     MATRIX CONTAINING THE SINGULAR VALUES, AND V IS N BY MIN(M,N)
C     CONTAINING THE RIGHT SINGULAR VECTORS.
C
C     THE ALGORITHM IMPLEMENTED IN THIS
C     ROUTINE HAS A HYBRID NATURE.  WHEN M IS APPROXIMATELY EQUAL TO N,
C     THE GOLUB-REINSCH ALGORITHM IS USED, BUT WHEN EITHER OF THE RATIOS
C     M/N OR N/M IS GREATER THAN ABOUT 2,
C     A MODIFIED VERSION OF THE GOLUB-REINSCH
C     ALGORITHM IS USED.  THIS MODIFIED ALGORITHM FIRST TRANSFORMS A
C                                                                T
C     INTO UPPER TRIANGULAR FORM BY HOUSEHOLDER TRANSFORMATIONS L
C     AND THEN USES THE GOLUB-REINSCH ALGORITHM TO FIND THE SINGULAR
C     VALUE DECOMPOSITION OF THE RESULTING UPPER TRIANGULAR MATRIX R.
C     WHEN U IS NEEDED EXPLICITLY IN THE CASE M.GE.N (OR V IN THE CASE
C     M.LT.N), AN EXTRA ARRAY Z (OF SIZE AT LEAST
C     MIN(M,N)**2) IS NEEDED, BUT OTHERWISE Z IS NOT REFERENCED
C     AND NO EXTRA STORAGE IS REQUIRED.  THIS HYBRID METHOD
C     SHOULD BE MORE EFFICIENT THAN THE GOLUB-REINSCH ALGORITHM WHEN
C     M/N OR N/M IS LARGE.  FOR DETAILS, SEE (2).
C
C     WHEN M .GE. N,
C     HYBSVD CAN ALSO BE USED TO COMPUTE THE MINIMAL LENGTH LEAST
C     SQUARES SOLUTION TO THE OVERDETERMINED LINEAR SYSTEM A*X=B.
C     IF M .LT. N (I.E. FOR UNDERDETERMINED SYSTEMS), THE RHS B
C     IS NOT PROCESSED.
C
C     NOTICE THAT THE SINGULAR VALUE DECOMPOSITION OF A MATRIX
C     IS UNIQUE ONLY UP TO THE SIGN OF THE CORRESPONDING COLUMNS
C     OF U AND V.
C
C     THIS ROUTINE HAS BEEN CHECKED BY THE PFORT VERIFIER (3) FOR
C     ADHERENCE TO A LARGE, CAREFULLY DEFINED, PORTABLE SUBSET OF
C     AMERICAN NATIONAL STANDARD FORTRAN CALLED PFORT.
C
C     REFERENCES:
C
C     (1) GOLUB,G.H. AND REINSCH,C. (1970) 'SINGULAR VALUE
C         DECOMPOSITION AND LEAST SQUARES SOLUTIONS,'
C         NUMER. MATH. 14,403-420, 1970.
C
C     (2) CHAN,T.F. (1982) 'AN IMPROVED ALGORITHM FOR COMPUTING
C         THE SINGULAR VALUE DECOMPOSITION,' ACM TOMS, VOL.8,
C         NO. 1, MARCH, 1982.
C
C     (3) RYDER,B.G. (1974) 'THE PFORT VERIFIER,' SOFTWARE -
C         PRACTICE AND EXPERIENCE, VOL.4, 359-377, 1974.
C
C     ON INPUT:
C
C        NA MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER A AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.  NOTE THAT NA MUST BE AT LEAST
C          AS LARGE AS M.
C
C        NU MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY U AS DECLARED IN THE CALLING PROGRAM DIMENSION
C          STATEMENT. NU MUST BE AT LEAST AS LARGE AS M.
C
C        NV MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER V AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT. NV MUST BE AT LEAST AS LARGE AS N.
C
C        NZ MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER Z AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.  NOTE THAT NZ MUST BE AT LEAST
C          AS LARGE AS MIN(M,N).
C
C        NB MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER B AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT. NB MUST BE AT LEAST AS LARGE AS M.
C
C        M IS THE NUMBER OF ROWS OF A (AND U).
C
C        N IS THE NUMBER OF COLUMNS OF A (AND NUMBER OF ROWS OF V).
C
C        A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED.
C
C        B CONTAINS THE IRHS RIGHT-HAND-SIDES OF THE OVERDETERMINED
C         LINEAR SYSTEM A*X=B. IF IRHS .GT. 0 AND M .GE. N,
C         THEN ON OUTPUT, THE FIRST N COMPONENTS OF THESE IRHS COLUMNS
C                       T
C         WILL CONTAIN U B. THUS, TO COMPUTE THE MINIMAL LENGTH LEAST
C                                               +
C         SQUARES SOLUTION, ONE MUST COMPUTE V*W  TIMES THE COLUMNS OF
C                   +                        +
C         B, WHERE W  IS A DIAGONAL MATRIX, W (I)=0 IF W(I) IS
C         NEGLIGIBLE, OTHERWISE IS 1/W(I). IF IRHS=0 OR M.LT.N,
C         B IS NOT REFERENCED.
C
C        IRHS IS THE NUMBER OF RIGHT-HAND-SIDES OF THE OVERDETERMINED
C         SYSTEM A*X=B. IRHS SHOULD BE SET TO ZERO IF ONLY THE SINGULAR
C         VALUE DECOMPOSITION OF A IS DESIRED.
C
C        MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE
C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
C
C        MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE
C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
C
C        WHEN HYBSVD IS USED TO COMPUTE THE MINIMAL LENGTH LEAST
C        SQUARES SOLUTION TO AN OVERDETERMINED SYSTEM, MATU SHOULD
C        BE SET TO .FALSE. , AND MATV SHOULD BE SET TO .TRUE.  .
C
C     ON OUTPUT:
C
C        A IS UNALTERED (UNLESS OVERWRITTEN BY U OR V).
C
C        W CONTAINS THE (NON-NEGATIVE) SINGULAR VALUES OF A (THE
C          DIAGONAL ELEMENTS OF W).  THEY ARE SORTED IN DESCENDING
C          ORDER.  IF AN ERROR EXIT IS MADE, THE SINGULAR VALUES
C          SHOULD BE CORRECT AND SORTED FOR INDICES IERR+1,...,MIN(M,N).
C
C        U CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE
C          DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE.  IF MATU IS
C          FALSE, THEN U IS EITHER USED AS A TEMPORARY STORAGE (IF
C          M .GE. N) OR NOT REFERENCED (IF M .LT. N).
C          U MAY COINCIDE WITH A IN THE CALLING SEQUENCE.
C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING
C          TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT.
C
C        V CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF
C          MATV HAS BEEN SET TO .TRUE.  IF MATV IS
C          FALSE, THEN V IS EITHER USED AS A TEMPORARY STORAGE (IF
C          M .LT. N) OR NOT REFERENCED (IF M .GE. N).
C          IF M .GE. N, V MAY ALSO COINCIDE WITH A.  IF AN ERROR
C          EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO INDICES OF
C          CORRECT SINGULAR VALUES SHOULD BE CORRECT.
C
C        Z CONTAINS THE MATRIX X IN THE SINGULAR VALUE DECOMPOSITION
C                  T
C          OF R=XSY,  IF THE MODIFIED ALGORITHM IS USED. IF THE
C          GOLUB-REINSCH PROCEDURE IS USED, THEN IT IS NOT REFERENCED.
C          IF MATU HAS BEEN SET TO .FALSE. IN THE CASE M.GE.N (OR
C          MATV SET TO .FALSE. IN THE CASE M.LT.N), THEN Z IS NOT
C          REFERENCED AND NO EXTRA STORAGE IS REQUIRED.
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C          -1         IF IRHS .LT. 0 .
C          -2         IF M .LT. 1 .OR. N .LT. 1
C          -3         IF NA .LT. M .OR. NU .LT. M .OR. NB .LT. M.
C          -4         IF NV .LT. N .
C          -5         IF NZ .LT. MIN(M,N).
C
C        RV1 IS A TEMPORARY STORAGE ARRAY OF LENGTH AT LEAST MIN(M,N).
C
C     PROGRAMMED BY : TONY CHAN
C                     BOX 2158, YALE STATION,
C                     COMPUTER SCIENCE DEPT, YALE UNIV.,
C                     NEW HAVEN, CT 06520.
C     LAST MODIFIED : JANUARY, 1982.
C
C     HYBSVD USES THE FOLLOWING FUNCTIONS AND SUBROUTINES.
C       INTERNAL  GRSVD, MGNSVD, SRELPR
C       FORTRAN   MIN0,DABS,DSQRT,DFLOAT,DSIGN,DMAX1
C       BLAS      DSWAP
C
C     -----------------------------------------------------------------
C     ERROR CHECK.
C
      ierr = 0
      IF (irhs.GE.0) GO TO 10
      ierr = -1
      RETURN
   10 IF (m.GE.1 .AND. n.GE.1) GO TO 20
      ierr = -2
      RETURN
   20 IF (na.GE.m .AND. nu.GE.m .AND. nb.GE.m) GO TO 30
      ierr = -3
      RETURN
   30 IF (nv.GE.n) GO TO 40
      ierr = -4
      RETURN
   40 IF (nz.GE.min0(m,n)) GO TO 50
      ierr = -5
      RETURN
   50 CONTINUE
C
C     FIRST COPIES A INTO EITHER U OR V ACCORDING TO WHETHER
C     M .GE. N OR M .LT. N, AND THEN CALLS SUBROUTINE MGNSVD
C     WHICH ASSUMES THAT NUMBER OF ROWS .GE. NUMBER OF COLUMNS.
C
      IF (m.LT.n) GO TO 80
C
C       M .GE. N  CASE.
C
      DO 70 i=1,m
        DO 60 j=1,n
          u(i,j) = a(i,j)
   60   CONTINUE
   70 CONTINUE
C
      CALL mgnsvd(nu, nv, nz, nb, m, n, w, matu, u, matv, v, z, b,
     * irhs, ierr, rv1)
      RETURN
C
   80 CONTINUE
C                              T
C       M .LT. N CASE. COPIES A  INTO V.
C
      DO 100 i=1,m
        DO 90 j=1,n
          v(j,i) = a(i,j)
   90   CONTINUE
  100 CONTINUE
      CALL mgnsvd(nv, nu, nz, nb, n, m, w, matv, v, matu, u, z, b, 0,
     * ierr, rv1)
      RETURN

mgnsvd()

subroutine mgnsvd	(	integer	NU,
		integer	NV,
		integer	NZ,
			NB,
		integer	M,
		integer	N,
		real*8, dimension(1)	W,
		logical	MATU,
		real*8, dimension(nu,1)	U,
		logical	MATV,
		real*8, dimension(nv,1)	V,
		real*8, dimension(nz,1)	Z
	)

     * b, irhs, ierr, rv1)
c
c     the description of SUBROUTINE mgnsvd is almost identical
c     to that for SUBROUTINE hybsvd above, with the exception
c     that mgnsvd assumes m .GE. n.
c     it also assumes that a copy of the matrix a is in the array u.
c
      INTEGER nu, nv, nz, m, n, irhs, ierr, ip1, i, j, k, im1, iback
      REAL*8 w(1), u(nu,1), v(nv,1), z(nz,1), b(nb,irhs), rv1(1)
      REAL*8 xovrpt, c, r, g, scale, dsign, dabs, dsqrt, f, s, h
      REAL*8 dfloat
      LOGICAL matu, matv
c
c     set VALUE for c. the VALUE for c depends on the relative
c     efficiency of floating point multiplications, floating point
c     additions and two-dimensional array indexings on the
c     computer WHERE this SUBROUTINE is to be run.  c should
c     usually be between 2 and 4.  for details on choosing c, see
c(2).  the algorithm is not sensitive to the VALUE of c
c     actually used as long as c is between 2 and 4.
c
      c = 4.d0
c
c     determine cross-over point
c
      IF (matu .AND. matv) xovrpt = (c+7.d0/3.d0)/c
      IF (matu .AND. .NOT.matv) xovrpt = (c+7.d0/3.d0)/c
      IF (.NOT.matu .AND. matv) xovrpt = 5.d0/3.d0
      IF (.NOT.matu .AND. .NOT.matv) xovrpt = 5.d0/3.d0
c
c     determine whether to USE golub-reinsch or the modified
c     algorithm.
c
      r = dfloat(m)/dfloat(n)
      IF (r.GE.xovrpt) GO TO 10
c
c     USE golub-reinsch procedure
c
      CALL grsvd(nu, nv, nb, m, n, w, matu, u, matv, v, b, irhs, ierr,
     * rv1)
      GO TO 330
c
c     USE modified algorithm
c
   10 CONTINUE
c
c     triangularize u by householder transformations, using
c     w and rv1 as temporary storage.
c
      DO 110 i=1,n
        g = 0.d0
        s = 0.d0
        scale = 0.d0
c
c         perform scaling of columns to avoid unnecssary overflow
c         or underflow
c
        DO 20 k=i,m
          scale = scale + dabs(u(k,i))
   20   CONTINUE
        IF (scale.EQ.0.d0) GO TO 110
        DO 30 k=i,m
          u(k,i) = u(k,i)/scale
          s = s + u(k,i)*u(k,i)
   30   CONTINUE
c
c         the vector e of the householder transformation i + ee'/H
C         WILL BE STORED IN COLUMN I OF U. THE TRANSFORMED ELEMENT
C         U(I,I) WILL BE STORED IN W(I) AND THE SCALAR H IN
C         RV1(I).
C
        F = U(I,I)
        G = -DSIGN(DSQRT(S),F)
        H = F*G - S
        U(I,I) = F - G
        RV1(I) = H
        W(I) = SCALE*G
C
.EQ.        IF (IN) GO TO 70
C
C         APPLY TRANSFORMATIONS TO REMAINING COLUMNS OF A
C
        IP1 = I + 1
        DO 60 J=IP1,N
          S = 0.d0
          DO 40 K=I,M
            S = S + U(K,I)*U(K,J)
   40     CONTINUE
          F = S/H
          DO 50 K=I,M
            U(K,J) = U(K,J) + F*U(K,I)
   50     CONTINUE
   60   CONTINUE
C
.GT.C         APPLY TRANSFORMATIONS TO COLUMNS OF B IF IRHS  0
C
.EQ.   70   IF (IRHS0) GO TO 110
        DO 100 J=1,IRHS
          S = 0.d0
          DO 80 K=I,M
            S = S + U(K,I)*B(K,J)
   80     CONTINUE
          F = S/H
          DO 90 K=I,M
            B(K,J) = B(K,J) + F*U(K,I)
   90     CONTINUE
  100   CONTINUE
  110 CONTINUE
C
C     COPY R INTO Z IF MATU = .TRUE.
C
.NOT.      IF (MATU) GO TO 290
      DO 130 I=1,N
        DO 120 J=I,N
          Z(J,I) = 0.d0
          Z(I,J) = U(I,J)
  120   CONTINUE
        Z(I,I) = W(I)
  130 CONTINUE
C
C     ACCUMULATE HOUSEHOLDER TRANSFORMATIONS IN U
C
      DO 240 IBACK=1,N
        I = N - IBACK + 1
        IP1 = I + 1
        G = W(I)
        H = RV1(I)
.EQ.        IF (IN) GO TO 150
C
        DO 140 J=IP1,N
          U(I,J) = 0.d0
  140   CONTINUE
C
.EQ.  150   IF (H0.d0) GO TO 210
.EQ.        IF (IN) GO TO 190
C
        DO 180 J=IP1,N
          S = 0.d0
          DO 160 K=IP1,M
            S = S + U(K,I)*U(K,J)
  160     CONTINUE
          F = S/H
          DO 170 K=I,M
            U(K,J) = U(K,J) + F*U(K,I)
  170     CONTINUE
  180   CONTINUE
C
  190   S = U(I,I)/H
        DO 200 J=I,M
          U(J,I) = U(J,I)*S
  200   CONTINUE
        GO TO 230
C
  210   DO 220 J=I,M
          U(J,I) = 0.d0
  220   CONTINUE
  230   U(I,I) = U(I,I) + 1.d0
  240 CONTINUE
C
C     COMPUTE SVD OF R (WHICH IS STORED IN Z)
C
      CALL GRSVD(NZ, NV, NB, N, N, W, MATU, Z, MATV, V, B, IRHS, IERR,
     * RV1)
C
C                                      T
C     FORM L*X TO OBTAIN U (WHERE R=XWY ). X IS RETURNED IN Z
C     BY GRSVD. THE MATRIX MULTIPLY IS DONE ONE ROW AT A TIME,
C     USING RV1 AS SCRATCH SPACE.
C
      DO 280 I=1,M
        DO 260 J=1,N
          S = 0.d0
          DO 250 K=1,N
            S = S + U(I,K)*Z(K,J)
  250     CONTINUE
          RV1(J) = S
  260   CONTINUE
        DO 270 J=1,N
          U(I,J) = RV1(J)
  270   CONTINUE
  280 CONTINUE
      GO TO 330
C
C     FORM R IN U BY ZEROING THE LOWER TRIANGULAR PART OF R IN U
C
.EQ.  290 IF (N1) GO TO 320
      DO 310 I=2,N
        IM1 = I - 1
        DO 300 J=1,IM1
          U(I,J) = 0.d0
  300   CONTINUE
        U(I,I) = W(I)
  310 CONTINUE
  320 U(1,1) = W(1)
C
      CALL GRSVD(NU, NV, NB, N, N, W, MATU, U, MATV, V, B, IRHS, IERR,
     * RV1)
  330 CONTINUE
      IERRP1 = IERR + 1
.LT..OR..LE..OR..EQ.      IF (IERR0  N1  IERRP1N) RETURN
C
C     SORT SINGULAR VALUES AND EXCHANGE COLUMNS OF U AND V ACCORDINGLY.
C     SELECTION SORT MINIMIZES SWAPPING OF U AND V.
C
      NM1 = N - 1
      DO 360 I=IERRP1,NM1
C...    FIND INDEX OF MAXIMUM SINGULAR VALUE
        ID = I
        IP1 = I + 1
        DO 340 J=IP1,N
.GT.          IF (W(J)W(ID)) ID = J
  340   CONTINUE
.EQ.        IF (IDI) GO TO 360
C...    SWAP SINGULAR VALUES AND VECTORS
        T = W(I)
        W(I) = W(ID)
        W(ID) = T
        IF (MATV) CALL DSWAP(N, V(1,I), 1, V(1,ID), 1)
        IF (MATU) CALL DSWAP(M, U(1,I), 1, U(1,ID), 1)
.LT.        IF (IRHS1) GO TO 360
        DO 350 KRHS=1,IRHS
          T = B(I,KRHS)
          B(I,KRHS) = B(ID,KRHS)
          B(ID,KRHS) = T
  350   CONTINUE
  360 CONTINUE
      RETURN
C     ************** LAST CARD OF HYBSVD *****************

srelpr()

real*8 function srelpr

(

real*8

DUMMY

)

      REAL*8 dummy
C
C     SRELPR COMPUTES THE RELATIVE PRECISION OF THE FLOATING POINT
C     ARITHMETIC OF THE MACHINE.
C
C     IF TROUBLE WITH AUTOMATIC COMPUTATION OF THESE QUANTITIES,
C     THEY CAN BE SET BY DIRECT ASSIGNMENT STATEMENTS.
C     ASSUME THE COMPUTER HAS
C
C        B = BASE OF ARITHMETIC
C        T = NUMBER OF BASE  B  DIGITS
C
C     THEN
C
C        SRELPR = B**(1-T)
C
      REAL*8 s
C
      srelpr = 1.d0
   10 srelpr = srelpr/2.d0
      s = 1.d0 + srelpr
      IF (s.GT.1.d0) GO TO 10
      srelpr = 2.d0*srelpr
      RETURN