LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slasd8()

subroutine slasd8 ( integer  ICOMPQ,
integer  K,
real, dimension( * )  D,
real, dimension( * )  Z,
real, dimension( * )  VF,
real, dimension( * )  VL,
real, dimension( * )  DIFL,
real, dimension( lddifr, * )  DIFR,
integer  LDDIFR,
real, dimension( * )  DSIGMA,
real, dimension( * )  WORK,
integer  INFO 
)

SLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. Used by sbdsdc.

Download SLASD8 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SLASD8 finds the square roots of the roots of the secular equation,
 as defined by the values in DSIGMA and Z. It makes the appropriate
 calls to SLASD4, and stores, for each  element in D, the distance
 to its two nearest poles (elements in DSIGMA). It also updates
 the arrays VF and VL, the first and last components of all the
 right singular vectors of the original bidiagonal matrix.

 SLASD8 is called from SLASD6.
Parameters
[in]ICOMPQ
          ICOMPQ is INTEGER
          Specifies whether singular vectors are to be computed in
          factored form in the calling routine:
          = 0: Compute singular values only.
          = 1: Compute singular vectors in factored form as well.
[in]K
          K is INTEGER
          The number of terms in the rational function to be solved
          by SLASD4.  K >= 1.
[out]D
          D is REAL array, dimension ( K )
          On output, D contains the updated singular values.
[in,out]Z
          Z is REAL array, dimension ( K )
          On entry, the first K elements of this array contain the
          components of the deflation-adjusted updating row vector.
          On exit, Z is updated.
[in,out]VF
          VF is REAL array, dimension ( K )
          On entry, VF contains  information passed through DBEDE8.
          On exit, VF contains the first K components of the first
          components of all right singular vectors of the bidiagonal
          matrix.
[in,out]VL
          VL is REAL array, dimension ( K )
          On entry, VL contains  information passed through DBEDE8.
          On exit, VL contains the first K components of the last
          components of all right singular vectors of the bidiagonal
          matrix.
[out]DIFL
          DIFL is REAL array, dimension ( K )
          On exit, DIFL(I) = D(I) - DSIGMA(I).
[out]DIFR
          DIFR is REAL array,
                   dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
                   dimension ( K ) if ICOMPQ = 0.
          On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
          defined and will not be referenced.

          If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
          normalizing factors for the right singular vector matrix.
[in]LDDIFR
          LDDIFR is INTEGER
          The leading dimension of DIFR, must be at least K.
[in,out]DSIGMA
          DSIGMA is REAL array, dimension ( K )
          On entry, the first K elements of this array contain the old
          roots of the deflated updating problem.  These are the poles
          of the secular equation.
          On exit, the elements of DSIGMA may be very slightly altered
          in value.
[out]WORK
          WORK is REAL array, dimension (3*K)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  if INFO = 1, a singular value did not converge
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 164 of file slasd8.f.

166 *
167 * -- LAPACK auxiliary routine --
168 * -- LAPACK is a software package provided by Univ. of Tennessee, --
169 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
170 *
171 * .. Scalar Arguments ..
172  INTEGER ICOMPQ, INFO, K, LDDIFR
173 * ..
174 * .. Array Arguments ..
175  REAL D( * ), DIFL( * ), DIFR( LDDIFR, * ),
176  $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
177  $ Z( * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  REAL ONE
184  parameter( one = 1.0e+0 )
185 * ..
186 * .. Local Scalars ..
187  INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
188  REAL DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
189 * ..
190 * .. External Subroutines ..
191  EXTERNAL scopy, slascl, slasd4, slaset, xerbla
192 * ..
193 * .. External Functions ..
194  REAL SDOT, SLAMC3, SNRM2
195  EXTERNAL sdot, slamc3, snrm2
196 * ..
197 * .. Intrinsic Functions ..
198  INTRINSIC abs, sign, sqrt
199 * ..
200 * .. Executable Statements ..
201 *
202 * Test the input parameters.
203 *
204  info = 0
205 *
206  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
207  info = -1
208  ELSE IF( k.LT.1 ) THEN
209  info = -2
210  ELSE IF( lddifr.LT.k ) THEN
211  info = -9
212  END IF
213  IF( info.NE.0 ) THEN
214  CALL xerbla( 'SLASD8', -info )
215  RETURN
216  END IF
217 *
218 * Quick return if possible
219 *
220  IF( k.EQ.1 ) THEN
221  d( 1 ) = abs( z( 1 ) )
222  difl( 1 ) = d( 1 )
223  IF( icompq.EQ.1 ) THEN
224  difl( 2 ) = one
225  difr( 1, 2 ) = one
226  END IF
227  RETURN
228  END IF
229 *
230 * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
231 * be computed with high relative accuracy (barring over/underflow).
232 * This is a problem on machines without a guard digit in
233 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
234 * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
235 * which on any of these machines zeros out the bottommost
236 * bit of DSIGMA(I) if it is 1; this makes the subsequent
237 * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
238 * occurs. On binary machines with a guard digit (almost all
239 * machines) it does not change DSIGMA(I) at all. On hexadecimal
240 * and decimal machines with a guard digit, it slightly
241 * changes the bottommost bits of DSIGMA(I). It does not account
242 * for hexadecimal or decimal machines without guard digits
243 * (we know of none). We use a subroutine call to compute
244 * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
245 * this code.
246 *
247  DO 10 i = 1, k
248  dsigma( i ) = slamc3( dsigma( i ), dsigma( i ) ) - dsigma( i )
249  10 CONTINUE
250 *
251 * Book keeping.
252 *
253  iwk1 = 1
254  iwk2 = iwk1 + k
255  iwk3 = iwk2 + k
256  iwk2i = iwk2 - 1
257  iwk3i = iwk3 - 1
258 *
259 * Normalize Z.
260 *
261  rho = snrm2( k, z, 1 )
262  CALL slascl( 'G', 0, 0, rho, one, k, 1, z, k, info )
263  rho = rho*rho
264 *
265 * Initialize WORK(IWK3).
266 *
267  CALL slaset( 'A', k, 1, one, one, work( iwk3 ), k )
268 *
269 * Compute the updated singular values, the arrays DIFL, DIFR,
270 * and the updated Z.
271 *
272  DO 40 j = 1, k
273  CALL slasd4( k, j, dsigma, z, work( iwk1 ), rho, d( j ),
274  $ work( iwk2 ), info )
275 *
276 * If the root finder fails, report the convergence failure.
277 *
278  IF( info.NE.0 ) THEN
279  RETURN
280  END IF
281  work( iwk3i+j ) = work( iwk3i+j )*work( j )*work( iwk2i+j )
282  difl( j ) = -work( j )
283  difr( j, 1 ) = -work( j+1 )
284  DO 20 i = 1, j - 1
285  work( iwk3i+i ) = work( iwk3i+i )*work( i )*
286  $ work( iwk2i+i ) / ( dsigma( i )-
287  $ dsigma( j ) ) / ( dsigma( i )+
288  $ dsigma( j ) )
289  20 CONTINUE
290  DO 30 i = j + 1, k
291  work( iwk3i+i ) = work( iwk3i+i )*work( i )*
292  $ work( iwk2i+i ) / ( dsigma( i )-
293  $ dsigma( j ) ) / ( dsigma( i )+
294  $ dsigma( j ) )
295  30 CONTINUE
296  40 CONTINUE
297 *
298 * Compute updated Z.
299 *
300  DO 50 i = 1, k
301  z( i ) = sign( sqrt( abs( work( iwk3i+i ) ) ), z( i ) )
302  50 CONTINUE
303 *
304 * Update VF and VL.
305 *
306  DO 80 j = 1, k
307  diflj = difl( j )
308  dj = d( j )
309  dsigj = -dsigma( j )
310  IF( j.LT.k ) THEN
311  difrj = -difr( j, 1 )
312  dsigjp = -dsigma( j+1 )
313  END IF
314  work( j ) = -z( j ) / diflj / ( dsigma( j )+dj )
315  DO 60 i = 1, j - 1
316  work( i ) = z( i ) / ( slamc3( dsigma( i ), dsigj )-diflj )
317  $ / ( dsigma( i )+dj )
318  60 CONTINUE
319  DO 70 i = j + 1, k
320  work( i ) = z( i ) / ( slamc3( dsigma( i ), dsigjp )+difrj )
321  $ / ( dsigma( i )+dj )
322  70 CONTINUE
323  temp = snrm2( k, work, 1 )
324  work( iwk2i+j ) = sdot( k, work, 1, vf, 1 ) / temp
325  work( iwk3i+j ) = sdot( k, work, 1, vl, 1 ) / temp
326  IF( icompq.EQ.1 ) THEN
327  difr( j, 2 ) = temp
328  END IF
329  80 CONTINUE
330 *
331  CALL scopy( k, work( iwk2 ), 1, vf, 1 )
332  CALL scopy( k, work( iwk3 ), 1, vl, 1 )
333 *
334  RETURN
335 *
336 * End of SLASD8
337 *
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slasd4(N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO)
SLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modif...
Definition: slasd4.f:153
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
real(wp) function snrm2(n, x, incx)
SNRM2
Definition: snrm2.f90:89
real function slamc3(A, B)
SLAMC3
Definition: slamch.f:169
Here is the call graph for this function:
Here is the caller graph for this function: