Annotation of rpl/lapack/lapack/dbdsdc.f, revision 1.10
1.10 ! bertrand 1: *> \brief \b DBDSDC
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DBDSDC + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsdc.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsdc.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsdc.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
! 22: * WORK, IWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER COMPQ, UPLO
! 26: * INTEGER INFO, LDU, LDVT, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER IQ( * ), IWORK( * )
! 30: * DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
! 31: * $ VT( LDVT, * ), WORK( * )
! 32: * ..
! 33: *
! 34: *
! 35: *> \par Purpose:
! 36: * =============
! 37: *>
! 38: *> \verbatim
! 39: *>
! 40: *> DBDSDC computes the singular value decomposition (SVD) of a real
! 41: *> N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
! 42: *> using a divide and conquer method, where S is a diagonal matrix
! 43: *> with non-negative diagonal elements (the singular values of B), and
! 44: *> U and VT are orthogonal matrices of left and right singular vectors,
! 45: *> respectively. DBDSDC can be used to compute all singular values,
! 46: *> and optionally, singular vectors or singular vectors in compact form.
! 47: *>
! 48: *> This code makes very mild assumptions about floating point
! 49: *> arithmetic. It will work on machines with a guard digit in
! 50: *> add/subtract, or on those binary machines without guard digits
! 51: *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
! 52: *> It could conceivably fail on hexadecimal or decimal machines
! 53: *> without guard digits, but we know of none. See DLASD3 for details.
! 54: *>
! 55: *> The code currently calls DLASDQ if singular values only are desired.
! 56: *> However, it can be slightly modified to compute singular values
! 57: *> using the divide and conquer method.
! 58: *> \endverbatim
! 59: *
! 60: * Arguments:
! 61: * ==========
! 62: *
! 63: *> \param[in] UPLO
! 64: *> \verbatim
! 65: *> UPLO is CHARACTER*1
! 66: *> = 'U': B is upper bidiagonal.
! 67: *> = 'L': B is lower bidiagonal.
! 68: *> \endverbatim
! 69: *>
! 70: *> \param[in] COMPQ
! 71: *> \verbatim
! 72: *> COMPQ is CHARACTER*1
! 73: *> Specifies whether singular vectors are to be computed
! 74: *> as follows:
! 75: *> = 'N': Compute singular values only;
! 76: *> = 'P': Compute singular values and compute singular
! 77: *> vectors in compact form;
! 78: *> = 'I': Compute singular values and singular vectors.
! 79: *> \endverbatim
! 80: *>
! 81: *> \param[in] N
! 82: *> \verbatim
! 83: *> N is INTEGER
! 84: *> The order of the matrix B. N >= 0.
! 85: *> \endverbatim
! 86: *>
! 87: *> \param[in,out] D
! 88: *> \verbatim
! 89: *> D is DOUBLE PRECISION array, dimension (N)
! 90: *> On entry, the n diagonal elements of the bidiagonal matrix B.
! 91: *> On exit, if INFO=0, the singular values of B.
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in,out] E
! 95: *> \verbatim
! 96: *> E is DOUBLE PRECISION array, dimension (N-1)
! 97: *> On entry, the elements of E contain the offdiagonal
! 98: *> elements of the bidiagonal matrix whose SVD is desired.
! 99: *> On exit, E has been destroyed.
! 100: *> \endverbatim
! 101: *>
! 102: *> \param[out] U
! 103: *> \verbatim
! 104: *> U is DOUBLE PRECISION array, dimension (LDU,N)
! 105: *> If COMPQ = 'I', then:
! 106: *> On exit, if INFO = 0, U contains the left singular vectors
! 107: *> of the bidiagonal matrix.
! 108: *> For other values of COMPQ, U is not referenced.
! 109: *> \endverbatim
! 110: *>
! 111: *> \param[in] LDU
! 112: *> \verbatim
! 113: *> LDU is INTEGER
! 114: *> The leading dimension of the array U. LDU >= 1.
! 115: *> If singular vectors are desired, then LDU >= max( 1, N ).
! 116: *> \endverbatim
! 117: *>
! 118: *> \param[out] VT
! 119: *> \verbatim
! 120: *> VT is DOUBLE PRECISION array, dimension (LDVT,N)
! 121: *> If COMPQ = 'I', then:
! 122: *> On exit, if INFO = 0, VT**T contains the right singular
! 123: *> vectors of the bidiagonal matrix.
! 124: *> For other values of COMPQ, VT is not referenced.
! 125: *> \endverbatim
! 126: *>
! 127: *> \param[in] LDVT
! 128: *> \verbatim
! 129: *> LDVT is INTEGER
! 130: *> The leading dimension of the array VT. LDVT >= 1.
! 131: *> If singular vectors are desired, then LDVT >= max( 1, N ).
! 132: *> \endverbatim
! 133: *>
! 134: *> \param[out] Q
! 135: *> \verbatim
! 136: *> Q is DOUBLE PRECISION array, dimension (LDQ)
! 137: *> If COMPQ = 'P', then:
! 138: *> On exit, if INFO = 0, Q and IQ contain the left
! 139: *> and right singular vectors in a compact form,
! 140: *> requiring O(N log N) space instead of 2*N**2.
! 141: *> In particular, Q contains all the DOUBLE PRECISION data in
! 142: *> LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
! 143: *> words of memory, where SMLSIZ is returned by ILAENV and
! 144: *> is equal to the maximum size of the subproblems at the
! 145: *> bottom of the computation tree (usually about 25).
! 146: *> For other values of COMPQ, Q is not referenced.
! 147: *> \endverbatim
! 148: *>
! 149: *> \param[out] IQ
! 150: *> \verbatim
! 151: *> IQ is INTEGER array, dimension (LDIQ)
! 152: *> If COMPQ = 'P', then:
! 153: *> On exit, if INFO = 0, Q and IQ contain the left
! 154: *> and right singular vectors in a compact form,
! 155: *> requiring O(N log N) space instead of 2*N**2.
! 156: *> In particular, IQ contains all INTEGER data in
! 157: *> LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
! 158: *> words of memory, where SMLSIZ is returned by ILAENV and
! 159: *> is equal to the maximum size of the subproblems at the
! 160: *> bottom of the computation tree (usually about 25).
! 161: *> For other values of COMPQ, IQ is not referenced.
! 162: *> \endverbatim
! 163: *>
! 164: *> \param[out] WORK
! 165: *> \verbatim
! 166: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 167: *> If COMPQ = 'N' then LWORK >= (4 * N).
! 168: *> If COMPQ = 'P' then LWORK >= (6 * N).
! 169: *> If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
! 170: *> \endverbatim
! 171: *>
! 172: *> \param[out] IWORK
! 173: *> \verbatim
! 174: *> IWORK is INTEGER array, dimension (8*N)
! 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: The algorithm failed to compute a singular value.
! 183: *> The update process of divide and conquer failed.
! 184: *> \endverbatim
! 185: *
! 186: * Authors:
! 187: * ========
! 188: *
! 189: *> \author Univ. of Tennessee
! 190: *> \author Univ. of California Berkeley
! 191: *> \author Univ. of Colorado Denver
! 192: *> \author NAG Ltd.
! 193: *
! 194: *> \date November 2011
! 195: *
! 196: *> \ingroup auxOTHERcomputational
! 197: *
! 198: *> \par Contributors:
! 199: * ==================
! 200: *>
! 201: *> Ming Gu and Huan Ren, Computer Science Division, University of
! 202: *> California at Berkeley, USA
! 203: *>
! 204: * =====================================================================
1.1 bertrand 205: SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
206: $ WORK, IWORK, INFO )
207: *
1.10 ! bertrand 208: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 209: * -- LAPACK is a software package provided by Univ. of Tennessee, --
210: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.10 ! bertrand 211: * November 2011
1.1 bertrand 212: *
213: * .. Scalar Arguments ..
214: CHARACTER COMPQ, UPLO
215: INTEGER INFO, LDU, LDVT, N
216: * ..
217: * .. Array Arguments ..
218: INTEGER IQ( * ), IWORK( * )
219: DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
220: $ VT( LDVT, * ), WORK( * )
221: * ..
222: *
223: * =====================================================================
224: * Changed dimension statement in comment describing E from (N) to
225: * (N-1). Sven, 17 Feb 05.
226: * =====================================================================
227: *
228: * .. Parameters ..
229: DOUBLE PRECISION ZERO, ONE, TWO
230: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
231: * ..
232: * .. Local Scalars ..
233: INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
234: $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
235: $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
236: $ SMLSZP, SQRE, START, WSTART, Z
237: DOUBLE PRECISION CS, EPS, ORGNRM, P, R, SN
238: * ..
239: * .. External Functions ..
240: LOGICAL LSAME
241: INTEGER ILAENV
242: DOUBLE PRECISION DLAMCH, DLANST
243: EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
244: * ..
245: * .. External Subroutines ..
246: EXTERNAL DCOPY, DLARTG, DLASCL, DLASD0, DLASDA, DLASDQ,
247: $ DLASET, DLASR, DSWAP, XERBLA
248: * ..
249: * .. Intrinsic Functions ..
250: INTRINSIC ABS, DBLE, INT, LOG, SIGN
251: * ..
252: * .. Executable Statements ..
253: *
254: * Test the input parameters.
255: *
256: INFO = 0
257: *
258: IUPLO = 0
259: IF( LSAME( UPLO, 'U' ) )
260: $ IUPLO = 1
261: IF( LSAME( UPLO, 'L' ) )
262: $ IUPLO = 2
263: IF( LSAME( COMPQ, 'N' ) ) THEN
264: ICOMPQ = 0
265: ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
266: ICOMPQ = 1
267: ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
268: ICOMPQ = 2
269: ELSE
270: ICOMPQ = -1
271: END IF
272: IF( IUPLO.EQ.0 ) THEN
273: INFO = -1
274: ELSE IF( ICOMPQ.LT.0 ) THEN
275: INFO = -2
276: ELSE IF( N.LT.0 ) THEN
277: INFO = -3
278: ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
279: $ N ) ) ) THEN
280: INFO = -7
281: ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
282: $ N ) ) ) THEN
283: INFO = -9
284: END IF
285: IF( INFO.NE.0 ) THEN
286: CALL XERBLA( 'DBDSDC', -INFO )
287: RETURN
288: END IF
289: *
290: * Quick return if possible
291: *
292: IF( N.EQ.0 )
293: $ RETURN
294: SMLSIZ = ILAENV( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
295: IF( N.EQ.1 ) THEN
296: IF( ICOMPQ.EQ.1 ) THEN
297: Q( 1 ) = SIGN( ONE, D( 1 ) )
298: Q( 1+SMLSIZ*N ) = ONE
299: ELSE IF( ICOMPQ.EQ.2 ) THEN
300: U( 1, 1 ) = SIGN( ONE, D( 1 ) )
301: VT( 1, 1 ) = ONE
302: END IF
303: D( 1 ) = ABS( D( 1 ) )
304: RETURN
305: END IF
306: NM1 = N - 1
307: *
308: * If matrix lower bidiagonal, rotate to be upper bidiagonal
309: * by applying Givens rotations on the left
310: *
311: WSTART = 1
312: QSTART = 3
313: IF( ICOMPQ.EQ.1 ) THEN
314: CALL DCOPY( N, D, 1, Q( 1 ), 1 )
315: CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
316: END IF
317: IF( IUPLO.EQ.2 ) THEN
318: QSTART = 5
319: WSTART = 2*N - 1
320: DO 10 I = 1, N - 1
321: CALL DLARTG( D( I ), E( I ), CS, SN, R )
322: D( I ) = R
323: E( I ) = SN*D( I+1 )
324: D( I+1 ) = CS*D( I+1 )
325: IF( ICOMPQ.EQ.1 ) THEN
326: Q( I+2*N ) = CS
327: Q( I+3*N ) = SN
328: ELSE IF( ICOMPQ.EQ.2 ) THEN
329: WORK( I ) = CS
330: WORK( NM1+I ) = -SN
331: END IF
332: 10 CONTINUE
333: END IF
334: *
335: * If ICOMPQ = 0, use DLASDQ to compute the singular values.
336: *
337: IF( ICOMPQ.EQ.0 ) THEN
338: CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
339: $ LDU, WORK( WSTART ), INFO )
340: GO TO 40
341: END IF
342: *
343: * If N is smaller than the minimum divide size SMLSIZ, then solve
344: * the problem with another solver.
345: *
346: IF( N.LE.SMLSIZ ) THEN
347: IF( ICOMPQ.EQ.2 ) THEN
348: CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
349: CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
350: CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
351: $ LDU, WORK( WSTART ), INFO )
352: ELSE IF( ICOMPQ.EQ.1 ) THEN
353: IU = 1
354: IVT = IU + N
355: CALL DLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
356: $ N )
357: CALL DLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
358: $ N )
359: CALL DLASDQ( 'U', 0, N, N, N, 0, D, E,
360: $ Q( IVT+( QSTART-1 )*N ), N,
361: $ Q( IU+( QSTART-1 )*N ), N,
362: $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
363: $ INFO )
364: END IF
365: GO TO 40
366: END IF
367: *
368: IF( ICOMPQ.EQ.2 ) THEN
369: CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
370: CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
371: END IF
372: *
373: * Scale.
374: *
375: ORGNRM = DLANST( 'M', N, D, E )
376: IF( ORGNRM.EQ.ZERO )
377: $ RETURN
378: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
379: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
380: *
1.9 bertrand 381: EPS = (0.9D+0)*DLAMCH( 'Epsilon' )
1.1 bertrand 382: *
383: MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
384: SMLSZP = SMLSIZ + 1
385: *
386: IF( ICOMPQ.EQ.1 ) THEN
387: IU = 1
388: IVT = 1 + SMLSIZ
389: DIFL = IVT + SMLSZP
390: DIFR = DIFL + MLVL
391: Z = DIFR + MLVL*2
392: IC = Z + MLVL
393: IS = IC + 1
394: POLES = IS + 1
395: GIVNUM = POLES + 2*MLVL
396: *
397: K = 1
398: GIVPTR = 2
399: PERM = 3
400: GIVCOL = PERM + MLVL
401: END IF
402: *
403: DO 20 I = 1, N
404: IF( ABS( D( I ) ).LT.EPS ) THEN
405: D( I ) = SIGN( EPS, D( I ) )
406: END IF
407: 20 CONTINUE
408: *
409: START = 1
410: SQRE = 0
411: *
412: DO 30 I = 1, NM1
413: IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
414: *
415: * Subproblem found. First determine its size and then
416: * apply divide and conquer on it.
417: *
418: IF( I.LT.NM1 ) THEN
419: *
420: * A subproblem with E(I) small for I < NM1.
421: *
422: NSIZE = I - START + 1
423: ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
424: *
425: * A subproblem with E(NM1) not too small but I = NM1.
426: *
427: NSIZE = N - START + 1
428: ELSE
429: *
430: * A subproblem with E(NM1) small. This implies an
431: * 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
432: * first.
433: *
434: NSIZE = I - START + 1
435: IF( ICOMPQ.EQ.2 ) THEN
436: U( N, N ) = SIGN( ONE, D( N ) )
437: VT( N, N ) = ONE
438: ELSE IF( ICOMPQ.EQ.1 ) THEN
439: Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
440: Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
441: END IF
442: D( N ) = ABS( D( N ) )
443: END IF
444: IF( ICOMPQ.EQ.2 ) THEN
445: CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
446: $ U( START, START ), LDU, VT( START, START ),
447: $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
448: ELSE
449: CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
450: $ E( START ), Q( START+( IU+QSTART-2 )*N ), N,
451: $ Q( START+( IVT+QSTART-2 )*N ),
452: $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
453: $ N ), Q( START+( DIFR+QSTART-2 )*N ),
454: $ Q( START+( Z+QSTART-2 )*N ),
455: $ Q( START+( POLES+QSTART-2 )*N ),
456: $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
457: $ N, IQ( START+PERM*N ),
458: $ Q( START+( GIVNUM+QSTART-2 )*N ),
459: $ Q( START+( IC+QSTART-2 )*N ),
460: $ Q( START+( IS+QSTART-2 )*N ),
461: $ WORK( WSTART ), IWORK, INFO )
1.5 bertrand 462: END IF
463: IF( INFO.NE.0 ) THEN
464: RETURN
1.1 bertrand 465: END IF
466: START = I + 1
467: END IF
468: 30 CONTINUE
469: *
470: * Unscale
471: *
472: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
473: 40 CONTINUE
474: *
475: * Use Selection Sort to minimize swaps of singular vectors
476: *
477: DO 60 II = 2, N
478: I = II - 1
479: KK = I
480: P = D( I )
481: DO 50 J = II, N
482: IF( D( J ).GT.P ) THEN
483: KK = J
484: P = D( J )
485: END IF
486: 50 CONTINUE
487: IF( KK.NE.I ) THEN
488: D( KK ) = D( I )
489: D( I ) = P
490: IF( ICOMPQ.EQ.1 ) THEN
491: IQ( I ) = KK
492: ELSE IF( ICOMPQ.EQ.2 ) THEN
493: CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
494: CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
495: END IF
496: ELSE IF( ICOMPQ.EQ.1 ) THEN
497: IQ( I ) = I
498: END IF
499: 60 CONTINUE
500: *
501: * If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
502: *
503: IF( ICOMPQ.EQ.1 ) THEN
504: IF( IUPLO.EQ.1 ) THEN
505: IQ( N ) = 1
506: ELSE
507: IQ( N ) = 0
508: END IF
509: END IF
510: *
511: * If B is lower bidiagonal, update U by those Givens rotations
512: * which rotated B to be upper bidiagonal
513: *
514: IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
515: $ CALL DLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )
516: *
517: RETURN
518: *
519: * End of DBDSDC
520: *
521: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dbdsdc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack