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