LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slasd3()

subroutine slasd3 ( integer  NL,
integer  NR,
integer  SQRE,
integer  K,
real, dimension( * )  D,
real, dimension( ldq, * )  Q,
integer  LDQ,
real, dimension( * )  DSIGMA,
real, dimension( ldu, * )  U,
integer  LDU,
real, dimension( ldu2, * )  U2,
integer  LDU2,
real, dimension( ldvt, * )  VT,
integer  LDVT,
real, dimension( ldvt2, * )  VT2,
integer  LDVT2,
integer, dimension( * )  IDXC,
integer, dimension( * )  CTOT,
real, dimension( * )  Z,
integer  INFO 
)

SLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and Z, and then updates the singular vectors by matrix multiplication. Used by sbdsdc.

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

Purpose:
 SLASD3 finds all the square roots of the roots of the secular
 equation, as defined by the values in D and Z.  It makes the
 appropriate calls to SLASD4 and then updates the singular
 vectors by matrix multiplication.

 This code makes very mild assumptions about floating point
 arithmetic. It will work on machines with a guard digit in
 add/subtract, or on those binary machines without guard digits
 which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
 It could conceivably fail on hexadecimal or decimal machines
 without guard digits, but we know of none.

 SLASD3 is called from SLASD1.
Parameters
[in]NL
          NL is INTEGER
         The row dimension of the upper block.  NL >= 1.
[in]NR
          NR is INTEGER
         The row dimension of the lower block.  NR >= 1.
[in]SQRE
          SQRE is INTEGER
         = 0: the lower block is an NR-by-NR square matrix.
         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

         The bidiagonal matrix has N = NL + NR + 1 rows and
         M = N + SQRE >= N columns.
[in]K
          K is INTEGER
         The size of the secular equation, 1 =< K = < N.
[out]D
          D is REAL array, dimension(K)
         On exit the square roots of the roots of the secular equation,
         in ascending order.
[out]Q
          Q is REAL array, dimension (LDQ,K)
[in]LDQ
          LDQ is INTEGER
         The leading dimension of the array Q.  LDQ >= K.
[in,out]DSIGMA
          DSIGMA is REAL array, dimension(K)
         The first K elements of this array contain the old roots
         of the deflated updating problem.  These are the poles
         of the secular equation.
[out]U
          U is REAL array, dimension (LDU, N)
         The last N - K columns of this matrix contain the deflated
         left singular vectors.
[in]LDU
          LDU is INTEGER
         The leading dimension of the array U.  LDU >= N.
[in]U2
          U2 is REAL array, dimension (LDU2, N)
         The first K columns of this matrix contain the non-deflated
         left singular vectors for the split problem.
[in]LDU2
          LDU2 is INTEGER
         The leading dimension of the array U2.  LDU2 >= N.
[out]VT
          VT is REAL array, dimension (LDVT, M)
         The last M - K columns of VT**T contain the deflated
         right singular vectors.
[in]LDVT
          LDVT is INTEGER
         The leading dimension of the array VT.  LDVT >= N.
[in,out]VT2
          VT2 is REAL array, dimension (LDVT2, N)
         The first K columns of VT2**T contain the non-deflated
         right singular vectors for the split problem.
[in]LDVT2
          LDVT2 is INTEGER
         The leading dimension of the array VT2.  LDVT2 >= N.
[in]IDXC
          IDXC is INTEGER array, dimension (N)
         The permutation used to arrange the columns of U (and rows of
         VT) into three groups:  the first group contains non-zero
         entries only at and above (or before) NL +1; the second
         contains non-zero entries only at and below (or after) NL+2;
         and the third is dense. The first column of U and the row of
         VT are treated separately, however.

         The rows of the singular vectors found by SLASD4
         must be likewise permuted before the matrix multiplies can
         take place.
[in]CTOT
          CTOT is INTEGER array, dimension (4)
         A count of the total number of the various types of columns
         in U (or rows in VT), as described in IDXC. The fourth column
         type is any column which has been deflated.
[in,out]Z
          Z is REAL array, dimension (K)
         The first K elements of this array contain the components
         of the deflation-adjusted updating row vector.
