![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, 2: $ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, 3: $ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, 4: $ IWORK, INFO ) 5: * 6: * -- LAPACK auxiliary routine (version 3.3.0) -- 7: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 9: * November 2010 10: * 11: * .. Scalar Arguments .. 12: INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 13: $ NR, SQRE 14: DOUBLE PRECISION ALPHA, BETA, C, S 15: * .. 16: * .. Array Arguments .. 17: INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ), 18: $ PERM( * ) 19: DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ), 20: $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), 21: $ VF( * ), VL( * ), WORK( * ), Z( * ) 22: * .. 23: * 24: * Purpose 25: * ======= 26: * 27: * DLASD6 computes the SVD of an updated upper bidiagonal matrix B 28: * obtained by merging two smaller ones by appending a row. This 29: * routine is used only for the problem which requires all singular 30: * values and optionally singular vector matrices in factored form. 31: * B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE. 32: * A related subroutine, DLASD1, handles the case in which all singular 33: * values and singular vectors of the bidiagonal matrix are desired. 34: * 35: * DLASD6 computes the SVD as follows: 36: * 37: * ( D1(in) 0 0 0 ) 38: * B = U(in) * ( Z1' a Z2' b ) * VT(in) 39: * ( 0 0 D2(in) 0 ) 40: * 41: * = U(out) * ( D(out) 0) * VT(out) 42: * 43: * where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M 44: * with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros 45: * elsewhere; and the entry b is empty if SQRE = 0. 46: * 47: * The singular values of B can be computed using D1, D2, the first 48: * components of all the right singular vectors of the lower block, and 49: * the last components of all the right singular vectors of the upper 50: * block. These components are stored and updated in VF and VL, 51: * respectively, in DLASD6. Hence U and VT are not explicitly 52: * referenced. 53: * 54: * The singular values are stored in D. The algorithm consists of two 55: * stages: 56: * 57: * The first stage consists of deflating the size of the problem 58: * when there are multiple singular values or if there is a zero 59: * in the Z vector. For each such occurence the dimension of the 60: * secular equation problem is reduced by one. This stage is 61: * performed by the routine DLASD7. 62: * 63: * The second stage consists of calculating the updated 64: * singular values. This is done by finding the roots of the 65: * secular equation via the routine DLASD4 (as called by DLASD8). 66: * This routine also updates VF and VL and computes the distances 67: * between the updated singular values and the old singular 68: * values. 69: * 70: * DLASD6 is called from DLASDA. 71: * 72: * Arguments 73: * ========= 74: * 75: * ICOMPQ (input) INTEGER 76: * Specifies whether singular vectors are to be computed in 77: * factored form: 78: * = 0: Compute singular values only. 79: * = 1: Compute singular vectors in factored form as well. 80: * 81: * NL (input) INTEGER 82: * The row dimension of the upper block. NL >= 1. 83: * 84: * NR (input) INTEGER 85: * The row dimension of the lower block. NR >= 1. 86: * 87: * SQRE (input) INTEGER 88: * = 0: the lower block is an NR-by-NR square matrix. 89: * = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 90: * 91: * The bidiagonal matrix has row dimension N = NL + NR + 1, 92: * and column dimension M = N + SQRE. 93: * 94: * D (input/output) DOUBLE PRECISION array, dimension ( NL+NR+1 ). 95: * On entry D(1:NL,1:NL) contains the singular values of the 96: * upper block, and D(NL+2:N) contains the singular values 97: * of the lower block. On exit D(1:N) contains the singular 98: * values of the modified matrix. 99: * 100: * VF (input/output) DOUBLE PRECISION array, dimension ( M ) 101: * On entry, VF(1:NL+1) contains the first components of all 102: * right singular vectors of the upper block; and VF(NL+2:M) 103: * contains the first components of all right singular vectors 104: * of the lower block. On exit, VF contains the first components 105: * of all right singular vectors of the bidiagonal matrix. 106: * 107: * VL (input/output) DOUBLE PRECISION array, dimension ( M ) 108: * On entry, VL(1:NL+1) contains the last components of all 109: * right singular vectors of the upper block; and VL(NL+2:M) 110: * contains the last components of all right singular vectors of 111: * the lower block. On exit, VL contains the last components of 112: * all right singular vectors of the bidiagonal matrix. 113: * 114: * ALPHA (input/output) DOUBLE PRECISION 115: * Contains the diagonal element associated with the added row. 116: * 117: * BETA (input/output) DOUBLE PRECISION 118: * Contains the off-diagonal element associated with the added 119: * row. 120: * 121: * IDXQ (output) INTEGER array, dimension ( N ) 122: * This contains the permutation which will reintegrate the 123: * subproblem just solved back into sorted order, i.e. 124: * D( IDXQ( I = 1, N ) ) will be in ascending order. 125: * 126: * PERM (output) INTEGER array, dimension ( N ) 127: * The permutations (from deflation and sorting) to be applied 128: * to each block. Not referenced if ICOMPQ = 0. 129: * 130: * GIVPTR (output) INTEGER 131: * The number of Givens rotations which took place in this 132: * subproblem. Not referenced if ICOMPQ = 0. 133: * 134: * GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 ) 135: * Each pair of numbers indicates a pair of columns to take place 136: * in a Givens rotation. Not referenced if ICOMPQ = 0. 137: * 138: * LDGCOL (input) INTEGER 139: * leading dimension of GIVCOL, must be at least N. 140: * 141: * GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 142: * Each number indicates the C or S value to be used in the 143: * corresponding Givens rotation. Not referenced if ICOMPQ = 0. 144: * 145: * LDGNUM (input) INTEGER 146: * The leading dimension of GIVNUM and POLES, must be at least N. 147: * 148: * POLES (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 149: * On exit, POLES(1,*) is an array containing the new singular 150: * values obtained from solving the secular equation, and 151: * POLES(2,*) is an array containing the poles in the secular 152: * equation. Not referenced if ICOMPQ = 0. 153: * 154: * DIFL (output) DOUBLE PRECISION array, dimension ( N ) 155: * On exit, DIFL(I) is the distance between I-th updated 156: * (undeflated) singular value and the I-th (undeflated) old 157: * singular value. 158: * 159: * DIFR (output) DOUBLE PRECISION array, 160: * dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and 161: * dimension ( N ) if ICOMPQ = 0. 162: * On exit, DIFR(I, 1) is the distance between I-th updated 163: * (undeflated) singular value and the I+1-th (undeflated) old 164: * singular value. 165: * 166: * If ICOMPQ = 1, DIFR(1:K,2) is an array containing the 167: * normalizing factors for the right singular vector matrix. 168: * 169: * See DLASD8 for details on DIFL and DIFR. 170: * 171: * Z (output) DOUBLE PRECISION array, dimension ( M ) 172: * The first elements of this array contain the components 173: * of the deflation-adjusted updating row vector. 174: * 175: * K (output) INTEGER 176: * Contains the dimension of the non-deflated matrix, 177: * This is the order of the related secular equation. 1 <= K <=N. 178: * 179: * C (output) DOUBLE PRECISION 180: * C contains garbage if SQRE =0 and the C-value of a Givens 181: * rotation related to the right null space if SQRE = 1. 182: * 183: * S (output) DOUBLE PRECISION 184: * S contains garbage if SQRE =0 and the S-value of a Givens 185: * rotation related to the right null space if SQRE = 1. 186: * 187: * WORK (workspace) DOUBLE PRECISION array, dimension ( 4 * M ) 188: * 189: * IWORK (workspace) INTEGER array, dimension ( 3 * N ) 190: * 191: * INFO (output) INTEGER 192: * = 0: successful exit. 193: * < 0: if INFO = -i, the i-th argument had an illegal value. 194: * > 0: if INFO = 1, a singular value did not converge 195: * 196: * Further Details 197: * =============== 198: * 199: * Based on contributions by 200: * Ming Gu and Huan Ren, Computer Science Division, University of 201: * California at Berkeley, USA 202: * 203: * ===================================================================== 204: * 205: * .. Parameters .. 206: DOUBLE PRECISION ONE, ZERO 207: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 208: * .. 209: * .. Local Scalars .. 210: INTEGER I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M, 211: $ N, N1, N2 212: DOUBLE PRECISION ORGNRM 213: * .. 214: * .. External Subroutines .. 215: EXTERNAL DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA 216: * .. 217: * .. Intrinsic Functions .. 218: INTRINSIC ABS, MAX 219: * .. 220: * .. Executable Statements .. 221: * 222: * Test the input parameters. 223: * 224: INFO = 0 225: N = NL + NR + 1 226: M = N + SQRE 227: * 228: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 229: INFO = -1 230: ELSE IF( NL.LT.1 ) THEN 231: INFO = -2 232: ELSE IF( NR.LT.1 ) THEN 233: INFO = -3 234: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 235: INFO = -4 236: ELSE IF( LDGCOL.LT.N ) THEN 237: INFO = -14 238: ELSE IF( LDGNUM.LT.N ) THEN 239: INFO = -16 240: END IF 241: IF( INFO.NE.0 ) THEN 242: CALL XERBLA( 'DLASD6', -INFO ) 243: RETURN 244: END IF 245: * 246: * The following values are for bookkeeping purposes only. They are 247: * integer pointers which indicate the portion of the workspace 248: * used by a particular array in DLASD7 and DLASD8. 249: * 250: ISIGMA = 1 251: IW = ISIGMA + N 252: IVFW = IW + M 253: IVLW = IVFW + M 254: * 255: IDX = 1 256: IDXC = IDX + N 257: IDXP = IDXC + N 258: * 259: * Scale. 260: * 261: ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 262: D( NL+1 ) = ZERO 263: DO 10 I = 1, N 264: IF( ABS( D( I ) ).GT.ORGNRM ) THEN 265: ORGNRM = ABS( D( I ) ) 266: END IF 267: 10 CONTINUE 268: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 269: ALPHA = ALPHA / ORGNRM 270: BETA = BETA / ORGNRM 271: * 272: * Sort and Deflate singular values. 273: * 274: CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF, 275: $ WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA, 276: $ WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ, 277: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, 278: $ INFO ) 279: * 280: * Solve Secular Equation, compute DIFL, DIFR, and update VF, VL. 281: * 282: CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM, 283: $ WORK( ISIGMA ), WORK( IW ), INFO ) 284: * 285: * Handle error returned 286: * 287: IF( INFO.NE.0 ) THEN 288: CALL XERBLA( 'DLASD8', -INFO ) 289: RETURN 290: END IF 291: * 292: * Save the poles if ICOMPQ = 1. 293: * 294: IF( ICOMPQ.EQ.1 ) THEN 295: CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 ) 296: CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 ) 297: END IF 298: * 299: * Unscale. 300: * 301: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 302: * 303: * Prepare the IDXQ sorting permutation. 304: * 305: N1 = K 306: N2 = N - K 307: CALL DLAMRG( N1, N2, D, 1, -1, IDXQ ) 308: * 309: RETURN 310: * 311: * End of DLASD6 312: * 313: END