LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlasd3.f
Go to the documentation of this file.
1 *> \brief \b 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.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLASD3 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd3.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd3.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd3.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLASD3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2,
22 * LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z,
23 * INFO )
24 *
25 * .. Scalar Arguments ..
26 * INTEGER INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR,
27 * $ SQRE
28 * ..
29 * .. Array Arguments ..
30 * INTEGER CTOT( * ), IDXC( * )
31 * DOUBLE PRECISION D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ),
32 * $ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ),
33 * $ Z( * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> DLASD3 finds all the square roots of the roots of the secular
43 *> equation, as defined by the values in D and Z. It makes the
44 *> appropriate calls to DLASD4 and then updates the singular
45 *> vectors by matrix multiplication.
46 *>
47 *> This code makes very mild assumptions about floating point
48 *> arithmetic. It will work on machines with a guard digit in
49 *> add/subtract, or on those binary machines without guard digits
50 *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
51 *> It could conceivably fail on hexadecimal or decimal machines
52 *> without guard digits, but we know of none.
53 *>
54 *> DLASD3 is called from DLASD1.
55 *> \endverbatim
56 *
57 * Arguments:
58 * ==========
59 *
60 *> \param[in] NL
61 *> \verbatim
62 *> NL is INTEGER
63 *> The row dimension of the upper block. NL >= 1.
64 *> \endverbatim
65 *>
66 *> \param[in] NR
67 *> \verbatim
68 *> NR is INTEGER
69 *> The row dimension of the lower block. NR >= 1.
70 *> \endverbatim
71 *>
72 *> \param[in] SQRE
73 *> \verbatim
74 *> SQRE is INTEGER
75 *> = 0: the lower block is an NR-by-NR square matrix.
76 *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
77 *>
78 *> The bidiagonal matrix has N = NL + NR + 1 rows and
79 *> M = N + SQRE >= N columns.
80 *> \endverbatim
81 *>
82 *> \param[in] K
83 *> \verbatim
84 *> K is INTEGER
85 *> The size of the secular equation, 1 =< K = < N.
86 *> \endverbatim
87 *>
88 *> \param[out] D
89 *> \verbatim
90 *> D is DOUBLE PRECISION array, dimension(K)
91 *> On exit the square roots of the roots of the secular equation,
92 *> in ascending order.
93 *> \endverbatim
94 *>
95 *> \param[out] Q
96 *> \verbatim
97 *> Q is DOUBLE PRECISION array, dimension (LDQ,K)
98 *> \endverbatim
99 *>
100 *> \param[in] LDQ
101 *> \verbatim
102 *> LDQ is INTEGER
103 *> The leading dimension of the array Q. LDQ >= K.
104 *> \endverbatim
105 *>
106 *> \param[in,out] DSIGMA
107 *> \verbatim
108 *> DSIGMA is DOUBLE PRECISION array, dimension(K)
109 *> The first K elements of this array contain the old roots
110 *> of the deflated updating problem. These are the poles
111 *> of the secular equation.
112 *> \endverbatim
113 *>
114 *> \param[out] U
115 *> \verbatim
116 *> U is DOUBLE PRECISION array, dimension (LDU, N)
117 *> The last N - K columns of this matrix contain the deflated
118 *> left singular vectors.
119 *> \endverbatim
120 *>
121 *> \param[in] LDU
122 *> \verbatim
123 *> LDU is INTEGER
124 *> The leading dimension of the array U. LDU >= N.
125 *> \endverbatim
126 *>
127 *> \param[in] U2
128 *> \verbatim
129 *> U2 is DOUBLE PRECISION array, dimension (LDU2, N)
130 *> The first K columns of this matrix contain the non-deflated
131 *> left singular vectors for the split problem.
132 *> \endverbatim
133 *>
134 *> \param[in] LDU2
135 *> \verbatim
136 *> LDU2 is INTEGER
137 *> The leading dimension of the array U2. LDU2 >= N.
138 *> \endverbatim
139 *>
140 *> \param[out] VT
141 *> \verbatim
142 *> VT is DOUBLE PRECISION array, dimension (LDVT, M)
143 *> The last M - K columns of VT**T contain the deflated
144 *> right singular vectors.
145 *> \endverbatim
146 *>
147 *> \param[in] LDVT
148 *> \verbatim
149 *> LDVT is INTEGER
150 *> The leading dimension of the array VT. LDVT >= N.
151 *> \endverbatim
152 *>
153 *> \param[in,out] VT2
154 *> \verbatim
155 *> VT2 is DOUBLE PRECISION array, dimension (LDVT2, N)
156 *> The first K columns of VT2**T contain the non-deflated
157 *> right singular vectors for the split problem.
158 *> \endverbatim
159 *>
160 *> \param[in] LDVT2
161 *> \verbatim
162 *> LDVT2 is INTEGER
163 *> The leading dimension of the array VT2. LDVT2 >= N.
164 *> \endverbatim
165 *>
166 *> \param[in] IDXC
167 *> \verbatim
168 *> IDXC is INTEGER array, dimension ( N )
169 *> The permutation used to arrange the columns of U (and rows of
170 *> VT) into three groups: the first group contains non-zero
171 *> entries only at and above (or before) NL +1; the second
172 *> contains non-zero entries only at and below (or after) NL+2;
173 *> and the third is dense. The first column of U and the row of
174 *> VT are treated separately, however.
175 *>
176 *> The rows of the singular vectors found by DLASD4
177 *> must be likewise permuted before the matrix multiplies can
178 *> take place.
179 *> \endverbatim
180 *>
181 *> \param[in] CTOT
182 *> \verbatim
183 *> CTOT is INTEGER array, dimension ( 4 )
184 *> A count of the total number of the various types of columns
185 *> in U (or rows in VT), as described in IDXC. The fourth column
186 *> type is any column which has been deflated.
187 *> \endverbatim
188 *>
189 *> \param[in,out] Z
190 *> \verbatim
191 *> Z is DOUBLE PRECISION array, dimension (K)
192 *> The first K elements of this array contain the components
193 *> of the deflation-adjusted updating row vector.
194 *> \endverbatim
195 *>
196 *> \param[out] INFO
197 *> \verbatim
198 *> INFO is INTEGER
199 *> = 0: successful exit.
200 *> < 0: if INFO = -i, the i-th argument had an illegal value.
201 *> > 0: if INFO = 1, a singular value did not converge
202 *> \endverbatim
203 *
204 * Authors:
205 * ========
206 *
207 *> \author Univ. of Tennessee
208 *> \author Univ. of California Berkeley
209 *> \author Univ. of Colorado Denver
210 *> \author NAG Ltd.
211 *
212 *> \ingroup OTHERauxiliary
213 *
214 *> \par Contributors:
215 * ==================
216 *>
217 *> Ming Gu and Huan Ren, Computer Science Division, University of
218 *> California at Berkeley, USA
219 *>
220 * =====================================================================
221  SUBROUTINE dlasd3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2,
222  $ LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z,
223  $ INFO )
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 *
466  END
subroutine dlasd3(NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2, LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z, INFO)
DLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and...
Definition: dlasd3.f:224
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