LAPACK 3.12.0
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 270 of file dlasda.f.

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 DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE
294 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION ALPHA, BETA
302* ..
303* .. External Subroutines ..
304 EXTERNAL dcopy, dlasd6, dlasdq, dlasdt, dlaset, 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, ldu,
337 $ u, ldu, work, info )
338 ELSE
339 CALL dlasdq( '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 dlasdt( 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 DLASDQ.
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 dlaset( 'A', nlp1, nlp1, zero, one, work( nwork1 ),
390 $ smlszp )
391 CALL dlasdq( '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 dcopy( nlp1, work( nwork1 ), 1, work( vfi ), 1 )
397 CALL dcopy( nlp1, work( itemp ), 1, work( vli ), 1 )
398 ELSE
399 CALL dlaset( 'A', nl, nl, zero, one, u( nlf, 1 ), ldu )
400 CALL dlaset( 'A', nlp1, nlp1, zero, one, vt( nlf, 1 ), ldu )
401 CALL dlasdq( '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 dcopy( nlp1, vt( nlf, 1 ), 1, work( vfi ), 1 )
405 CALL dcopy( 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 dlaset( 'A', nrp1, nrp1, zero, one, work( nwork1 ),
424 $ smlszp )
425 CALL dlasdq( '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 dcopy( nrp1, work( nwork1 ), 1, work( vfi ), 1 )
431 CALL dcopy( nrp1, work( itemp ), 1, work( vli ), 1 )
432 ELSE
433 CALL dlaset( 'A', nr, nr, zero, one, u( nrf, 1 ), ldu )
434 CALL dlaset( 'A', nrp1, nrp1, zero, one, vt( nrf, 1 ), ldu )
435 CALL dlasdq( '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 dcopy( nrp1, vt( nrf, 1 ), 1, work( vfi ), 1 )
439 CALL dcopy( 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 dlasd6( 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 dlasd6( 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 DLASDA
510*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
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:313
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:211
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:105
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:110
Here is the call graph for this function:
Here is the caller graph for this function: