LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slasdq.f
Go to the documentation of this file.
00001 *> \brief \b SLASDQ
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLASDQ + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasdq.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasdq.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasdq.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
00022 *                          U, LDU, C, LDC, WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          UPLO
00026 *       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       REAL               C( LDC, * ), D( * ), E( * ), U( LDU, * ),
00030 *      $                   VT( LDVT, * ), WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> SLASDQ computes the singular value decomposition (SVD) of a real
00040 *> (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
00041 *> E, accumulating the transformations if desired. Letting B denote
00042 *> the input bidiagonal matrix, the algorithm computes orthogonal
00043 *> matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose
00044 *> of P). The singular values S are overwritten on D.
00045 *>
00046 *> The input matrix U  is changed to U  * Q  if desired.
00047 *> The input matrix VT is changed to P**T * VT if desired.
00048 *> The input matrix C  is changed to Q**T * C  if desired.
00049 *>
00050 *> See "Computing  Small Singular Values of Bidiagonal Matrices With
00051 *> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
00052 *> LAPACK Working Note #3, for a detailed description of the algorithm.
00053 *> \endverbatim
00054 *
00055 *  Arguments:
00056 *  ==========
00057 *
00058 *> \param[in] UPLO
00059 *> \verbatim
00060 *>          UPLO is CHARACTER*1
00061 *>        On entry, UPLO specifies whether the input bidiagonal matrix
00062 *>        is upper or lower bidiagonal, and wether it is square are
00063 *>        not.
00064 *>           UPLO = 'U' or 'u'   B is upper bidiagonal.
00065 *>           UPLO = 'L' or 'l'   B is lower bidiagonal.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] SQRE
00069 *> \verbatim
00070 *>          SQRE is INTEGER
00071 *>        = 0: then the input matrix is N-by-N.
00072 *>        = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
00073 *>             (N+1)-by-N if UPLU = 'L'.
00074 *>
00075 *>        The bidiagonal matrix has
00076 *>        N = NL + NR + 1 rows and
00077 *>        M = N + SQRE >= N columns.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] N
00081 *> \verbatim
00082 *>          N is INTEGER
00083 *>        On entry, N specifies the number of rows and columns
00084 *>        in the matrix. N must be at least 0.
00085 *> \endverbatim
00086 *>
00087 *> \param[in] NCVT
00088 *> \verbatim
00089 *>          NCVT is INTEGER
00090 *>        On entry, NCVT specifies the number of columns of
00091 *>        the matrix VT. NCVT must be at least 0.
00092 *> \endverbatim
00093 *>
00094 *> \param[in] NRU
00095 *> \verbatim
00096 *>          NRU is INTEGER
00097 *>        On entry, NRU specifies the number of rows of
00098 *>        the matrix U. NRU must be at least 0.
00099 *> \endverbatim
00100 *>
00101 *> \param[in] NCC
00102 *> \verbatim
00103 *>          NCC is INTEGER
00104 *>        On entry, NCC specifies the number of columns of
00105 *>        the matrix C. NCC must be at least 0.
00106 *> \endverbatim
00107 *>
00108 *> \param[in,out] D
00109 *> \verbatim
00110 *>          D is REAL array, dimension (N)
00111 *>        On entry, D contains the diagonal entries of the
00112 *>        bidiagonal matrix whose SVD is desired. On normal exit,
00113 *>        D contains the singular values in ascending order.
00114 *> \endverbatim
00115 *>
00116 *> \param[in,out] E
00117 *> \verbatim
00118 *>          E is REAL array.
00119 *>        dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
00120 *>        On entry, the entries of E contain the offdiagonal entries
00121 *>        of the bidiagonal matrix whose SVD is desired. On normal
00122 *>        exit, E will contain 0. If the algorithm does not converge,
00123 *>        D and E will contain the diagonal and superdiagonal entries
00124 *>        of a bidiagonal matrix orthogonally equivalent to the one
00125 *>        given as input.
00126 *> \endverbatim
00127 *>
00128 *> \param[in,out] VT
00129 *> \verbatim
00130 *>          VT is REAL array, dimension (LDVT, NCVT)
00131 *>        On entry, contains a matrix which on exit has been
00132 *>        premultiplied by P**T, dimension N-by-NCVT if SQRE = 0
00133 *>        and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
00134 *> \endverbatim
00135 *>
00136 *> \param[in] LDVT
00137 *> \verbatim
00138 *>          LDVT is INTEGER
00139 *>        On entry, LDVT specifies the leading dimension of VT as
00140 *>        declared in the calling (sub) program. LDVT must be at
00141 *>        least 1. If NCVT is nonzero LDVT must also be at least N.
00142 *> \endverbatim
00143 *>
00144 *> \param[in,out] U
00145 *> \verbatim
00146 *>          U is REAL array, dimension (LDU, N)
00147 *>        On entry, contains a  matrix which on exit has been
00148 *>        postmultiplied by Q, dimension NRU-by-N if SQRE = 0
00149 *>        and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
00150 *> \endverbatim
00151 *>
00152 *> \param[in] LDU
00153 *> \verbatim
00154 *>          LDU is INTEGER
00155 *>        On entry, LDU  specifies the leading dimension of U as
00156 *>        declared in the calling (sub) program. LDU must be at
00157 *>        least max( 1, NRU ) .
00158 *> \endverbatim
00159 *>
00160 *> \param[in,out] C
00161 *> \verbatim
00162 *>          C is REAL array, dimension (LDC, NCC)
00163 *>        On entry, contains an N-by-NCC matrix which on exit
00164 *>        has been premultiplied by Q**T  dimension N-by-NCC if SQRE = 0
00165 *>        and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
00166 *> \endverbatim
00167 *>
00168 *> \param[in] LDC
00169 *> \verbatim
00170 *>          LDC is INTEGER
00171 *>        On entry, LDC  specifies the leading dimension of C as
00172 *>        declared in the calling (sub) program. LDC must be at
00173 *>        least 1. If NCC is nonzero, LDC must also be at least N.
00174 *> \endverbatim
00175 *>
00176 *> \param[out] WORK
00177 *> \verbatim
00178 *>          WORK is REAL array, dimension (4*N)
00179 *>        Workspace. Only referenced if one of NCVT, NRU, or NCC is
00180 *>        nonzero, and if N is at least 2.
00181 *> \endverbatim
00182 *>
00183 *> \param[out] INFO
00184 *> \verbatim
00185 *>          INFO is INTEGER
00186 *>        On exit, a value of 0 indicates a successful exit.
00187 *>        If INFO < 0, argument number -INFO is illegal.
00188 *>        If INFO > 0, the algorithm did not converge, and INFO
00189 *>        specifies how many superdiagonals did not converge.
00190 *> \endverbatim
00191 *
00192 *  Authors:
00193 *  ========
00194 *
00195 *> \author Univ. of Tennessee 
00196 *> \author Univ. of California Berkeley 
00197 *> \author Univ. of Colorado Denver 
00198 *> \author NAG Ltd. 
00199 *
00200 *> \date November 2011
00201 *
00202 *> \ingroup auxOTHERauxiliary
00203 *
00204 *> \par Contributors:
00205 *  ==================
00206 *>
00207 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00208 *>     California at Berkeley, USA
00209 *>
00210 *  =====================================================================
00211       SUBROUTINE SLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
00212      $                   U, LDU, C, LDC, WORK, INFO )
00213 *
00214 *  -- LAPACK auxiliary routine (version 3.4.0) --
00215 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00216 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00217 *     November 2011
00218 *
00219 *     .. Scalar Arguments ..
00220       CHARACTER          UPLO
00221       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
00222 *     ..
00223 *     .. Array Arguments ..
00224       REAL               C( LDC, * ), D( * ), E( * ), U( LDU, * ),
00225      $                   VT( LDVT, * ), WORK( * )
00226 *     ..
00227 *
00228 *  =====================================================================
00229 *
00230 *     .. Parameters ..
00231       REAL               ZERO
00232       PARAMETER          ( ZERO = 0.0E+0 )
00233 *     ..
00234 *     .. Local Scalars ..
00235       LOGICAL            ROTATE
00236       INTEGER            I, ISUB, IUPLO, J, NP1, SQRE1
00237       REAL               CS, R, SMIN, SN
00238 *     ..
00239 *     .. External Subroutines ..
00240       EXTERNAL           SBDSQR, SLARTG, SLASR, SSWAP, XERBLA
00241 *     ..
00242 *     .. External Functions ..
00243       LOGICAL            LSAME
00244       EXTERNAL           LSAME
00245 *     ..
00246 *     .. Intrinsic Functions ..
00247       INTRINSIC          MAX
00248 *     ..
00249 *     .. Executable Statements ..
00250 *
00251 *     Test the input parameters.
00252 *
00253       INFO = 0
00254       IUPLO = 0
00255       IF( LSAME( UPLO, 'U' ) )
00256      $   IUPLO = 1
00257       IF( LSAME( UPLO, 'L' ) )
00258      $   IUPLO = 2
00259       IF( IUPLO.EQ.0 ) THEN
00260          INFO = -1
00261       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
00262          INFO = -2
00263       ELSE IF( N.LT.0 ) THEN
00264          INFO = -3
00265       ELSE IF( NCVT.LT.0 ) THEN
00266          INFO = -4
00267       ELSE IF( NRU.LT.0 ) THEN
00268          INFO = -5
00269       ELSE IF( NCC.LT.0 ) THEN
00270          INFO = -6
00271       ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
00272      $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
00273          INFO = -10
00274       ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
00275          INFO = -12
00276       ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
00277      $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
00278          INFO = -14
00279       END IF
00280       IF( INFO.NE.0 ) THEN
00281          CALL XERBLA( 'SLASDQ', -INFO )
00282          RETURN
00283       END IF
00284       IF( N.EQ.0 )
00285      $   RETURN
00286 *
00287 *     ROTATE is true if any singular vectors desired, false otherwise
00288 *
00289       ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
00290       NP1 = N + 1
00291       SQRE1 = SQRE
00292 *
00293 *     If matrix non-square upper bidiagonal, rotate to be lower
00294 *     bidiagonal.  The rotations are on the right.
00295 *
00296       IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
00297          DO 10 I = 1, N - 1
00298             CALL SLARTG( D( I ), E( I ), CS, SN, R )
00299             D( I ) = R
00300             E( I ) = SN*D( I+1 )
00301             D( I+1 ) = CS*D( I+1 )
00302             IF( ROTATE ) THEN
00303                WORK( I ) = CS
00304                WORK( N+I ) = SN
00305             END IF
00306    10    CONTINUE
00307          CALL SLARTG( D( N ), E( N ), CS, SN, R )
00308          D( N ) = R
00309          E( N ) = ZERO
00310          IF( ROTATE ) THEN
00311             WORK( N ) = CS
00312             WORK( N+N ) = SN
00313          END IF
00314          IUPLO = 2
00315          SQRE1 = 0
00316 *
00317 *        Update singular vectors if desired.
00318 *
00319          IF( NCVT.GT.0 )
00320      $      CALL SLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
00321      $                  WORK( NP1 ), VT, LDVT )
00322       END IF
00323 *
00324 *     If matrix lower bidiagonal, rotate to be upper bidiagonal
00325 *     by applying Givens rotations on the left.
00326 *
00327       IF( IUPLO.EQ.2 ) THEN
00328          DO 20 I = 1, N - 1
00329             CALL SLARTG( D( I ), E( I ), CS, SN, R )
00330             D( I ) = R
00331             E( I ) = SN*D( I+1 )
00332             D( I+1 ) = CS*D( I+1 )
00333             IF( ROTATE ) THEN
00334                WORK( I ) = CS
00335                WORK( N+I ) = SN
00336             END IF
00337    20    CONTINUE
00338 *
00339 *        If matrix (N+1)-by-N lower bidiagonal, one additional
00340 *        rotation is needed.
00341 *
00342          IF( SQRE1.EQ.1 ) THEN
00343             CALL SLARTG( D( N ), E( N ), CS, SN, R )
00344             D( N ) = R
00345             IF( ROTATE ) THEN
00346                WORK( N ) = CS
00347                WORK( N+N ) = SN
00348             END IF
00349          END IF
00350 *
00351 *        Update singular vectors if desired.
00352 *
00353          IF( NRU.GT.0 ) THEN
00354             IF( SQRE1.EQ.0 ) THEN
00355                CALL SLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
00356      $                     WORK( NP1 ), U, LDU )
00357             ELSE
00358                CALL SLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
00359      $                     WORK( NP1 ), U, LDU )
00360             END IF
00361          END IF
00362          IF( NCC.GT.0 ) THEN
00363             IF( SQRE1.EQ.0 ) THEN
00364                CALL SLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
00365      $                     WORK( NP1 ), C, LDC )
00366             ELSE
00367                CALL SLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
00368      $                     WORK( NP1 ), C, LDC )
00369             END IF
00370          END IF
00371       END IF
00372 *
00373 *     Call SBDSQR to compute the SVD of the reduced real
00374 *     N-by-N upper bidiagonal matrix.
00375 *
00376       CALL SBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
00377      $             LDC, WORK, INFO )
00378 *
00379 *     Sort the singular values into ascending order (insertion sort on
00380 *     singular values, but only one transposition per singular vector)
00381 *
00382       DO 40 I = 1, N
00383 *
00384 *        Scan for smallest D(I).
00385 *
00386          ISUB = I
00387          SMIN = D( I )
00388          DO 30 J = I + 1, N
00389             IF( D( J ).LT.SMIN ) THEN
00390                ISUB = J
00391                SMIN = D( J )
00392             END IF
00393    30    CONTINUE
00394          IF( ISUB.NE.I ) THEN
00395 *
00396 *           Swap singular values and vectors.
00397 *
00398             D( ISUB ) = D( I )
00399             D( I ) = SMIN
00400             IF( NCVT.GT.0 )
00401      $         CALL SSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
00402             IF( NRU.GT.0 )
00403      $         CALL SSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
00404             IF( NCC.GT.0 )
00405      $         CALL SSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
00406          END IF
00407    40 CONTINUE
00408 *
00409       RETURN
00410 *
00411 *     End of SLASDQ
00412 *
00413       END
 All Files Functions