text:
      SUBROUTINE SG03BD( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q,
     $                   LDQ, Z, LDZ, B, LDB, SCALE, ALPHAR, ALPHAI,
     $                   BETA, DWORK, LDWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute the Cholesky factor U of the matrix X,
C
C                 T
C        X = op(U)  * op(U),
C
C     which is the solution of either the generalized
C     c-stable continuous-time Lyapunov equation
C
C             T                    T
C        op(A)  * X * op(E) + op(E)  * X * op(A)
C
C                 2        T
C        = - SCALE  * op(B)  * op(B),                                (1)
C
C     or the generalized d-stable discrete-time Lyapunov equation
C
C             T                    T
C        op(A)  * X * op(A) - op(E)  * X * op(E)
C
C                 2        T
C        = - SCALE  * op(B)  * op(B),                                (2)
C
C     without first finding X and without the need to form the matrix
C     op(B)**T * op(B).
C
C     op(K) is either K or K**T for K = A, B, E, U. A and E are N-by-N
C     matrices, op(B) is an M-by-N matrix. The resulting matrix U is an
C     N-by-N upper triangular matrix with non-negative entries on its
C     main diagonal. SCALE is an output scale factor set to avoid
C     overflow in U.
C
C     In the continuous-time case (1) the pencil A - lambda * E must be
C     c-stable (that is, all eigenvalues must have negative real parts).
C     In the discrete-time case (2) the pencil A - lambda * E must be
C     d-stable (that is, the moduli of all eigenvalues must be smaller
C     than one).
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     DICO    CHARACTER*1
C             Specifies which type of the equation is considered:
C             = 'C':  Continuous-time equation (1);
C             = 'D':  Discrete-time equation (2).
C
C     FACT    CHARACTER*1
C             Specifies whether the generalized real Schur
C             factorization of the pencil A - lambda * E is supplied
C             on entry or not:
C             = 'N':  Factorization is not supplied;
C             = 'F':  Factorization is supplied.
C
C     TRANS   CHARACTER*1
C             Specifies whether the transposed equation is to be solved
C             or not:
C             = 'N':  op(A) = A,    op(E) = E;
C             = 'T':  op(A) = A**T, op(E) = E**T.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the matrix A.  N >= 0.
C
C     M       (input) INTEGER
C             The number of rows in the matrix op(B).  M >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, if FACT = 'F', then the leading N-by-N upper
C             Hessenberg part of this array must contain the
C             generalized Schur factor A_s of the matrix A (see
C             definition (3) in section METHOD). A_s must be an upper
C             quasitriangular matrix. The elements below the upper
C             Hessenberg part of the array A are not referenced.
C             If FACT = 'N', then the leading N-by-N part of this
C             array must contain the matrix A.
C             On exit, the leading N-by-N part of this array contains
C             the generalized Schur factor A_s of the matrix A. (A_s is
C             an upper quasitriangular matrix.)
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= MAX(1,N).
C
C     E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
C             On entry, if FACT = 'F', then the leading N-by-N upper
C             triangular part of this array must contain the
C             generalized Schur factor E_s of the matrix E (see
C             definition (4) in section METHOD). The elements below the
C             upper triangular part of the array E are not referenced.
C             If FACT = 'N', then the leading N-by-N part of this
C             array must contain the coefficient matrix E of the
C             equation.
C             On exit, the leading N-by-N part of this array contains
C             the generalized Schur factor E_s of the matrix E. (E_s is
C             an upper triangular matrix.)
C
C     LDE     INTEGER
C             The leading dimension of the array E.  LDE >= MAX(1,N).
C
C     Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
C             On entry, if FACT = 'F', then the leading N-by-N part of
C             this array must contain the orthogonal matrix Q from
C             the generalized Schur factorization (see definitions (3)
C             and (4) in section METHOD).
C             If FACT = 'N', Q need not be set on entry.
C             On exit, the leading N-by-N part of this array contains
C             the orthogonal matrix Q from the generalized Schur
C             factorization.
C
C     LDQ     INTEGER
C             The leading dimension of the array Q.  LDQ >= MAX(1,N).
C
C     Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
C             On entry, if FACT = 'F', then the leading N-by-N part of
C             this array must contain the orthogonal matrix Z from
C             the generalized Schur factorization (see definitions (3)
C             and (4) in section METHOD).
C             If FACT = 'N', Z need not be set on entry.
C             On exit, the leading N-by-N part of this array contains
C             the orthogonal matrix Z from the generalized Schur
C             factorization.
C
C     LDZ     INTEGER
C             The leading dimension of the array Z.  LDZ >= MAX(1,N).
C
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,N1)
C             On entry, if TRANS = 'T', the leading N-by-M part of this
C             array must contain the matrix B and N1 >= MAX(M,N).
C             If TRANS = 'N', the leading M-by-N part of this array
C             must contain the matrix B and N1 >= N.
C             On exit, the leading N-by-N part of this array contains
C             the Cholesky factor U of the solution matrix X of the
C             problem, X = op(U)**T * op(U).
C             If M = 0 and N > 0, then U is set to zero.
C
C     LDB     INTEGER
C             The leading dimension of the array B.
C             If TRANS = 'T', LDB >= MAX(1,N).
C             If TRANS = 'N', LDB >= MAX(1,M,N).
C
C     SCALE   (output) DOUBLE PRECISION
C             The scale factor set to avoid overflow in U.
C             0 < SCALE <= 1.
C
C     ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
C     ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
C     BETA    (output) DOUBLE PRECISION array, dimension (N)
C             If INFO = 0, 3, 5, 6, or 7, then
C             (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, are the
C             eigenvalues of the matrix pencil A - lambda * E.
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK.
C
C     LDWORK  INTEGER
C             The dimension of the array DWORK.
C             LDWORK >= MAX(1,4*N,6*N-6),  if FACT = 'N';
C             LDWORK >= MAX(1,2*N,6*N-6),  if FACT = 'F'.
C             For good performance, LDWORK should be larger.
C
C             If LDWORK = -1, then a workspace query is assumed; the
C             routine only calculates the optimal size of the DWORK
C             array, returns this value as the first entry of the DWORK
C             array, and no error message related to LDWORK is issued by
C             XERBLA.
C
C     Error indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             = 1:  the pencil A - lambda * E is (nearly) singular;
C                   perturbed values were used to solve the equation
C                   (but the reduced (quasi)triangular matrices A and E
C                   are unchanged);
C             = 2:  FACT = 'F' and the matrix contained in the upper
C                   Hessenberg part of the array A is not in upper
C                   quasitriangular form;
C             = 3:  FACT = 'F' and there is a 2-by-2 block on the main
C                   diagonal of the pencil A_s - lambda * E_s whose
C                   eigenvalues are not conjugate complex;
C             = 4:  FACT = 'N' and the pencil A - lambda * E cannot be
C                   reduced to generalized Schur form: LAPACK routine
C                   DGEGS (or DGGES) has failed to converge;
C             = 5:  DICO = 'C' and the pencil A - lambda * E is not
C                   c-stable;
C             = 6:  DICO = 'D' and the pencil A - lambda * E is not
C                   d-stable;
C             = 7:  the LAPACK routine DSYEVX utilized to factorize M3
C                   failed to converge in the discrete-time case (see
C                   section METHOD for SLICOT Library routine SG03BU).
C                   This error is unlikely to occur.
C
C     METHOD
C
C     An extension [2] of Hammarling's method [1] to generalized
C     Lyapunov equations is utilized to solve (1) or (2).
C
C     First the pencil A - lambda * E is reduced to real generalized
C     Schur form A_s - lambda * E_s by means of orthogonal
C     transformations (QZ-algorithm):
C
C        A_s = Q**T * A * Z   (upper quasitriangular)                (3)
C
C        E_s = Q**T * E * Z   (upper triangular).                    (4)
C
C     If the pencil A - lambda * E has already been factorized prior to
C     calling the routine however, then the factors A_s, E_s, Q and Z
C     may be supplied and the initial factorization omitted.
C
C     Depending on the parameters TRANS and M the N-by-N upper
C     triangular matrix B_s is defined as follows. In any case Q_B is
C     an M-by-M orthogonal matrix, which need not be accumulated.
C
C     1. If TRANS = 'N' and M < N, B_s is the upper triangular matrix
C        from the QR-factorization
C
C           ( Q_B  O )           ( B * Z )
C           (        ) * B_s  =  (       ),
C           (  O   I )           (   O   )
C
C        where the O's are zero matrices of proper size and I is the
C        identity matrix of order N-M.
C
C     2. If TRANS = 'N' and M >= N, B_s is the upper triangular matrix
C        from the (rectangular) QR-factorization
C
C                 ( B_s )
C           Q_B * (     )  =  B * Z,
C                 (  O  )
C
C        where O is the (M-N)-by-N zero matrix.
C
C     3. If TRANS = 'T' and M < N, B_s is the upper triangular matrix
C        from the RQ-factorization
C
C                       ( Q_B  O )
C           (B_s  O ) * (        )  =  ( Q**T * B   O ).
C                       (  O   I )
C
C     4. If TRANS = 'T' and M >= N, B_s is the upper triangular matrix
C        from the (rectangular) RQ-factorization
C
C           ( B_s   O ) * Q_B  =  Q**T * B,
C
C        where O is the N-by-(M-N) zero matrix.
C
C     Assuming SCALE = 1, the transformation of A, E and B described
C     above leads to the reduced continuous-time equation
C
C                 T        T
C          op(A_s)  op(U_s)  op(U_s) op(E_s)
C
C                 T        T
C        + op(E_s)  op(U_s)  op(U_s) op(A_s)
C
C                    T
C        =  - op(B_s)  op(B_s)                                       (5)
C
C     or to the reduced discrete-time equation
C
C                 T        T
C          op(A_s)  op(U_s)  op(U_s) op(A_s)
C
C                 T        T
C        - op(E_s)  op(U_s)  op(U_s) op(E_s)
C
C                    T
C        =  - op(B_s)  op(B_s).                                      (6)
C
C     For brevity we restrict ourself to equation (5) and the case
C     TRANS = 'N'. The other three cases can be treated in a similar
C     fashion.
C
C     We use the following partitioning for the matrices A_s, E_s, B_s
C     and U_s
C
C                 ( A11   A12 )          ( E11   E12 )
C           A_s = (           ),   E_s = (           ),
C                 (   0   A22 )          (   0   E22 )
C
C                 ( B11   B12 )          ( U11   U12 )
C           B_s = (           ),   U_s = (           ).              (7)
C                 (   0   B22 )          (   0   U22 )
C
C     The size of the (1,1)-blocks is 1-by-1 (iff A_s(2,1) = 0.0) or
C     2-by-2.
C
C     We compute U11 and U12**T in three steps.
C
C     Step I:
C
C        From (5) and (7) we get the 1-by-1 or 2-by-2 equation
C
C                T      T                   T      T
C             A11  * U11  * U11 * E11  + E11  * U11  * U11 * A11
C
C                    T
C             = - B11  * B11.
C
C        For brevity, details are omitted here. See [2]. The technique
C        for computing U11 is similar to those applied to standard
C        Lyapunov equations in Hammarling's algorithm ([1], section 6).
C
C        Furthermore, the auxiliary matrices M1 and M2 defined as
C        follows
C
C                               -1      -1
C           M1 = U11 * A11 * E11   * U11
C
C                         -1      -1
C           M2 = B11 * E11   * U11
C
C        are computed in a numerically reliable way.
C
C     Step II:
C
C        The generalized Sylvester equation
C
C              T      T      T      T
C           A22  * U12  + E22  * U12  * M1  =
C
C                T           T      T      T      T
C           - B12  * M2 - A12  * U11  - E12  * U11  * M1
C
C        is solved for U12**T.
C
C     Step III:
C
C        It can be shown that
C
C              T      T                  T      T
C           A22  * U22  * U22 * E22 + E22  * U22  * U22 * A22  =
C
C                T              T
C           - B22  * B22 - y * y                                     (8)
C
C        holds, where y is defined as
C
C                  T        T      T      T      T       T
C           y = B12  - ( E12  * U11  + E22  * U12  ) * M2 .
C
C        If B22_tilde is the square triangular matrix arising from the
C        (rectangular) QR-factorization
C
C                       ( B22_tilde )     ( B22  )
C           Q_B_tilde * (           )  =  (      ),
C                       (     O     )     ( y**T )
C
C        where Q_B_tilde is an orthogonal matrix of order N, then
C
C                T              T                T
C           - B22  * B22 - y * y   =  - B22_tilde  * B22_tilde.
C
C        Replacing the right hand side in (8) by the term
C        - B22_tilde**T * B22_tilde leads to a reduced generalized
C        Lyapunov equation of lower dimension compared to (5).
C
C     The recursive application of the steps I to III yields the
C     solution U_s of the equation (5).
C
C     It remains to compute the solution matrix U of the original
C     problem (1) or (2) from the matrix U_s. To this end we transform
C     the solution back (with respect to the transformation that led
C     from (1) to (5) (from (2) to (6)) and apply the QR-factorization
C     (RQ-factorization). The upper triangular solution matrix U is
C     obtained by
C
C        Q_U * U  =  U_s * Q**T     (if TRANS = 'N')
C
C     or
C
C        U * Q_U  =  Z * U_s        (if TRANS = 'T')
C
C     where Q_U is an N-by-N orthogonal matrix. Again, the orthogonal
C     matrix Q_U need not be accumulated.
C
C     REFERENCES
C
C     [1] Hammarling, S.J.
C         Numerical solution of the stable, non-negative definite
C         Lyapunov equation.
C         IMA J. Num. Anal., 2, pp. 303-323, 1982.
C
C     [2] Penzl, T.
C         Numerical solution of generalized Lyapunov equations.
C         Advances in Comp. Math., vol. 8, pp. 33-48, 1998.
C
C     NUMERICAL ASPECTS
C
C     The number of flops required by the routine is given by the
C     following table. Note that we count a single floating point
C     arithmetic operation as one flop.
C
C                 |           FACT = 'F'                  FACT = 'N'
C        ---------+--------------------------------------------------
C         M <= N  |     (13*N**3+6*M*N**2         (211*N**3+6*M*N**2
C                 |   +6*M**2*N-2*M**3)/3        +6*M**2*N-2*M**3)/3
C                 |
C          M > N  | (11*N**3+12*M*N**2)/3     (209*N**3+12*M*N**2)/3
C
C     FURTHER COMMENTS
C
C     The Lyapunov equation may be very ill-conditioned. In particular,
C     if DICO = 'D' and the pencil A - lambda * E has a pair of almost
C     reciprocal eigenvalues, or DICO = 'C' and the pencil has an almost
C     degenerate pair of eigenvalues, then the Lyapunov equation will be
C     ill-conditioned. Perturbed values were used to solve the equation.
C     A condition estimate can be obtained from the routine SG03AD.
C     When setting the error indicator INFO, the routine does not test
C     for near instability in the equation but only for exact
C     instability.
C
C     CONTRIBUTOR
C
C     T. Penzl, Technical University Chemnitz, Germany, Aug. 1998.
C
C     REVISIONS
C
C     Sep. 1998 (V. Sima).
C     May 1999 (V. Sima).
C     March 2002 (A. Varga).
C     Feb. 2004 (V. Sima).
C     July 2011 (V. Sima).
C
C     KEYWORDS
C
C     Lyapunov equation
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  MONE, ONE, ZERO
      PARAMETER         ( MONE = -1.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
C     .. Scalar Arguments ..
      DOUBLE PRECISION  SCALE
      INTEGER           INFO, LDA, LDB, LDE, LDQ, LDWORK, LDZ, M, N
      CHARACTER         DICO, FACT, TRANS
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*),
     $                  BETA(*), DWORK(*), E(LDE,*), Q(LDQ,*), Z(LDZ,*)
C     .. Local Scalars ..
      DOUBLE PRECISION  S1, S2, SAFMIN, WI, WR1, WR2
      INTEGER           I, INFO1, MAXMN, MINGG, MINMN, MINWRK, OPTWRK
      LOGICAL           ISDISC, ISFACT, ISTRAN, LQUERY
C     .. Local Arrays ..
      DOUBLE PRECISION  E1(2,2)
      LOGICAL           BWORK(1)
C     .. External Functions ..
      DOUBLE PRECISION  DLAMCH, DLAPY2
      LOGICAL           DELCTG, LSAME
      EXTERNAL          DELCTG, DLAMCH, DLAPY2, LSAME
C     .. External Subroutines ..
      EXTERNAL          DCOPY, DGEGS, DGEMM, DGEMV, DGEQRF, DGERQF,
     $                  DGGES, DLACPY, DLAG2, DLASET, DSCAL, DTRMM,
     $                  SG03BU, SG03BV, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, DBLE, INT, MAX, MIN, SIGN
C     .. Executable Statements ..
C
C     Decode input parameters.
C
      ISDISC = LSAME( DICO,  'D' )
      ISFACT = LSAME( FACT,  'F' )
      ISTRAN = LSAME( TRANS, 'T' )
      LQUERY = LDWORK.EQ.-1
C
C     Check the scalar input parameters.
C
      INFO = 0
      IF ( .NOT.( ISDISC .OR. LSAME( DICO, 'C' ) ) ) THEN
         INFO = -1
      ELSEIF ( .NOT.( ISFACT .OR. LSAME( FACT,  'N' ) ) ) THEN
         INFO = -2
      ELSEIF ( .NOT.( ISTRAN .OR. LSAME( TRANS, 'N' ) ) ) THEN
         INFO = -3
      ELSEIF ( N .LT. 0 ) THEN
         INFO = -4
      ELSEIF ( M .LT. 0 ) THEN
         INFO = -5
      ELSEIF ( LDA .LT. MAX( 1, N ) ) THEN
         INFO = -7
      ELSEIF ( LDE .LT. MAX( 1, N ) ) THEN
         INFO = -9
      ELSEIF ( LDQ .LT. MAX( 1, N ) ) THEN
         INFO = -11
      ELSEIF ( LDZ .LT. MAX( 1, N ) ) THEN
         INFO = -13
      ELSEIF ( ( ISTRAN .AND. ( LDB .LT. MAX( 1, N ) ) ) .OR.
     $    ( .NOT.ISTRAN .AND. ( LDB .LT. MAX( 1, M, N ) ) ) ) THEN
         INFO = -15
      ELSE
C
C        Compute minimal and optimal workspace.
C
         IF (ISFACT ) THEN
            MINWRK = MAX( 1, 2*N, 6*N-6 )
         ELSE
            MINWRK = MAX( 1, 4*N, 6*N-6 )
         END IF
         MINGG = MAX( MINWRK, 8*N + 16 )
         IF( LQUERY ) THEN
            CALL DGGES( 'Vectors', 'Vectors', 'Not ordered', DELCTG, N,
     $                  A, LDA, E, LDE, I, ALPHAR, ALPHAI, BETA, Q, LDQ,
     $                  Z, LDZ, DWORK, -1, BWORK, INFO1 )
            MAXMN  = MAX( M, N )
            OPTWRK = MAX( MINGG, INT( DWORK(1) ), N*MAXMN )
            IF ( .NOT.ISTRAN ) THEN
               CALL DGEQRF( MAXMN, N, B, LDB, DWORK, DWORK, -1, INFO1 )
            ELSE
               CALL DGERQF( N, MAXMN, B, LDB, DWORK, DWORK, -1, INFO1 )
            END IF
            OPTWRK = MAX( OPTWRK, INT( DWORK(1) ) + N )
         ELSEIF ( LDWORK .LT. MINWRK ) THEN
            INFO = -21
         END IF
      END IF
C
      IF ( INFO .NE. 0 ) THEN
         CALL XERBLA( 'SG03BD', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         DWORK(1) = OPTWRK
         RETURN
      END IF
C
      SCALE = ONE
C
C     Quick return if possible.
C
      MINMN = MIN( M, N )
      IF ( MINMN .EQ. 0 ) THEN
         IF ( N.GT.0 )
     $      CALL DLASET( 'Full', N, N, ZERO, ZERO, B, LDB )
         DWORK(1) = ONE
         RETURN
      ENDIF
C
      IF ( ISFACT ) THEN
C
C        Make sure the upper Hessenberg part of A is quasitriangular.
C
         DO 20 I = 1, N-2
            IF ( A(I+1,I).NE.ZERO .AND. A(I+2,I+1).NE.ZERO ) THEN
               INFO = 2
               RETURN
            END IF
   20    CONTINUE
      END IF
C
      IF ( .NOT.ISFACT ) THEN
C
C        Reduce the pencil A - lambda * E to generalized Schur form.
C
C           A := Q**T * A * Z   (upper quasitriangular)
C           E := Q**T * E * Z   (upper triangular)
C
         IF ( LDWORK .LT. MINGG ) THEN
C
C           Use DGEGS for backward compatibilty with LDWORK value.
C           ( Workspace: >= MAX(1,4*N) )
C
            CALL DGEGS( 'Vectors', 'Vectors', N, A, LDA, E, LDE, ALPHAR,
     $                  ALPHAI, BETA, Q, LDQ, Z, LDZ, DWORK, LDWORK,
     $                  INFO1 )
         ELSE
C
C           Use DGGES. The workspace is increased to avoid an error
C           return, while it should not really be larger than above.
C           ( Workspace: >= MAX(1,8*N+16) )
C
            CALL DGGES( 'Vectors', 'Vectors', 'Not ordered', DELCTG, N,
     $                  A, LDA, E, LDE, I, ALPHAR, ALPHAI, BETA, Q, LDQ,
     $                  Z, LDZ, DWORK, LDWORK, BWORK, INFO1 )
         END IF
         IF ( INFO1 .NE. 0 ) THEN
            INFO = 4
            RETURN
         END IF
         OPTWRK = INT( DWORK(1) )
      ELSE
         OPTWRK = MINWRK
      END IF
C
      IF ( ISFACT ) THEN
C
C        If the matrix pencil A - lambda * E has been in generalized
C        Schur form on entry, compute its eigenvalues.
C
         SAFMIN = DLAMCH( 'Safe minimum' )
         E1(2,1) = ZERO
         I = 1
C        WHILE ( I .LE. N ) DO
   30    IF ( I .LE. N ) THEN
            IF ( ( I.EQ.N ) .OR. ( A(MIN( I+1, N ),I).EQ.ZERO ) ) THEN
               ALPHAR(I) = A(I,I)
               ALPHAI(I) = ZERO
               BETA(I) = E(I,I)
               I = I+1
            ELSE
               E1(1,1) = E(I,I)
               E1(1,2) = E(I,I+1)
               E1(2,2) = E(I+1,I+1)
               CALL DLAG2( A(I,I), LDA, E1, 2, SAFMIN, S1, S2, WR1, WR2,
     $                     WI )
               IF ( WI .EQ. ZERO ) INFO = 3
               ALPHAR(I) = WR1
               ALPHAI(I) = WI
               BETA(I) = S1
               ALPHAR(I+1) = WR2
               ALPHAI(I+1) = -WI
               BETA(I+1) = S2
               I = I+2
            END IF
         GOTO 30
         END IF
C        END WHILE 30
         IF ( INFO.NE.0 ) RETURN
      END IF
C
C     Check on the stability of the matrix pencil A - lambda * E.
C
      DO 40 I = 1, N
         IF ( ISDISC ) THEN
            IF ( DLAPY2( ALPHAR(I), ALPHAI(I) ) .GE. ABS( BETA(I) ) )
     $         THEN
               INFO = 6
               RETURN
            END IF
         ELSE
            IF ( ( ALPHAR(I).EQ.ZERO ) .OR. ( BETA(I).EQ.ZERO ) .OR.
     $         ( SIGN( ONE,ALPHAR(I) )*SIGN( ONE, BETA(I) ) .GE. ZERO) )
     $         THEN
               INFO = 5
               RETURN
            END IF
         END IF
   40 CONTINUE
C
C     Transformation of the right hand side.
C
C        B := B * Z  or  B := Q**T * B
C
C     Use BLAS 3 if there is enough workspace. Otherwise, use BLAS 2.
C
C     ( Workspace: max(1,N) )
C
      IF ( .NOT.ISTRAN ) THEN
         IF ( LDWORK .GE. N*M ) THEN
            CALL DGEMM(  'NoTranspose', 'NoTranspose', M, N, N, ONE, B,
     $                   LDB, Z, LDZ, ZERO, DWORK, M )
            CALL DLACPY( 'All', M, N, DWORK, M, B, LDB )
         ELSE
            DO 60 I = 1, M
               CALL DCOPY( N, B(I,1), LDB, DWORK, 1 )
               CALL DGEMV( 'Transpose', N, N, ONE, Z, LDZ, DWORK, 1,
     $                     ZERO, B(I,1), LDB )
 60         CONTINUE
         END IF
         IF ( M .LT. N )
     $      CALL DLASET( 'All', N-M, N, ZERO, ZERO, B(M+1,1), LDB )
      ELSE
         IF ( LDWORK .GE. N*M ) THEN
            CALL DLACPY( 'All', N, M, B, LDB, DWORK, N )
            CALL DGEMM(  'Transpose', 'NoTranspose', N, M, N, ONE, Q,
     $                   LDQ, DWORK, N, ZERO, B, LDB )
         ELSE
            DO 80 I = 1, M
               CALL DCOPY( N, B(1,I), 1, DWORK, 1 )
               CALL DGEMV( 'Transpose', N, N, ONE, Q, LDQ, DWORK, 1,
     $                     ZERO, B(1,I), 1 )
 80         CONTINUE
         END IF
         IF ( M .LT. N )
     $      CALL DLASET( 'All', N, N-M, ZERO, ZERO, B(1,M+1), LDB )
      END IF
      OPTWRK = MAX( OPTWRK, N*M )
C
C     Overwrite B with the triangular matrix of its QR-factorization
C     or its RQ-factorization.
C     (The entries on the main diagonal are non-negative.)
C
C     ( Workspace: >= max(1,2*N) )
C
      IF ( .NOT.ISTRAN ) THEN
         IF ( M .GE. 2 ) THEN
            CALL DGEQRF( M, N, B, LDB, DWORK, DWORK(N+1), LDWORK-N,
     $                   INFO1 )
            CALL DLASET( 'Lower', MAX( M, N )-1, MIN( M, N ), ZERO,
     $                   ZERO, B(2,1), LDB )
         END IF
         DO 100 I = 1, MINMN
            IF ( B(I,I) .LT. ZERO )
     $         CALL DSCAL( N+1-I, MONE, B(I,I), LDB )
  100    CONTINUE
      ELSE
         IF ( M .GE. 2 ) THEN
            CALL DGERQF( N, M, B, LDB, DWORK, DWORK(N+1), LDWORK-N,
     $                   INFO1 )
            IF ( N .GE. M ) THEN
               CALL DLASET( 'Lower', M-1, M-1, ZERO, ZERO, B(N-M+2,1),
     $                      LDB )
               IF ( N .GT. M ) THEN
                  DO 120 I = M, 1, -1
                     CALL DCOPY( N, B(1,I), 1, B(1,I+N-M), 1 )
  120             CONTINUE
                  CALL DLASET( 'All', N, N-M, ZERO, ZERO, B(1,1), LDB )
               END IF
            ELSE
               IF ( N .GT. 1 )
     $            CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO,
     $                         B(2,M-N+1), LDB )
               DO 140 I = 1, N
                  CALL DCOPY( N, B(1,M-N+I), 1, B(1,I), 1 )
  140          CONTINUE
               CALL DLASET( 'All', N, M-N, ZERO, ZERO, B(1,N+1), LDB )
            END IF
         ELSE
            IF ( N .NE. 1 ) THEN
               CALL DCOPY( N, B(1,1), 1, B(1,N), 1 )
               CALL DLASET( 'All', N, 1, ZERO, ZERO, B(1,1), LDB )
            END IF
         END IF
         DO 160 I = N - MINMN + 1, N
            IF ( B(I,I) .LT. ZERO )
     $         CALL DSCAL( I, MONE, B(1,I), 1 )
  160    CONTINUE
      END IF
      OPTWRK = MAX( OPTWRK, INT( DWORK(N+1) ) + N )
C
C     Solve the reduced generalized Lyapunov equation.
C
C     ( Workspace: 6*N-6 )
C
      IF ( ISDISC ) THEN
         CALL SG03BU( TRANS, N, A, LDA, E, LDE, B, LDB, SCALE, DWORK,
     $                INFO1 )
         IF ( INFO1 .NE. 0 ) THEN
            IF ( INFO1 .EQ. 1 ) INFO = 1
            IF ( INFO1 .EQ. 2 ) INFO = 3
            IF ( INFO1 .EQ. 3 ) INFO = 6
            IF ( INFO1 .EQ. 4 ) INFO = 7
            IF ( INFO  .NE. 1 )
     $         RETURN
         END IF
      ELSE
         CALL SG03BV( TRANS, N, A, LDA, E, LDE, B, LDB, SCALE, DWORK,
     $                INFO1 )
         IF ( INFO1 .NE. 0 ) THEN
            IF ( INFO1 .EQ. 1 ) INFO = 1
            IF ( INFO1 .GE. 2 ) INFO = 3
            IF ( INFO1 .EQ. 3 ) INFO = 5
            IF ( INFO  .NE. 1 )
     $         RETURN
         END IF
      END IF
C
C     Transform the solution matrix back.
C
C        U := U * Q**T   or   U := Z * U
C
C     Use BLAS 3 if there is enough workspace. Otherwise, use BLAS 2.
C
C     ( Workspace: max(1,N) )
C
      IF ( .NOT.ISTRAN ) THEN
         IF ( LDWORK .GE. N*N ) THEN
            CALL DLACPY( 'All', N, N, Q, LDQ, DWORK, N )
            CALL DTRMM( 'Right', 'Upper', 'Transpose', 'NonUnit', N, N,
     $                  ONE, B, LDB, DWORK, N)
            DO 170 I = 1, N
               CALL DCOPY( N, DWORK(N*(I-1)+1), 1, B(I,1), LDB )
  170       CONTINUE
         ELSE
            DO 180 I = 1, N
               CALL DCOPY( N-I+1, B(I,I), LDB, DWORK, 1 )
               CALL DGEMV( 'NoTranspose', N, N-I+1, ONE, Q(1,I), LDQ,
     $                     DWORK, 1, ZERO, B(I,1), LDB )
  180       CONTINUE
         END IF
      ELSE
         IF ( LDWORK .GE. N*N ) THEN
            CALL DLACPY( 'All', N, N, Z, LDZ, DWORK, N )
            CALL DTRMM(  'Right', 'Upper', 'NoTranspose', 'NonUnit', N,
     $                   N, ONE, B, LDB, DWORK, N )
            CALL DLACPY( 'All', N, N, DWORK, N, B, LDB )
         ELSE
            DO 200 I = 1, N
               CALL DCOPY( I, B(1,I), 1, DWORK, 1 )
               CALL DGEMV( 'NoTranspose', N, I, ONE, Z, LDZ, DWORK, 1,
     $                     ZERO, B(1,I), 1 )
 200        CONTINUE
         END IF
      END IF
      OPTWRK = MAX( OPTWRK, N*N )
C
C     Overwrite U with the triangular matrix of its QR-factorization
C     or its RQ-factorization.
C     (The entries on the main diagonal are non-negative.)
C
C     ( Workspace: >= max(1,2*N) )
C
      IF ( .NOT.ISTRAN ) THEN
         CALL DGEQRF( N, N, B, LDB, DWORK, DWORK(N+1), LDWORK-N, INFO1 )
         IF ( N .GT. 1 )
     $      CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO, B(2,1), LDB )
         DO 220 I = 1, N
            IF ( B(I,I) .LT. ZERO )
     $         CALL DSCAL( N+1-I, MONE, B(I,I), LDB )
  220    CONTINUE
      ELSE
         CALL DGERQF( N, N, B, LDB, DWORK, DWORK(N+1), LDWORK-N, INFO1 )
         IF ( N .GT. 1 )
     $      CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO, B(2,1), LDB )
         DO 240 I = 1, N
            IF ( B(I,I) .LT. ZERO )
     $         CALL DSCAL( I, MONE, B(1,I), 1 )
  240    CONTINUE
      END IF
      OPTWRK = MAX( OPTWRK, INT( DWORK(N+1) ) + N )
C
      DWORK(1) = DBLE( MAX( OPTWRK, MINWRK ) )
      RETURN
C *** Last line of SG03BD ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
      double complex function Fcc_qpgpgpqm(j1,j2,j3,j4,j5,j6,za,zb)
      implicit none
      integer j1,j2,j3,j4,j5,j6
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      double complex L0,Lsm1,Lsm1_2me
      double precision t

      Fcc_qpgpgpqm=
     .-(((-Lsm1(-s(j1,j2),-t(j1,j2,j3),-s(j2,j3),-t(j1,j2,j3))-
     .Lsm1(-s(j2,j3),-t(j2,j3,j4),-s(j3,j4),-t(j2,j3,j4))-
     .Lsm1_2me(t(j1,j2,j3),t(j2,j3,j4),s(j2,j3),s(j5,j6)))*za(j4,j5)**2)
     ./
     .(za(j1,j2)*za(j2,j3)*za(j3,j4)*za(j5,j6)))-
     .(2d0*L0(-t(j2,j3,j4),-s(j3,j4))*za(j2,j5)*za(j4,j5)*zb(j2,j3))/
     .(s(j3,j4)*za(j1,j2)*za(j2,j3)*za(j5,j6))-
     .(2d0*L0(-s(j5,j6),-t(j2,j3,j4))*za(j1,j5)*za(j4,j5)*
     .(za(j2,j4)*zb(j1,j2)+za(j3,j4)*zb(j1,j3)))/
     .(za(j1,j2)*za(j2,j3)*za(j3,j4)*za(j5,j6)*t(j2,j3,j4))

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
Tau Kappa Epsilon (often shortened to TKE, pronounced like Teke) is a fraternity at UC Davis.  http://www.tke.org Tau Kappa Epsilon is the largest college social fraternity by number of chapters worldwide, with chapters across the US and Canada. The Davis chapter of TKE was founded in 1989 and reaches around 60 active members every year. As a fraternity, TKE has strong involvement at UC Davis by participating in several philanthropies and performing community service as well as hosting their annual philanthropy, TKE Sweetheart, which benefits St. Jude Childrens Research Hospital.

Beyond participating in the fraternity, Tekes are also involved with other organizations such as ASUCD, IFC, Clinica Tepati, BloodSource, and various sports teams including soccer, tennis, and rugby. Additionally, TKE members have majors in and pursue careers in all fields of study, with recent graduates having received acceptances from UC Davis Medical, Johns Hopkins Medical, Stanford School of Business, Georgetown Law, USC law, and UCLA Law.

Recent Chapter Awards:
1st Place  2014 EAS Grilloff & Volleyball Tournament
3rd Place  2014 Kappa Alpha Thetas Mr. Casanova
2014 IFC Award of Excellence: Most Outstanding Member
2nd Place  2013 Delta Gammas Anchorsplash
Most Creative Performance  2013 Pi Beta Phis Arrow Jam
1st Place  2013 Kappa Alpha Thetas Mr. Casanova
1st Place  2012 Alpha Chi Omega Runway
2012 IFC Award of Excellence: New Member Evaluation and Recruitment
3rd Place  2012 Kappa Alpha Thetas Mr. Casanova
2nd Place  2010 Kappa Alpha Thetas Mr. Casanova

ive never heard anything about their parties or philanthropy, but i do see their posters featuring Ronnie Reagan everywhere. also, was the above photo taken at the bottom of the ocean? can we do any better than that?  Users/ArlenAbraham