[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 221 of file slasd3.f.

224 *
225 * -- LAPACK auxiliary routine --
226 * -- LAPACK is a software package provided by Univ. of Tennessee, --
227 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
228 *
229 * .. Scalar Arguments ..
230  INTEGER INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR,
231  $ SQRE
232 * ..
233 * .. Array Arguments ..
234  INTEGER CTOT( * ), IDXC( * )
235  REAL D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ),
236  $ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ),
237  $ Z( * )
238 * ..
239 *
240 * =====================================================================
241 *
242 * .. Parameters ..
243  REAL ONE, ZERO, NEGONE
244  parameter( one = 1.0e+0, zero = 0.0e+0,
245  $ negone = -1.0e+0 )
246 * ..
247 * .. Local Scalars ..
248  INTEGER CTEMP, I, J, JC, KTEMP, M, N, NLP1, NLP2, NRP1
249  REAL RHO, TEMP
250 * ..
251 * .. External Functions ..
252  REAL SLAMC3, SNRM2
253  EXTERNAL slamc3, snrm2
254 * ..
255 * .. External Subroutines ..
256  EXTERNAL scopy, sgemm, slacpy, slascl, slasd4, xerbla
257 * ..
258 * .. Intrinsic Functions ..
259  INTRINSIC abs, sign, sqrt
260 * ..
261 * .. Executable Statements ..
262 *
263 * Test the input parameters.
264 *
265  info = 0
266 *
267  IF( nl.LT.1 ) THEN
268  info = -1
269  ELSE IF( nr.LT.1 ) THEN
270  info = -2
271  ELSE IF( ( sqre.NE.1 ) .AND. ( sqre.NE.0 ) ) THEN
272  info = -3
273  END IF
274 *
275  n = nl + nr + 1
276  m = n + sqre
277  nlp1 = nl + 1
278  nlp2 = nl + 2
279 *
280  IF( ( k.LT.1 ) .OR. ( k.GT.n ) ) THEN
281  info = -4
282  ELSE IF( ldq.LT.k ) THEN
283  info = -7
284  ELSE IF( ldu.LT.n ) THEN
285  info = -10
286  ELSE IF( ldu2.LT.n ) THEN
287  info = -12
288  ELSE IF( ldvt.LT.m ) THEN
289  info = -14
290  ELSE IF( ldvt2.LT.m ) THEN
291  info = -16
292  END IF
293  IF( info.NE.0 ) THEN
294  CALL xerbla( 'SLASD3', -info )
295  RETURN
296  END IF
297 *
298 * Quick return if possible
299 *
300  IF( k.EQ.1 ) THEN
301  d( 1 ) = abs( z( 1 ) )
302  CALL scopy( m, vt2( 1, 1 ), ldvt2, vt( 1, 1 ), ldvt )
303  IF( z( 1 ).GT.zero ) THEN
304  CALL scopy( n, u2( 1, 1 ), 1, u( 1, 1 ), 1 )
305  ELSE
306  DO 10 i = 1, n
307  u( i, 1 ) = -u2( i, 1 )
308  10 CONTINUE
309  END IF
310  RETURN
311  END IF
312 *
313 * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
314 * be computed with high relative accuracy (barring over/underflow).
315 * This is a problem on machines without a guard digit in
316 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
317 * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
318 * which on any of these machines zeros out the bottommost
319 * bit of DSIGMA(I) if it is 1; this makes the subsequent
320 * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
321 * occurs. On binary machines with a guard digit (almost all
322 * machines) it does not change DSIGMA(I) at all. On hexadecimal
323 * and decimal machines with a guard digit, it slightly
324 * changes the bottommost bits of DSIGMA(I). It does not account
325 * for hexadecimal or decimal machines without guard digits
326 * (we know of none). We use a subroutine call to compute
327 * 2*DSIGMA(I) to prevent optimizing compilers from eliminating
328 * this code.
329 *
330  DO 20 i = 1, k
331  dsigma( i ) = slamc3( dsigma( i ), dsigma( i ) ) - dsigma( i )
332  20 CONTINUE
333 *
334 * Keep a copy of Z.
335 *
336  CALL scopy( k, z, 1, q, 1 )
337 *
338 * Normalize Z.
339 *
340  rho = snrm2( k, z, 1 )
341  CALL slascl( 'G', 0, 0, rho, one, k, 1, z, k, info )
342  rho = rho*rho
343 *
344 * Find the new singular values.
345 *
346  DO 30 j = 1, k
347  CALL slasd4( k, j, dsigma, z, u( 1, j ), rho, d( j ),
348  $ vt( 1, j ), info )
349 *
350 * If the zero finder fails, report the convergence failure.
351 *
352  IF( info.NE.0 ) THEN
353  RETURN
354  END IF
355  30 CONTINUE
356 *
357 * Compute updated Z.
358 *
359  DO 60 i = 1, k
360  z( i ) = u( i, k )*vt( i, k )
361  DO 40 j = 1, i - 1
362  z( i ) = z( i )*( u( i, j )*vt( i, j ) /
363  $ ( dsigma( i )-dsigma( j ) ) /
364  $ ( dsigma( i )+dsigma( j ) ) )
365  40 CONTINUE
366  DO 50 j = i, k - 1
367  z( i ) = z( i )*( u( i, j )*vt( i, j ) /
368  $ ( dsigma( i )-dsigma( j+1 ) ) /
369  $ ( dsigma( i )+dsigma( j+1 ) ) )
370  50 CONTINUE
371  z( i ) = sign( sqrt( abs( z( i ) ) ), q( i, 1 ) )
372  60 CONTINUE
373 *
374 * Compute left singular vectors of the modified diagonal matrix,
375 * and store related information for the right singular vectors.
376 *
377  DO 90 i = 1, k
378  vt( 1, i ) = z( 1 ) / u( 1, i ) / vt( 1, i )
379  u( 1, i ) = negone
380  DO 70 j = 2, k
381  vt( j, i ) = z( j ) / u( j, i ) / vt( j, i )
382  u( j, i ) = dsigma( j )*vt( j, i )
383  70 CONTINUE
384  temp = snrm2( k, u( 1, i ), 1 )
385  q( 1, i ) = u( 1, i ) / temp
386  DO 80 j = 2, k
387  jc = idxc( j )
388  q( j, i ) = u( jc, i ) / temp
389  80 CONTINUE
390  90 CONTINUE
391 *
392 * Update the left singular vector matrix.
393 *
394  IF( k.EQ.2 ) THEN
395  CALL sgemm( 'N', 'N', n, k, k, one, u2, ldu2, q, ldq, zero, u,
396  $ ldu )
397  GO TO 100
398  END IF
399  IF( ctot( 1 ).GT.0 ) THEN
400  CALL sgemm( 'N', 'N', nl, k, ctot( 1 ), one, u2( 1, 2 ), ldu2,
401  $ q( 2, 1 ), ldq, zero, u( 1, 1 ), ldu )
402  IF( ctot( 3 ).GT.0 ) THEN
403  ktemp = 2 + ctot( 1 ) + ctot( 2 )
404  CALL sgemm( 'N', 'N', nl, k, ctot( 3 ), one, u2( 1, ktemp ),
405  $ ldu2, q( ktemp, 1 ), ldq, one, u( 1, 1 ), ldu )
406  END IF
407  ELSE IF( ctot( 3 ).GT.0 ) THEN
408  ktemp = 2 + ctot( 1 ) + ctot( 2 )
409  CALL sgemm( 'N', 'N', nl, k, ctot( 3 ), one, u2( 1, ktemp ),
410  $ ldu2, q( ktemp, 1 ), ldq, zero, u( 1, 1 ), ldu )
411  ELSE
412  CALL slacpy( 'F', nl, k, u2, ldu2, u, ldu )
413  END IF
414  CALL scopy( k, q( 1, 1 ), ldq, u( nlp1, 1 ), ldu )
415  ktemp = 2 + ctot( 1 )
416  ctemp = ctot( 2 ) + ctot( 3 )
417  CALL sgemm( 'N', 'N', nr, k, ctemp, one, u2( nlp2, ktemp ), ldu2,
418  $ q( ktemp, 1 ), ldq, zero, u( nlp2, 1 ), ldu )
419 *
420 * Generate the right singular vectors.
421 *
422  100 CONTINUE
423  DO 120 i = 1, k
424  temp = snrm2( k, vt( 1, i ), 1 )
425  q( i, 1 ) = vt( 1, i ) / temp
426  DO 110 j = 2, k
427  jc = idxc( j )
428  q( i, j ) = vt( jc, i ) / temp
429  110 CONTINUE
430  120 CONTINUE
431 *
432 * Update the right singular vector matrix.
433 *
434  IF( k.EQ.2 ) THEN
435  CALL sgemm( 'N', 'N', k, m, k, one, q, ldq, vt2, ldvt2, zero,
436  $ vt, ldvt )
437  RETURN
438  END IF
439  ktemp = 1 + ctot( 1 )
440  CALL sgemm( 'N', 'N', k, nlp1, ktemp, one, q( 1, 1 ), ldq,
441  $ vt2( 1, 1 ), ldvt2, zero, vt( 1, 1 ), ldvt )
442  ktemp = 2 + ctot( 1 ) + ctot( 2 )
443  IF( ktemp.LE.ldvt2 )
444  $ CALL sgemm( 'N', 'N', k, nlp1, ctot( 3 ), one, q( 1, ktemp ),
445  $ ldq, vt2( ktemp, 1 ), ldvt2, one, vt( 1, 1 ),
446  $ ldvt )
447 *
448  ktemp = ctot( 1 ) + 1
449  nrp1 = nr + sqre
450  IF( ktemp.GT.1 ) THEN
451  DO 130 i = 1, k
452  q( i, ktemp ) = q( i, 1 )
453  130 CONTINUE
454  DO 140 i = nlp2, m
455  vt2( ktemp, i ) = vt2( 1, i )
456  140 CONTINUE
457  END IF
458  ctemp = 1 + ctot( 2 ) + ctot( 3 )
459  CALL sgemm( 'N', 'N', k, nrp1, ctemp, one, q( 1, ktemp ), ldq,
460  $ vt2( ktemp, nlp2 ), ldvt2, zero, vt( 1, nlp2 ), ldvt )
461 *
462  RETURN
463 *
464 * End of SLASD3
465 *
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 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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
real(wp) function snrm2(n, x, incx)
SNRM2
Definition: snrm2.f90:89
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
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: