LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ dlasd3()

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

DLASD3 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.

Purpose:
``` DLASD3 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 DLASD4 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.

DLASD3 is called from DLASD1.```
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 DOUBLE PRECISION array, dimension(K) On exit the square roots of the roots of the secular equation, in ascending order.``` [out] Q ` Q is DOUBLE PRECISION array, dimension (LDQ,K)` [in] LDQ ``` LDQ is INTEGER The leading dimension of the array Q. LDQ >= K.``` [in,out] DSIGMA ``` DSIGMA is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLASD4 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 DOUBLE PRECISION 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```
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 221 of file dlasd3.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  DOUBLE PRECISION D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ),
236  \$ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ),
237  \$ Z( * )
238 * ..
239 *
240 * =====================================================================
241 *
242 * .. Parameters ..
243  DOUBLE PRECISION ONE, ZERO, NEGONE
244  parameter( one = 1.0d+0, zero = 0.0d+0,
245  \$ negone = -1.0d+0 )
246 * ..
247 * .. Local Scalars ..
248  INTEGER CTEMP, I, J, JC, KTEMP, M, N, NLP1, NLP2, NRP1
249  DOUBLE PRECISION RHO, TEMP
250 * ..
251 * .. External Functions ..
252  DOUBLE PRECISION DLAMC3, DNRM2
253  EXTERNAL dlamc3, dnrm2
254 * ..
255 * .. External Subroutines ..
256  EXTERNAL dcopy, dgemm, dlacpy, dlascl, dlasd4, 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( 'DLASD3', -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 dcopy( m, vt2( 1, 1 ), ldvt2, vt( 1, 1 ), ldvt )
303  IF( z( 1 ).GT.zero ) THEN
304  CALL dcopy( 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 ) = dlamc3( dsigma( i ), dsigma( i ) ) - dsigma( i )
332  20 CONTINUE
333 *
334 * Keep a copy of Z.
335 *
336  CALL dcopy( k, z, 1, q, 1 )
337 *
338 * Normalize Z.
339 *
340  rho = dnrm2( k, z, 1 )
341  CALL dlascl( '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 dlasd4( 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 = dnrm2( 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 dgemm( '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 dgemm( '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 dgemm( '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 dgemm( '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 dlacpy( 'F', nl, k, u2, ldu2, u, ldu )
413  END IF
414  CALL dcopy( k, q( 1, 1 ), ldq, u( nlp1, 1 ), ldu )
415  ktemp = 2 + ctot( 1 )
416  ctemp = ctot( 2 ) + ctot( 3 )
417  CALL dgemm( '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 = dnrm2( 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 dgemm( '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 dgemm( '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 dgemm( '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 dgemm( '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 DLASD3
465 *
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:169
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlasd4(N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO)
DLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modif...
Definition: dlasd4.f:153
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition: dnrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: