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

◆ dlasda()

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

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

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

Purpose:
!>
!> Using a divide and conquer approach, DLASDA 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, DLASD0, 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension ( M-1 )
!>         Contains the subdiagonal entries of the bidiagonal matrix.
!>         On exit, E has been destroyed.
!>
[out]U
!>          U is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension ( LDU, NLVL ),
!>         where NLVL = floor(log_2 (N/SMLSIZ))).
!>
[out]DIFR
!>          DIFR is DOUBLE PRECISION 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 DLASD8 for details.
!>
[out]Z
!>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dlasda.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 DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE
293 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION ALPHA, BETA
301* ..
302* .. External Subroutines ..
303 EXTERNAL dcopy, dlasd6, dlasdq, dlasdt, dlaset,
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( 'DLASDA', -info )
327 RETURN
328 END IF
329*
330 m = n + sqre
331*
332* If the input matrix is too small, call DLASDQ to find the SVD.
333*
334 IF( n.LE.smlsiz ) THEN
335 IF( icompq.EQ.0 ) THEN
336 CALL dlasdq( 'U', sqre, n, 0, 0, 0, d, e, vt, ldu, u,
337 $ ldu,
338 $ u, ldu, work, info )
339 ELSE
340 CALL dlasdq( '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 dlasdt( 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 DLASDQ.
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 dlaset( 'A', nlp1, nlp1, zero, one, work( nwork1 ),
392 $ smlszp )
393 CALL dlasdq( '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 dcopy( nlp1, work( nwork1 ), 1, work( vfi ), 1 )
399 CALL dcopy( nlp1, work( itemp ), 1, work( vli ), 1 )
400 ELSE
401 CALL dlaset( 'A', nl, nl, zero, one, u( nlf, 1 ), ldu )
402 CALL dlaset( 'A', nlp1, nlp1, zero, one, vt( nlf, 1 ),
403 $ ldu )
404 CALL dlasdq( '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 dcopy( nlp1, vt( nlf, 1 ), 1, work( vfi ), 1 )
408 CALL dcopy( 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 dlaset( 'A', nrp1, nrp1, zero, one, work( nwork1 ),
427 $ smlszp )
428 CALL dlasdq( '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 dcopy( nrp1, work( nwork1 ), 1, work( vfi ), 1 )
434 CALL dcopy( nrp1, work( itemp ), 1, work( vli ), 1 )
435 ELSE
436 CALL dlaset( 'A', nr, nr, zero, one, u( nrf, 1 ), ldu )
437 CALL dlaset( 'A', nrp1, nrp1, zero, one, vt( nrf, 1 ),
438 $ ldu )
439 CALL dlasdq( '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 dcopy( nrp1, vt( nrf, 1 ), 1, work( vfi ), 1 )
443 CALL dcopy( 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 dlasd6( 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 dlasd6( 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 DLASDA
514*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlasd6
subroutine dlasd6(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)
DLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by...
Definition dlasd6.f:312
dlasdq
subroutine dlasdq(uplo, sqre, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e....
Definition dlasdq.f:210
dlasdt
subroutine dlasdt(n, lvl, nd, inode, ndiml, ndimr, msub)
DLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.
Definition dlasdt.f:103
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:108
Here is the call graph for this function:
Here is the caller graph for this function: