LAPACK 3.3.1
Linear Algebra PACKage

dlasdq.f

Go to the documentation of this file.
00001       SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
00002      $                   U, LDU, C, LDC, WORK, INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          UPLO
00011       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
00015      $                   VT( LDVT, * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DLASDQ computes the singular value decomposition (SVD) of a real
00022 *  (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
00023 *  E, accumulating the transformations if desired. Letting B denote
00024 *  the input bidiagonal matrix, the algorithm computes orthogonal
00025 *  matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose
00026 *  of P). The singular values S are overwritten on D.
00027 *
00028 *  The input matrix U  is changed to U  * Q  if desired.
00029 *  The input matrix VT is changed to P**T * VT if desired.
00030 *  The input matrix C  is changed to Q**T * C  if desired.
00031 *
00032 *  See "Computing  Small Singular Values of Bidiagonal Matrices With
00033 *  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
00034 *  LAPACK Working Note #3, for a detailed description of the algorithm.
00035 *
00036 *  Arguments
00037 *  =========
00038 *
00039 *  UPLO  (input) CHARACTER*1
00040 *        On entry, UPLO specifies whether the input bidiagonal matrix
00041 *        is upper or lower bidiagonal, and wether it is square are
00042 *        not.
00043 *           UPLO = 'U' or 'u'   B is upper bidiagonal.
00044 *           UPLO = 'L' or 'l'   B is lower bidiagonal.
00045 *
00046 *  SQRE  (input) INTEGER
00047 *        = 0: then the input matrix is N-by-N.
00048 *        = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
00049 *             (N+1)-by-N if UPLU = 'L'.
00050 *
00051 *        The bidiagonal matrix has
00052 *        N = NL + NR + 1 rows and
00053 *        M = N + SQRE >= N columns.
00054 *
00055 *  N     (input) INTEGER
00056 *        On entry, N specifies the number of rows and columns
00057 *        in the matrix. N must be at least 0.
00058 *
00059 *  NCVT  (input) INTEGER
00060 *        On entry, NCVT specifies the number of columns of
00061 *        the matrix VT. NCVT must be at least 0.
00062 *
00063 *  NRU   (input) INTEGER
00064 *        On entry, NRU specifies the number of rows of
00065 *        the matrix U. NRU must be at least 0.
00066 *
00067 *  NCC   (input) INTEGER
00068 *        On entry, NCC specifies the number of columns of
00069 *        the matrix C. NCC must be at least 0.
00070 *
00071 *  D     (input/output) DOUBLE PRECISION array, dimension (N)
00072 *        On entry, D contains the diagonal entries of the
00073 *        bidiagonal matrix whose SVD is desired. On normal exit,
00074 *        D contains the singular values in ascending order.
00075 *
00076 *  E     (input/output) DOUBLE PRECISION array.
00077 *        dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
00078 *        On entry, the entries of E contain the offdiagonal entries
00079 *        of the bidiagonal matrix whose SVD is desired. On normal
00080 *        exit, E will contain 0. If the algorithm does not converge,
00081 *        D and E will contain the diagonal and superdiagonal entries
00082 *        of a bidiagonal matrix orthogonally equivalent to the one
00083 *        given as input.
00084 *
00085 *  VT    (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
00086 *        On entry, contains a matrix which on exit has been
00087 *        premultiplied by P**T, dimension N-by-NCVT if SQRE = 0
00088 *        and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
00089 *
00090 *  LDVT  (input) INTEGER
00091 *        On entry, LDVT specifies the leading dimension of VT as
00092 *        declared in the calling (sub) program. LDVT must be at
00093 *        least 1. If NCVT is nonzero LDVT must also be at least N.
00094 *
00095 *  U     (input/output) DOUBLE PRECISION array, dimension (LDU, N)
00096 *        On entry, contains a  matrix which on exit has been
00097 *        postmultiplied by Q, dimension NRU-by-N if SQRE = 0
00098 *        and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
00099 *
00100 *  LDU   (input) INTEGER
00101 *        On entry, LDU  specifies the leading dimension of U as
00102 *        declared in the calling (sub) program. LDU must be at
00103 *        least max( 1, NRU ) .
00104 *
00105 *  C     (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
00106 *        On entry, contains an N-by-NCC matrix which on exit
00107 *        has been premultiplied by Q**T  dimension N-by-NCC if SQRE = 0
00108 *        and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
00109 *
00110 *  LDC   (input) INTEGER
00111 *        On entry, LDC  specifies the leading dimension of C as
00112 *        declared in the calling (sub) program. LDC must be at
00113 *        least 1. If NCC is nonzero, LDC must also be at least N.
00114 *
00115 *  WORK  (workspace) DOUBLE PRECISION array, dimension (4*N)
00116 *        Workspace. Only referenced if one of NCVT, NRU, or NCC is
00117 *        nonzero, and if N is at least 2.
00118 *
00119 *  INFO  (output) INTEGER
00120 *        On exit, a value of 0 indicates a successful exit.
00121 *        If INFO < 0, argument number -INFO is illegal.
00122 *        If INFO > 0, the algorithm did not converge, and INFO
00123 *        specifies how many superdiagonals did not converge.
00124 *
00125 *  Further Details
00126 *  ===============
00127 *
00128 *  Based on contributions by
00129 *     Ming Gu and Huan Ren, Computer Science Division, University of
00130 *     California at Berkeley, USA
00131 *
00132 *  =====================================================================
00133 *
00134 *     .. Parameters ..
00135       DOUBLE PRECISION   ZERO
00136       PARAMETER          ( ZERO = 0.0D+0 )
00137 *     ..
00138 *     .. Local Scalars ..
00139       LOGICAL            ROTATE
00140       INTEGER            I, ISUB, IUPLO, J, NP1, SQRE1
00141       DOUBLE PRECISION   CS, R, SMIN, SN
00142 *     ..
00143 *     .. External Subroutines ..
00144       EXTERNAL           DBDSQR, DLARTG, DLASR, DSWAP, XERBLA
00145 *     ..
00146 *     .. External Functions ..
00147       LOGICAL            LSAME
00148       EXTERNAL           LSAME
00149 *     ..
00150 *     .. Intrinsic Functions ..
00151       INTRINSIC          MAX
00152 *     ..
00153 *     .. Executable Statements ..
00154 *
00155 *     Test the input parameters.
00156 *
00157       INFO = 0
00158       IUPLO = 0
00159       IF( LSAME( UPLO, 'U' ) )
00160      $   IUPLO = 1
00161       IF( LSAME( UPLO, 'L' ) )
00162      $   IUPLO = 2
00163       IF( IUPLO.EQ.0 ) THEN
00164          INFO = -1
00165       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
00166          INFO = -2
00167       ELSE IF( N.LT.0 ) THEN
00168          INFO = -3
00169       ELSE IF( NCVT.LT.0 ) THEN
00170          INFO = -4
00171       ELSE IF( NRU.LT.0 ) THEN
00172          INFO = -5
00173       ELSE IF( NCC.LT.0 ) THEN
00174          INFO = -6
00175       ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
00176      $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
00177          INFO = -10
00178       ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
00179          INFO = -12
00180       ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
00181      $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
00182          INFO = -14
00183       END IF
00184       IF( INFO.NE.0 ) THEN
00185          CALL XERBLA( 'DLASDQ', -INFO )
00186          RETURN
00187       END IF
00188       IF( N.EQ.0 )
00189      $   RETURN
00190 *
00191 *     ROTATE is true if any singular vectors desired, false otherwise
00192 *
00193       ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
00194       NP1 = N + 1
00195       SQRE1 = SQRE
00196 *
00197 *     If matrix non-square upper bidiagonal, rotate to be lower
00198 *     bidiagonal.  The rotations are on the right.
00199 *
00200       IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
00201          DO 10 I = 1, N - 1
00202             CALL DLARTG( D( I ), E( I ), CS, SN, R )
00203             D( I ) = R
00204             E( I ) = SN*D( I+1 )
00205             D( I+1 ) = CS*D( I+1 )
00206             IF( ROTATE ) THEN
00207                WORK( I ) = CS
00208                WORK( N+I ) = SN
00209             END IF
00210    10    CONTINUE
00211          CALL DLARTG( D( N ), E( N ), CS, SN, R )
00212          D( N ) = R
00213          E( N ) = ZERO
00214          IF( ROTATE ) THEN
00215             WORK( N ) = CS
00216             WORK( N+N ) = SN
00217          END IF
00218          IUPLO = 2
00219          SQRE1 = 0
00220 *
00221 *        Update singular vectors if desired.
00222 *
00223          IF( NCVT.GT.0 )
00224      $      CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
00225      $                  WORK( NP1 ), VT, LDVT )
00226       END IF
00227 *
00228 *     If matrix lower bidiagonal, rotate to be upper bidiagonal
00229 *     by applying Givens rotations on the left.
00230 *
00231       IF( IUPLO.EQ.2 ) THEN
00232          DO 20 I = 1, N - 1
00233             CALL DLARTG( D( I ), E( I ), CS, SN, R )
00234             D( I ) = R
00235             E( I ) = SN*D( I+1 )
00236             D( I+1 ) = CS*D( I+1 )
00237             IF( ROTATE ) THEN
00238                WORK( I ) = CS
00239                WORK( N+I ) = SN
00240             END IF
00241    20    CONTINUE
00242 *
00243 *        If matrix (N+1)-by-N lower bidiagonal, one additional
00244 *        rotation is needed.
00245 *
00246          IF( SQRE1.EQ.1 ) THEN
00247             CALL DLARTG( D( N ), E( N ), CS, SN, R )
00248             D( N ) = R
00249             IF( ROTATE ) THEN
00250                WORK( N ) = CS
00251                WORK( N+N ) = SN
00252             END IF
00253          END IF
00254 *
00255 *        Update singular vectors if desired.
00256 *
00257          IF( NRU.GT.0 ) THEN
00258             IF( SQRE1.EQ.0 ) THEN
00259                CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
00260      $                     WORK( NP1 ), U, LDU )
00261             ELSE
00262                CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
00263      $                     WORK( NP1 ), U, LDU )
00264             END IF
00265          END IF
00266          IF( NCC.GT.0 ) THEN
00267             IF( SQRE1.EQ.0 ) THEN
00268                CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
00269      $                     WORK( NP1 ), C, LDC )
00270             ELSE
00271                CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
00272      $                     WORK( NP1 ), C, LDC )
00273             END IF
00274          END IF
00275       END IF
00276 *
00277 *     Call DBDSQR to compute the SVD of the reduced real
00278 *     N-by-N upper bidiagonal matrix.
00279 *
00280       CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
00281      $             LDC, WORK, INFO )
00282 *
00283 *     Sort the singular values into ascending order (insertion sort on
00284 *     singular values, but only one transposition per singular vector)
00285 *
00286       DO 40 I = 1, N
00287 *
00288 *        Scan for smallest D(I).
00289 *
00290          ISUB = I
00291          SMIN = D( I )
00292          DO 30 J = I + 1, N
00293             IF( D( J ).LT.SMIN ) THEN
00294                ISUB = J
00295                SMIN = D( J )
00296             END IF
00297    30    CONTINUE
00298          IF( ISUB.NE.I ) THEN
00299 *
00300 *           Swap singular values and vectors.
00301 *
00302             D( ISUB ) = D( I )
00303             D( I ) = SMIN
00304             IF( NCVT.GT.0 )
00305      $         CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
00306             IF( NRU.GT.0 )
00307      $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
00308             IF( NCC.GT.0 )
00309      $         CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
00310          END IF
00311    40 CONTINUE
00312 *
00313       RETURN
00314 *
00315 *     End of DLASDQ
00316 *
00317       END
 All Files Functions