Annotation of rpl/lapack/lapack/dlasd1.f, revision 1.10
1.10 ! bertrand 1: *> \brief \b DLASD1
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLASD1 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd1.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd1.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd1.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
! 22: * IDXQ, IWORK, WORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * INTEGER INFO, LDU, LDVT, NL, NR, SQRE
! 26: * DOUBLE PRECISION ALPHA, BETA
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER IDXQ( * ), IWORK( * )
! 30: * DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
! 40: *> where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
! 41: *>
! 42: *> A related subroutine DLASD7 handles the case in which the singular
! 43: *> values (and the singular vectors in factored form) are desired.
! 44: *>
! 45: *> DLASD1 computes the SVD as follows:
! 46: *>
! 47: *> ( D1(in) 0 0 0 )
! 48: *> B = U(in) * ( Z1**T a Z2**T b ) * VT(in)
! 49: *> ( 0 0 D2(in) 0 )
! 50: *>
! 51: *> = U(out) * ( D(out) 0) * VT(out)
! 52: *>
! 53: *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
! 54: *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
! 55: *> elsewhere; and the entry b is empty if SQRE = 0.
! 56: *>
! 57: *> The left singular vectors of the original matrix are stored in U, and
! 58: *> the transpose of the right singular vectors are stored in VT, and the
! 59: *> singular values are in D. The algorithm consists of three stages:
! 60: *>
! 61: *> The first stage consists of deflating the size of the problem
! 62: *> when there are multiple singular values or when there are zeros in
! 63: *> the Z vector. For each such occurence the dimension of the
! 64: *> secular equation problem is reduced by one. This stage is
! 65: *> performed by the routine DLASD2.
! 66: *>
! 67: *> The second stage consists of calculating the updated
! 68: *> singular values. This is done by finding the square roots of the
! 69: *> roots of the secular equation via the routine DLASD4 (as called
! 70: *> by DLASD3). This routine also calculates the singular vectors of
! 71: *> the current problem.
! 72: *>
! 73: *> The final stage consists of computing the updated singular vectors
! 74: *> directly using the updated singular values. The singular vectors
! 75: *> for the current problem are multiplied with the singular vectors
! 76: *> from the overall problem.
! 77: *> \endverbatim
! 78: *
! 79: * Arguments:
! 80: * ==========
! 81: *
! 82: *> \param[in] NL
! 83: *> \verbatim
! 84: *> NL is INTEGER
! 85: *> The row dimension of the upper block. NL >= 1.
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[in] NR
! 89: *> \verbatim
! 90: *> NR is INTEGER
! 91: *> The row dimension of the lower block. NR >= 1.
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in] SQRE
! 95: *> \verbatim
! 96: *> SQRE is INTEGER
! 97: *> = 0: the lower block is an NR-by-NR square matrix.
! 98: *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
! 99: *>
! 100: *> The bidiagonal matrix has row dimension N = NL + NR + 1,
! 101: *> and column dimension M = N + SQRE.
! 102: *> \endverbatim
! 103: *>
! 104: *> \param[in,out] D
! 105: *> \verbatim
! 106: *> D is DOUBLE PRECISION array,
! 107: *> dimension (N = NL+NR+1).
! 108: *> On entry D(1:NL,1:NL) contains the singular values of the
! 109: *> upper block; and D(NL+2:N) contains the singular values of
! 110: *> the lower block. On exit D(1:N) contains the singular values
! 111: *> of the modified matrix.
! 112: *> \endverbatim
! 113: *>
! 114: *> \param[in,out] ALPHA
! 115: *> \verbatim
! 116: *> ALPHA is DOUBLE PRECISION
! 117: *> Contains the diagonal element associated with the added row.
! 118: *> \endverbatim
! 119: *>
! 120: *> \param[in,out] BETA
! 121: *> \verbatim
! 122: *> BETA is DOUBLE PRECISION
! 123: *> Contains the off-diagonal element associated with the added
! 124: *> row.
! 125: *> \endverbatim
! 126: *>
! 127: *> \param[in,out] U
! 128: *> \verbatim
! 129: *> U is DOUBLE PRECISION array, dimension(LDU,N)
! 130: *> On entry U(1:NL, 1:NL) contains the left singular vectors of
! 131: *> the upper block; U(NL+2:N, NL+2:N) contains the left singular
! 132: *> vectors of the lower block. On exit U contains the left
! 133: *> singular vectors of the bidiagonal matrix.
! 134: *> \endverbatim
! 135: *>
! 136: *> \param[in] LDU
! 137: *> \verbatim
! 138: *> LDU is INTEGER
! 139: *> The leading dimension of the array U. LDU >= max( 1, N ).
! 140: *> \endverbatim
! 141: *>
! 142: *> \param[in,out] VT
! 143: *> \verbatim
! 144: *> VT is DOUBLE PRECISION array, dimension(LDVT,M)
! 145: *> where M = N + SQRE.
! 146: *> On entry VT(1:NL+1, 1:NL+1)**T contains the right singular
! 147: *> vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains
! 148: *> the right singular vectors of the lower block. On exit
! 149: *> VT**T contains the right singular vectors of the
! 150: *> bidiagonal matrix.
! 151: *> \endverbatim
! 152: *>
! 153: *> \param[in] LDVT
! 154: *> \verbatim
! 155: *> LDVT is INTEGER
! 156: *> The leading dimension of the array VT. LDVT >= max( 1, M ).
! 157: *> \endverbatim
! 158: *>
! 159: *> \param[out] IDXQ
! 160: *> \verbatim
! 161: *> IDXQ is INTEGER array, dimension(N)
! 162: *> This contains the permutation which will reintegrate the
! 163: *> subproblem just solved back into sorted order, i.e.
! 164: *> D( IDXQ( I = 1, N ) ) will be in ascending order.
! 165: *> \endverbatim
! 166: *>
! 167: *> \param[out] IWORK
! 168: *> \verbatim
! 169: *> IWORK is INTEGER array, dimension( 4 * N )
! 170: *> \endverbatim
! 171: *>
! 172: *> \param[out] WORK
! 173: *> \verbatim
! 174: *> WORK is DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
! 175: *> \endverbatim
! 176: *>
! 177: *> \param[out] INFO
! 178: *> \verbatim
! 179: *> INFO is INTEGER
! 180: *> = 0: successful exit.
! 181: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 182: *> > 0: if INFO = 1, a singular value did not converge
! 183: *> \endverbatim
! 184: *
! 185: * Authors:
! 186: * ========
! 187: *
! 188: *> \author Univ. of Tennessee
! 189: *> \author Univ. of California Berkeley
! 190: *> \author Univ. of Colorado Denver
! 191: *> \author NAG Ltd.
! 192: *
! 193: *> \date November 2011
! 194: *
! 195: *> \ingroup auxOTHERauxiliary
! 196: *
! 197: *> \par Contributors:
! 198: * ==================
! 199: *>
! 200: *> Ming Gu and Huan Ren, Computer Science Division, University of
! 201: *> California at Berkeley, USA
! 202: *>
! 203: * =====================================================================
1.1 bertrand 204: SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
205: $ IDXQ, IWORK, WORK, INFO )
206: *
1.10 ! bertrand 207: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 208: * -- LAPACK is a software package provided by Univ. of Tennessee, --
209: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.10 ! bertrand 210: * November 2011
1.1 bertrand 211: *
212: * .. Scalar Arguments ..
213: INTEGER INFO, LDU, LDVT, NL, NR, SQRE
214: DOUBLE PRECISION ALPHA, BETA
215: * ..
216: * .. Array Arguments ..
217: INTEGER IDXQ( * ), IWORK( * )
218: DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
219: * ..
220: *
221: * =====================================================================
222: *
223: * .. Parameters ..
224: *
225: DOUBLE PRECISION ONE, ZERO
226: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
227: * ..
228: * .. Local Scalars ..
229: INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
230: $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
231: DOUBLE PRECISION ORGNRM
232: * ..
233: * .. External Subroutines ..
234: EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
235: * ..
236: * .. Intrinsic Functions ..
237: INTRINSIC ABS, MAX
238: * ..
239: * .. Executable Statements ..
240: *
241: * Test the input parameters.
242: *
243: INFO = 0
244: *
245: IF( NL.LT.1 ) THEN
246: INFO = -1
247: ELSE IF( NR.LT.1 ) THEN
248: INFO = -2
249: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
250: INFO = -3
251: END IF
252: IF( INFO.NE.0 ) THEN
253: CALL XERBLA( 'DLASD1', -INFO )
254: RETURN
255: END IF
256: *
257: N = NL + NR + 1
258: M = N + SQRE
259: *
260: * The following values are for bookkeeping purposes only. They are
261: * integer pointers which indicate the portion of the workspace
262: * used by a particular array in DLASD2 and DLASD3.
263: *
264: LDU2 = N
265: LDVT2 = M
266: *
267: IZ = 1
268: ISIGMA = IZ + M
269: IU2 = ISIGMA + N
270: IVT2 = IU2 + LDU2*N
271: IQ = IVT2 + LDVT2*M
272: *
273: IDX = 1
274: IDXC = IDX + N
275: COLTYP = IDXC + N
276: IDXP = COLTYP + N
277: *
278: * Scale.
279: *
280: ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
281: D( NL+1 ) = ZERO
282: DO 10 I = 1, N
283: IF( ABS( D( I ) ).GT.ORGNRM ) THEN
284: ORGNRM = ABS( D( I ) )
285: END IF
286: 10 CONTINUE
287: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
288: ALPHA = ALPHA / ORGNRM
289: BETA = BETA / ORGNRM
290: *
291: * Deflate singular values.
292: *
293: CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
294: $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
295: $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
296: $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
297: *
298: * Solve Secular Equation and update singular vectors.
299: *
300: LDQ = K
301: CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
302: $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
303: $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
304: $ INFO )
305: IF( INFO.NE.0 ) THEN
306: RETURN
307: END IF
308: *
309: * Unscale.
310: *
311: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
312: *
313: * Prepare the IDXQ sorting permutation.
314: *
315: N1 = K
316: N2 = N - K
317: CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
318: *
319: RETURN
320: *
321: * End of DLASD1
322: *
323: END
CVSweb interface <joel.bertrand@systella.fr>