LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slasda()

subroutine slasda ( integer icompq,
integer smlsiz,
integer n,
integer sqre,
real, dimension( * ) d,
real, dimension( * ) e,
real, dimension( ldu, * ) u,
integer ldu,
real, dimension( ldu, * ) vt,
integer, dimension( * ) k,
real, dimension( ldu, * ) difl,
real, dimension( ldu, * ) difr,
real, dimension( ldu, * ) z,
real, dimension( ldu, * ) poles,
integer, dimension( * ) givptr,
integer, dimension( ldgcol, * ) givcol,
integer ldgcol,
integer, dimension( ldgcol, * ) perm,
real, dimension( ldu, * ) givnum,
real, dimension( * ) c,
real, dimension( * ) s,
real, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc.

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

Purpose:
!>
!> Using a divide and conquer approach, SLASDA computes the singular
!> value decomposition (SVD) of a real upper bidiagonal N-by-M matrix
!> B with diagonal D and offdiagonal E, where M = N + SQRE. The
!> algorithm computes the singular values in the SVD B = U * S * VT.
!> The orthogonal matrices U and VT are optionally computed in
!> compact form.
!>
!> A related subroutine, SLASD0, computes the singular values and
!> the singular vectors in explicit form.
!>
Parameters
[in]ICOMPQ
!>          ICOMPQ is INTEGER
!>         Specifies whether singular vectors are to be computed
!>         in compact form, as follows
!>         = 0: Compute singular values only.
!>         = 1: Compute singular vectors of upper bidiagonal
!>              matrix in compact form.
!>
[in]SMLSIZ
!>          SMLSIZ is INTEGER
!>         The maximum size of the subproblems at the bottom of the
!>         computation tree.
!>
[in]N
!>          N is INTEGER
!>         The row dimension of the upper bidiagonal matrix. This is
!>         also the dimension of the main diagonal array D.
!>
[in]SQRE
!>          SQRE is INTEGER
!>         Specifies the column dimension of the bidiagonal matrix.
!>         = 0: The bidiagonal matrix has column dimension M = N;
!>         = 1: The bidiagonal matrix has column dimension M = N + 1.
!>
[in,out]D
!>          D is REAL array, dimension ( N )
!>         On entry D contains the main diagonal of the bidiagonal
!>         matrix. On exit D, if INFO = 0, contains its singular values.
!>
[in]E
!>          E is REAL array, dimension ( M-1 )
!>         Contains the subdiagonal entries of the bidiagonal matrix.
!>         On exit, E has been destroyed.
!>
[out]U
!>          U is REAL array,
!>         dimension ( LDU, SMLSIZ ) if ICOMPQ = 1, and not referenced
!>         if ICOMPQ = 0. If ICOMPQ = 1, on exit, U contains the left
!>         singular vector matrices of all subproblems at the bottom
!>         level.
!>
[in]LDU
!>          LDU is INTEGER, LDU = > N.
!>         The leading dimension of arrays U, VT, DIFL, DIFR, POLES,
!>         GIVNUM, and Z.
!>
[out]VT
!>          VT is REAL array,
!>         dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced
!>         if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT**T contains the right
!>         singular vector matrices of all subproblems at the bottom
!>         level.
!>
[out]K
!>          K is INTEGER array, dimension ( N )
!>         if ICOMPQ = 1 and dimension 1 if ICOMPQ = 0.
!>         If ICOMPQ = 1, on exit, K(I) is the dimension of the I-th
!>         secular equation on the computation tree.
!>
[out]DIFL
!>          DIFL is REAL array, dimension ( LDU, NLVL ),
!>         where NLVL = floor(log_2 (N/SMLSIZ))).
!>
[out]DIFR
!>          DIFR is REAL array,
!>                  dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1 and
!>                  dimension ( N ) if ICOMPQ = 0.
!>         If ICOMPQ = 1, on exit, DIFL(1:N, I) and DIFR(1:N, 2 * I - 1)
!>         record distances between singular values on the I-th
!>         level and singular values on the (I -1)-th level, and
!>         DIFR(1:N, 2 * I ) contains the normalizing factors for
!>         the right singular vector matrix. See SLASD8 for details.
!>
[out]Z
!>          Z is REAL array,
!>                  dimension ( LDU, NLVL ) if ICOMPQ = 1 and
!>                  dimension ( N ) if ICOMPQ = 0.
!>         The first K elements of Z(1, I) contain the components of
!>         the deflation-adjusted updating row vector for subproblems
!>         on the I-th level.
!>
[out]POLES
!>          POLES is REAL array,
!>         dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced
!>         if ICOMPQ = 0. If ICOMPQ = 1, on exit, POLES(1, 2*I - 1) and
!>         POLES(1, 2*I) contain  the new and old singular values
!>         involved in the secular equations on the I-th level.
!>
[out]GIVPTR
!>          GIVPTR is INTEGER array,
!>         dimension ( N ) if ICOMPQ = 1, and not referenced if
!>         ICOMPQ = 0. If ICOMPQ = 1, on exit, GIVPTR( I ) records
!>         the number of Givens rotations performed on the I-th
!>         problem on the computation tree.
!>
[out]GIVCOL
!>          GIVCOL is INTEGER array,
!>         dimension ( LDGCOL, 2 * NLVL ) if ICOMPQ = 1, and not
!>         referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
!>         GIVCOL(1, 2 *I - 1) and GIVCOL(1, 2 *I) record the locations
!>         of Givens rotations performed on the I-th level on the
!>         computation tree.
!>
[in]LDGCOL
!>          LDGCOL is INTEGER, LDGCOL = > N.
!>         The leading dimension of arrays GIVCOL and PERM.
!>
[out]PERM
!>          PERM is INTEGER array, dimension ( LDGCOL, NLVL )
!>         if ICOMPQ = 1, and not referenced
!>         if ICOMPQ = 0. If ICOMPQ = 1, on exit, PERM(1, I) records
!>         permutations done on the I-th level of the computation tree.
!>
[out]GIVNUM
!>          GIVNUM is REAL array,
!>         dimension ( LDU,  2 * NLVL ) if ICOMPQ = 1, and not
!>         referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
!>         GIVNUM(1, 2 *I - 1) and GIVNUM(1, 2 *I) record the C- and S-
!>         values of Givens rotations performed on the I-th level on
!>         the computation tree.
!>
[out]C
!>          C is REAL array,
!>         dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0.
!>         If ICOMPQ = 1 and the I-th subproblem is not square, on exit,
!>         C( I ) contains the C-value of a Givens rotation related to
!>         the right null space of the I-th subproblem.
!>
[out]S
!>          S is REAL array, dimension ( N ) if
!>         ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1
!>         and the I-th subproblem is not square, on exit, S( I )
!>         contains the S-value of a Givens rotation related to
!>         the right null space of the I-th subproblem.
!>
[out]WORK
!>          WORK is REAL array, dimension
!>         (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)).
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (7*N).
!>
[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 268 of file slasda.f.

272*
273* -- LAPACK auxiliary routine --
274* -- LAPACK is a software package provided by Univ. of Tennessee, --
275* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276*
277* .. Scalar Arguments ..
278 INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE
279* ..
280* .. Array Arguments ..
281 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
282 $ K( * ), PERM( LDGCOL, * )
283 REAL C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ),
284 $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ),
285 $ S( * ), U( LDU, * ), VT( LDU, * ), WORK( * ),
286 $ Z( LDU, * )
287* ..
288*
289* =====================================================================
290*
291* .. Parameters ..
292 REAL ZERO, ONE
293 parameter( zero = 0.0e+0, one = 1.0e+0 )
294* ..
295* .. Local Scalars ..
296 INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK,
297 $ J, LF, LL, LVL, LVL2, M, NCC, ND, NDB1, NDIML,
298 $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, NRU,
299 $ NWORK1, NWORK2, SMLSZP, SQREI, VF, VFI, VL, VLI
300 REAL ALPHA, BETA
301* ..
302* .. External Subroutines ..
303 EXTERNAL scopy, slasd6, slasdq, slasdt, slaset,
304 $ 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,
337 $ ldu,
338 $ u, ldu, work, info )
339 ELSE
340 CALL slasdq( 'U', sqre, n, m, n, 0, d, e, vt, ldu, u,
341 $ ldu,
342 $ u, ldu, work, info )
343 END IF
344 RETURN
345 END IF
346*
347* Book-keeping and set up the computation tree.
348*
349 inode = 1
350 ndiml = inode + n
351 ndimr = ndiml + n
352 idxq = ndimr + n
353 iwk = idxq + n
354*
355 ncc = 0
356 nru = 0
357*
358 smlszp = smlsiz + 1
359 vf = 1
360 vl = vf + m
361 nwork1 = vl + m
362 nwork2 = nwork1 + smlszp*smlszp
363*
364 CALL slasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
365 $ iwork( ndimr ), smlsiz )
366*
367* for the nodes on bottom level of the tree, solve
368* their subproblems by SLASDQ.
369*
370 ndb1 = ( nd+1 ) / 2
371 DO 30 i = ndb1, nd
372*
373* IC : center row of each node
374* NL : number of rows of left subproblem
375* NR : number of rows of right subproblem
376* NLF: starting row of the left subproblem
377* NRF: starting row of the right subproblem
378*
379 i1 = i - 1
380 ic = iwork( inode+i1 )
381 nl = iwork( ndiml+i1 )
382 nlp1 = nl + 1
383 nr = iwork( ndimr+i1 )
384 nlf = ic - nl
385 nrf = ic + 1
386 idxqi = idxq + nlf - 2
387 vfi = vf + nlf - 1
388 vli = vl + nlf - 1
389 sqrei = 1
390 IF( icompq.EQ.0 ) THEN
391 CALL slaset( 'A', nlp1, nlp1, zero, one, work( nwork1 ),
392 $ smlszp )
393 CALL slasdq( 'U', sqrei, nl, nlp1, nru, ncc, d( nlf ),
394 $ e( nlf ), work( nwork1 ), smlszp,
395 $ work( nwork2 ), nl, work( nwork2 ), nl,
396 $ work( nwork2 ), info )
397 itemp = nwork1 + nl*smlszp
398 CALL scopy( nlp1, work( nwork1 ), 1, work( vfi ), 1 )
399 CALL scopy( nlp1, work( itemp ), 1, work( vli ), 1 )
400 ELSE
401 CALL slaset( 'A', nl, nl, zero, one, u( nlf, 1 ), ldu )
402 CALL slaset( 'A', nlp1, nlp1, zero, one, vt( nlf, 1 ),
403 $ ldu )
404 CALL slasdq( 'U', sqrei, nl, nlp1, nl, ncc, d( nlf ),
405 $ e( nlf ), vt( nlf, 1 ), ldu, u( nlf, 1 ), ldu,
406 $ u( nlf, 1 ), ldu, work( nwork1 ), info )
407 CALL scopy( nlp1, vt( nlf, 1 ), 1, work( vfi ), 1 )
408 CALL scopy( nlp1, vt( nlf, nlp1 ), 1, work( vli ), 1 )
409 END IF
410 IF( info.NE.0 ) THEN
411 RETURN
412 END IF
413 DO 10 j = 1, nl
414 iwork( idxqi+j ) = j
415 10 CONTINUE
416 IF( ( i.EQ.nd ) .AND. ( sqre.EQ.0 ) ) THEN
417 sqrei = 0
418 ELSE
419 sqrei = 1
420 END IF
421 idxqi = idxqi + nlp1
422 vfi = vfi + nlp1
423 vli = vli + nlp1
424 nrp1 = nr + sqrei
425 IF( icompq.EQ.0 ) THEN
426 CALL slaset( 'A', nrp1, nrp1, zero, one, work( nwork1 ),
427 $ smlszp )
428 CALL slasdq( 'U', sqrei, nr, nrp1, nru, ncc, d( nrf ),
429 $ e( nrf ), work( nwork1 ), smlszp,
430 $ work( nwork2 ), nr, work( nwork2 ), nr,
431 $ work( nwork2 ), info )
432 itemp = nwork1 + ( nrp1-1 )*smlszp
433 CALL scopy( nrp1, work( nwork1 ), 1, work( vfi ), 1 )
434 CALL scopy( nrp1, work( itemp ), 1, work( vli ), 1 )
435 ELSE
436 CALL slaset( 'A', nr, nr, zero, one, u( nrf, 1 ), ldu )
437 CALL slaset( 'A', nrp1, nrp1, zero, one, vt( nrf, 1 ),
438 $ ldu )
439 CALL slasdq( 'U', sqrei, nr, nrp1, nr, ncc, d( nrf ),
440 $ e( nrf ), vt( nrf, 1 ), ldu, u( nrf, 1 ), ldu,
441 $ u( nrf, 1 ), ldu, work( nwork1 ), info )
442 CALL scopy( nrp1, vt( nrf, 1 ), 1, work( vfi ), 1 )
443 CALL scopy( nrp1, vt( nrf, nrp1 ), 1, work( vli ), 1 )
444 END IF
445 IF( info.NE.0 ) THEN
446 RETURN
447 END IF
448 DO 20 j = 1, nr
449 iwork( idxqi+j ) = j
450 20 CONTINUE
451 30 CONTINUE
452*
453* Now conquer each subproblem bottom-up.
454*
455 j = 2**nlvl
456 DO 50 lvl = nlvl, 1, -1
457 lvl2 = lvl*2 - 1
458*
459* Find the first node LF and last node LL on
460* the current level LVL.
461*
462 IF( lvl.EQ.1 ) THEN
463 lf = 1
464 ll = 1
465 ELSE
466 lf = 2**( lvl-1 )
467 ll = 2*lf - 1
468 END IF
469 DO 40 i = lf, ll
470 im1 = i - 1
471 ic = iwork( inode+im1 )
472 nl = iwork( ndiml+im1 )
473 nr = iwork( ndimr+im1 )
474 nlf = ic - nl
475 nrf = ic + 1
476 IF( i.EQ.ll ) THEN
477 sqrei = sqre
478 ELSE
479 sqrei = 1
480 END IF
481 vfi = vf + nlf - 1
482 vli = vl + nlf - 1
483 idxqi = idxq + nlf - 1
484 alpha = d( ic )
485 beta = e( ic )
486 IF( icompq.EQ.0 ) THEN
487 CALL slasd6( icompq, nl, nr, sqrei, d( nlf ),
488 $ work( vfi ), work( vli ), alpha, beta,
489 $ iwork( idxqi ), perm, givptr( 1 ), givcol,
490 $ ldgcol, givnum, ldu, poles, difl, difr, z,
491 $ k( 1 ), c( 1 ), s( 1 ), work( nwork1 ),
492 $ iwork( iwk ), info )
493 ELSE
494 j = j - 1
495 CALL slasd6( icompq, nl, nr, sqrei, d( nlf ),
496 $ work( vfi ), work( vli ), alpha, beta,
497 $ iwork( idxqi ), perm( nlf, lvl ),
498 $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
499 $ givnum( nlf, lvl2 ), ldu,
500 $ poles( nlf, lvl2 ), difl( nlf, lvl ),
501 $ difr( nlf, lvl2 ), z( nlf, lvl ), k( j ),
502 $ c( j ), s( j ), work( nwork1 ),
503 $ iwork( iwk ), info )
504 END IF
505 IF( info.NE.0 ) THEN
506 RETURN
507 END IF
508 40 CONTINUE
509 50 CONTINUE
510*
511 RETURN
512*
513* End of SLASDA
514*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
slasd6
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:312
slasdq
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:210
slasdt
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:103
slaset
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:108
Here is the call graph for this function:
Here is the caller graph for this function: