LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slasda.f
Go to the documentation of this file.
1*> \brief \b SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagonal d and off-diagonal e. 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 SLASDA + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasda.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasda.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasda.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLASDA( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K,
22* DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL,
23* PERM, GIVNUM, C, S, WORK, IWORK, INFO )
24*
25* .. Scalar Arguments ..
26* INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE
27* ..
28* .. Array Arguments ..
29* INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
30* $ K( * ), PERM( LDGCOL, * )
31* REAL C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ),
32* $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ),
33* $ S( * ), U( LDU, * ), VT( LDU, * ), WORK( * ),
34* $ Z( LDU, * )
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*> \verbatim
42*>
43*> Using a divide and conquer approach, SLASDA computes the singular
44*> value decomposition (SVD) of a real upper bidiagonal N-by-M matrix
45*> B with diagonal D and offdiagonal E, where M = N + SQRE. The
46*> algorithm computes the singular values in the SVD B = U * S * VT.
47*> The orthogonal matrices U and VT are optionally computed in
48*> compact form.
49*>
50*> A related subroutine, SLASD0, computes the singular values and
51*> the singular vectors in explicit form.
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] ICOMPQ
58*> \verbatim
59*> ICOMPQ is INTEGER
60*> Specifies whether singular vectors are to be computed
61*> in compact form, as follows
62*> = 0: Compute singular values only.
63*> = 1: Compute singular vectors of upper bidiagonal
64*> matrix in compact form.
65*> \endverbatim
66*>
67*> \param[in] SMLSIZ
68*> \verbatim
69*> SMLSIZ is INTEGER
70*> The maximum size of the subproblems at the bottom of the
71*> computation tree.
72*> \endverbatim
73*>
74*> \param[in] N
75*> \verbatim
76*> N is INTEGER
77*> The row dimension of the upper bidiagonal matrix. This is
78*> also the dimension of the main diagonal array D.
79*> \endverbatim
80*>
81*> \param[in] SQRE
82*> \verbatim
83*> SQRE is INTEGER
84*> Specifies the column dimension of the bidiagonal matrix.
85*> = 0: The bidiagonal matrix has column dimension M = N;
86*> = 1: The bidiagonal matrix has column dimension M = N + 1.
87*> \endverbatim
88*>
89*> \param[in,out] D
90*> \verbatim
91*> D is REAL array, dimension ( N )
92*> On entry D contains the main diagonal of the bidiagonal
93*> matrix. On exit D, if INFO = 0, contains its singular values.
94*> \endverbatim
95*>
96*> \param[in] E
97*> \verbatim
98*> E is REAL array, dimension ( M-1 )
99*> Contains the subdiagonal entries of the bidiagonal matrix.
100*> On exit, E has been destroyed.
101*> \endverbatim
102*>
103*> \param[out] U
104*> \verbatim
105*> U is REAL array,
106*> dimension ( LDU, SMLSIZ ) if ICOMPQ = 1, and not referenced
107*> if ICOMPQ = 0. If ICOMPQ = 1, on exit, U contains the left
108*> singular vector matrices of all subproblems at the bottom
109*> level.
110*> \endverbatim
111*>
112*> \param[in] LDU
113*> \verbatim
114*> LDU is INTEGER, LDU = > N.
115*> The leading dimension of arrays U, VT, DIFL, DIFR, POLES,
116*> GIVNUM, and Z.
117*> \endverbatim
118*>
119*> \param[out] VT
120*> \verbatim
121*> VT is REAL array,
122*> dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced
123*> if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT**T contains the right
124*> singular vector matrices of all subproblems at the bottom
125*> level.
126*> \endverbatim
127*>
128*> \param[out] K
129*> \verbatim
130*> K is INTEGER array, dimension ( N )
131*> if ICOMPQ = 1 and dimension 1 if ICOMPQ = 0.
132*> If ICOMPQ = 1, on exit, K(I) is the dimension of the I-th
133*> secular equation on the computation tree.
134*> \endverbatim
135*>
136*> \param[out] DIFL
137*> \verbatim
138*> DIFL is REAL array, dimension ( LDU, NLVL ),
139*> where NLVL = floor(log_2 (N/SMLSIZ))).
140*> \endverbatim
141*>
142*> \param[out] DIFR
143*> \verbatim
144*> DIFR is REAL array,
145*> dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1 and
146*> dimension ( N ) if ICOMPQ = 0.
147*> If ICOMPQ = 1, on exit, DIFL(1:N, I) and DIFR(1:N, 2 * I - 1)
148*> record distances between singular values on the I-th
149*> level and singular values on the (I -1)-th level, and
150*> DIFR(1:N, 2 * I ) contains the normalizing factors for
151*> the right singular vector matrix. See SLASD8 for details.
152*> \endverbatim
153*>
154*> \param[out] Z
155*> \verbatim
156*> Z is REAL array,
157*> dimension ( LDU, NLVL ) if ICOMPQ = 1 and
158*> dimension ( N ) if ICOMPQ = 0.
159*> The first K elements of Z(1, I) contain the components of
160*> the deflation-adjusted updating row vector for subproblems
161*> on the I-th level.
162*> \endverbatim
163*>
164*> \param[out] POLES
165*> \verbatim
166*> POLES is REAL array,
167*> dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced
168*> if ICOMPQ = 0. If ICOMPQ = 1, on exit, POLES(1, 2*I - 1) and
169*> POLES(1, 2*I) contain the new and old singular values
170*> involved in the secular equations on the I-th level.
171*> \endverbatim
172*>
173*> \param[out] GIVPTR
174*> \verbatim
175*> GIVPTR is INTEGER array,
176*> dimension ( N ) if ICOMPQ = 1, and not referenced if
177*> ICOMPQ = 0. If ICOMPQ = 1, on exit, GIVPTR( I ) records
178*> the number of Givens rotations performed on the I-th
179*> problem on the computation tree.
180*> \endverbatim
181*>
182*> \param[out] GIVCOL
183*> \verbatim
184*> GIVCOL is INTEGER array,
185*> dimension ( LDGCOL, 2 * NLVL ) if ICOMPQ = 1, and not
186*> referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
187*> GIVCOL(1, 2 *I - 1) and GIVCOL(1, 2 *I) record the locations
188*> of Givens rotations performed on the I-th level on the
189*> computation tree.
190*> \endverbatim
191*>
192*> \param[in] LDGCOL
193*> \verbatim
194*> LDGCOL is INTEGER, LDGCOL = > N.
195*> The leading dimension of arrays GIVCOL and PERM.
196*> \endverbatim
197*>
198*> \param[out] PERM
199*> \verbatim
200*> PERM is INTEGER array, dimension ( LDGCOL, NLVL )
201*> if ICOMPQ = 1, and not referenced
202*> if ICOMPQ = 0. If ICOMPQ = 1, on exit, PERM(1, I) records
203*> permutations done on the I-th level of the computation tree.
204*> \endverbatim
205*>
206*> \param[out] GIVNUM
207*> \verbatim
208*> GIVNUM is REAL array,
209*> dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not
210*> referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
211*> GIVNUM(1, 2 *I - 1) and GIVNUM(1, 2 *I) record the C- and S-
212*> values of Givens rotations performed on the I-th level on
213*> the computation tree.
214*> \endverbatim
215*>
216*> \param[out] C
217*> \verbatim
218*> C is REAL array,
219*> dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0.
220*> If ICOMPQ = 1 and the I-th subproblem is not square, on exit,
221*> C( I ) contains the C-value of a Givens rotation related to
222*> the right null space of the I-th subproblem.
223*> \endverbatim
224*>
225*> \param[out] S
226*> \verbatim
227*> S is REAL array, dimension ( N ) if
228*> ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1
229*> and the I-th subproblem is not square, on exit, S( I )
230*> contains the S-value of a Givens rotation related to
231*> the right null space of the I-th subproblem.
232*> \endverbatim
233*>
234*> \param[out] WORK
235*> \verbatim
236*> WORK is REAL array, dimension
237*> (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)).
238*> \endverbatim
239*>
240*> \param[out] IWORK
241*> \verbatim
242*> IWORK is INTEGER array, dimension (7*N).
243*> \endverbatim
244*>
245*> \param[out] INFO
246*> \verbatim
247*> INFO is INTEGER
248*> = 0: successful exit.
249*> < 0: if INFO = -i, the i-th argument had an illegal value.
250*> > 0: if INFO = 1, a singular value did not converge
251*> \endverbatim
252*
253* Authors:
254* ========
255*
256*> \author Univ. of Tennessee
257*> \author Univ. of California Berkeley
258*> \author Univ. of Colorado Denver
259*> \author NAG Ltd.
260*
261*> \ingroup OTHERauxiliary
262*
263*> \par Contributors:
264* ==================
265*>
266*> Ming Gu and Huan Ren, Computer Science Division, University of
267*> California at Berkeley, USA
268*>
269* =====================================================================
270 SUBROUTINE slasda( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K,
271 $ DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL,
272 $ PERM, GIVNUM, C, S, WORK, IWORK, INFO )
273*
274* -- LAPACK auxiliary routine --
275* -- LAPACK is a software package provided by Univ. of Tennessee, --
276* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
277*
278* .. Scalar Arguments ..
279 INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE
280* ..
281* .. Array Arguments ..
282 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
283 $ K( * ), PERM( LDGCOL, * )
284 REAL C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ),
285 $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ),
286 $ s( * ), u( ldu, * ), vt( ldu, * ), work( * ),
287 $ z( ldu, * )
288* ..
289*
290* =====================================================================
291*
292* .. Parameters ..
293 REAL ZERO, ONE
294 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
295* ..
296* .. Local Scalars ..
297 INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK,
298 $ J, LF, LL, LVL, LVL2, M, NCC, ND, NDB1, NDIML,
299 $ ndimr, nl, nlf, nlp1, nlvl, nr, nrf, nrp1, nru,
300 $ nwork1, nwork2, smlszp, sqrei, vf, vfi, vl, vli
301 REAL ALPHA, BETA
302* ..
303* .. External Subroutines ..
304 EXTERNAL scopy, slasd6, slasdq, slasdt, slaset, xerbla
305* ..
306* .. Executable Statements ..
307*
308* Test the input parameters.
309*
310 info = 0
311*
312 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
313 info = -1
314 ELSE IF( smlsiz.LT.3 ) THEN
315 info = -2
316 ELSE IF( n.LT.0 ) THEN
317 info = -3
318 ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
319 info = -4
320 ELSE IF( ldu.LT.( n+sqre ) ) THEN
321 info = -8
322 ELSE IF( ldgcol.LT.n ) THEN
323 info = -17
324 END IF
325 IF( info.NE.0 ) THEN
326 CALL xerbla( 'SLASDA', -info )
327 RETURN
328 END IF
329*
330 m = n + sqre
331*
332* If the input matrix is too small, call SLASDQ to find the SVD.
333*
334 IF( n.LE.smlsiz ) THEN
335 IF( icompq.EQ.0 ) THEN
336 CALL slasdq( 'U', sqre, n, 0, 0, 0, d, e, vt, ldu, u, ldu,
337 $ u, ldu, work, info )
338 ELSE
339 CALL slasdq( 'U', sqre, n, m, n, 0, d, e, vt, ldu, u, ldu,
340 $ u, ldu, work, info )
341 END IF
342 RETURN
343 END IF
344*
345* Book-keeping and set up the computation tree.
346*
347 inode = 1
348 ndiml = inode + n
349 ndimr = ndiml + n
350 idxq = ndimr + n
351 iwk = idxq + n
352*
353 ncc = 0
354 nru = 0
355*
356 smlszp = smlsiz + 1
357 vf = 1
358 vl = vf + m
359 nwork1 = vl + m
360 nwork2 = nwork1 + smlszp*smlszp
361*
362 CALL slasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
363 $ iwork( ndimr ), smlsiz )
364*
365* for the nodes on bottom level of the tree, solve
366* their subproblems by SLASDQ.
367*
368 ndb1 = ( nd+1 ) / 2
369 DO 30 i = ndb1, nd
370*
371* IC : center row of each node
372* NL : number of rows of left subproblem
373* NR : number of rows of right subproblem
374* NLF: starting row of the left subproblem
375* NRF: starting row of the right subproblem
376*
377 i1 = i - 1
378 ic = iwork( inode+i1 )
379 nl = iwork( ndiml+i1 )
380 nlp1 = nl + 1
381 nr = iwork( ndimr+i1 )
382 nlf = ic - nl
383 nrf = ic + 1
384 idxqi = idxq + nlf - 2
385 vfi = vf + nlf - 1
386 vli = vl + nlf - 1
387 sqrei = 1
388 IF( icompq.EQ.0 ) THEN
389 CALL slaset( 'A', nlp1, nlp1, zero, one, work( nwork1 ),
390 $ smlszp )
391 CALL slasdq( 'U', sqrei, nl, nlp1, nru, ncc, d( nlf ),
392 $ e( nlf ), work( nwork1 ), smlszp,
393 $ work( nwork2 ), nl, work( nwork2 ), nl,
394 $ work( nwork2 ), info )
395 itemp = nwork1 + nl*smlszp
396 CALL scopy( nlp1, work( nwork1 ), 1, work( vfi ), 1 )
397 CALL scopy( nlp1, work( itemp ), 1, work( vli ), 1 )
398 ELSE
399 CALL slaset( 'A', nl, nl, zero, one, u( nlf, 1 ), ldu )
400 CALL slaset( 'A', nlp1, nlp1, zero, one, vt( nlf, 1 ), ldu )
401 CALL slasdq( 'U', sqrei, nl, nlp1, nl, ncc, d( nlf ),
402 $ e( nlf ), vt( nlf, 1 ), ldu, u( nlf, 1 ), ldu,
403 $ u( nlf, 1 ), ldu, work( nwork1 ), info )
404 CALL scopy( nlp1, vt( nlf, 1 ), 1, work( vfi ), 1 )
405 CALL scopy( nlp1, vt( nlf, nlp1 ), 1, work( vli ), 1 )
406 END IF
407 IF( info.NE.0 ) THEN
408 RETURN
409 END IF
410 DO 10 j = 1, nl
411 iwork( idxqi+j ) = j
412 10 CONTINUE
413 IF( ( i.EQ.nd ) .AND. ( sqre.EQ.0 ) ) THEN
414 sqrei = 0
415 ELSE
416 sqrei = 1
417 END IF
418 idxqi = idxqi + nlp1
419 vfi = vfi + nlp1
420 vli = vli + nlp1
421 nrp1 = nr + sqrei
422 IF( icompq.EQ.0 ) THEN
423 CALL slaset( 'A', nrp1, nrp1, zero, one, work( nwork1 ),
424 $ smlszp )
425 CALL slasdq( 'U', sqrei, nr, nrp1, nru, ncc, d( nrf ),
426 $ e( nrf ), work( nwork1 ), smlszp,
427 $ work( nwork2 ), nr, work( nwork2 ), nr,
428 $ work( nwork2 ), info )
429 itemp = nwork1 + ( nrp1-1 )*smlszp
430 CALL scopy( nrp1, work( nwork1 ), 1, work( vfi ), 1 )
431 CALL scopy( nrp1, work( itemp ), 1, work( vli ), 1 )
432 ELSE
433 CALL slaset( 'A', nr, nr, zero, one, u( nrf, 1 ), ldu )
434 CALL slaset( 'A', nrp1, nrp1, zero, one, vt( nrf, 1 ), ldu )
435 CALL slasdq( 'U', sqrei, nr, nrp1, nr, ncc, d( nrf ),
436 $ e( nrf ), vt( nrf, 1 ), ldu, u( nrf, 1 ), ldu,
437 $ u( nrf, 1 ), ldu, work( nwork1 ), info )
438 CALL scopy( nrp1, vt( nrf, 1 ), 1, work( vfi ), 1 )
439 CALL scopy( nrp1, vt( nrf, nrp1 ), 1, work( vli ), 1 )
440 END IF
441 IF( info.NE.0 ) THEN
442 RETURN
443 END IF
444 DO 20 j = 1, nr
445 iwork( idxqi+j ) = j
446 20 CONTINUE
447 30 CONTINUE
448*
449* Now conquer each subproblem bottom-up.
450*
451 j = 2**nlvl
452 DO 50 lvl = nlvl, 1, -1
453 lvl2 = lvl*2 - 1
454*
455* Find the first node LF and last node LL on
456* the current level LVL.
457*
458 IF( lvl.EQ.1 ) THEN
459 lf = 1
460 ll = 1
461 ELSE
462 lf = 2**( lvl-1 )
463 ll = 2*lf - 1
464 END IF
465 DO 40 i = lf, ll
466 im1 = i - 1
467 ic = iwork( inode+im1 )
468 nl = iwork( ndiml+im1 )
469 nr = iwork( ndimr+im1 )
470 nlf = ic - nl
471 nrf = ic + 1
472 IF( i.EQ.ll ) THEN
473 sqrei = sqre
474 ELSE
475 sqrei = 1
476 END IF
477 vfi = vf + nlf - 1
478 vli = vl + nlf - 1
479 idxqi = idxq + nlf - 1
480 alpha = d( ic )
481 beta = e( ic )
482 IF( icompq.EQ.0 ) THEN
483 CALL slasd6( icompq, nl, nr, sqrei, d( nlf ),
484 $ work( vfi ), work( vli ), alpha, beta,
485 $ iwork( idxqi ), perm, givptr( 1 ), givcol,
486 $ ldgcol, givnum, ldu, poles, difl, difr, z,
487 $ k( 1 ), c( 1 ), s( 1 ), work( nwork1 ),
488 $ iwork( iwk ), info )
489 ELSE
490 j = j - 1
491 CALL slasd6( icompq, nl, nr, sqrei, d( nlf ),
492 $ work( vfi ), work( vli ), alpha, beta,
493 $ iwork( idxqi ), perm( nlf, lvl ),
494 $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
495 $ givnum( nlf, lvl2 ), ldu,
496 $ poles( nlf, lvl2 ), difl( nlf, lvl ),
497 $ difr( nlf, lvl2 ), z( nlf, lvl ), k( j ),
498 $ c( j ), s( j ), work( nwork1 ),
499 $ iwork( iwk ), info )
500 END IF
501 IF( info.NE.0 ) THEN
502 RETURN
503 END IF
504 40 CONTINUE
505 50 CONTINUE
506*
507 RETURN
508*
509* End of SLASDA
510*
511 END
subroutine slasd6(ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, IWORK, INFO)
SLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by...
Definition: slasd6.f:313
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 slasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e....
Definition: slasdq.f:211
subroutine slasdt(N, LVL, ND, INODE, NDIML, NDIMR, MSUB)
SLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.
Definition: slasdt.f:105
subroutine slasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
Definition: slasda.f:273
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82