![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.2.2.
1: SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, 2: $ IDXQ, IWORK, 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, NL, NR, SQRE 11: DOUBLE PRECISION ALPHA, BETA 12: * .. 13: * .. Array Arguments .. 14: INTEGER IDXQ( * ), IWORK( * ) 15: DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, 22: * where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0. 23: * 24: * A related subroutine DLASD7 handles the case in which the singular 25: * values (and the singular vectors in factored form) are desired. 26: * 27: * DLASD1 computes the SVD as follows: 28: * 29: * ( D1(in) 0 0 0 ) 30: * B = U(in) * ( Z1' a Z2' b ) * VT(in) 31: * ( 0 0 D2(in) 0 ) 32: * 33: * = U(out) * ( D(out) 0) * VT(out) 34: * 35: * where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M 36: * with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros 37: * elsewhere; and the entry b is empty if SQRE = 0. 38: * 39: * The left singular vectors of the original matrix are stored in U, and 40: * the transpose of the right singular vectors are stored in VT, and the 41: * singular values are in D. The algorithm consists of three stages: 42: * 43: * The first stage consists of deflating the size of the problem 44: * when there are multiple singular values or when there are zeros in 45: * the Z vector. For each such occurence the dimension of the 46: * secular equation problem is reduced by one. This stage is 47: * performed by the routine DLASD2. 48: * 49: * The second stage consists of calculating the updated 50: * singular values. This is done by finding the square roots of the 51: * roots of the secular equation via the routine DLASD4 (as called 52: * by DLASD3). This routine also calculates the singular vectors of 53: * the current problem. 54: * 55: * The final stage consists of computing the updated singular vectors 56: * directly using the updated singular values. The singular vectors 57: * for the current problem are multiplied with the singular vectors 58: * from the overall problem. 59: * 60: * Arguments 61: * ========= 62: * 63: * NL (input) INTEGER 64: * The row dimension of the upper block. NL >= 1. 65: * 66: * NR (input) INTEGER 67: * The row dimension of the lower block. NR >= 1. 68: * 69: * SQRE (input) INTEGER 70: * = 0: the lower block is an NR-by-NR square matrix. 71: * = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 72: * 73: * The bidiagonal matrix has row dimension N = NL + NR + 1, 74: * and column dimension M = N + SQRE. 75: * 76: * D (input/output) DOUBLE PRECISION array, 77: * dimension (N = NL+NR+1). 78: * On entry D(1:NL,1:NL) contains the singular values of the 79: * upper block; and D(NL+2:N) contains the singular values of 80: * the lower block. On exit D(1:N) contains the singular values 81: * of the modified matrix. 82: * 83: * ALPHA (input/output) DOUBLE PRECISION 84: * Contains the diagonal element associated with the added row. 85: * 86: * BETA (input/output) DOUBLE PRECISION 87: * Contains the off-diagonal element associated with the added 88: * row. 89: * 90: * U (input/output) DOUBLE PRECISION array, dimension(LDU,N) 91: * On entry U(1:NL, 1:NL) contains the left singular vectors of 92: * the upper block; U(NL+2:N, NL+2:N) contains the left singular 93: * vectors of the lower block. On exit U contains the left 94: * singular vectors of the bidiagonal matrix. 95: * 96: * LDU (input) INTEGER 97: * The leading dimension of the array U. LDU >= max( 1, N ). 98: * 99: * VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M) 100: * where M = N + SQRE. 101: * On entry VT(1:NL+1, 1:NL+1)' contains the right singular 102: * vectors of the upper block; VT(NL+2:M, NL+2:M)' contains 103: * the right singular vectors of the lower block. On exit 104: * VT' contains the right singular vectors of the 105: * bidiagonal matrix. 106: * 107: * LDVT (input) INTEGER 108: * The leading dimension of the array VT. LDVT >= max( 1, M ). 109: * 110: * IDXQ (output) INTEGER array, dimension(N) 111: * This contains the permutation which will reintegrate the 112: * subproblem just solved back into sorted order, i.e. 113: * D( IDXQ( I = 1, N ) ) will be in ascending order. 114: * 115: * IWORK (workspace) INTEGER array, dimension( 4 * N ) 116: * 117: * WORK (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M ) 118: * 119: * INFO (output) INTEGER 120: * = 0: successful exit. 121: * < 0: if INFO = -i, the i-th argument had an illegal value. 122: * > 0: if INFO = 1, a singular value did not converge 123: * 124: * Further Details 125: * =============== 126: * 127: * Based on contributions by 128: * Ming Gu and Huan Ren, Computer Science Division, University of 129: * California at Berkeley, USA 130: * 131: * ===================================================================== 132: * 133: * .. Parameters .. 134: * 135: DOUBLE PRECISION ONE, ZERO 136: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 137: * .. 138: * .. Local Scalars .. 139: INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2, 140: $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2 141: DOUBLE PRECISION ORGNRM 142: * .. 143: * .. External Subroutines .. 144: EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA 145: * .. 146: * .. Intrinsic Functions .. 147: INTRINSIC ABS, MAX 148: * .. 149: * .. Executable Statements .. 150: * 151: * Test the input parameters. 152: * 153: INFO = 0 154: * 155: IF( NL.LT.1 ) THEN 156: INFO = -1 157: ELSE IF( NR.LT.1 ) THEN 158: INFO = -2 159: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 160: INFO = -3 161: END IF 162: IF( INFO.NE.0 ) THEN 163: CALL XERBLA( 'DLASD1', -INFO ) 164: RETURN 165: END IF 166: * 167: N = NL + NR + 1 168: M = N + SQRE 169: * 170: * The following values are for bookkeeping purposes only. They are 171: * integer pointers which indicate the portion of the workspace 172: * used by a particular array in DLASD2 and DLASD3. 173: * 174: LDU2 = N 175: LDVT2 = M 176: * 177: IZ = 1 178: ISIGMA = IZ + M 179: IU2 = ISIGMA + N 180: IVT2 = IU2 + LDU2*N 181: IQ = IVT2 + LDVT2*M 182: * 183: IDX = 1 184: IDXC = IDX + N 185: COLTYP = IDXC + N 186: IDXP = COLTYP + N 187: * 188: * Scale. 189: * 190: ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 191: D( NL+1 ) = ZERO 192: DO 10 I = 1, N 193: IF( ABS( D( I ) ).GT.ORGNRM ) THEN 194: ORGNRM = ABS( D( I ) ) 195: END IF 196: 10 CONTINUE 197: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 198: ALPHA = ALPHA / ORGNRM 199: BETA = BETA / ORGNRM 200: * 201: * Deflate singular values. 202: * 203: CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU, 204: $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2, 205: $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ), 206: $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO ) 207: * 208: * Solve Secular Equation and update singular vectors. 209: * 210: LDQ = K 211: CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ), 212: $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ), 213: $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ), 214: $ INFO ) 215: IF( INFO.NE.0 ) THEN 216: RETURN 217: END IF 218: * 219: * Unscale. 220: * 221: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 222: * 223: * Prepare the IDXQ sorting permutation. 224: * 225: N1 = K 226: N2 = N - K 227: CALL DLAMRG( N1, N2, D, 1, -1, IDXQ ) 228: * 229: RETURN 230: * 231: * End of DLASD1 232: * 233: END