Annotation of rpl/lapack/lapack/dlasd1.f, revision 1.1
1.1 ! bertrand 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) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 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, an 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack