Chimie Théorique » Opt'n Path

      SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,

     $                   U, LDU, C, LDC, WORK, INFO )

*

*  -- LAPACK auxiliary routine (version 3.2) --

*  -- LAPACK is a software package provided by Univ. of Tennessee,    --

*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

*     November 2006

*

*     .. Scalar Arguments ..

      CHARACTER          UPLO

      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE

*     ..

*     .. Array Arguments ..

      DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),

     $                   VT( LDVT, * ), WORK( * )

*     ..

*

*  Purpose

*  =======

*

*  DLASDQ computes the singular value decomposition (SVD) of a real

*  (upper or lower) bidiagonal matrix with diagonal D and offdiagonal

*  E, accumulating the transformations if desired. Letting B denote

*  the input bidiagonal matrix, the algorithm computes orthogonal

*  matrices Q and P such that B = Q * S * P' (P' denotes the transpose

*  of P). The singular values S are overwritten on D.

*

*  The input matrix U  is changed to U  * Q  if desired.

*  The input matrix VT is changed to P' * VT if desired.

*  The input matrix C  is changed to Q' * C  if desired.

*

*  See "Computing  Small Singular Values of Bidiagonal Matrices With

*  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,

*  LAPACK Working Note #3, for a detailed description of the algorithm.

*

*  Arguments

*  =========

*

*  UPLO  (input) CHARACTER*1

*        On entry, UPLO specifies whether the input bidiagonal matrix

*        is upper or lower bidiagonal, and wether it is square are

*        not.

*           UPLO = 'U' or 'u'   B is upper bidiagonal.

*           UPLO = 'L' or 'l'   B is lower bidiagonal.

*

*  SQRE  (input) INTEGER

*        = 0: then the input matrix is N-by-N.

*        = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and

*             (N+1)-by-N if UPLU = 'L'.

*

*        The bidiagonal matrix has

*        N = NL + NR + 1 rows and

*        M = N + SQRE >= N columns.

*

*  N     (input) INTEGER

*        On entry, N specifies the number of rows and columns

*        in the matrix. N must be at least 0.

*

*  NCVT  (input) INTEGER

*        On entry, NCVT specifies the number of columns of

*        the matrix VT. NCVT must be at least 0.

*

*  NRU   (input) INTEGER

*        On entry, NRU specifies the number of rows of

*        the matrix U. NRU must be at least 0.

*

*  NCC   (input) INTEGER

*        On entry, NCC specifies the number of columns of

*        the matrix C. NCC must be at least 0.

*

*  D     (input/output) DOUBLE PRECISION array, dimension (N)

*        On entry, D contains the diagonal entries of the

*        bidiagonal matrix whose SVD is desired. On normal exit,

*        D contains the singular values in ascending order.

*

*  E     (input/output) DOUBLE PRECISION array.

*        dimension is (N-1) if SQRE = 0 and N if SQRE = 1.

*        On entry, the entries of E contain the offdiagonal entries

*        of the bidiagonal matrix whose SVD is desired. On normal

*        exit, E will contain 0. If the algorithm does not converge,

*        D and E will contain the diagonal and superdiagonal entries

*        of a bidiagonal matrix orthogonally equivalent to the one

*        given as input.

*

*  VT    (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)

*        On entry, contains a matrix which on exit has been

*        premultiplied by P', dimension N-by-NCVT if SQRE = 0

*        and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).

*

*  LDVT  (input) INTEGER

*        On entry, LDVT specifies the leading dimension of VT as

*        declared in the calling (sub) program. LDVT must be at

*        least 1. If NCVT is nonzero LDVT must also be at least N.

*

*  U     (input/output) DOUBLE PRECISION array, dimension (LDU, N)

*        On entry, contains a  matrix which on exit has been

*        postmultiplied by Q, dimension NRU-by-N if SQRE = 0

*        and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).

*

*  LDU   (input) INTEGER

*        On entry, LDU  specifies the leading dimension of U as

*        declared in the calling (sub) program. LDU must be at

*        least max( 1, NRU ) .

*

*  C     (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)

*        On entry, contains an N-by-NCC matrix which on exit

*        has been premultiplied by Q'  dimension N-by-NCC if SQRE = 0

*        and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).

*

*  LDC   (input) INTEGER

*        On entry, LDC  specifies the leading dimension of C as

*        declared in the calling (sub) program. LDC must be at

*        least 1. If NCC is nonzero, LDC must also be at least N.

*

*  WORK  (workspace) DOUBLE PRECISION array, dimension (4*N)

*        Workspace. Only referenced if one of NCVT, NRU, or NCC is

*        nonzero, and if N is at least 2.

*

*  INFO  (output) INTEGER

*        On exit, a value of 0 indicates a successful exit.

*        If INFO < 0, argument number -INFO is illegal.

*        If INFO > 0, the algorithm did not converge, and INFO

*        specifies how many superdiagonals did not converge.

*

*  Further Details

*  ===============

*

*  Based on contributions by

*     Ming Gu and Huan Ren, Computer Science Division, University of

*     California at Berkeley, USA

*

*  =====================================================================

*

*     .. Parameters ..

      DOUBLE PRECISION   ZERO

      PARAMETER          ( ZERO = 0.0D+0 )

*     ..