Their iceblocking parties were a frequent event when I was a freshman and had to walk to parties. Imagine 100 strangers putting their mouths on the same indentation on a block of ice and the TKEs pouring Kamikazes down the channel. Console your fears of salivaborne diseases with the ideas that alcohol kills everything and that hopefully the ice has melted sufficiently since the last person put their mouth there that youve got fresh ice. Also dont expect to feel your mouth again for a while.    Users/DanielMedinaCleghorn

20060430 02:35:08 nbsp Im glad my old party house (502 Oeste) is being put to good use.  Users/JaimeRaba

20060430 02:47:16 nbsp My brother is a TKE (albeit on the opposite side of the country).  Ive heard good things.  Anyone who makes Michael Morris black out running a tread mill is ok in my book. Users/ZacMorris

20120317 19:18:43 nbsp i do know they have a really awesome philanthropy with TKE Sweetheart and most of the people i know have nothing but good things to say about them. Seem pretty down to earth too. theyve been exploding in recent years on campus, from what ive seen.  Users/shaikhshake

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      PROGRAM TPWRZS
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
C OPEN GKS, OPEN WORKSTATION OF TYPE 1, ACTIVATE WORKSTATION
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C INVOKE DEMO DRIVER
C
      CALL TPWRZS1(IERR)
C
C DEACTIVATE AND CLOSE WORKSTATION, CLOSE GKS.
C
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS
C
      STOP
      END
C
      SUBROUTINE TPWRZS1 (IERROR)
C
C PURPOSE                To provide a simple demonstration of
C                        entry PWRZS with the SRFACE utility.
C
C USAGE                  CALL TPWRZS1 (IERROR)
C
C ARGUMENTS
C
C ON OUTPUT              IERROR
C                          An integer variable
C                          = 0, if the test was successful,
C                          = 1, otherwise
C
C I/O                    If the test is successful, the message
C
C               PWRZS TEST EXECUTED--SEE PLOTS TO CERTIFY
C
C                        is printed on unit 6.  In addition, 1
C                        frame is produced on the machine graphics
C                        device.  In order to determine if the test
C                        was successful, it is necessary to examine
C                        the plot.
C
C PRECISION              Single
C
C REQUIRED ROUTINES      PWRZS, SRFACE
C
C REQUIRED GKS LEVEL     0A
C
C LANGUAGE               FORTRAN 77
C
C ALGORITHM              A function of 2 variables is defined and the
C                        values of the function on a 2-D rectangular
C                        grid are stored in an array.  This routine
C                        calls SRFACE to draw a surface representation
C                        of the array values.  PWRZS is then called 3
C                        times to label the front, side, and back of
C                        the picture.
C
      DIMENSION       Z(20,30)   ,X(20)      ,Y(30)      ,MM(20,30,2),
     1                S(6)
C
C Load the SRFACE common block needed to supress a NEWFM call.
C
      COMMON /SRFIP1/ IFR        ,ISTP       ,IROTS      ,IDRX       ,
     1                IDRY       ,IDRZ       ,IUPPER     ,ISKIRT     ,
     2                NCLA       ,THETA      ,HSKIRT     ,CHI        ,
     3                CLO        ,CINC       ,ISPVAL
C
C Define the center of a plot title string on a square grid of size
C 0. to 1.
C
      DATA TX/0.4375/, TY/0.9667/
C
C Specify grid loop indices and a line of sight.
C
      DATA M/20/, N/30/
      DATA S/4.,5.,3.,0.,0.,0./
C
C Initial the error parameter.
C
      IERROR = 1
C
C Define the function values and store them in the Z array.
C
      DO  10 I=1,M
         X(I) = -1.+REAL(I-1)/REAL(M-1)*2.
   10 CONTINUE
      DO  20 J=1,N
         Y(J) = -1.+REAL(J-1)/REAL(N-1)*2.
   20 CONTINUE
      DO  40 J=1,N
         DO  30 I=1,M
            Z(I,J) = EXP(-2.*SQRT(X(I)**2+Y(J)**2))
   30    CONTINUE
   40 CONTINUE
C
C Set SRFACE parameters to supress the FRAME call and draw contours.
C
      IFR = 0
      IDRZ = 1
C
C Select normalization trans number 0.
C
      CALL GSELNT (0)
C
C Label the plot.
C
      CALL PLCHLQ (TX,TY,'DEMONSTRATION PLOT FOR PWRZS',16.,0.,0.)
C
C Draw the surface plot.
C
      CALL SRFACE (X,Y,Z,MM,M,M,N,S,0.)
C
C Put the PWRZS labels on the picture.
C
      ISIZE = 35
      CALL PWRZS (0.,1.1,0.,'FRONT',5,ISIZE,-1,3,0)
      CALL PWRZS (1.1,0.,0.,'SIDE',4,ISIZE,2,-1,0)
      CALL PWRZS (0.,-1.1,.2,' BACK BACK BACK BACK BACK',25,ISIZE,-1,
     1            3,0)
      CALL FRAME
C
      IERROR = 0
      WRITE (6,1001)
C
C Restore the SRFACE parameters to their default values.
C
      IFR = 1
      IDRZ = 0
C
      RETURN
C
 1001 FORMAT (' PWRZS TEST EXECUTED--SEE PLOT TO CERTIFY')
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


program test_get_output_item_count

  use bmif_2_0, only: BMI_FAILURE
  use bmiprmsstreamflow
  use fixtures, only: status

  implicit none

  integer, parameter :: expected = 22
  type (bmi_prms_streamflow) :: m
  integer :: count

  status = m%get_output_item_count(count)

  if (count /= expected) then
     stop BMI_FAILURE
  end if
end program test_get_output_item_count


text:
 Wiki mapping
Pages in the wiki can be mapped.  See wiki:wikispot:Help with Maps on the Wiki Spot hub for instructions on mapping pages.  When a page has been given a map location a Map icon will appear on the page.  Click it to bring up a map of the pages location along with locations of nearby places.

 General Davis map information

The city of Davis has maintained a rather extensive collection of useful and interesting maps over at their http://www.cityofdavis.org/gis/ Geographic Information System site.  The http://www.cityofdavis.org/gis/library/ downloadable map layers page is especially useful for the wiki.  You can download .zip files containing layer files, which you can view with:
  http://www.esri.com/software/arcexplorer/download.html free GIS software from ESRI  Windows only :(
  http://qgis.org Free open source GIS software (like the wiki)  Windows, Mac & Linux
These files contain a lot of useful data.  For instance: the parks file contains the locations/shapes of all the parks, size of all the parks, and their names.

Among the most interesting are the 2005 aerial photos, which you can http://www2.dcn.org/orgs/orthophotos/ download individually in high resolution or http://www.cityofdavis.org/gis/2003aerial browse seamlessly using flash.

The http://www.cityofdavis.org/gis/unitrans Unitrans route map is likely the most useful to students, though.

The city also has a http://www.city.davis.ca.us/gis/bikemap.pdf bike map available in PDF form. This map shows Bike Paths bike paths in addition to city streets.

OpenStreetMap has very good coverage of Davis, and you can improve it too.

20070611 18:57:32 nbsp So what page does one go to view an overview map of all wiki locations? Users/AlexMandel

   Right now you cant.  But its wiki:wikispot:Feature Requests planned.

20070616 20:11:36 nbsp exactly? If we are using google maps now... then.... how do we use it? Can we link it to each page like it used to be. This was one of my favorite features of the wiki. Users/CarlosBarahona

20070616 20:41:07 nbsp Maybe a kml file could be available for download? Users/JasonAller

20071105 21:11:46 nbsp Im currently working to integrate a new solution. I have a plan and all, anyone want to help me out? Users/AlexMandel

20090429 12:52:58 nbsp I miss the old Dwiki map. It was not the best one, but at least it was something. Why did it go? Users/JoePomidor

20120604 13:55:24 nbsp When will we get the bees knees? https://trianglewiki.org/map/#zoom12&lat35.78754&lon78.67784&layersBTT Users/BillKendrick


text:
Pam Wilson is a certified holistic health practitioner, certified massage therapist with 1000+ hours training, has a masters degree in counseling psychology, and has trained and practiced in aromatherapy, biofeedback, meditation, hypnotherapy, breathwork and emotional release work. Pams intuition facilitates energy movement through an escalen approach (hypnotic method using slow rhythmic movement, bringing about a sense of well being through deep states of consciousness). The method focuses on mindbodyspirit as a whole.

for more information, see Vitality for Life

20080706 11:55:25 nbsp Welcome to the Wiki. The reason the comment you removed was restored to the Vitality for Life page is that as a Wiki Community community the wiki has decided on certain behaviors. It might be worth the time to read Welcome to the Wiki/Business Owner; it will help explain how you can make the wiki a positive experience for you without clashing with established wiki social norms. Users/JasonAller


text:
!> -------------------------------------------
!! Supply regression utilities to GARD
!!
!! Provides least squares regression and logistic regression
!! includes output of both the coefficients and the residual error, optionally computed with supplied weights
!!
!! --------------------------------------------
module  regression_mod
interface
    module function compute_regression(x, training_x, training_y, coefficients, y_test, error, weights, used_vars) result(y)
        implicit none
        real,    intent(in),    dimension(:)   :: x
        real,    intent(in),    dimension(:,:) :: training_x
        real,    intent(in),    dimension(:)   :: training_y
        real(8), intent(inout), dimension(:)   :: coefficients
        real,    intent(inout),                             optional :: error
        real,    intent(in),    dimension(:),               optional :: y_test
        real,    intent(inout), dimension(:),  allocatable, optional :: weights
        integer, intent(inout), dimension(:),               optional :: used_vars
        real :: y
    end function compute_regression

    module function compute_logistic_regression(x, training_x, binary_values, coefficients) result(output)
        implicit none
        real,    intent(in) :: x(:)
        real,    intent(in) :: training_x(:,:)
        real,    intent(in) :: binary_values(:)
        real(8), intent(out):: coefficients(:)
        real :: output
    end function compute_logistic_regression

end interface

end module regression_mod

max_stars_repo_name: benlazarine/GARD

max_issues_repo_name: benlazarine/GARD

max_forks_repo_name: benlazarine/GARD


text:
*> \brief \b DGESDD
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DGESDD + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
*                          LWORK, IWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          JOBZ
*       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
*       ..
*       .. Array Arguments ..
*       INTEGER            IWORK( * )
*       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
*      $                   VT( LDVT, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGESDD computes the singular value decomposition (SVD) of a real
*> M-by-N matrix A, optionally computing the left and right singular
*> vectors.  If singular vectors are desired, it uses a
*> divide-and-conquer algorithm.
*>
*> The SVD is written
*>
*>      A = U * SIGMA * transpose(V)
*>
*> where SIGMA is an M-by-N matrix which is zero except for its
*> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
*> V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
*> are the singular values of A; they are real and non-negative, and
*> are returned in descending order.  The first min(m,n) columns of
*> U and V are the left and right singular vectors of A.
*>
*> Note that the routine returns VT = V**T, not V.
*>
*> The divide and conquer algorithm makes very mild assumptions about
*> floating point arithmetic. It will work on machines with a guard
*> digit in add/subtract, or on those binary machines without guard
*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
*> without guard digits, but we know of none.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOBZ
*> \verbatim
*>          JOBZ is CHARACTER*1
*>          Specifies options for computing all or part of the matrix U:
*>          = 'A':  all M columns of U and all N rows of V**T are
*>                  returned in the arrays U and VT;
*>          = 'S':  the first min(M,N) columns of U and the first
*>                  min(M,N) rows of V**T are returned in the arrays U
*>                  and VT;
*>          = 'O':  If M >= N, the first N columns of U are overwritten
*>                  on the array A and all rows of V**T are returned in
*>                  the array VT;
*>                  otherwise, all columns of U are returned in the
*>                  array U and the first M rows of V**T are overwritten
*>                  in the array A;
*>          = 'N':  no columns of U or rows of V**T are computed.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the input matrix A.  M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the input matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA,N)
*>          On entry, the M-by-N matrix A.
*>          On exit,
*>          if JOBZ = 'O',  A is overwritten with the first N columns
*>                          of U (the left singular vectors, stored
*>                          columnwise) if M >= N;
*>                          A is overwritten with the first M rows
*>                          of V**T (the right singular vectors, stored
*>                          rowwise) otherwise.
*>          if JOBZ .ne. 'O', the contents of A are destroyed.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,M).
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*>          S is DOUBLE PRECISION array, dimension (min(M,N))
*>          The singular values of A, sorted so that S(i) >= S(i+1).
*> \endverbatim
*>
*> \param[out] U
*> \verbatim
*>          U is DOUBLE PRECISION array, dimension (LDU,UCOL)
*>          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
*>          UCOL = min(M,N) if JOBZ = 'S'.
*>          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
*>          orthogonal matrix U;
*>          if JOBZ = 'S', U contains the first min(M,N) columns of U
*>          (the left singular vectors, stored columnwise);
*>          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER
*>          The leading dimension of the array U.  LDU >= 1; if
*>          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
*> \endverbatim
*>
*> \param[out] VT
*> \verbatim
*>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
*>          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
*>          N-by-N orthogonal matrix V**T;
*>          if JOBZ = 'S', VT contains the first min(M,N) rows of
*>          V**T (the right singular vectors, stored rowwise);
*>          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
*> \endverbatim
*>
*> \param[in] LDVT
*> \verbatim
*>          LDVT is INTEGER
*>          The leading dimension of the array VT.  LDVT >= 1; if
*>          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
*>          if JOBZ = 'S', LDVT >= min(M,N).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The dimension of the array WORK. LWORK >= 1.
*>          If JOBZ = 'N',
*>            LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)).
*>          If JOBZ = 'O',
*>            LWORK >= 3*min(M,N) +
*>                     max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
*>          If JOBZ = 'S' or 'A'
*>            LWORK >= 3*min(M,N) +
*>                     max(max(M,N),4*min(M,N)*min(M,N)+3*min(M,N)+max(M,N)).
*>          For good performance, LWORK should generally be larger.
*>          If LWORK = -1 but other input arguments are legal, WORK(1)
*>          returns the optimal LWORK.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (8*min(M,N))
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit.
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*>          > 0:  DBDSDC did not converge, updating process failed.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date September 2012
*
*> \ingroup doubleGEsing
*
*> \par Contributors:
*  ==================
*>
*>     Ming Gu and Huan Ren, Computer Science Division, University of
*>     California at Berkeley, USA
*>
*  =====================================================================
      SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
     $                   LWORK, IWORK, INFO )
*
*  -- LAPACK driver routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      CHARACTER          JOBZ
      INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
     $                   VT( LDVT, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
      INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
     $                   IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
     $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
     $                   MNTHR, NWORK, WRKBL
      DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
*     ..
*     .. Local Arrays ..
      INTEGER            IDUM( 1 )
      DOUBLE PRECISION   DUM( 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY,
     $                   DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR,
     $                   XERBLA
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      DOUBLE PRECISION   DLAMCH, DLANGE
      EXTERNAL           DLAMCH, DLANGE, ILAENV, LSAME
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          INT, MAX, MIN, SQRT
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments
*
      INFO = 0
      MINMN = MIN( M, N )
      WNTQA = LSAME( JOBZ, 'A' )
      WNTQS = LSAME( JOBZ, 'S' )
      WNTQAS = WNTQA .OR. WNTQS
      WNTQO = LSAME( JOBZ, 'O' )
      WNTQN = LSAME( JOBZ, 'N' )
      LQUERY = ( LWORK.EQ.-1 )
*
      IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
         INFO = -1
      ELSE IF( M.LT.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -5
      ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
     $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
         INFO = -8
      ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
     $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
     $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
         INFO = -10
      END IF
*
*     Compute workspace
*      (Note: Comments in the code beginning "Workspace:" describe the
*       minimal amount of workspace needed at that point in the code,
*       as well as the preferred amount for good performance.
*       NB refers to the optimal block size for the immediately
*       following subroutine, as returned by ILAENV.)
*
      IF( INFO.EQ.0 ) THEN
         MINWRK = 1
         MAXWRK = 1
         IF( M.GE.N .AND. MINMN.GT.0 ) THEN
*
*           Compute space needed for DBDSDC
*
            MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
            IF( WNTQN ) THEN
               BDSPAC = 7*N
            ELSE
               BDSPAC = 3*N*N + 4*N
            END IF
            IF( M.GE.MNTHR ) THEN
               IF( WNTQN ) THEN
*
*                 Path 1 (M much larger than N, JOBZ='N')
*
                  WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,
     $                    -1 )
                  WRKBL = MAX( WRKBL, 3*N+2*N*
     $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
                  MAXWRK = MAX( WRKBL, BDSPAC+N )
                  MINWRK = BDSPAC + N
               ELSE IF( WNTQO ) THEN
*
*                 Path 2 (M much larger than N, JOBZ='O')
*
                  WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
                  WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
     $                    N, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+2*N*
     $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
                  MAXWRK = WRKBL + 2*N*N
                  MINWRK = BDSPAC + 2*N*N + 3*N
               ELSE IF( WNTQS ) THEN
*
*                 Path 3 (M much larger than N, JOBZ='S')
*
                  WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
                  WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
     $                    N, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+2*N*
     $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
                  MAXWRK = WRKBL + N*N
                  MINWRK = BDSPAC + N*N + 3*N
               ELSE IF( WNTQA ) THEN
*
*                 Path 4 (M much larger than N, JOBZ='A')
*
                  WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
                  WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
     $                    M, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+2*N*
     $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
                  MAXWRK = WRKBL + N*N
                  MINWRK = BDSPAC + N*N + 2*N + M
               END IF
            ELSE
*
*              Path 5 (M at least N, but not much larger)
*
               WRKBL = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
     $                 -1 )
               IF( WNTQN ) THEN
                  MAXWRK = MAX( WRKBL, BDSPAC+3*N )
                  MINWRK = 3*N + MAX( M, BDSPAC )
               ELSE IF( WNTQO ) THEN
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
                  MAXWRK = WRKBL + M*N
                  MINWRK = 3*N + MAX( M, N*N+BDSPAC )
               ELSE IF( WNTQS ) THEN
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
                  MAXWRK = MAX( WRKBL, BDSPAC+3*N )
                  MINWRK = 3*N + MAX( M, BDSPAC )
               ELSE IF( WNTQA ) THEN
                  WRKBL = MAX( WRKBL, 3*N+M*
     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*N+N*
     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
                  MAXWRK = MAX( MAXWRK, BDSPAC+3*N )
                  MINWRK = 3*N + MAX( M, BDSPAC )
               END IF
            END IF
         ELSE IF( MINMN.GT.0 ) THEN
*
*           Compute space needed for DBDSDC
*
            MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
            IF( WNTQN ) THEN
               BDSPAC = 7*M
            ELSE
               BDSPAC = 3*M*M + 4*M
            END IF
            IF( N.GE.MNTHR ) THEN
               IF( WNTQN ) THEN
*
*                 Path 1t (N much larger than M, JOBZ='N')
*
                  WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,
     $                    -1 )
                  WRKBL = MAX( WRKBL, 3*M+2*M*
     $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
                  MAXWRK = MAX( WRKBL, BDSPAC+M )
                  MINWRK = BDSPAC + M
               ELSE IF( WNTQO ) THEN
*
*                 Path 2t (N much larger than M, JOBZ='O')
*
                  WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
                  WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
     $                    N, M, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+2*M*
     $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
                  MAXWRK = WRKBL + 2*M*M
                  MINWRK = BDSPAC + 2*M*M + 3*M
               ELSE IF( WNTQS ) THEN
*
*                 Path 3t (N much larger than M, JOBZ='S')
*
                  WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
                  WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
     $                    N, M, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+2*M*
     $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
                  MAXWRK = WRKBL + M*M
                  MINWRK = BDSPAC + M*M + 3*M
               ELSE IF( WNTQA ) THEN
*
*                 Path 4t (N much larger than M, JOBZ='A')
*
                  WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
                  WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
     $                    N, M, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+2*M*
     $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
                  MAXWRK = WRKBL + M*M
                  MINWRK = BDSPAC + M*M + 3*M
               END IF
            ELSE
*
*              Path 5t (N greater than M, but not much larger)
*
               WRKBL = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
     $                 -1 )
               IF( WNTQN ) THEN
                  MAXWRK = MAX( WRKBL, BDSPAC+3*M )
                  MINWRK = 3*M + MAX( N, BDSPAC )
               ELSE IF( WNTQO ) THEN
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
                  MAXWRK = WRKBL + M*N
                  MINWRK = 3*M + MAX( N, M*M+BDSPAC )
               ELSE IF( WNTQS ) THEN
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
                  MAXWRK = MAX( WRKBL, BDSPAC+3*M )
                  MINWRK = 3*M + MAX( N, BDSPAC )
               ELSE IF( WNTQA ) THEN
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
                  WRKBL = MAX( WRKBL, 3*M+M*
     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) )
                  MAXWRK = MAX( WRKBL, BDSPAC+3*M )
                  MINWRK = 3*M + MAX( N, BDSPAC )
               END IF
            END IF
         END IF
         MAXWRK = MAX( MAXWRK, MINWRK )
         WORK( 1 ) = MAXWRK
*
         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
            INFO = -12
         END IF
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGESDD', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
         RETURN
      END IF
*
*     Get machine constants
*
      EPS = DLAMCH( 'P' )
      SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
      BIGNUM = ONE / SMLNUM
*
*     Scale A if max element outside range [SMLNUM,BIGNUM]
*
      ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
      ISCL = 0
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
         ISCL = 1
         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
      ELSE IF( ANRM.GT.BIGNUM ) THEN
         ISCL = 1
         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
      END IF
*
      IF( M.GE.N ) THEN
*
*        A has at least as many rows as columns. If A has sufficiently
*        more rows than columns, first reduce using the QR
*        decomposition (if sufficient workspace available)
*
         IF( M.GE.MNTHR ) THEN
*
            IF( WNTQN ) THEN
*
*              Path 1 (M much larger than N, JOBZ='N')
*              No singular vectors to be computed
*
               ITAU = 1
               NWORK = ITAU + N
*
*              Compute A=Q*R
*              (Workspace: need 2*N, prefer N+N*NB)
*
               CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Zero out below R
*
               CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
               IE = 1
               ITAUQ = IE + N
               ITAUP = ITAUQ + N
               NWORK = ITAUP + N
*
*              Bidiagonalize R in A
*              (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
               CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                      IERR )
               NWORK = IE + N
*
*              Perform bidiagonal SVD, computing singular values only
*              (Workspace: need N+BDSPAC)
*
               CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
*
            ELSE IF( WNTQO ) THEN
*
*              Path 2 (M much larger than N, JOBZ = 'O')
*              N left singular vectors to be overwritten on A and
*              N right singular vectors to be computed in VT
*
               IR = 1
*
*              WORK(IR) is LDWRKR by N
*
               IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN
                  LDWRKR = LDA
               ELSE
                  LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / N
               END IF
               ITAU = IR + LDWRKR*N
               NWORK = ITAU + N
*
*              Compute A=Q*R
*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
               CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy R to WORK(IR), zeroing out below it
*
               CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
               CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
     $                      LDWRKR )
*
*              Generate Q in A
*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
               CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = ITAU
               ITAUQ = IE + N
               ITAUP = ITAUQ + N
               NWORK = ITAUP + N
*
*              Bidiagonalize R in VT, copying result to WORK(IR)
*              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
               CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              WORK(IU) is N by N
*
               IU = NWORK
               NWORK = IU + N*N
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in WORK(IU) and computing right
*              singular vectors of bidiagonal matrix in VT
*              (Workspace: need N+N*N+BDSPAC)
*
               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
     $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
*
*              Overwrite WORK(IU) by left singular vectors of R
*              and VT by right singular vectors of R
*              (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
     $                      WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Multiply Q in A by left singular vectors of R in
*              WORK(IU), storing result in WORK(IR) and copying to A
*              (Workspace: need 2*N*N, prefer N*N+M*N)
*
               DO 10 I = 1, M, LDWRKR
                  CHUNK = MIN( M-I+1, LDWRKR )
                  CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
     $                        LDA, WORK( IU ), N, ZERO, WORK( IR ),
     $                        LDWRKR )
                  CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
     $                         A( I, 1 ), LDA )
   10          CONTINUE
*
            ELSE IF( WNTQS ) THEN
*
*              Path 3 (M much larger than N, JOBZ='S')
*              N left singular vectors to be computed in U and
*              N right singular vectors to be computed in VT
*
               IR = 1
*
*              WORK(IR) is N by N
*
               LDWRKR = N
               ITAU = IR + LDWRKR*N
               NWORK = ITAU + N
*
*              Compute A=Q*R
*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
               CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy R to WORK(IR), zeroing out below it
*
               CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
               CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
     $                      LDWRKR )
*
*              Generate Q in A
*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
               CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = ITAU
               ITAUQ = IE + N
               ITAUP = ITAUQ + N
               NWORK = ITAUP + N
*
*              Bidiagonalize R in WORK(IR)
*              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
               CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagoal matrix in U and computing right singular
*              vectors of bidiagonal matrix in VT
*              (Workspace: need N+BDSPAC)
*
               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
*
*              Overwrite U by left singular vectors of R and VT
*              by right singular vectors of R
*              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
               CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Multiply Q in A by left singular vectors of R in
*              WORK(IR), storing result in U
*              (Workspace: need N*N)
*
               CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
               CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
     $                     LDWRKR, ZERO, U, LDU )
*
            ELSE IF( WNTQA ) THEN
*
*              Path 4 (M much larger than N, JOBZ='A')
*              M left singular vectors to be computed in U and
*              N right singular vectors to be computed in VT
*
               IU = 1
*
*              WORK(IU) is N by N
*
               LDWRKU = N
               ITAU = IU + LDWRKU*N
               NWORK = ITAU + N
*
*              Compute A=Q*R, copying result to U
*              (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
*
               CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
*
*              Generate Q in U
*              (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
               CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*              Produce R in A, zeroing out other entries
*
               CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
               IE = ITAU
               ITAUQ = IE + N
               ITAUP = ITAUQ + N
               NWORK = ITAUP + N
*
*              Bidiagonalize R in A
*              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
               CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                      IERR )
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in WORK(IU) and computing right
*              singular vectors of bidiagonal matrix in VT
*              (Workspace: need N+N*N+BDSPAC)
*
               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
     $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
*
*              Overwrite WORK(IU) by left singular vectors of R and VT
*              by right singular vectors of R
*              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
     $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Multiply Q in U by left singular vectors of R in
*              WORK(IU), storing result in A
*              (Workspace: need N*N)
*
               CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
     $                     LDWRKU, ZERO, A, LDA )
*
*              Copy left singular vectors of A from A to U
*
               CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
*
            END IF
*
         ELSE
*
*           M .LT. MNTHR
*
*           Path 5 (M at least N, but not much larger)
*           Reduce to bidiagonal form without QR decomposition
*
            IE = 1
            ITAUQ = IE + N
            ITAUP = ITAUQ + N
            NWORK = ITAUP + N
*
*           Bidiagonalize A
*           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
*
            CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
     $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                   IERR )
            IF( WNTQN ) THEN
*
*              Perform bidiagonal SVD, only computing singular values
*              (Workspace: need N+BDSPAC)
*
               CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
            ELSE IF( WNTQO ) THEN
               IU = NWORK
               IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
*
*                 WORK( IU ) is M by N
*
                  LDWRKU = M
                  NWORK = IU + LDWRKU*N
                  CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
     $                         LDWRKU )
               ELSE
*
*                 WORK( IU ) is N by N
*
                  LDWRKU = N
                  NWORK = IU + LDWRKU*N
*
*                 WORK(IR) is LDWRKR by N
*
                  IR = NWORK
                  LDWRKR = ( LWORK-N*N-3*N ) / N
               END IF
               NWORK = IU + LDWRKU*N
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in WORK(IU) and computing right
*              singular vectors of bidiagonal matrix in VT
*              (Workspace: need N+N*N+BDSPAC)
*
               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
     $                      LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
     $                      IWORK, INFO )
*
*              Overwrite VT by right singular vectors of A
*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
               CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
               IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
*
*                 Overwrite WORK(IU) by left singular vectors of A
*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                  CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
     $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*                 Copy left singular vectors of A from WORK(IU) to A
*
                  CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
               ELSE
*
*                 Generate Q in A
*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
                  CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*                 Multiply Q in A by left singular vectors of
*                 bidiagonal matrix in WORK(IU), storing result in
*                 WORK(IR) and copying to A
*                 (Workspace: need 2*N*N, prefer N*N+M*N)
*
                  DO 20 I = 1, M, LDWRKR
                     CHUNK = MIN( M-I+1, LDWRKR )
                     CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
     $                           LDA, WORK( IU ), LDWRKU, ZERO,
     $                           WORK( IR ), LDWRKR )
                     CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
     $                            A( I, 1 ), LDA )
   20             CONTINUE
               END IF
*
            ELSE IF( WNTQS ) THEN
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in U and computing right singular
*              vectors of bidiagonal matrix in VT
*              (Workspace: need N+BDSPAC)
*
               CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
*
*              Overwrite U by left singular vectors of A and VT
*              by right singular vectors of A
*              (Workspace: need 3*N, prefer 2*N+N*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
            ELSE IF( WNTQA ) THEN
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in U and computing right singular
*              vectors of bidiagonal matrix in VT
*              (Workspace: need N+BDSPAC)
*
               CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
*
*              Set the right corner of U to identity matrix
*
               IF( M.GT.N ) THEN
                  CALL DLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ),
     $                         LDU )
               END IF
*
*              Overwrite U by left singular vectors of A and VT
*              by right singular vectors of A
*              (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
            END IF
*
         END IF
*
      ELSE
*
*        A has more columns than rows. If A has sufficiently more
*        columns than rows, first reduce using the LQ decomposition (if
*        sufficient workspace available)
*
         IF( N.GE.MNTHR ) THEN
*
            IF( WNTQN ) THEN
*
*              Path 1t (N much larger than M, JOBZ='N')
*              No singular vectors to be computed
*
               ITAU = 1
               NWORK = ITAU + M
*
*              Compute A=L*Q
*              (Workspace: need 2*M, prefer M+M*NB)
*
               CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Zero out above L
*
               CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
               IE = 1
               ITAUQ = IE + M
               ITAUP = ITAUQ + M
               NWORK = ITAUP + M
*
*              Bidiagonalize L in A
*              (Workspace: need 4*M, prefer 3*M+2*M*NB)
*
               CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                      IERR )
               NWORK = IE + M
*
*              Perform bidiagonal SVD, computing singular values only
*              (Workspace: need M+BDSPAC)
*
               CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
*
            ELSE IF( WNTQO ) THEN
*
*              Path 2t (N much larger than M, JOBZ='O')
*              M right singular vectors to be overwritten on A and
*              M left singular vectors to be computed in U
*
               IVT = 1
*
*              IVT is M by M
*
               IL = IVT + M*M
               IF( LWORK.GE.M*N+M*M+3*M+BDSPAC ) THEN
*
*                 WORK(IL) is M by N
*
                  LDWRKL = M
                  CHUNK = N
               ELSE
                  LDWRKL = M
                  CHUNK = ( LWORK-M*M ) / M
               END IF
               ITAU = IL + LDWRKL*M
               NWORK = ITAU + M
*
*              Compute A=L*Q
*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
               CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy L to WORK(IL), zeroing about above it
*
               CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
               CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
     $                      WORK( IL+LDWRKL ), LDWRKL )
*
*              Generate Q in A
*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
               CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = ITAU
               ITAUQ = IE + M
               ITAUP = ITAUQ + M
               NWORK = ITAUP + M
*
*              Bidiagonalize L in WORK(IL)
*              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
*
               CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in U, and computing right singular
*              vectors of bidiagonal matrix in WORK(IVT)
*              (Workspace: need M+M*M+BDSPAC)
*
               CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
     $                      WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
     $                      IWORK, INFO )
*
*              Overwrite U by left singular vectors of L and WORK(IVT)
*              by right singular vectors of L
*              (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
     $                      WORK( ITAUP ), WORK( IVT ), M,
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*              Multiply right singular vectors of L in WORK(IVT) by Q
*              in A, storing result in WORK(IL) and copying to A
*              (Workspace: need 2*M*M, prefer M*M+M*N)
*
               DO 30 I = 1, N, CHUNK
                  BLK = MIN( N-I+1, CHUNK )
                  CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
     $                        A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
                  CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
     $                         A( 1, I ), LDA )
   30          CONTINUE
*
            ELSE IF( WNTQS ) THEN
*
*              Path 3t (N much larger than M, JOBZ='S')
*              M right singular vectors to be computed in VT and
*              M left singular vectors to be computed in U
*
               IL = 1
*
*              WORK(IL) is M by M
*
               LDWRKL = M
               ITAU = IL + LDWRKL*M
               NWORK = ITAU + M
*
*              Compute A=L*Q
*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
               CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy L to WORK(IL), zeroing out above it
*
               CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
               CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
     $                      WORK( IL+LDWRKL ), LDWRKL )
*
*              Generate Q in A
*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
               CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = ITAU
               ITAUQ = IE + M
               ITAUP = ITAUQ + M
               NWORK = ITAUP + M
*
*              Bidiagonalize L in WORK(IU), copying result to U
*              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
*
               CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in U and computing right singular
*              vectors of bidiagonal matrix in VT
*              (Workspace: need M+BDSPAC)
*
               CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
*
*              Overwrite U by left singular vectors of L and VT
*              by right singular vectors of L
*              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Multiply right singular vectors of L in WORK(IL) by
*              Q in A, storing result in VT
*              (Workspace: need M*M)
*
               CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
               CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
     $                     A, LDA, ZERO, VT, LDVT )
*
            ELSE IF( WNTQA ) THEN
*
*              Path 4t (N much larger than M, JOBZ='A')
*              N right singular vectors to be computed in VT and
*              M left singular vectors to be computed in U
*
               IVT = 1
*
*              WORK(IVT) is M by M
*
               LDWKVT = M
               ITAU = IVT + LDWKVT*M
               NWORK = ITAU + M
*
*              Compute A=L*Q, copying result to VT
*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
               CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
*              Generate Q in VT
*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
               CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*              Produce L in A, zeroing out other entries
*
               CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
               IE = ITAU
               ITAUQ = IE + M
               ITAUP = ITAUQ + M
               NWORK = ITAUP + M
*
*              Bidiagonalize L in A
*              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
*
               CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                      IERR )
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in U and computing right singular
*              vectors of bidiagonal matrix in WORK(IVT)
*              (Workspace: need M+M*M+BDSPAC)
*
               CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
     $                      WORK( IVT ), LDWKVT, DUM, IDUM,
     $                      WORK( NWORK ), IWORK, INFO )
*
*              Overwrite U by left singular vectors of L and WORK(IVT)
*              by right singular vectors of L
*              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
     $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*              Multiply right singular vectors of L in WORK(IVT) by
*              Q in VT, storing result in A
*              (Workspace: need M*M)
*
               CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
     $                     VT, LDVT, ZERO, A, LDA )
*
*              Copy right singular vectors of A from A to VT
*
               CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
*
            END IF
*
         ELSE
*
*           N .LT. MNTHR
*
*           Path 5t (N greater than M, but not much larger)
*           Reduce to bidiagonal form without LQ decomposition
*
            IE = 1
            ITAUQ = IE + M
            ITAUP = ITAUQ + M
            NWORK = ITAUP + M
*
*           Bidiagonalize A
*           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
*
            CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
     $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                   IERR )
            IF( WNTQN ) THEN
*
*              Perform bidiagonal SVD, only computing singular values
*              (Workspace: need M+BDSPAC)
*
               CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
            ELSE IF( WNTQO ) THEN
               LDWKVT = M
               IVT = NWORK
               IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
*
*                 WORK( IVT ) is M by N
*
                  CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
     $                         LDWKVT )
                  NWORK = IVT + LDWKVT*N
               ELSE
*
*                 WORK( IVT ) is M by M
*
                  NWORK = IVT + LDWKVT*M
                  IL = NWORK
*
*                 WORK(IL) is M by CHUNK
*
                  CHUNK = ( LWORK-M*M-3*M ) / M
               END IF
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in U and computing right singular
*              vectors of bidiagonal matrix in WORK(IVT)
*              (Workspace: need M*M+BDSPAC)
*
               CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
     $                      WORK( IVT ), LDWKVT, DUM, IDUM,
     $                      WORK( NWORK ), IWORK, INFO )
*
*              Overwrite U by left singular vectors of A
*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
               IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
*
*                 Overwrite WORK(IVT) by left singular vectors of A
*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                  CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
     $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*                 Copy right singular vectors of A from WORK(IVT) to A
*
                  CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
               ELSE
*
*                 Generate P**T in A
*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
*
                  CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*                 Multiply Q in A by right singular vectors of
*                 bidiagonal matrix in WORK(IVT), storing result in
*                 WORK(IL) and copying to A
*                 (Workspace: need 2*M*M, prefer M*M+M*N)
*
                  DO 40 I = 1, N, CHUNK
                     BLK = MIN( N-I+1, CHUNK )
                     CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
     $                           LDWKVT, A( 1, I ), LDA, ZERO,
     $                           WORK( IL ), M )
                     CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
     $                            LDA )
   40             CONTINUE
               END IF
            ELSE IF( WNTQS ) THEN
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in U and computing right singular
*              vectors of bidiagonal matrix in VT
*              (Workspace: need M+BDSPAC)
*
               CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
               CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
*
*              Overwrite U by left singular vectors of A and VT
*              by right singular vectors of A
*              (Workspace: need 3*M, prefer 2*M+M*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
            ELSE IF( WNTQA ) THEN
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in U and computing right singular
*              vectors of bidiagonal matrix in VT
*              (Workspace: need M+BDSPAC)
*
               CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
               CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
*
*              Set the right corner of VT to identity matrix
*
               IF( N.GT.M ) THEN
                  CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ),
     $                         LDVT )
               END IF
*
*              Overwrite U by left singular vectors of A and VT
*              by right singular vectors of A
*              (Workspace: need 2*M+N, prefer 2*M+N*NB)
*
               CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
            END IF
*
         END IF
*
      END IF
*
*     Undo scaling if necessary
*
      IF( ISCL.EQ.1 ) THEN
         IF( ANRM.GT.BIGNUM )
     $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
     $                   IERR )
         IF( ANRM.LT.SMLNUM )
     $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
     $                   IERR )
      END IF
*
*     Return optimal workspace in WORK(1)
*
      WORK( 1 ) = MAXWRK
*
      RETURN
*
*     End of DGESDD
*
      END

max_issues_repo_name: WatsonZhouAnda/Cheetah-Software

max_forks_repo_name: WatsonZhouAnda/Cheetah-Software