Annotation of rpl/lapack/lapack/dbdsvdx.f, revision 1.3
1.1 bertrand 1: *> \brief \b DBDSVDX
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DBDSVDX + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsvdx.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsvdx.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsvdx.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
22: * $ NS, S, Z, LDZ, WORK, IWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER JOBZ, RANGE, UPLO
26: * INTEGER IL, INFO, IU, LDZ, N, NS
27: * DOUBLE PRECISION VL, VU
28: * ..
29: * .. Array Arguments ..
30: * INTEGER IWORK( * )
31: * DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ),
32: * Z( LDZ, * )
33: * ..
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DBDSVDX computes the singular value decomposition (SVD) of a real
41: *> N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT,
42: *> where S is a diagonal matrix with non-negative diagonal elements
43: *> (the singular values of B), and U and VT are orthogonal matrices
44: *> of left and right singular vectors, respectively.
45: *>
46: *> Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ]
47: *> and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the
48: *> singular value decompositon of B through the eigenvalues and
49: *> eigenvectors of the N*2-by-N*2 tridiagonal matrix
50: *>
51: *> | 0 d_1 |
52: *> | d_1 0 e_1 |
53: *> TGK = | e_1 0 d_2 |
54: *> | d_2 . . |
55: *> | . . . |
56: *>
57: *> If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then
58: *> (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) /
59: *> sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and
60: *> P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ].
61: *>
62: *> Given a TGK matrix, one can either a) compute -s,-v and change signs
63: *> so that the singular values (and corresponding vectors) are already in
64: *> descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder
65: *> the values (and corresponding vectors). DBDSVDX implements a) by
66: *> calling DSTEVX (bisection plus inverse iteration, to be replaced
67: *> with a version of the Multiple Relative Robust Representation
68: *> algorithm. (See P. Willems and B. Lang, A framework for the MR^3
69: *> algorithm: theory and implementation, SIAM J. Sci. Comput.,
70: *> 35:740-766, 2013.)
71: *> \endverbatim
72: *
73: * Arguments:
74: * ==========
75: *
76: *> \param[in] UPLO
77: *> \verbatim
78: *> UPLO is CHARACTER*1
79: *> = 'U': B is upper bidiagonal;
80: *> = 'L': B is lower bidiagonal.
81: *> \endverbatim
82: *>
1.2 bertrand 83: *> \param[in] JOBZ
1.1 bertrand 84: *> \verbatim
85: *> JOBZ is CHARACTER*1
86: *> = 'N': Compute singular values only;
87: *> = 'V': Compute singular values and singular vectors.
88: *> \endverbatim
89: *>
90: *> \param[in] RANGE
91: *> \verbatim
92: *> RANGE is CHARACTER*1
93: *> = 'A': all singular values will be found.
94: *> = 'V': all singular values in the half-open interval [VL,VU)
95: *> will be found.
96: *> = 'I': the IL-th through IU-th singular values will be found.
97: *> \endverbatim
98: *>
99: *> \param[in] N
100: *> \verbatim
101: *> N is INTEGER
102: *> The order of the bidiagonal matrix. N >= 0.
103: *> \endverbatim
104: *>
105: *> \param[in] D
106: *> \verbatim
107: *> D is DOUBLE PRECISION array, dimension (N)
108: *> The n diagonal elements of the bidiagonal matrix B.
109: *> \endverbatim
110: *>
111: *> \param[in] E
112: *> \verbatim
113: *> E is DOUBLE PRECISION array, dimension (max(1,N-1))
114: *> The (n-1) superdiagonal elements of the bidiagonal matrix
115: *> B in elements 1 to N-1.
116: *> \endverbatim
117: *>
118: *> \param[in] VL
119: *> \verbatim
1.2 bertrand 120: *> VL is DOUBLE PRECISION
121: *> If RANGE='V', the lower bound of the interval to
122: *> be searched for singular values. VU > VL.
123: *> Not referenced if RANGE = 'A' or 'I'.
1.1 bertrand 124: *> \endverbatim
125: *>
126: *> \param[in] VU
127: *> \verbatim
128: *> VU is DOUBLE PRECISION
1.2 bertrand 129: *> If RANGE='V', the upper bound of the interval to
1.1 bertrand 130: *> be searched for singular values. VU > VL.
131: *> Not referenced if RANGE = 'A' or 'I'.
132: *> \endverbatim
133: *>
134: *> \param[in] IL
135: *> \verbatim
136: *> IL is INTEGER
1.2 bertrand 137: *> If RANGE='I', the index of the
138: *> smallest singular value to be returned.
139: *> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
140: *> Not referenced if RANGE = 'A' or 'V'.
1.1 bertrand 141: *> \endverbatim
142: *>
143: *> \param[in] IU
144: *> \verbatim
145: *> IU is INTEGER
1.2 bertrand 146: *> If RANGE='I', the index of the
147: *> largest singular value to be returned.
1.1 bertrand 148: *> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
149: *> Not referenced if RANGE = 'A' or 'V'.
150: *> \endverbatim
151: *>
152: *> \param[out] NS
153: *> \verbatim
154: *> NS is INTEGER
155: *> The total number of singular values found. 0 <= NS <= N.
156: *> If RANGE = 'A', NS = N, and if RANGE = 'I', NS = IU-IL+1.
157: *> \endverbatim
158: *>
159: *> \param[out] S
160: *> \verbatim
161: *> S is DOUBLE PRECISION array, dimension (N)
162: *> The first NS elements contain the selected singular values in
163: *> ascending order.
164: *> \endverbatim
165: *>
166: *> \param[out] Z
167: *> \verbatim
168: *> Z is DOUBLE PRECISION array, dimension (2*N,K) )
169: *> If JOBZ = 'V', then if INFO = 0 the first NS columns of Z
170: *> contain the singular vectors of the matrix B corresponding to
171: *> the selected singular values, with U in rows 1 to N and V
172: *> in rows N+1 to N*2, i.e.
173: *> Z = [ U ]
174: *> [ V ]
175: *> If JOBZ = 'N', then Z is not referenced.
176: *> Note: The user must ensure that at least K = NS+1 columns are
177: *> supplied in the array Z; if RANGE = 'V', the exact value of
178: *> NS is not known in advance and an upper bound must be used.
179: *> \endverbatim
180: *>
181: *> \param[in] LDZ
182: *> \verbatim
183: *> LDZ is INTEGER
184: *> The leading dimension of the array Z. LDZ >= 1, and if
185: *> JOBZ = 'V', LDZ >= max(2,N*2).
186: *> \endverbatim
187: *>
188: *> \param[out] WORK
189: *> \verbatim
190: *> WORK is DOUBLE PRECISION array, dimension (14*N)
191: *> \endverbatim
192: *>
193: *> \param[out] IWORK
194: *> \verbatim
195: *> IWORK is INTEGER array, dimension (12*N)
196: *> If JOBZ = 'V', then if INFO = 0, the first NS elements of
197: *> IWORK are zero. If INFO > 0, then IWORK contains the indices
198: *> of the eigenvectors that failed to converge in DSTEVX.
1.2 bertrand 199: *> \endverbatim
1.1 bertrand 200: *>
1.2 bertrand 201: *> \param[out] INFO
202: *> \verbatim
1.1 bertrand 203: *> INFO is INTEGER
204: *> = 0: successful exit
205: *> < 0: if INFO = -i, the i-th argument had an illegal value
206: *> > 0: if INFO = i, then i eigenvectors failed to converge
207: *> in DSTEVX. The indices of the eigenvectors
208: *> (as returned by DSTEVX) are stored in the
209: *> array IWORK.
210: *> if INFO = N*2 + 1, an internal error occurred.
211: *> \endverbatim
212: *
213: * Authors:
214: * ========
215: *
216: *> \author Univ. of Tennessee
217: *> \author Univ. of California Berkeley
218: *> \author Univ. of Colorado Denver
219: *> \author NAG Ltd.
220: *
1.2 bertrand 221: *> \date June 2016
1.1 bertrand 222: *
223: *> \ingroup doubleOTHEReigen
224: *
225: * =====================================================================
226: SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
227: $ NS, S, Z, LDZ, WORK, IWORK, INFO)
228: *
1.2 bertrand 229: * -- LAPACK driver routine (version 3.6.1) --
1.1 bertrand 230: * -- LAPACK is a software package provided by Univ. of Tennessee, --
231: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232: * November 2016
233: *
234: * .. Scalar Arguments ..
235: CHARACTER JOBZ, RANGE, UPLO
236: INTEGER IL, INFO, IU, LDZ, N, NS
237: DOUBLE PRECISION VL, VU
238: * ..
239: * .. Array Arguments ..
240: INTEGER IWORK( * )
241: DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ),
242: $ Z( LDZ, * )
243: * ..
244: *
245: * =====================================================================
246: *
247: * .. Parameters ..
248: DOUBLE PRECISION ZERO, ONE, TEN, HNDRD, MEIGTH
249: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0,
250: $ HNDRD = 100.0D0, MEIGTH = -0.1250D0 )
251: DOUBLE PRECISION FUDGE
252: PARAMETER ( FUDGE = 2.0D0 )
253: * ..
254: * .. Local Scalars ..
255: CHARACTER RNGVX
256: LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ
257: INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR,
258: $ IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV,
259: $ IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K,
260: $ NTGK, NRU, NRV, NSL
261: DOUBLE PRECISION ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX,
262: $ SMIN, SQRT2, THRESH, TOL, ULP,
263: $ VLTGK, VUTGK, ZJTJI
264: * ..
265: * .. External Functions ..
266: LOGICAL LSAME
267: INTEGER IDAMAX
268: DOUBLE PRECISION DDOT, DLAMCH, DNRM2
269: EXTERNAL IDAMAX, LSAME, DAXPY, DDOT, DLAMCH, DNRM2
270: * ..
271: * .. External Subroutines ..
272: EXTERNAL DCOPY, DLASET, DSCAL, DSWAP
273: * ..
274: * .. Intrinsic Functions ..
275: INTRINSIC ABS, DBLE, SIGN, SQRT
276: * ..
277: * .. Executable Statements ..
278: *
279: * Test the input parameters.
280: *
281: ALLSV = LSAME( RANGE, 'A' )
282: VALSV = LSAME( RANGE, 'V' )
283: INDSV = LSAME( RANGE, 'I' )
284: WANTZ = LSAME( JOBZ, 'V' )
285: LOWER = LSAME( UPLO, 'L' )
286: *
287: INFO = 0
288: IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
289: INFO = -1
290: ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
291: INFO = -2
292: ELSE IF( .NOT.( ALLSV .OR. VALSV .OR. INDSV ) ) THEN
293: INFO = -3
294: ELSE IF( N.LT.0 ) THEN
295: INFO = -4
296: ELSE IF( N.GT.0 ) THEN
297: IF( VALSV ) THEN
298: IF( VL.LT.ZERO ) THEN
299: INFO = -7
300: ELSE IF( VU.LE.VL ) THEN
301: INFO = -8
302: END IF
303: ELSE IF( INDSV ) THEN
304: IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
305: INFO = -9
306: ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
307: INFO = -10
308: END IF
309: END IF
310: END IF
311: IF( INFO.EQ.0 ) THEN
312: IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N*2 ) ) INFO = -14
313: END IF
314: *
315: IF( INFO.NE.0 ) THEN
316: CALL XERBLA( 'DBDSVDX', -INFO )
317: RETURN
318: END IF
319: *
320: * Quick return if possible (N.LE.1)
321: *
322: NS = 0
323: IF( N.EQ.0 ) RETURN
324: *
325: IF( N.EQ.1 ) THEN
326: IF( ALLSV .OR. INDSV ) THEN
327: NS = 1
328: S( 1 ) = ABS( D( 1 ) )
329: ELSE
330: IF( VL.LT.ABS( D( 1 ) ) .AND. VU.GE.ABS( D( 1 ) ) ) THEN
331: NS = 1
332: S( 1 ) = ABS( D( 1 ) )
333: END IF
334: END IF
335: IF( WANTZ ) THEN
336: Z( 1, 1 ) = SIGN( ONE, D( 1 ) )
337: Z( 2, 1 ) = ONE
338: END IF
339: RETURN
340: END IF
341: *
342: ABSTOL = 2*DLAMCH( 'Safe Minimum' )
343: ULP = DLAMCH( 'Precision' )
344: EPS = DLAMCH( 'Epsilon' )
345: SQRT2 = SQRT( 2.0D0 )
346: ORTOL = SQRT( ULP )
347: *
348: * Criterion for splitting is taken from DBDSQR when singular
349: * values are computed to relative accuracy TOL. (See J. Demmel and
350: * W. Kahan, Accurate singular values of bidiagonal matrices, SIAM
351: * J. Sci. and Stat. Comput., 11:873–912, 1990.)
352: *
353: TOL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )*EPS
354: *
355: * Compute approximate maximum, minimum singular values.
356: *
357: I = IDAMAX( N, D, 1 )
358: SMAX = ABS( D( I ) )
359: I = IDAMAX( N-1, E, 1 )
360: SMAX = MAX( SMAX, ABS( E( I ) ) )
361: *
362: * Compute threshold for neglecting D's and E's.
363: *
364: SMIN = ABS( D( 1 ) )
365: IF( SMIN.NE.ZERO ) THEN
366: MU = SMIN
367: DO I = 2, N
368: MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
369: SMIN = MIN( SMIN, MU )
370: IF( SMIN.EQ.ZERO ) EXIT
371: END DO
372: END IF
373: SMIN = SMIN / SQRT( DBLE( N ) )
374: THRESH = TOL*SMIN
375: *
376: * Check for zeros in D and E (splits), i.e. submatrices.
377: *
378: DO I = 1, N-1
379: IF( ABS( D( I ) ).LE.THRESH ) D( I ) = ZERO
380: IF( ABS( E( I ) ).LE.THRESH ) E( I ) = ZERO
381: END DO
382: IF( ABS( D( N ) ).LE.THRESH ) D( N ) = ZERO
383: *
384: * Pointers for arrays used by DSTEVX.
385: *
386: IDTGK = 1
387: IETGK = IDTGK + N*2
388: ITEMP = IETGK + N*2
389: IIFAIL = 1
390: IIWORK = IIFAIL + N*2
391: *
392: * Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode.
393: * VL,VU or IL,IU are redefined to conform to implementation a)
394: * described in the leading comments.
395: *
396: ILTGK = 0
397: IUTGK = 0
398: VLTGK = ZERO
399: VUTGK = ZERO
400: *
401: IF( ALLSV ) THEN
402: *
403: * All singular values will be found. We aim at -s (see
404: * leading comments) with RNGVX = 'I'. IL and IU are set
405: * later (as ILTGK and IUTGK) according to the dimension
406: * of the active submatrix.
407: *
408: RNGVX = 'I'
1.2 bertrand 409: IF( WANTZ ) CALL DLASET( 'F', N*2, N+1, ZERO, ZERO, Z, LDZ )
1.1 bertrand 410: ELSE IF( VALSV ) THEN
411: *
412: * Find singular values in a half-open interval. We aim
413: * at -s (see leading comments) and we swap VL and VU
414: * (as VUTGK and VLTGK), changing their signs.
415: *
416: RNGVX = 'V'
417: VLTGK = -VU
418: VUTGK = -VL
419: WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
420: CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
421: CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
422: CALL DSTEVX( 'N', 'V', N*2, WORK( IDTGK ), WORK( IETGK ),
423: $ VLTGK, VUTGK, ILTGK, ILTGK, ABSTOL, NS, S,
424: $ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),
425: $ IWORK( IIFAIL ), INFO )
426: IF( NS.EQ.0 ) THEN
427: RETURN
428: ELSE
1.2 bertrand 429: IF( WANTZ ) CALL DLASET( 'F', N*2, NS, ZERO, ZERO, Z, LDZ )
1.1 bertrand 430: END IF
431: ELSE IF( INDSV ) THEN
432: *
433: * Find the IL-th through the IU-th singular values. We aim
434: * at -s (see leading comments) and indices are mapped into
435: * values, therefore mimicking DSTEBZ, where
436: *
437: * GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
438: * GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
439: *
440: ILTGK = IL
441: IUTGK = IU
442: RNGVX = 'V'
443: WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
444: CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
445: CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
446: CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ),
447: $ VLTGK, VLTGK, ILTGK, ILTGK, ABSTOL, NS, S,
448: $ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),
449: $ IWORK( IIFAIL ), INFO )
450: VLTGK = S( 1 ) - FUDGE*SMAX*ULP*N
451: WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
452: CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
453: CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
454: CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ),
455: $ VUTGK, VUTGK, IUTGK, IUTGK, ABSTOL, NS, S,
456: $ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ),
457: $ IWORK( IIFAIL ), INFO )
458: VUTGK = S( 1 ) + FUDGE*SMAX*ULP*N
459: VUTGK = MIN( VUTGK, ZERO )
460: *
461: * If VLTGK=VUTGK, DSTEVX returns an error message,
462: * so if needed we change VUTGK slightly.
463: *
464: IF( VLTGK.EQ.VUTGK ) VLTGK = VLTGK - TOL
465: *
1.2 bertrand 466: IF( WANTZ ) CALL DLASET( 'F', N*2, IU-IL+1, ZERO, ZERO, Z, LDZ)
1.1 bertrand 467: END IF
468: *
469: * Initialize variables and pointers for S, Z, and WORK.
470: *
471: * NRU, NRV: number of rows in U and V for the active submatrix
472: * IDBEG, ISBEG: offsets for the entries of D and S
473: * IROWZ, ICOLZ: offsets for the rows and columns of Z
474: * IROWU, IROWV: offsets for the rows of U and V
475: *
476: NS = 0
477: NRU = 0
478: NRV = 0
479: IDBEG = 1
480: ISBEG = 1
481: IROWZ = 1
482: ICOLZ = 1
483: IROWU = 2
484: IROWV = 1
485: SPLIT = .FALSE.
486: SVEQ0 = .FALSE.
487: *
488: * Form the tridiagonal TGK matrix.
489: *
490: S( 1:N ) = ZERO
491: WORK( IETGK+2*N-1 ) = ZERO
492: WORK( IDTGK:IDTGK+2*N-1 ) = ZERO
493: CALL DCOPY( N, D, 1, WORK( IETGK ), 2 )
494: CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 )
495: *
496: *
497: * Check for splits in two levels, outer level
498: * in E and inner level in D.
499: *
500: DO IEPTR = 2, N*2, 2
501: IF( WORK( IETGK+IEPTR-1 ).EQ.ZERO ) THEN
502: *
503: * Split in E (this piece of B is square) or bottom
504: * of the (input bidiagonal) matrix.
505: *
506: ISPLT = IDBEG
507: IDEND = IEPTR - 1
508: DO IDPTR = IDBEG, IDEND, 2
509: IF( WORK( IETGK+IDPTR-1 ).EQ.ZERO ) THEN
510: *
511: * Split in D (rectangular submatrix). Set the number
512: * of rows in U and V (NRU and NRV) accordingly.
513: *
514: IF( IDPTR.EQ.IDBEG ) THEN
515: *
516: * D=0 at the top.
517: *
518: SVEQ0 = .TRUE.
519: IF( IDBEG.EQ.IDEND) THEN
520: NRU = 1
521: NRV = 1
522: END IF
523: ELSE IF( IDPTR.EQ.IDEND ) THEN
524: *
525: * D=0 at the bottom.
526: *
527: SVEQ0 = .TRUE.
528: NRU = (IDEND-ISPLT)/2 + 1
529: NRV = NRU
530: IF( ISPLT.NE.IDBEG ) THEN
531: NRU = NRU + 1
532: END IF
533: ELSE
534: IF( ISPLT.EQ.IDBEG ) THEN
535: *
536: * Split: top rectangular submatrix.
537: *
538: NRU = (IDPTR-IDBEG)/2
539: NRV = NRU + 1
540: ELSE
541: *
542: * Split: middle square submatrix.
543: *
544: NRU = (IDPTR-ISPLT)/2 + 1
545: NRV = NRU
546: END IF
547: END IF
548: ELSE IF( IDPTR.EQ.IDEND ) THEN
549: *
550: * Last entry of D in the active submatrix.
551: *
552: IF( ISPLT.EQ.IDBEG ) THEN
553: *
554: * No split (trivial case).
555: *
556: NRU = (IDEND-IDBEG)/2 + 1
557: NRV = NRU
558: ELSE
559: *
560: * Split: bottom rectangular submatrix.
561: *
562: NRV = (IDEND-ISPLT)/2 + 1
563: NRU = NRV + 1
564: END IF
565: END IF
566: *
567: NTGK = NRU + NRV
568: *
569: IF( NTGK.GT.0 ) THEN
570: *
571: * Compute eigenvalues/vectors of the active
572: * submatrix according to RANGE:
573: * if RANGE='A' (ALLSV) then RNGVX = 'I'
574: * if RANGE='V' (VALSV) then RNGVX = 'V'
575: * if RANGE='I' (INDSV) then RNGVX = 'V'
576: *
577: ILTGK = 1
578: IUTGK = NTGK / 2
579: IF( ALLSV .OR. VUTGK.EQ.ZERO ) THEN
580: IF( SVEQ0 .OR.
581: $ SMIN.LT.EPS .OR.
582: $ MOD(NTGK,2).GT.0 ) THEN
583: * Special case: eigenvalue equal to zero or very
584: * small, additional eigenvector is needed.
585: IUTGK = IUTGK + 1
586: END IF
587: END IF
588: *
589: * Workspace needed by DSTEVX:
590: * WORK( ITEMP: ): 2*5*NTGK
591: * IWORK( 1: ): 2*6*NTGK
592: *
593: CALL DSTEVX( JOBZ, RNGVX, NTGK, WORK( IDTGK+ISPLT-1 ),
594: $ WORK( IETGK+ISPLT-1 ), VLTGK, VUTGK,
595: $ ILTGK, IUTGK, ABSTOL, NSL, S( ISBEG ),
596: $ Z( IROWZ,ICOLZ ), LDZ, WORK( ITEMP ),
597: $ IWORK( IIWORK ), IWORK( IIFAIL ),
598: $ INFO )
599: IF( INFO.NE.0 ) THEN
600: * Exit with the error code from DSTEVX.
601: RETURN
602: END IF
603: EMIN = ABS( MAXVAL( S( ISBEG:ISBEG+NSL-1 ) ) )
604: *
605: IF( NSL.GT.0 .AND. WANTZ ) THEN
606: *
607: * Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:),
608: * changing the sign of v as discussed in the leading
609: * comments. The norms of u and v may be (slightly)
610: * different from 1/sqrt(2) if the corresponding
611: * eigenvalues are very small or too close. We check
612: * those norms and, if needed, reorthogonalize the
613: * vectors.
614: *
615: IF( NSL.GT.1 .AND.
616: $ VUTGK.EQ.ZERO .AND.
617: $ MOD(NTGK,2).EQ.0 .AND.
618: $ EMIN.EQ.0 .AND. .NOT.SPLIT ) THEN
619: *
620: * D=0 at the top or bottom of the active submatrix:
621: * one eigenvalue is equal to zero; concatenate the
622: * eigenvectors corresponding to the two smallest
623: * eigenvalues.
624: *
625: Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) =
626: $ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) +
627: $ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 )
628: Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) =
629: $ ZERO
630: * IF( IUTGK*2.GT.NTGK ) THEN
631: * Eigenvalue equal to zero or very small.
632: * NSL = NSL - 1
633: * END IF
634: END IF
635: *
636: DO I = 0, MIN( NSL-1, NRU-1 )
637: NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 )
638: IF( NRMU.EQ.ZERO ) THEN
639: INFO = N*2 + 1
640: RETURN
641: END IF
642: CALL DSCAL( NRU, ONE/NRMU,
643: $ Z( IROWU,ICOLZ+I ), 2 )
644: IF( NRMU.NE.ONE .AND.
645: $ ABS( NRMU-ORTOL )*SQRT2.GT.ONE )
646: $ THEN
647: DO J = 0, I-1
648: ZJTJI = -DDOT( NRU, Z( IROWU, ICOLZ+J ),
649: $ 2, Z( IROWU, ICOLZ+I ), 2 )
650: CALL DAXPY( NRU, ZJTJI,
651: $ Z( IROWU, ICOLZ+J ), 2,
652: $ Z( IROWU, ICOLZ+I ), 2 )
653: END DO
654: NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 )
655: CALL DSCAL( NRU, ONE/NRMU,
656: $ Z( IROWU,ICOLZ+I ), 2 )
657: END IF
658: END DO
659: DO I = 0, MIN( NSL-1, NRV-1 )
660: NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 )
661: IF( NRMV.EQ.ZERO ) THEN
662: INFO = N*2 + 1
663: RETURN
664: END IF
665: CALL DSCAL( NRV, -ONE/NRMV,
666: $ Z( IROWV,ICOLZ+I ), 2 )
667: IF( NRMV.NE.ONE .AND.
668: $ ABS( NRMV-ORTOL )*SQRT2.GT.ONE )
669: $ THEN
670: DO J = 0, I-1
671: ZJTJI = -DDOT( NRV, Z( IROWV, ICOLZ+J ),
672: $ 2, Z( IROWV, ICOLZ+I ), 2 )
673: CALL DAXPY( NRU, ZJTJI,
674: $ Z( IROWV, ICOLZ+J ), 2,
675: $ Z( IROWV, ICOLZ+I ), 2 )
676: END DO
677: NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 )
678: CALL DSCAL( NRV, ONE/NRMV,
679: $ Z( IROWV,ICOLZ+I ), 2 )
680: END IF
681: END DO
682: IF( VUTGK.EQ.ZERO .AND.
683: $ IDPTR.LT.IDEND .AND.
684: $ MOD(NTGK,2).GT.0 ) THEN
685: *
686: * D=0 in the middle of the active submatrix (one
687: * eigenvalue is equal to zero): save the corresponding
688: * eigenvector for later use (when bottom of the
689: * active submatrix is reached).
690: *
691: SPLIT = .TRUE.
692: Z( IROWZ:IROWZ+NTGK-1,N+1 ) =
693: $ Z( IROWZ:IROWZ+NTGK-1,NS+NSL )
694: Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) =
695: $ ZERO
696: END IF
697: END IF !** WANTZ **!
698: *
699: NSL = MIN( NSL, NRU )
700: SVEQ0 = .FALSE.
701: *
702: * Absolute values of the eigenvalues of TGK.
703: *
704: DO I = 0, NSL-1
705: S( ISBEG+I ) = ABS( S( ISBEG+I ) )
706: END DO
707: *
708: * Update pointers for TGK, S and Z.
709: *
710: ISBEG = ISBEG + NSL
711: IROWZ = IROWZ + NTGK
712: ICOLZ = ICOLZ + NSL
713: IROWU = IROWZ
714: IROWV = IROWZ + 1
715: ISPLT = IDPTR + 1
716: NS = NS + NSL
717: NRU = 0
718: NRV = 0
719: END IF !** NTGK.GT.0 **!
1.2 bertrand 720: IF( IROWZ.LT.N*2 .AND. WANTZ ) THEN
721: Z( 1:IROWZ-1, ICOLZ ) = ZERO
722: END IF
1.1 bertrand 723: END DO !** IDPTR loop **!
1.2 bertrand 724: IF( SPLIT .AND. WANTZ ) THEN
1.1 bertrand 725: *
726: * Bring back eigenvector corresponding
727: * to eigenvalue equal to zero.
728: *
729: Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) =
730: $ Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) +
731: $ Z( IDBEG:IDEND-NTGK+1,N+1 )
732: Z( IDBEG:IDEND-NTGK+1,N+1 ) = 0
733: END IF
734: IROWV = IROWV - 1
735: IROWU = IROWU + 1
736: IDBEG = IEPTR + 1
737: SVEQ0 = .FALSE.
738: SPLIT = .FALSE.
739: END IF !** Check for split in E **!
740: END DO !** IEPTR loop **!
741: *
742: * Sort the singular values into decreasing order (insertion sort on
743: * singular values, but only one transposition per singular vector)
744: *
745: DO I = 1, NS-1
746: K = 1
747: SMIN = S( 1 )
748: DO J = 2, NS + 1 - I
749: IF( S( J ).LE.SMIN ) THEN
750: K = J
751: SMIN = S( J )
752: END IF
753: END DO
754: IF( K.NE.NS+1-I ) THEN
755: S( K ) = S( NS+1-I )
756: S( NS+1-I ) = SMIN
1.2 bertrand 757: IF( WANTZ ) CALL DSWAP( N*2, Z( 1,K ), 1, Z( 1,NS+1-I ), 1 )
1.1 bertrand 758: END IF
759: END DO
760: *
761: * If RANGE=I, check for singular values/vectors to be discarded.
762: *
763: IF( INDSV ) THEN
764: K = IU - IL + 1
765: IF( K.LT.NS ) THEN
766: S( K+1:NS ) = ZERO
1.2 bertrand 767: IF( WANTZ ) Z( 1:N*2,K+1:NS ) = ZERO
1.1 bertrand 768: NS = K
769: END IF
770: END IF
771: *
772: * Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ).
773: * If B is a lower diagonal, swap U and V.
774: *
1.2 bertrand 775: IF( WANTZ ) THEN
1.1 bertrand 776: DO I = 1, NS
777: CALL DCOPY( N*2, Z( 1,I ), 1, WORK, 1 )
778: IF( LOWER ) THEN
779: CALL DCOPY( N, WORK( 2 ), 2, Z( N+1,I ), 1 )
780: CALL DCOPY( N, WORK( 1 ), 2, Z( 1 ,I ), 1 )
781: ELSE
782: CALL DCOPY( N, WORK( 2 ), 2, Z( 1 ,I ), 1 )
783: CALL DCOPY( N, WORK( 1 ), 2, Z( N+1,I ), 1 )
784: END IF
785: END DO
1.2 bertrand 786: END IF
1.1 bertrand 787: *
788: RETURN
789: *
790: * End of DBDSVDX
791: *
792: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dbdsvdx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack