Annotation of rpl/lapack/lapack/dlasd6.f, revision 1.11
1.11 ! bertrand 1: *> \brief \b DLASD6
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLASD6 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd6.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd6.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd6.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
! 22: * IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
! 23: * LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
! 24: * IWORK, INFO )
! 25: *
! 26: * .. Scalar Arguments ..
! 27: * INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
! 28: * $ NR, SQRE
! 29: * DOUBLE PRECISION ALPHA, BETA, C, S
! 30: * ..
! 31: * .. Array Arguments ..
! 32: * INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
! 33: * $ PERM( * )
! 34: * DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ),
! 35: * $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
! 36: * $ VF( * ), VL( * ), WORK( * ), Z( * )
! 37: * ..
! 38: *
! 39: *
! 40: *> \par Purpose:
! 41: * =============
! 42: *>
! 43: *> \verbatim
! 44: *>
! 45: *> DLASD6 computes the SVD of an updated upper bidiagonal matrix B
! 46: *> obtained by merging two smaller ones by appending a row. This
! 47: *> routine is used only for the problem which requires all singular
! 48: *> values and optionally singular vector matrices in factored form.
! 49: *> B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE.
! 50: *> A related subroutine, DLASD1, handles the case in which all singular
! 51: *> values and singular vectors of the bidiagonal matrix are desired.
! 52: *>
! 53: *> DLASD6 computes the SVD as follows:
! 54: *>
! 55: *> ( D1(in) 0 0 0 )
! 56: *> B = U(in) * ( Z1**T a Z2**T b ) * VT(in)
! 57: *> ( 0 0 D2(in) 0 )
! 58: *>
! 59: *> = U(out) * ( D(out) 0) * VT(out)
! 60: *>
! 61: *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
! 62: *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
! 63: *> elsewhere; and the entry b is empty if SQRE = 0.
! 64: *>
! 65: *> The singular values of B can be computed using D1, D2, the first
! 66: *> components of all the right singular vectors of the lower block, and
! 67: *> the last components of all the right singular vectors of the upper
! 68: *> block. These components are stored and updated in VF and VL,
! 69: *> respectively, in DLASD6. Hence U and VT are not explicitly
! 70: *> referenced.
! 71: *>
! 72: *> The singular values are stored in D. The algorithm consists of two
! 73: *> stages:
! 74: *>
! 75: *> The first stage consists of deflating the size of the problem
! 76: *> when there are multiple singular values or if there is a zero
! 77: *> in the Z vector. For each such occurence the dimension of the
! 78: *> secular equation problem is reduced by one. This stage is
! 79: *> performed by the routine DLASD7.
! 80: *>
! 81: *> The second stage consists of calculating the updated
! 82: *> singular values. This is done by finding the roots of the
! 83: *> secular equation via the routine DLASD4 (as called by DLASD8).
! 84: *> This routine also updates VF and VL and computes the distances
! 85: *> between the updated singular values and the old singular
! 86: *> values.
! 87: *>
! 88: *> DLASD6 is called from DLASDA.
! 89: *> \endverbatim
! 90: *
! 91: * Arguments:
! 92: * ==========
! 93: *
! 94: *> \param[in] ICOMPQ
! 95: *> \verbatim
! 96: *> ICOMPQ is INTEGER
! 97: *> Specifies whether singular vectors are to be computed in
! 98: *> factored form:
! 99: *> = 0: Compute singular values only.
! 100: *> = 1: Compute singular vectors in factored form as well.
! 101: *> \endverbatim
! 102: *>
! 103: *> \param[in] NL
! 104: *> \verbatim
! 105: *> NL is INTEGER
! 106: *> The row dimension of the upper block. NL >= 1.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[in] NR
! 110: *> \verbatim
! 111: *> NR is INTEGER
! 112: *> The row dimension of the lower block. NR >= 1.
! 113: *> \endverbatim
! 114: *>
! 115: *> \param[in] SQRE
! 116: *> \verbatim
! 117: *> SQRE is INTEGER
! 118: *> = 0: the lower block is an NR-by-NR square matrix.
! 119: *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
! 120: *>
! 121: *> The bidiagonal matrix has row dimension N = NL + NR + 1,
! 122: *> and column dimension M = N + SQRE.
! 123: *> \endverbatim
! 124: *>
! 125: *> \param[in,out] D
! 126: *> \verbatim
! 127: *> D is DOUBLE PRECISION array, dimension ( NL+NR+1 ).
! 128: *> On entry D(1:NL,1:NL) contains the singular values of the
! 129: *> upper block, and D(NL+2:N) contains the singular values
! 130: *> of the lower block. On exit D(1:N) contains the singular
! 131: *> values of the modified matrix.
! 132: *> \endverbatim
! 133: *>
! 134: *> \param[in,out] VF
! 135: *> \verbatim
! 136: *> VF is DOUBLE PRECISION array, dimension ( M )
! 137: *> On entry, VF(1:NL+1) contains the first components of all
! 138: *> right singular vectors of the upper block; and VF(NL+2:M)
! 139: *> contains the first components of all right singular vectors
! 140: *> of the lower block. On exit, VF contains the first components
! 141: *> of all right singular vectors of the bidiagonal matrix.
! 142: *> \endverbatim
! 143: *>
! 144: *> \param[in,out] VL
! 145: *> \verbatim
! 146: *> VL is DOUBLE PRECISION array, dimension ( M )
! 147: *> On entry, VL(1:NL+1) contains the last components of all
! 148: *> right singular vectors of the upper block; and VL(NL+2:M)
! 149: *> contains the last components of all right singular vectors of
! 150: *> the lower block. On exit, VL contains the last components of
! 151: *> all right singular vectors of the bidiagonal matrix.
! 152: *> \endverbatim
! 153: *>
! 154: *> \param[in,out] ALPHA
! 155: *> \verbatim
! 156: *> ALPHA is DOUBLE PRECISION
! 157: *> Contains the diagonal element associated with the added row.
! 158: *> \endverbatim
! 159: *>
! 160: *> \param[in,out] BETA
! 161: *> \verbatim
! 162: *> BETA is DOUBLE PRECISION
! 163: *> Contains the off-diagonal element associated with the added
! 164: *> row.
! 165: *> \endverbatim
! 166: *>
! 167: *> \param[out] IDXQ
! 168: *> \verbatim
! 169: *> IDXQ is INTEGER array, dimension ( N )
! 170: *> This contains the permutation which will reintegrate the
! 171: *> subproblem just solved back into sorted order, i.e.
! 172: *> D( IDXQ( I = 1, N ) ) will be in ascending order.
! 173: *> \endverbatim
! 174: *>
! 175: *> \param[out] PERM
! 176: *> \verbatim
! 177: *> PERM is INTEGER array, dimension ( N )
! 178: *> The permutations (from deflation and sorting) to be applied
! 179: *> to each block. Not referenced if ICOMPQ = 0.
! 180: *> \endverbatim
! 181: *>
! 182: *> \param[out] GIVPTR
! 183: *> \verbatim
! 184: *> GIVPTR is INTEGER
! 185: *> The number of Givens rotations which took place in this
! 186: *> subproblem. Not referenced if ICOMPQ = 0.
! 187: *> \endverbatim
! 188: *>
! 189: *> \param[out] GIVCOL
! 190: *> \verbatim
! 191: *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
! 192: *> Each pair of numbers indicates a pair of columns to take place
! 193: *> in a Givens rotation. Not referenced if ICOMPQ = 0.
! 194: *> \endverbatim
! 195: *>
! 196: *> \param[in] LDGCOL
! 197: *> \verbatim
! 198: *> LDGCOL is INTEGER
! 199: *> leading dimension of GIVCOL, must be at least N.
! 200: *> \endverbatim
! 201: *>
! 202: *> \param[out] GIVNUM
! 203: *> \verbatim
! 204: *> GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
! 205: *> Each number indicates the C or S value to be used in the
! 206: *> corresponding Givens rotation. Not referenced if ICOMPQ = 0.
! 207: *> \endverbatim
! 208: *>
! 209: *> \param[in] LDGNUM
! 210: *> \verbatim
! 211: *> LDGNUM is INTEGER
! 212: *> The leading dimension of GIVNUM and POLES, must be at least N.
! 213: *> \endverbatim
! 214: *>
! 215: *> \param[out] POLES
! 216: *> \verbatim
! 217: *> POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
! 218: *> On exit, POLES(1,*) is an array containing the new singular
! 219: *> values obtained from solving the secular equation, and
! 220: *> POLES(2,*) is an array containing the poles in the secular
! 221: *> equation. Not referenced if ICOMPQ = 0.
! 222: *> \endverbatim
! 223: *>
! 224: *> \param[out] DIFL
! 225: *> \verbatim
! 226: *> DIFL is DOUBLE PRECISION array, dimension ( N )
! 227: *> On exit, DIFL(I) is the distance between I-th updated
! 228: *> (undeflated) singular value and the I-th (undeflated) old
! 229: *> singular value.
! 230: *> \endverbatim
! 231: *>
! 232: *> \param[out] DIFR
! 233: *> \verbatim
! 234: *> DIFR is DOUBLE PRECISION array,
! 235: *> dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and
! 236: *> dimension ( N ) if ICOMPQ = 0.
! 237: *> On exit, DIFR(I, 1) is the distance between I-th updated
! 238: *> (undeflated) singular value and the I+1-th (undeflated) old
! 239: *> singular value.
! 240: *>
! 241: *> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
! 242: *> normalizing factors for the right singular vector matrix.
! 243: *>
! 244: *> See DLASD8 for details on DIFL and DIFR.
! 245: *> \endverbatim
! 246: *>
! 247: *> \param[out] Z
! 248: *> \verbatim
! 249: *> Z is DOUBLE PRECISION array, dimension ( M )
! 250: *> The first elements of this array contain the components
! 251: *> of the deflation-adjusted updating row vector.
! 252: *> \endverbatim
! 253: *>
! 254: *> \param[out] K
! 255: *> \verbatim
! 256: *> K is INTEGER
! 257: *> Contains the dimension of the non-deflated matrix,
! 258: *> This is the order of the related secular equation. 1 <= K <=N.
! 259: *> \endverbatim
! 260: *>
! 261: *> \param[out] C
! 262: *> \verbatim
! 263: *> C is DOUBLE PRECISION
! 264: *> C contains garbage if SQRE =0 and the C-value of a Givens
! 265: *> rotation related to the right null space if SQRE = 1.
! 266: *> \endverbatim
! 267: *>
! 268: *> \param[out] S
! 269: *> \verbatim
! 270: *> S is DOUBLE PRECISION
! 271: *> S contains garbage if SQRE =0 and the S-value of a Givens
! 272: *> rotation related to the right null space if SQRE = 1.
! 273: *> \endverbatim
! 274: *>
! 275: *> \param[out] WORK
! 276: *> \verbatim
! 277: *> WORK is DOUBLE PRECISION array, dimension ( 4 * M )
! 278: *> \endverbatim
! 279: *>
! 280: *> \param[out] IWORK
! 281: *> \verbatim
! 282: *> IWORK is INTEGER array, dimension ( 3 * N )
! 283: *> \endverbatim
! 284: *>
! 285: *> \param[out] INFO
! 286: *> \verbatim
! 287: *> INFO is INTEGER
! 288: *> = 0: successful exit.
! 289: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 290: *> > 0: if INFO = 1, a singular value did not converge
! 291: *> \endverbatim
! 292: *
! 293: * Authors:
! 294: * ========
! 295: *
! 296: *> \author Univ. of Tennessee
! 297: *> \author Univ. of California Berkeley
! 298: *> \author Univ. of Colorado Denver
! 299: *> \author NAG Ltd.
! 300: *
! 301: *> \date November 2011
! 302: *
! 303: *> \ingroup auxOTHERauxiliary
! 304: *
! 305: *> \par Contributors:
! 306: * ==================
! 307: *>
! 308: *> Ming Gu and Huan Ren, Computer Science Division, University of
! 309: *> California at Berkeley, USA
! 310: *>
! 311: * =====================================================================
1.1 bertrand 312: SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
313: $ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
314: $ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
315: $ IWORK, INFO )
316: *
1.11 ! bertrand 317: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 318: * -- LAPACK is a software package provided by Univ. of Tennessee, --
319: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11 ! bertrand 320: * November 2011
1.1 bertrand 321: *
322: * .. Scalar Arguments ..
323: INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
324: $ NR, SQRE
325: DOUBLE PRECISION ALPHA, BETA, C, S
326: * ..
327: * .. Array Arguments ..
328: INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
329: $ PERM( * )
330: DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ),
331: $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
332: $ VF( * ), VL( * ), WORK( * ), Z( * )
333: * ..
334: *
335: * =====================================================================
336: *
337: * .. Parameters ..
338: DOUBLE PRECISION ONE, ZERO
339: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
340: * ..
341: * .. Local Scalars ..
342: INTEGER I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M,
343: $ N, N1, N2
344: DOUBLE PRECISION ORGNRM
345: * ..
346: * .. External Subroutines ..
347: EXTERNAL DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA
348: * ..
349: * .. Intrinsic Functions ..
350: INTRINSIC ABS, MAX
351: * ..
352: * .. Executable Statements ..
353: *
354: * Test the input parameters.
355: *
356: INFO = 0
357: N = NL + NR + 1
358: M = N + SQRE
359: *
360: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
361: INFO = -1
362: ELSE IF( NL.LT.1 ) THEN
363: INFO = -2
364: ELSE IF( NR.LT.1 ) THEN
365: INFO = -3
366: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
367: INFO = -4
368: ELSE IF( LDGCOL.LT.N ) THEN
369: INFO = -14
370: ELSE IF( LDGNUM.LT.N ) THEN
371: INFO = -16
372: END IF
373: IF( INFO.NE.0 ) THEN
374: CALL XERBLA( 'DLASD6', -INFO )
375: RETURN
376: END IF
377: *
378: * The following values are for bookkeeping purposes only. They are
379: * integer pointers which indicate the portion of the workspace
380: * used by a particular array in DLASD7 and DLASD8.
381: *
382: ISIGMA = 1
383: IW = ISIGMA + N
384: IVFW = IW + M
385: IVLW = IVFW + M
386: *
387: IDX = 1
388: IDXC = IDX + N
389: IDXP = IDXC + N
390: *
391: * Scale.
392: *
393: ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
394: D( NL+1 ) = ZERO
395: DO 10 I = 1, N
396: IF( ABS( D( I ) ).GT.ORGNRM ) THEN
397: ORGNRM = ABS( D( I ) )
398: END IF
399: 10 CONTINUE
400: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
401: ALPHA = ALPHA / ORGNRM
402: BETA = BETA / ORGNRM
403: *
404: * Sort and Deflate singular values.
405: *
406: CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF,
407: $ WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA,
408: $ WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ,
409: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S,
410: $ INFO )
411: *
412: * Solve Secular Equation, compute DIFL, DIFR, and update VF, VL.
413: *
414: CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM,
415: $ WORK( ISIGMA ), WORK( IW ), INFO )
416: *
1.8 bertrand 417: * Handle error returned
418: *
419: IF( INFO.NE.0 ) THEN
420: CALL XERBLA( 'DLASD8', -INFO )
421: RETURN
422: END IF
423: *
1.1 bertrand 424: * Save the poles if ICOMPQ = 1.
425: *
426: IF( ICOMPQ.EQ.1 ) THEN
427: CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 )
428: CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 )
429: END IF
430: *
431: * Unscale.
432: *
433: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
434: *
435: * Prepare the IDXQ sorting permutation.
436: *
437: N1 = K
438: N2 = N - K
439: CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
440: *
441: RETURN
442: *
443: * End of DLASD6
444: *
445: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd6.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack