Annotation of rpl/lapack/lapack/dlasd7.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b DLASD7
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLASD7 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd7.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd7.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd7.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
! 22: * VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
! 23: * PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
! 24: * C, S, 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, * ), IDX( * ), IDXP( * ),
! 33: * $ IDXQ( * ), PERM( * )
! 34: * DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
! 35: * $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
! 36: * $ ZW( * )
! 37: * ..
! 38: *
! 39: *
! 40: *> \par Purpose:
! 41: * =============
! 42: *>
! 43: *> \verbatim
! 44: *>
! 45: *> DLASD7 merges the two sets of singular values together into a single
! 46: *> sorted set. Then it tries to deflate the size of the problem. There
! 47: *> are two ways in which deflation can occur: when two or more singular
! 48: *> values are close together or if there is a tiny entry in the Z
! 49: *> vector. For each such occurrence the order of the related
! 50: *> secular equation problem is reduced by one.
! 51: *>
! 52: *> DLASD7 is called from DLASD6.
! 53: *> \endverbatim
! 54: *
! 55: * Arguments:
! 56: * ==========
! 57: *
! 58: *> \param[in] ICOMPQ
! 59: *> \verbatim
! 60: *> ICOMPQ is INTEGER
! 61: *> Specifies whether singular vectors are to be computed
! 62: *> in compact form, as follows:
! 63: *> = 0: Compute singular values only.
! 64: *> = 1: Compute singular vectors of upper
! 65: *> bidiagonal matrix in compact form.
! 66: *> \endverbatim
! 67: *>
! 68: *> \param[in] NL
! 69: *> \verbatim
! 70: *> NL is INTEGER
! 71: *> The row dimension of the upper block. NL >= 1.
! 72: *> \endverbatim
! 73: *>
! 74: *> \param[in] NR
! 75: *> \verbatim
! 76: *> NR is INTEGER
! 77: *> The row dimension of the lower block. NR >= 1.
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[in] SQRE
! 81: *> \verbatim
! 82: *> SQRE is INTEGER
! 83: *> = 0: the lower block is an NR-by-NR square matrix.
! 84: *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
! 85: *>
! 86: *> The bidiagonal matrix has
! 87: *> N = NL + NR + 1 rows and
! 88: *> M = N + SQRE >= N columns.
! 89: *> \endverbatim
! 90: *>
! 91: *> \param[out] K
! 92: *> \verbatim
! 93: *> K is INTEGER
! 94: *> Contains the dimension of the non-deflated matrix, this is
! 95: *> the order of the related secular equation. 1 <= K <=N.
! 96: *> \endverbatim
! 97: *>
! 98: *> \param[in,out] D
! 99: *> \verbatim
! 100: *> D is DOUBLE PRECISION array, dimension ( N )
! 101: *> On entry D contains the singular values of the two submatrices
! 102: *> to be combined. On exit D contains the trailing (N-K) updated
! 103: *> singular values (those which were deflated) sorted into
! 104: *> increasing order.
! 105: *> \endverbatim
! 106: *>
! 107: *> \param[out] Z
! 108: *> \verbatim
! 109: *> Z is DOUBLE PRECISION array, dimension ( M )
! 110: *> On exit Z contains the updating row vector in the secular
! 111: *> equation.
! 112: *> \endverbatim
! 113: *>
! 114: *> \param[out] ZW
! 115: *> \verbatim
! 116: *> ZW is DOUBLE PRECISION array, dimension ( M )
! 117: *> Workspace for Z.
! 118: *> \endverbatim
! 119: *>
! 120: *> \param[in,out] VF
! 121: *> \verbatim
! 122: *> VF is DOUBLE PRECISION array, dimension ( M )
! 123: *> On entry, VF(1:NL+1) contains the first components of all
! 124: *> right singular vectors of the upper block; and VF(NL+2:M)
! 125: *> contains the first components of all right singular vectors
! 126: *> of the lower block. On exit, VF contains the first components
! 127: *> of all right singular vectors of the bidiagonal matrix.
! 128: *> \endverbatim
! 129: *>
! 130: *> \param[out] VFW
! 131: *> \verbatim
! 132: *> VFW is DOUBLE PRECISION array, dimension ( M )
! 133: *> Workspace for VF.
! 134: *> \endverbatim
! 135: *>
! 136: *> \param[in,out] VL
! 137: *> \verbatim
! 138: *> VL is DOUBLE PRECISION array, dimension ( M )
! 139: *> On entry, VL(1:NL+1) contains the last components of all
! 140: *> right singular vectors of the upper block; and VL(NL+2:M)
! 141: *> contains the last components of all right singular vectors
! 142: *> of the lower block. On exit, VL contains the last components
! 143: *> of all right singular vectors of the bidiagonal matrix.
! 144: *> \endverbatim
! 145: *>
! 146: *> \param[out] VLW
! 147: *> \verbatim
! 148: *> VLW is DOUBLE PRECISION array, dimension ( M )
! 149: *> Workspace for VL.
! 150: *> \endverbatim
! 151: *>
! 152: *> \param[in] ALPHA
! 153: *> \verbatim
! 154: *> ALPHA is DOUBLE PRECISION
! 155: *> Contains the diagonal element associated with the added row.
! 156: *> \endverbatim
! 157: *>
! 158: *> \param[in] BETA
! 159: *> \verbatim
! 160: *> BETA is DOUBLE PRECISION
! 161: *> Contains the off-diagonal element associated with the added
! 162: *> row.
! 163: *> \endverbatim
! 164: *>
! 165: *> \param[out] DSIGMA
! 166: *> \verbatim
! 167: *> DSIGMA is DOUBLE PRECISION array, dimension ( N )
! 168: *> Contains a copy of the diagonal elements (K-1 singular values
! 169: *> and one zero) in the secular equation.
! 170: *> \endverbatim
! 171: *>
! 172: *> \param[out] IDX
! 173: *> \verbatim
! 174: *> IDX is INTEGER array, dimension ( N )
! 175: *> This will contain the permutation used to sort the contents of
! 176: *> D into ascending order.
! 177: *> \endverbatim
! 178: *>
! 179: *> \param[out] IDXP
! 180: *> \verbatim
! 181: *> IDXP is INTEGER array, dimension ( N )
! 182: *> This will contain the permutation used to place deflated
! 183: *> values of D at the end of the array. On output IDXP(2:K)
! 184: *> points to the nondeflated D-values and IDXP(K+1:N)
! 185: *> points to the deflated singular values.
! 186: *> \endverbatim
! 187: *>
! 188: *> \param[in] IDXQ
! 189: *> \verbatim
! 190: *> IDXQ is INTEGER array, dimension ( N )
! 191: *> This contains the permutation which separately sorts the two
! 192: *> sub-problems in D into ascending order. Note that entries in
! 193: *> the first half of this permutation must first be moved one
! 194: *> position backward; and entries in the second half
! 195: *> must first have NL+1 added to their values.
! 196: *> \endverbatim
! 197: *>
! 198: *> \param[out] PERM
! 199: *> \verbatim
! 200: *> PERM is INTEGER array, dimension ( N )
! 201: *> The permutations (from deflation and sorting) to be applied
! 202: *> to each singular block. Not referenced if ICOMPQ = 0.
! 203: *> \endverbatim
! 204: *>
! 205: *> \param[out] GIVPTR
! 206: *> \verbatim
! 207: *> GIVPTR is INTEGER
! 208: *> The number of Givens rotations which took place in this
! 209: *> subproblem. Not referenced if ICOMPQ = 0.
! 210: *> \endverbatim
! 211: *>
! 212: *> \param[out] GIVCOL
! 213: *> \verbatim
! 214: *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
! 215: *> Each pair of numbers indicates a pair of columns to take place
! 216: *> in a Givens rotation. Not referenced if ICOMPQ = 0.
! 217: *> \endverbatim
! 218: *>
! 219: *> \param[in] LDGCOL
! 220: *> \verbatim
! 221: *> LDGCOL is INTEGER
! 222: *> The leading dimension of GIVCOL, must be at least N.
! 223: *> \endverbatim
! 224: *>
! 225: *> \param[out] GIVNUM
! 226: *> \verbatim
! 227: *> GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
! 228: *> Each number indicates the C or S value to be used in the
! 229: *> corresponding Givens rotation. Not referenced if ICOMPQ = 0.
! 230: *> \endverbatim
! 231: *>
! 232: *> \param[in] LDGNUM
! 233: *> \verbatim
! 234: *> LDGNUM is INTEGER
! 235: *> The leading dimension of GIVNUM, must be at least N.
! 236: *> \endverbatim
! 237: *>
! 238: *> \param[out] C
! 239: *> \verbatim
! 240: *> C is DOUBLE PRECISION
! 241: *> C contains garbage if SQRE =0 and the C-value of a Givens
! 242: *> rotation related to the right null space if SQRE = 1.
! 243: *> \endverbatim
! 244: *>
! 245: *> \param[out] S
! 246: *> \verbatim
! 247: *> S is DOUBLE PRECISION
! 248: *> S contains garbage if SQRE =0 and the S-value of a Givens
! 249: *> rotation related to the right null space if SQRE = 1.
! 250: *> \endverbatim
! 251: *>
! 252: *> \param[out] INFO
! 253: *> \verbatim
! 254: *> INFO is INTEGER
! 255: *> = 0: successful exit.
! 256: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 257: *> \endverbatim
! 258: *
! 259: * Authors:
! 260: * ========
! 261: *
! 262: *> \author Univ. of Tennessee
! 263: *> \author Univ. of California Berkeley
! 264: *> \author Univ. of Colorado Denver
! 265: *> \author NAG Ltd.
! 266: *
! 267: *> \date November 2011
! 268: *
! 269: *> \ingroup auxOTHERauxiliary
! 270: *
! 271: *> \par Contributors:
! 272: * ==================
! 273: *>
! 274: *> Ming Gu and Huan Ren, Computer Science Division, University of
! 275: *> California at Berkeley, USA
! 276: *>
! 277: * =====================================================================
1.1 bertrand 278: SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
279: $ VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
280: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
281: $ C, S, INFO )
282: *
1.8 ! bertrand 283: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 284: * -- LAPACK is a software package provided by Univ. of Tennessee, --
285: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 286: * November 2011
1.1 bertrand 287: *
288: * .. Scalar Arguments ..
289: INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
290: $ NR, SQRE
291: DOUBLE PRECISION ALPHA, BETA, C, S
292: * ..
293: * .. Array Arguments ..
294: INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
295: $ IDXQ( * ), PERM( * )
296: DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
297: $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
298: $ ZW( * )
299: * ..
300: *
301: * =====================================================================
302: *
303: * .. Parameters ..
304: DOUBLE PRECISION ZERO, ONE, TWO, EIGHT
305: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
306: $ EIGHT = 8.0D+0 )
307: * ..
308: * .. Local Scalars ..
309: *
310: INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
311: $ NLP1, NLP2
312: DOUBLE PRECISION EPS, HLFTOL, TAU, TOL, Z1
313: * ..
314: * .. External Subroutines ..
315: EXTERNAL DCOPY, DLAMRG, DROT, XERBLA
316: * ..
317: * .. External Functions ..
318: DOUBLE PRECISION DLAMCH, DLAPY2
319: EXTERNAL DLAMCH, DLAPY2
320: * ..
321: * .. Intrinsic Functions ..
322: INTRINSIC ABS, MAX
323: * ..
324: * .. Executable Statements ..
325: *
326: * Test the input parameters.
327: *
328: INFO = 0
329: N = NL + NR + 1
330: M = N + SQRE
331: *
332: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
333: INFO = -1
334: ELSE IF( NL.LT.1 ) THEN
335: INFO = -2
336: ELSE IF( NR.LT.1 ) THEN
337: INFO = -3
338: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
339: INFO = -4
340: ELSE IF( LDGCOL.LT.N ) THEN
341: INFO = -22
342: ELSE IF( LDGNUM.LT.N ) THEN
343: INFO = -24
344: END IF
345: IF( INFO.NE.0 ) THEN
346: CALL XERBLA( 'DLASD7', -INFO )
347: RETURN
348: END IF
349: *
350: NLP1 = NL + 1
351: NLP2 = NL + 2
352: IF( ICOMPQ.EQ.1 ) THEN
353: GIVPTR = 0
354: END IF
355: *
356: * Generate the first part of the vector Z and move the singular
357: * values in the first part of D one position backward.
358: *
359: Z1 = ALPHA*VL( NLP1 )
360: VL( NLP1 ) = ZERO
361: TAU = VF( NLP1 )
362: DO 10 I = NL, 1, -1
363: Z( I+1 ) = ALPHA*VL( I )
364: VL( I ) = ZERO
365: VF( I+1 ) = VF( I )
366: D( I+1 ) = D( I )
367: IDXQ( I+1 ) = IDXQ( I ) + 1
368: 10 CONTINUE
369: VF( 1 ) = TAU
370: *
371: * Generate the second part of the vector Z.
372: *
373: DO 20 I = NLP2, M
374: Z( I ) = BETA*VF( I )
375: VF( I ) = ZERO
376: 20 CONTINUE
377: *
378: * Sort the singular values into increasing order
379: *
380: DO 30 I = NLP2, N
381: IDXQ( I ) = IDXQ( I ) + NLP1
382: 30 CONTINUE
383: *
384: * DSIGMA, IDXC, IDXC, and ZW are used as storage space.
385: *
386: DO 40 I = 2, N
387: DSIGMA( I ) = D( IDXQ( I ) )
388: ZW( I ) = Z( IDXQ( I ) )
389: VFW( I ) = VF( IDXQ( I ) )
390: VLW( I ) = VL( IDXQ( I ) )
391: 40 CONTINUE
392: *
393: CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
394: *
395: DO 50 I = 2, N
396: IDXI = 1 + IDX( I )
397: D( I ) = DSIGMA( IDXI )
398: Z( I ) = ZW( IDXI )
399: VF( I ) = VFW( IDXI )
400: VL( I ) = VLW( IDXI )
401: 50 CONTINUE
402: *
403: * Calculate the allowable deflation tolerence
404: *
405: EPS = DLAMCH( 'Epsilon' )
406: TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
407: TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
408: *
409: * There are 2 kinds of deflation -- first a value in the z-vector
410: * is small, second two (or more) singular values are very close
411: * together (their difference is small).
412: *
413: * If the value in the z-vector is small, we simply permute the
414: * array so that the corresponding singular value is moved to the
415: * end.
416: *
417: * If two values in the D-vector are close, we perform a two-sided
418: * rotation designed to make one of the corresponding z-vector
419: * entries zero, and then permute the array so that the deflated
420: * singular value is moved to the end.
421: *
422: * If there are multiple singular values then the problem deflates.
423: * Here the number of equal singular values are found. As each equal
424: * singular value is found, an elementary reflector is computed to
425: * rotate the corresponding singular subspace so that the
426: * corresponding components of Z are zero in this new basis.
427: *
428: K = 1
429: K2 = N + 1
430: DO 60 J = 2, N
431: IF( ABS( Z( J ) ).LE.TOL ) THEN
432: *
433: * Deflate due to small z component.
434: *
435: K2 = K2 - 1
436: IDXP( K2 ) = J
437: IF( J.EQ.N )
438: $ GO TO 100
439: ELSE
440: JPREV = J
441: GO TO 70
442: END IF
443: 60 CONTINUE
444: 70 CONTINUE
445: J = JPREV
446: 80 CONTINUE
447: J = J + 1
448: IF( J.GT.N )
449: $ GO TO 90
450: IF( ABS( Z( J ) ).LE.TOL ) THEN
451: *
452: * Deflate due to small z component.
453: *
454: K2 = K2 - 1
455: IDXP( K2 ) = J
456: ELSE
457: *
458: * Check if singular values are close enough to allow deflation.
459: *
460: IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
461: *
462: * Deflation is possible.
463: *
464: S = Z( JPREV )
465: C = Z( J )
466: *
467: * Find sqrt(a**2+b**2) without overflow or
468: * destructive underflow.
469: *
470: TAU = DLAPY2( C, S )
471: Z( J ) = TAU
472: Z( JPREV ) = ZERO
473: C = C / TAU
474: S = -S / TAU
475: *
476: * Record the appropriate Givens rotation
477: *
478: IF( ICOMPQ.EQ.1 ) THEN
479: GIVPTR = GIVPTR + 1
480: IDXJP = IDXQ( IDX( JPREV )+1 )
481: IDXJ = IDXQ( IDX( J )+1 )
482: IF( IDXJP.LE.NLP1 ) THEN
483: IDXJP = IDXJP - 1
484: END IF
485: IF( IDXJ.LE.NLP1 ) THEN
486: IDXJ = IDXJ - 1
487: END IF
488: GIVCOL( GIVPTR, 2 ) = IDXJP
489: GIVCOL( GIVPTR, 1 ) = IDXJ
490: GIVNUM( GIVPTR, 2 ) = C
491: GIVNUM( GIVPTR, 1 ) = S
492: END IF
493: CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
494: CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
495: K2 = K2 - 1
496: IDXP( K2 ) = JPREV
497: JPREV = J
498: ELSE
499: K = K + 1
500: ZW( K ) = Z( JPREV )
501: DSIGMA( K ) = D( JPREV )
502: IDXP( K ) = JPREV
503: JPREV = J
504: END IF
505: END IF
506: GO TO 80
507: 90 CONTINUE
508: *
509: * Record the last singular value.
510: *
511: K = K + 1
512: ZW( K ) = Z( JPREV )
513: DSIGMA( K ) = D( JPREV )
514: IDXP( K ) = JPREV
515: *
516: 100 CONTINUE
517: *
518: * Sort the singular values into DSIGMA. The singular values which
519: * were not deflated go into the first K slots of DSIGMA, except
520: * that DSIGMA(1) is treated separately.
521: *
522: DO 110 J = 2, N
523: JP = IDXP( J )
524: DSIGMA( J ) = D( JP )
525: VFW( J ) = VF( JP )
526: VLW( J ) = VL( JP )
527: 110 CONTINUE
528: IF( ICOMPQ.EQ.1 ) THEN
529: DO 120 J = 2, N
530: JP = IDXP( J )
531: PERM( J ) = IDXQ( IDX( JP )+1 )
532: IF( PERM( J ).LE.NLP1 ) THEN
533: PERM( J ) = PERM( J ) - 1
534: END IF
535: 120 CONTINUE
536: END IF
537: *
538: * The deflated singular values go back into the last N - K slots of
539: * D.
540: *
541: CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
542: *
543: * Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
544: * VL(M).
545: *
546: DSIGMA( 1 ) = ZERO
547: HLFTOL = TOL / TWO
548: IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
549: $ DSIGMA( 2 ) = HLFTOL
550: IF( M.GT.N ) THEN
551: Z( 1 ) = DLAPY2( Z1, Z( M ) )
552: IF( Z( 1 ).LE.TOL ) THEN
553: C = ONE
554: S = ZERO
555: Z( 1 ) = TOL
556: ELSE
557: C = Z1 / Z( 1 )
558: S = -Z( M ) / Z( 1 )
559: END IF
560: CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
561: CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
562: ELSE
563: IF( ABS( Z1 ).LE.TOL ) THEN
564: Z( 1 ) = TOL
565: ELSE
566: Z( 1 ) = Z1
567: END IF
568: END IF
569: *
570: * Restore Z, VF, and VL.
571: *
572: CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
573: CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
574: CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
575: *
576: RETURN
577: *
578: * End of DLASD7
579: *
580: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd7.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack