|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 |
SUBROUTINE SLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK,
$ WORK, INFO ) * * -- LAPACK auxiliary routine (version 3.2.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * June 2010 * * .. Scalar Arguments .. INTEGER INFO, LDU, LDVT, N, SMLSIZ, SQRE * .. * .. Array Arguments .. INTEGER IWORK( * ) REAL D( * ), E( * ), U( LDU, * ), VT( LDVT, * ), $ WORK( * ) * .. * * Purpose * ======= * * Using a divide and conquer approach, SLASD0 computes the singular * value decomposition (SVD) of a real upper bidiagonal N-by-M * matrix B with diagonal D and offdiagonal E, where M = N + SQRE. * The algorithm computes orthogonal matrices U and VT such that * B = U * S * VT. The singular values S are overwritten on D. * * A related subroutine, SLASDA, computes only the singular values, * and optionally, the singular vectors in compact form. * * Arguments * ========= * * N (input) INTEGER * On entry, the row dimension of the upper bidiagonal matrix. * This is also the dimension of the main diagonal array D. * * SQRE (input) INTEGER * Specifies the column dimension of the bidiagonal matrix. * = 0: The bidiagonal matrix has column dimension M = N; * = 1: The bidiagonal matrix has column dimension M = N+1; * * D (input/output) REAL array, dimension (N) * On entry D contains the main diagonal of the bidiagonal * matrix. * On exit D, if INFO = 0, contains its singular values. * * E (input) REAL array, dimension (M-1) * Contains the subdiagonal entries of the bidiagonal matrix. * On exit, E has been destroyed. * * U (output) REAL array, dimension at least (LDQ, N) * On exit, U contains the left singular vectors. * * LDU (input) INTEGER * On entry, leading dimension of U. * * VT (output) REAL array, dimension at least (LDVT, M) * On exit, VT**T contains the right singular vectors. * * LDVT (input) INTEGER * On entry, leading dimension of VT. * * SMLSIZ (input) INTEGER * On entry, maximum size of the subproblems at the * bottom of the computation tree. * * IWORK (workspace) INTEGER array, dimension (8*N) * * WORK (workspace) REAL array, dimension (3*M**2+2*M) * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: if INFO = 1, a singular value did not converge * * Further Details * =============== * * Based on contributions by * Ming Gu and Huan Ren, Computer Science Division, University of * California at Berkeley, USA * * ===================================================================== * * .. Local Scalars .. INTEGER I, I1, IC, IDXQ, IDXQC, IM1, INODE, ITEMP, IWK, $ J, LF, LL, LVL, M, NCC, ND, NDB1, NDIML, NDIMR, $ NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQREI REAL ALPHA, BETA * .. * .. External Subroutines .. EXTERNAL SLASD1, SLASDQ, SLASDT, XERBLA * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 * IF( N.LT.0 ) THEN INFO = -1 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN INFO = -2 END IF * M = N + SQRE * IF( LDU.LT.N ) THEN INFO = -6 ELSE IF( LDVT.LT.M ) THEN INFO = -8 ELSE IF( SMLSIZ.LT.3 ) THEN INFO = -9 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'SLASD0', -INFO ) RETURN END IF * * If the input matrix is too small, call SLASDQ to find the SVD. * IF( N.LE.SMLSIZ ) THEN CALL SLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDVT, U, LDU, U, $ LDU, WORK, INFO ) RETURN END IF * * Set up the computation tree. * INODE = 1 NDIML = INODE + N NDIMR = NDIML + N IDXQ = NDIMR + N IWK = IDXQ + N CALL SLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), $ IWORK( NDIMR ), SMLSIZ ) * * For the nodes on bottom level of the tree, solve * their subproblems by SLASDQ. * NDB1 = ( ND+1 ) / 2 NCC = 0 DO 30 I = NDB1, ND * * IC : center row of each node * NL : number of rows of left subproblem * NR : number of rows of right subproblem * NLF: starting row of the left subproblem * NRF: starting row of the right subproblem * I1 = I - 1 IC = IWORK( INODE+I1 ) NL = IWORK( NDIML+I1 ) NLP1 = NL + 1 NR = IWORK( NDIMR+I1 ) NRP1 = NR + 1 NLF = IC - NL NRF = IC + 1 SQREI = 1 CALL SLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), E( NLF ), $ VT( NLF, NLF ), LDVT, U( NLF, NLF ), LDU, $ U( NLF, NLF ), LDU, WORK, INFO ) IF( INFO.NE.0 ) THEN RETURN END IF ITEMP = IDXQ + NLF - 2 DO 10 J = 1, NL IWORK( ITEMP+J ) = J 10 CONTINUE IF( I.EQ.ND ) THEN SQREI = SQRE ELSE SQREI = 1 END IF NRP1 = NR + SQREI CALL SLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), E( NRF ), $ VT( NRF, NRF ), LDVT, U( NRF, NRF ), LDU, $ U( NRF, NRF ), LDU, WORK, INFO ) IF( INFO.NE.0 ) THEN RETURN END IF ITEMP = IDXQ + IC DO 20 J = 1, NR IWORK( ITEMP+J-1 ) = J 20 CONTINUE 30 CONTINUE * * Now conquer each subproblem bottom-up. * DO 50 LVL = NLVL, 1, -1 * * Find the first node LF and last node LL on the * current level LVL. * IF( LVL.EQ.1 ) THEN LF = 1 LL = 1 ELSE LF = 2**( LVL-1 ) LL = 2*LF - 1 END IF DO 40 I = LF, LL IM1 = I - 1 IC = IWORK( INODE+IM1 ) NL = IWORK( NDIML+IM1 ) NR = IWORK( NDIMR+IM1 ) NLF = IC - NL IF( ( SQRE.EQ.0 ) .AND. ( I.EQ.LL ) ) THEN SQREI = SQRE ELSE SQREI = 1 END IF IDXQC = IDXQ + NLF - 1 ALPHA = D( IC ) BETA = E( IC ) CALL SLASD1( NL, NR, SQREI, D( NLF ), ALPHA, BETA, $ U( NLF, NLF ), LDU, VT( NLF, NLF ), LDVT, $ IWORK( IDXQC ), IWORK( IWK ), WORK, INFO ) IF( INFO.NE.0 ) THEN RETURN END IF 40 CONTINUE 50 CONTINUE * RETURN * * End of SLASD0 * END |