 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK, 00002 $ WORK, INFO ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2.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 * June 2010 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDU, LDVT, N, SMLSIZ, SQRE 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IWORK( * ) 00014 DOUBLE PRECISION D( * ), E( * ), U( LDU, * ), VT( LDVT, * ), 00015 $ WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * Using a divide and conquer approach, DLASD0 computes the singular 00022 * value decomposition (SVD) of a real upper bidiagonal N-by-M 00023 * matrix B with diagonal D and offdiagonal E, where M = N + SQRE. 00024 * The algorithm computes orthogonal matrices U and VT such that 00025 * B = U * S * VT. The singular values S are overwritten on D. 00026 * 00027 * A related subroutine, DLASDA, computes only the singular values, 00028 * and optionally, the singular vectors in compact form. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * N (input) INTEGER 00034 * On entry, the row dimension of the upper bidiagonal matrix. 00035 * This is also the dimension of the main diagonal array D. 00036 * 00037 * SQRE (input) INTEGER 00038 * Specifies the column dimension of the bidiagonal matrix. 00039 * = 0: The bidiagonal matrix has column dimension M = N; 00040 * = 1: The bidiagonal matrix has column dimension M = N+1; 00041 * 00042 * D (input/output) DOUBLE PRECISION array, dimension (N) 00043 * On entry D contains the main diagonal of the bidiagonal 00044 * matrix. 00045 * On exit D, if INFO = 0, contains its singular values. 00046 * 00047 * E (input) DOUBLE PRECISION array, dimension (M-1) 00048 * Contains the subdiagonal entries of the bidiagonal matrix. 00049 * On exit, E has been destroyed. 00050 * 00051 * U (output) DOUBLE PRECISION array, dimension at least (LDQ, N) 00052 * On exit, U contains the left singular vectors. 00053 * 00054 * LDU (input) INTEGER 00055 * On entry, leading dimension of U. 00056 * 00057 * VT (output) DOUBLE PRECISION array, dimension at least (LDVT, M) 00058 * On exit, VT**T contains the right singular vectors. 00059 * 00060 * LDVT (input) INTEGER 00061 * On entry, leading dimension of VT. 00062 * 00063 * SMLSIZ (input) INTEGER 00064 * On entry, maximum size of the subproblems at the 00065 * bottom of the computation tree. 00066 * 00067 * IWORK (workspace) INTEGER work array. 00068 * Dimension must be at least (8 * N) 00069 * 00070 * WORK (workspace) DOUBLE PRECISION work array. 00071 * Dimension must be at least (3 * M**2 + 2 * M) 00072 * 00073 * INFO (output) INTEGER 00074 * = 0: successful exit. 00075 * < 0: if INFO = -i, the i-th argument had an illegal value. 00076 * > 0: if INFO = 1, a singular value did not converge 00077 * 00078 * Further Details 00079 * =============== 00080 * 00081 * Based on contributions by 00082 * Ming Gu and Huan Ren, Computer Science Division, University of 00083 * California at Berkeley, USA 00084 * 00085 * ===================================================================== 00086 * 00087 * .. Local Scalars .. 00088 INTEGER I, I1, IC, IDXQ, IDXQC, IM1, INODE, ITEMP, IWK, 00089 $ J, LF, LL, LVL, M, NCC, ND, NDB1, NDIML, NDIMR, 00090 $ NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQREI 00091 DOUBLE PRECISION ALPHA, BETA 00092 * .. 00093 * .. External Subroutines .. 00094 EXTERNAL DLASD1, DLASDQ, DLASDT, XERBLA 00095 * .. 00096 * .. Executable Statements .. 00097 * 00098 * Test the input parameters. 00099 * 00100 INFO = 0 00101 * 00102 IF( N.LT.0 ) THEN 00103 INFO = -1 00104 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00105 INFO = -2 00106 END IF 00107 * 00108 M = N + SQRE 00109 * 00110 IF( LDU.LT.N ) THEN 00111 INFO = -6 00112 ELSE IF( LDVT.LT.M ) THEN 00113 INFO = -8 00114 ELSE IF( SMLSIZ.LT.3 ) THEN 00115 INFO = -9 00116 END IF 00117 IF( INFO.NE.0 ) THEN 00118 CALL XERBLA( 'DLASD0', -INFO ) 00119 RETURN 00120 END IF 00121 * 00122 * If the input matrix is too small, call DLASDQ to find the SVD. 00123 * 00124 IF( N.LE.SMLSIZ ) THEN 00125 CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDVT, U, LDU, U, 00126 $ LDU, WORK, INFO ) 00127 RETURN 00128 END IF 00129 * 00130 * Set up the computation tree. 00131 * 00132 INODE = 1 00133 NDIML = INODE + N 00134 NDIMR = NDIML + N 00135 IDXQ = NDIMR + N 00136 IWK = IDXQ + N 00137 CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), 00138 $ IWORK( NDIMR ), SMLSIZ ) 00139 * 00140 * For the nodes on bottom level of the tree, solve 00141 * their subproblems by DLASDQ. 00142 * 00143 NDB1 = ( ND+1 ) / 2 00144 NCC = 0 00145 DO 30 I = NDB1, ND 00146 * 00147 * IC : center row of each node 00148 * NL : number of rows of left subproblem 00149 * NR : number of rows of right subproblem 00150 * NLF: starting row of the left subproblem 00151 * NRF: starting row of the right subproblem 00152 * 00153 I1 = I - 1 00154 IC = IWORK( INODE+I1 ) 00155 NL = IWORK( NDIML+I1 ) 00156 NLP1 = NL + 1 00157 NR = IWORK( NDIMR+I1 ) 00158 NRP1 = NR + 1 00159 NLF = IC - NL 00160 NRF = IC + 1 00161 SQREI = 1 00162 CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), E( NLF ), 00163 $ VT( NLF, NLF ), LDVT, U( NLF, NLF ), LDU, 00164 $ U( NLF, NLF ), LDU, WORK, INFO ) 00165 IF( INFO.NE.0 ) THEN 00166 RETURN 00167 END IF 00168 ITEMP = IDXQ + NLF - 2 00169 DO 10 J = 1, NL 00170 IWORK( ITEMP+J ) = J 00171 10 CONTINUE 00172 IF( I.EQ.ND ) THEN 00173 SQREI = SQRE 00174 ELSE 00175 SQREI = 1 00176 END IF 00177 NRP1 = NR + SQREI 00178 CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), E( NRF ), 00179 $ VT( NRF, NRF ), LDVT, U( NRF, NRF ), LDU, 00180 $ U( NRF, NRF ), LDU, WORK, INFO ) 00181 IF( INFO.NE.0 ) THEN 00182 RETURN 00183 END IF 00184 ITEMP = IDXQ + IC 00185 DO 20 J = 1, NR 00186 IWORK( ITEMP+J-1 ) = J 00187 20 CONTINUE 00188 30 CONTINUE 00189 * 00190 * Now conquer each subproblem bottom-up. 00191 * 00192 DO 50 LVL = NLVL, 1, -1 00193 * 00194 * Find the first node LF and last node LL on the 00195 * current level LVL. 00196 * 00197 IF( LVL.EQ.1 ) THEN 00198 LF = 1 00199 LL = 1 00200 ELSE 00201 LF = 2**( LVL-1 ) 00202 LL = 2*LF - 1 00203 END IF 00204 DO 40 I = LF, LL 00205 IM1 = I - 1 00206 IC = IWORK( INODE+IM1 ) 00207 NL = IWORK( NDIML+IM1 ) 00208 NR = IWORK( NDIMR+IM1 ) 00209 NLF = IC - NL 00210 IF( ( SQRE.EQ.0 ) .AND. ( I.EQ.LL ) ) THEN 00211 SQREI = SQRE 00212 ELSE 00213 SQREI = 1 00214 END IF 00215 IDXQC = IDXQ + NLF - 1 00216 ALPHA = D( IC ) 00217 BETA = E( IC ) 00218 CALL DLASD1( NL, NR, SQREI, D( NLF ), ALPHA, BETA, 00219 $ U( NLF, NLF ), LDU, VT( NLF, NLF ), LDVT, 00220 $ IWORK( IDXQC ), IWORK( IWK ), WORK, INFO ) 00221 IF( INFO.NE.0 ) THEN 00222 RETURN 00223 END IF 00224 40 CONTINUE 00225 50 CONTINUE 00226 * 00227 RETURN 00228 * 00229 * End of DLASD0 00230 * 00231 END
1.7.3 
