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