Annotation of rpl/lapack/lapack/dlasq1.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b DLASQ1
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLASQ1 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq1.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq1.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq1.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * INTEGER INFO, N
! 25: * ..
! 26: * .. Array Arguments ..
! 27: * DOUBLE PRECISION D( * ), E( * ), WORK( * )
! 28: * ..
! 29: *
! 30: *
! 31: *> \par Purpose:
! 32: * =============
! 33: *>
! 34: *> \verbatim
! 35: *>
! 36: *> DLASQ1 computes the singular values of a real N-by-N bidiagonal
! 37: *> matrix with diagonal D and off-diagonal E. The singular values
! 38: *> are computed to high relative accuracy, in the absence of
! 39: *> denormalization, underflow and overflow. The algorithm was first
! 40: *> presented in
! 41: *>
! 42: *> "Accurate singular values and differential qd algorithms" by K. V.
! 43: *> Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
! 44: *> 1994,
! 45: *>
! 46: *> and the present implementation is described in "An implementation of
! 47: *> the dqds Algorithm (Positive Case)", LAPACK Working Note.
! 48: *> \endverbatim
! 49: *
! 50: * Arguments:
! 51: * ==========
! 52: *
! 53: *> \param[in] N
! 54: *> \verbatim
! 55: *> N is INTEGER
! 56: *> The number of rows and columns in the matrix. N >= 0.
! 57: *> \endverbatim
! 58: *>
! 59: *> \param[in,out] D
! 60: *> \verbatim
! 61: *> D is DOUBLE PRECISION array, dimension (N)
! 62: *> On entry, D contains the diagonal elements of the
! 63: *> bidiagonal matrix whose SVD is desired. On normal exit,
! 64: *> D contains the singular values in decreasing order.
! 65: *> \endverbatim
! 66: *>
! 67: *> \param[in,out] E
! 68: *> \verbatim
! 69: *> E is DOUBLE PRECISION array, dimension (N)
! 70: *> On entry, elements E(1:N-1) contain the off-diagonal elements
! 71: *> of the bidiagonal matrix whose SVD is desired.
! 72: *> On exit, E is overwritten.
! 73: *> \endverbatim
! 74: *>
! 75: *> \param[out] WORK
! 76: *> \verbatim
! 77: *> WORK is DOUBLE PRECISION array, dimension (4*N)
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[out] INFO
! 81: *> \verbatim
! 82: *> INFO is INTEGER
! 83: *> = 0: successful exit
! 84: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 85: *> > 0: the algorithm failed
! 86: *> = 1, a split was marked by a positive value in E
! 87: *> = 2, current block of Z not diagonalized after 100*N
! 88: *> iterations (in inner while loop) On exit D and E
! 89: *> represent a matrix with the same singular values
! 90: *> which the calling subroutine could use to finish the
! 91: *> computation, or even feed back into DLASQ1
! 92: *> = 3, termination criterion of outer while loop not met
! 93: *> (program created more than N unreduced blocks)
! 94: *> \endverbatim
! 95: *
! 96: * Authors:
! 97: * ========
! 98: *
! 99: *> \author Univ. of Tennessee
! 100: *> \author Univ. of California Berkeley
! 101: *> \author Univ. of Colorado Denver
! 102: *> \author NAG Ltd.
! 103: *
! 104: *> \date November 2011
1.1 bertrand 105: *
1.8 ! bertrand 106: *> \ingroup auxOTHERcomputational
1.1 bertrand 107: *
1.8 ! bertrand 108: * =====================================================================
! 109: SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
1.1 bertrand 110: *
1.8 ! bertrand 111: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 112: * -- LAPACK is a software package provided by Univ. of Tennessee, --
113: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 114: * November 2011
1.1 bertrand 115: *
116: * .. Scalar Arguments ..
117: INTEGER INFO, N
118: * ..
119: * .. Array Arguments ..
120: DOUBLE PRECISION D( * ), E( * ), WORK( * )
121: * ..
122: *
123: * =====================================================================
124: *
125: * .. Parameters ..
126: DOUBLE PRECISION ZERO
127: PARAMETER ( ZERO = 0.0D0 )
128: * ..
129: * .. Local Scalars ..
130: INTEGER I, IINFO
131: DOUBLE PRECISION EPS, SCALE, SAFMIN, SIGMN, SIGMX
132: * ..
133: * .. External Subroutines ..
134: EXTERNAL DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
135: * ..
136: * .. External Functions ..
137: DOUBLE PRECISION DLAMCH
138: EXTERNAL DLAMCH
139: * ..
140: * .. Intrinsic Functions ..
141: INTRINSIC ABS, MAX, SQRT
142: * ..
143: * .. Executable Statements ..
144: *
145: INFO = 0
146: IF( N.LT.0 ) THEN
147: INFO = -2
148: CALL XERBLA( 'DLASQ1', -INFO )
149: RETURN
150: ELSE IF( N.EQ.0 ) THEN
151: RETURN
152: ELSE IF( N.EQ.1 ) THEN
153: D( 1 ) = ABS( D( 1 ) )
154: RETURN
155: ELSE IF( N.EQ.2 ) THEN
156: CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
157: D( 1 ) = SIGMX
158: D( 2 ) = SIGMN
159: RETURN
160: END IF
161: *
162: * Estimate the largest singular value.
163: *
164: SIGMX = ZERO
165: DO 10 I = 1, N - 1
166: D( I ) = ABS( D( I ) )
167: SIGMX = MAX( SIGMX, ABS( E( I ) ) )
168: 10 CONTINUE
169: D( N ) = ABS( D( N ) )
170: *
171: * Early return if SIGMX is zero (matrix is already diagonal).
172: *
173: IF( SIGMX.EQ.ZERO ) THEN
174: CALL DLASRT( 'D', N, D, IINFO )
175: RETURN
176: END IF
177: *
178: DO 20 I = 1, N
179: SIGMX = MAX( SIGMX, D( I ) )
180: 20 CONTINUE
181: *
182: * Copy D and E into WORK (in the Z format) and scale (squaring the
183: * input data makes scaling by a power of the radix pointless).
184: *
185: EPS = DLAMCH( 'Precision' )
186: SAFMIN = DLAMCH( 'Safe minimum' )
187: SCALE = SQRT( EPS / SAFMIN )
188: CALL DCOPY( N, D, 1, WORK( 1 ), 2 )
189: CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 )
190: CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
191: $ IINFO )
192: *
193: * Compute the q's and e's.
194: *
195: DO 30 I = 1, 2*N - 1
196: WORK( I ) = WORK( I )**2
197: 30 CONTINUE
198: WORK( 2*N ) = ZERO
199: *
200: CALL DLASQ2( N, WORK, INFO )
201: *
202: IF( INFO.EQ.0 ) THEN
203: DO 40 I = 1, N
204: D( I ) = SQRT( WORK( I ) )
205: 40 CONTINUE
206: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
1.8 ! bertrand 207: ELSE IF( INFO.EQ.2 ) THEN
! 208: *
! 209: * Maximum number of iterations exceeded. Move data from WORK
! 210: * into D and E so the calling subroutine can try to finish
! 211: *
! 212: DO I = 1, N
! 213: D( I ) = SQRT( WORK( 2*I-1 ) )
! 214: E( I ) = SQRT( WORK( 2*I ) )
! 215: END DO
! 216: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
! 217: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, E, N, IINFO )
1.1 bertrand 218: END IF
219: *
220: RETURN
221: *
222: * End of DLASQ1
223: *
224: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack