![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLASD0 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLASD0 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd0.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd0.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd0.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK, 00022 * WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDU, LDVT, N, SMLSIZ, SQRE 00026 * .. 00027 * .. Array Arguments .. 00028 * INTEGER IWORK( * ) 00029 * DOUBLE PRECISION D( * ), E( * ), U( LDU, * ), VT( LDVT, * ), 00030 * $ WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> Using a divide and conquer approach, DLASD0 computes the singular 00040 *> value decomposition (SVD) of a real upper bidiagonal N-by-M 00041 *> matrix B with diagonal D and offdiagonal E, where M = N + SQRE. 00042 *> The algorithm computes orthogonal matrices U and VT such that 00043 *> B = U * S * VT. The singular values S are overwritten on D. 00044 *> 00045 *> A related subroutine, DLASDA, computes only the singular values, 00046 *> and optionally, the singular vectors in compact form. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] N 00053 *> \verbatim 00054 *> N is INTEGER 00055 *> On entry, the row dimension of the upper bidiagonal matrix. 00056 *> This is also the dimension of the main diagonal array D. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] SQRE 00060 *> \verbatim 00061 *> SQRE is INTEGER 00062 *> Specifies the column dimension of the bidiagonal matrix. 00063 *> = 0: The bidiagonal matrix has column dimension M = N; 00064 *> = 1: The bidiagonal matrix has column dimension M = N+1; 00065 *> \endverbatim 00066 *> 00067 *> \param[in,out] D 00068 *> \verbatim 00069 *> D is DOUBLE PRECISION array, dimension (N) 00070 *> On entry D contains the main diagonal of the bidiagonal 00071 *> matrix. 00072 *> On exit D, if INFO = 0, contains its singular values. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] E 00076 *> \verbatim 00077 *> E is DOUBLE PRECISION array, dimension (M-1) 00078 *> Contains the subdiagonal entries of the bidiagonal matrix. 00079 *> On exit, E has been destroyed. 00080 *> \endverbatim 00081 *> 00082 *> \param[out] U 00083 *> \verbatim 00084 *> U is DOUBLE PRECISION array, dimension at least (LDQ, N) 00085 *> On exit, U contains the left singular vectors. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] LDU 00089 *> \verbatim 00090 *> LDU is INTEGER 00091 *> On entry, leading dimension of U. 00092 *> \endverbatim 00093 *> 00094 *> \param[out] VT 00095 *> \verbatim 00096 *> VT is DOUBLE PRECISION array, dimension at least (LDVT, M) 00097 *> On exit, VT**T contains the right singular vectors. 00098 *> \endverbatim 00099 *> 00100 *> \param[in] LDVT 00101 *> \verbatim 00102 *> LDVT is INTEGER 00103 *> On entry, leading dimension of VT. 00104 *> \endverbatim 00105 *> 00106 *> \param[in] SMLSIZ 00107 *> \verbatim 00108 *> SMLSIZ is INTEGER 00109 *> On entry, maximum size of the subproblems at the 00110 *> bottom of the computation tree. 00111 *> \endverbatim 00112 *> 00113 *> \param[out] IWORK 00114 *> \verbatim 00115 *> IWORK is INTEGER work array. 00116 *> Dimension must be at least (8 * N) 00117 *> \endverbatim 00118 *> 00119 *> \param[out] WORK 00120 *> \verbatim 00121 *> WORK is DOUBLE PRECISION work array. 00122 *> Dimension must be at least (3 * M**2 + 2 * M) 00123 *> \endverbatim 00124 *> 00125 *> \param[out] INFO 00126 *> \verbatim 00127 *> INFO is INTEGER 00128 *> = 0: successful exit. 00129 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00130 *> > 0: if INFO = 1, a singular value did not converge 00131 *> \endverbatim 00132 * 00133 * Authors: 00134 * ======== 00135 * 00136 *> \author Univ. of Tennessee 00137 *> \author Univ. of California Berkeley 00138 *> \author Univ. of Colorado Denver 00139 *> \author NAG Ltd. 00140 * 00141 *> \date November 2011 00142 * 00143 *> \ingroup auxOTHERauxiliary 00144 * 00145 *> \par Contributors: 00146 * ================== 00147 *> 00148 *> Ming Gu and Huan Ren, Computer Science Division, University of 00149 *> California at Berkeley, USA 00150 *> 00151 * ===================================================================== 00152 SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK, 00153 $ WORK, INFO ) 00154 * 00155 * -- LAPACK auxiliary routine (version 3.4.0) -- 00156 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00157 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00158 * November 2011 00159 * 00160 * .. Scalar Arguments .. 00161 INTEGER INFO, LDU, LDVT, N, SMLSIZ, SQRE 00162 * .. 00163 * .. Array Arguments .. 00164 INTEGER IWORK( * ) 00165 DOUBLE PRECISION D( * ), E( * ), U( LDU, * ), VT( LDVT, * ), 00166 $ WORK( * ) 00167 * .. 00168 * 00169 * ===================================================================== 00170 * 00171 * .. Local Scalars .. 00172 INTEGER I, I1, IC, IDXQ, IDXQC, IM1, INODE, ITEMP, IWK, 00173 $ J, LF, LL, LVL, M, NCC, ND, NDB1, NDIML, NDIMR, 00174 $ NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQREI 00175 DOUBLE PRECISION ALPHA, BETA 00176 * .. 00177 * .. External Subroutines .. 00178 EXTERNAL DLASD1, DLASDQ, DLASDT, XERBLA 00179 * .. 00180 * .. Executable Statements .. 00181 * 00182 * Test the input parameters. 00183 * 00184 INFO = 0 00185 * 00186 IF( N.LT.0 ) THEN 00187 INFO = -1 00188 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00189 INFO = -2 00190 END IF 00191 * 00192 M = N + SQRE 00193 * 00194 IF( LDU.LT.N ) THEN 00195 INFO = -6 00196 ELSE IF( LDVT.LT.M ) THEN 00197 INFO = -8 00198 ELSE IF( SMLSIZ.LT.3 ) THEN 00199 INFO = -9 00200 END IF 00201 IF( INFO.NE.0 ) THEN 00202 CALL XERBLA( 'DLASD0', -INFO ) 00203 RETURN 00204 END IF 00205 * 00206 * If the input matrix is too small, call DLASDQ to find the SVD. 00207 * 00208 IF( N.LE.SMLSIZ ) THEN 00209 CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDVT, U, LDU, U, 00210 $ LDU, WORK, INFO ) 00211 RETURN 00212 END IF 00213 * 00214 * Set up the computation tree. 00215 * 00216 INODE = 1 00217 NDIML = INODE + N 00218 NDIMR = NDIML + N 00219 IDXQ = NDIMR + N 00220 IWK = IDXQ + N 00221 CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), 00222 $ IWORK( NDIMR ), SMLSIZ ) 00223 * 00224 * For the nodes on bottom level of the tree, solve 00225 * their subproblems by DLASDQ. 00226 * 00227 NDB1 = ( ND+1 ) / 2 00228 NCC = 0 00229 DO 30 I = NDB1, ND 00230 * 00231 * IC : center row of each node 00232 * NL : number of rows of left subproblem 00233 * NR : number of rows of right subproblem 00234 * NLF: starting row of the left subproblem 00235 * NRF: starting row of the right subproblem 00236 * 00237 I1 = I - 1 00238 IC = IWORK( INODE+I1 ) 00239 NL = IWORK( NDIML+I1 ) 00240 NLP1 = NL + 1 00241 NR = IWORK( NDIMR+I1 ) 00242 NRP1 = NR + 1 00243 NLF = IC - NL 00244 NRF = IC + 1 00245 SQREI = 1 00246 CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), E( NLF ), 00247 $ VT( NLF, NLF ), LDVT, U( NLF, NLF ), LDU, 00248 $ U( NLF, NLF ), LDU, WORK, INFO ) 00249 IF( INFO.NE.0 ) THEN 00250 RETURN 00251 END IF 00252 ITEMP = IDXQ + NLF - 2 00253 DO 10 J = 1, NL 00254 IWORK( ITEMP+J ) = J 00255 10 CONTINUE 00256 IF( I.EQ.ND ) THEN 00257 SQREI = SQRE 00258 ELSE 00259 SQREI = 1 00260 END IF 00261 NRP1 = NR + SQREI 00262 CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), E( NRF ), 00263 $ VT( NRF, NRF ), LDVT, U( NRF, NRF ), LDU, 00264 $ U( NRF, NRF ), LDU, WORK, INFO ) 00265 IF( INFO.NE.0 ) THEN 00266 RETURN 00267 END IF 00268 ITEMP = IDXQ + IC 00269 DO 20 J = 1, NR 00270 IWORK( ITEMP+J-1 ) = J 00271 20 CONTINUE 00272 30 CONTINUE 00273 * 00274 * Now conquer each subproblem bottom-up. 00275 * 00276 DO 50 LVL = NLVL, 1, -1 00277 * 00278 * Find the first node LF and last node LL on the 00279 * current level LVL. 00280 * 00281 IF( LVL.EQ.1 ) THEN 00282 LF = 1 00283 LL = 1 00284 ELSE 00285 LF = 2**( LVL-1 ) 00286 LL = 2*LF - 1 00287 END IF 00288 DO 40 I = LF, LL 00289 IM1 = I - 1 00290 IC = IWORK( INODE+IM1 ) 00291 NL = IWORK( NDIML+IM1 ) 00292 NR = IWORK( NDIMR+IM1 ) 00293 NLF = IC - NL 00294 IF( ( SQRE.EQ.0 ) .AND. ( I.EQ.LL ) ) THEN 00295 SQREI = SQRE 00296 ELSE 00297 SQREI = 1 00298 END IF 00299 IDXQC = IDXQ + NLF - 1 00300 ALPHA = D( IC ) 00301 BETA = E( IC ) 00302 CALL DLASD1( NL, NR, SQREI, D( NLF ), ALPHA, BETA, 00303 $ U( NLF, NLF ), LDU, VT( NLF, NLF ), LDVT, 00304 $ IWORK( IDXQC ), IWORK( IWK ), WORK, INFO ) 00305 IF( INFO.NE.0 ) THEN 00306 RETURN 00307 END IF 00308 40 CONTINUE 00309 50 CONTINUE 00310 * 00311 RETURN 00312 * 00313 * End of DLASD0 00314 * 00315 END