1*> \brief \b DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr. 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 December 2016 105* 106*> \ingroup auxOTHERcomputational 107* 108* ===================================================================== 109 SUBROUTINE DLASQ1( N, D, E, WORK, INFO ) 110* 111* -- LAPACK computational routine (version 3.7.0) -- 112* -- LAPACK is a software package provided by Univ. of Tennessee, -- 113* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 114* December 2016 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 = -1 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 ) 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 ) 218 END IF 219* 220 RETURN 221* 222* End of DLASQ1 223* 224 END 225