*     .. Local Scalars ..

      LOGICAL            ROTATE

      INTEGER            I, ISUB, IUPLO, J, NP1, SQRE1

      DOUBLE PRECISION   CS, R, SMIN, SN

*     ..

*     .. External Subroutines ..

      EXTERNAL           DBDSQR, DLARTG, DLASR, DSWAP, XERBLA

*     ..

*     .. External Functions ..

      LOGICAL            LSAME

      EXTERNAL           LSAME

*     ..

*     .. Intrinsic Functions ..

      INTRINSIC          MAX

*     ..

*     .. Executable Statements ..

*

*     Test the input parameters.

*

      INFO = 0

      IUPLO = 0

      IF( LSAME( UPLO, 'U' ) )

     $   IUPLO = 1

      IF( LSAME( UPLO, 'L' ) )

     $   IUPLO = 2

      IF( IUPLO.EQ.0 ) THEN

         INFO = -1

      ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN

         INFO = -2

      ELSE IF( N.LT.0 ) THEN

         INFO = -3

      ELSE IF( NCVT.LT.0 ) THEN

         INFO = -4

      ELSE IF( NRU.LT.0 ) THEN

         INFO = -5

      ELSE IF( NCC.LT.0 ) THEN

         INFO = -6

      ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.

     $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN

         INFO = -10

      ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN

         INFO = -12

      ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.

     $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN

         INFO = -14

      END IF

      IF( INFO.NE.0 ) THEN

         CALL XERBLA( 'DLASDQ', -INFO )

         RETURN

      END IF

      IF( N.EQ.0 )

     $   RETURN

*

*     ROTATE is true if any singular vectors desired, false otherwise

*

      ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )

      NP1 = N + 1

      SQRE1 = SQRE

*

*     If matrix non-square upper bidiagonal, rotate to be lower

*     bidiagonal.  The rotations are on the right.

*

      IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN

         DO 10 I = 1, N - 1

            CALL DLARTG( D( I ), E( I ), CS, SN, R )

            D( I ) = R

            E( I ) = SN*D( I+1 )

            D( I+1 ) = CS*D( I+1 )

            IF( ROTATE ) THEN

               WORK( I ) = CS

               WORK( N+I ) = SN

            END IF

   10    CONTINUE

         CALL DLARTG( D( N ), E( N ), CS, SN, R )

         D( N ) = R

         E( N ) = ZERO

         IF( ROTATE ) THEN

            WORK( N ) = CS

            WORK( N+N ) = SN

         END IF

         IUPLO = 2

         SQRE1 = 0

*

*        Update singular vectors if desired.

*

         IF( NCVT.GT.0 )

     $      CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),

     $                  WORK( NP1 ), VT, LDVT )

      END IF

*

*     If matrix lower bidiagonal, rotate to be upper bidiagonal

*     by applying Givens rotations on the left.

*

      IF( IUPLO.EQ.2 ) THEN

         DO 20 I = 1, N - 1

            CALL DLARTG( D( I ), E( I ), CS, SN, R )

            D( I ) = R

            E( I ) = SN*D( I+1 )

            D( I+1 ) = CS*D( I+1 )

            IF( ROTATE ) THEN

               WORK( I ) = CS

               WORK( N+I ) = SN

            END IF

   20    CONTINUE

*

*        If matrix (N+1)-by-N lower bidiagonal, one additional

*        rotation is needed.

*

         IF( SQRE1.EQ.1 ) THEN

            CALL DLARTG( D( N ), E( N ), CS, SN, R )

            D( N ) = R

            IF( ROTATE ) THEN

               WORK( N ) = CS

               WORK( N+N ) = SN

            END IF

         END IF

*

*        Update singular vectors if desired.

*

         IF( NRU.GT.0 ) THEN

            IF( SQRE1.EQ.0 ) THEN

               CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),

     $                     WORK( NP1 ), U, LDU )

            ELSE

               CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),

     $                     WORK( NP1 ), U, LDU )

            END IF

         END IF

         IF( NCC.GT.0 ) THEN

            IF( SQRE1.EQ.0 ) THEN

               CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),

     $                     WORK( NP1 ), C, LDC )

            ELSE

               CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),

     $                     WORK( NP1 ), C, LDC )

            END IF

         END IF

      END IF

*

*     Call DBDSQR to compute the SVD of the reduced real

*     N-by-N upper bidiagonal matrix.

*

      CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,

     $             LDC, WORK, INFO )

*

*     Sort the singular values into ascending order (insertion sort on

*     singular values, but only one transposition per singular vector)

*

      DO 40 I = 1, N

*

*        Scan for smallest D(I).

*

         ISUB = I

         SMIN = D( I )

         DO 30 J = I + 1, N

            IF( D( J ).LT.SMIN ) THEN

               ISUB = J

               SMIN = D( J )

            END IF

   30    CONTINUE

         IF( ISUB.NE.I ) THEN

*

*           Swap singular values and vectors.

*

            D( ISUB ) = D( I )

            D( I ) = SMIN

            IF( NCVT.GT.0 )

     $         CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )

            IF( NRU.GT.0 )

     $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )

            IF( NCC.GT.0 )

     $         CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )

         END IF

   40 CONTINUE

*

      RETURN

*

*     End of DLASDQ

*

      END