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

◆ zlalsa()

subroutine zlalsa ( integer  icompq,
integer  smlsiz,
integer  n,
integer  nrhs,
complex*16, dimension( ldb, * )  b,
integer  ldb,
complex*16, dimension( ldbx, * )  bx,
integer  ldbx,
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( * )  rwork,
integer, dimension( * )  iwork,
integer  info 
)

ZLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.

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

Purpose:
 ZLALSA is an intermediate step in solving the least squares problem
 by computing the SVD of the coefficient matrix in compact form (The
 singular vectors are computed as products of simple orthogonal
 matrices.).

 If ICOMPQ = 0, ZLALSA applies the inverse of the left singular vector
 matrix of an upper bidiagonal matrix to the right hand side; and if
 ICOMPQ = 1, ZLALSA applies the right singular vector matrix to the
 right hand side. The singular vector matrices were generated in
 compact form by ZLALSA.
Parameters
[in]ICOMPQ
          ICOMPQ is INTEGER
         Specifies whether the left or the right singular vector
         matrix is involved.
         = 0: Left singular vector matrix
         = 1: Right singular vector matrix
[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 and column dimensions of the upper bidiagonal matrix.
[in]NRHS
          NRHS is INTEGER
         The number of columns of B and BX. NRHS must be at least 1.
[in,out]B
          B is COMPLEX*16 array, dimension ( LDB, NRHS )
         On input, B contains the right hand sides of the least
         squares problem in rows 1 through M.
         On output, B contains the solution X in rows 1 through N.
[in]LDB
          LDB is INTEGER
         The leading dimension of B in the calling subprogram.
         LDB must be at least max(1,MAX( M, N ) ).
[out]BX
          BX is COMPLEX*16 array, dimension ( LDBX, NRHS )
         On exit, the result of applying the left or right singular
         vector matrix to B.
[in]LDBX
          LDBX is INTEGER
         The leading dimension of BX.
[in]U
          U is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
         On entry, 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.
[in]VT
          VT is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).
         On entry, VT**H contains the right singular vector matrices of
         all subproblems at the bottom level.
[in]K
          K is INTEGER array, dimension ( N ).
[in]DIFL
          DIFL is DOUBLE PRECISION array, dimension ( LDU, NLVL ).
         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
[in]DIFR
          DIFR is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
         distances between singular values on the I-th level and
         singular values on the (I -1)-th level, and DIFR(*, 2 * I)
         record the normalizing factors of the right singular vectors
         matrices of subproblems on I-th level.
[in]Z
          Z is DOUBLE PRECISION array, dimension ( LDU, NLVL ).
         On entry, Z(1, I) contains the components of the deflation-
         adjusted updating row vector for subproblems on the I-th
         level.
[in]POLES
          POLES is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
         On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
         singular values involved in the secular equations on the I-th
         level.
[in]GIVPTR
          GIVPTR is INTEGER array, dimension ( N ).
         On entry, GIVPTR( I ) records the number of Givens
         rotations performed on the I-th problem on the computation
         tree.
[in]GIVCOL
          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records 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.
[in]PERM
          PERM is INTEGER array, dimension ( LDGCOL, NLVL ).
         On entry, PERM(*, I) records permutations done on the I-th
         level of the computation tree.
[in]GIVNUM
          GIVNUM is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
         On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
         values of Givens rotations performed on the I-th level on the
         computation tree.
[in]C
          C is DOUBLE PRECISION array, dimension ( N ).
         On entry, if the I-th subproblem is not square,
         C( I ) contains the C-value of a Givens rotation related to
         the right null space of the I-th subproblem.
[in]S
          S is DOUBLE PRECISION array, dimension ( N ).
         On entry, if the I-th subproblem is not square,
         S( I ) contains the S-value of a Givens rotation related to
         the right null space of the I-th subproblem.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension at least
         MAX( (SMLSZ+1)*NRHS*3, N*(1+NRHS) + 2*NRHS ).
[out]IWORK
          IWORK is INTEGER array, dimension (3*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 263 of file zlalsa.f.

267*
268* -- LAPACK computational routine --
269* -- LAPACK is a software package provided by Univ. of Tennessee, --
270* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
271*
272* .. Scalar Arguments ..
273 INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
274 $ SMLSIZ
275* ..
276* .. Array Arguments ..
277 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
278 $ K( * ), PERM( LDGCOL, * )
279 DOUBLE PRECISION C( * ), DIFL( LDU, * ), DIFR( LDU, * ),
280 $ GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ),
281 $ S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * )
282 COMPLEX*16 B( LDB, * ), BX( LDBX, * )
283* ..
284*
285* =====================================================================
286*
287* .. Parameters ..
288 DOUBLE PRECISION ZERO, ONE
289 parameter( zero = 0.0d0, one = 1.0d0 )
290* ..
291* .. Local Scalars ..
292 INTEGER I, I1, IC, IM1, INODE, J, JCOL, JIMAG, JREAL,
293 $ JROW, LF, LL, LVL, LVL2, ND, NDB1, NDIML,
294 $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQRE
295* ..
296* .. External Subroutines ..
297 EXTERNAL dgemm, dlasdt, xerbla, zcopy, zlals0
298* ..
299* .. Intrinsic Functions ..
300 INTRINSIC dble, dcmplx, dimag
301* ..
302* .. Executable Statements ..
303*
304* Test the input parameters.
305*
306 info = 0
307*
308 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
309 info = -1
310 ELSE IF( smlsiz.LT.3 ) THEN
311 info = -2
312 ELSE IF( n.LT.smlsiz ) THEN
313 info = -3
314 ELSE IF( nrhs.LT.1 ) THEN
315 info = -4
316 ELSE IF( ldb.LT.n ) THEN
317 info = -6
318 ELSE IF( ldbx.LT.n ) THEN
319 info = -8
320 ELSE IF( ldu.LT.n ) THEN
321 info = -10
322 ELSE IF( ldgcol.LT.n ) THEN
323 info = -19
324 END IF
325 IF( info.NE.0 ) THEN
326 CALL xerbla( 'ZLALSA', -info )
327 RETURN
328 END IF
329*
330* Book-keeping and setting up the computation tree.
331*
332 inode = 1
333 ndiml = inode + n
334 ndimr = ndiml + n
335*
336 CALL dlasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
337 $ iwork( ndimr ), smlsiz )
338*
339* The following code applies back the left singular vector factors.
340* For applying back the right singular vector factors, go to 170.
341*
342 IF( icompq.EQ.1 ) THEN
343 GO TO 170
344 END IF
345*
346* The nodes on the bottom level of the tree were solved
347* by DLASDQ. The corresponding left and right singular vector
348* matrices are in explicit form. First apply back the left
349* singular vector matrices.
350*
351 ndb1 = ( nd+1 ) / 2
352 DO 130 i = ndb1, nd
353*
354* IC : center row of each node
355* NL : number of rows of left subproblem
356* NR : number of rows of right subproblem
357* NLF: starting row of the left subproblem
358* NRF: starting row of the right subproblem
359*
360 i1 = i - 1
361 ic = iwork( inode+i1 )
362 nl = iwork( ndiml+i1 )
363 nr = iwork( ndimr+i1 )
364 nlf = ic - nl
365 nrf = ic + 1
366*
367* Since B and BX are complex, the following call to DGEMM
368* is performed in two steps (real and imaginary parts).
369*
370* CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
371* $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
372*
373 j = nl*nrhs*2
374 DO 20 jcol = 1, nrhs
375 DO 10 jrow = nlf, nlf + nl - 1
376 j = j + 1
377 rwork( j ) = dble( b( jrow, jcol ) )
378 10 CONTINUE
379 20 CONTINUE
380 CALL dgemm( 'T', 'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
381 $ rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1 ), nl )
382 j = nl*nrhs*2
383 DO 40 jcol = 1, nrhs
384 DO 30 jrow = nlf, nlf + nl - 1
385 j = j + 1
386 rwork( j ) = dimag( b( jrow, jcol ) )
387 30 CONTINUE
388 40 CONTINUE
389 CALL dgemm( 'T', 'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
390 $ rwork( 1+nl*nrhs*2 ), nl, zero, rwork( 1+nl*nrhs ),
391 $ nl )
392 jreal = 0
393 jimag = nl*nrhs
394 DO 60 jcol = 1, nrhs
395 DO 50 jrow = nlf, nlf + nl - 1
396 jreal = jreal + 1
397 jimag = jimag + 1
398 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
399 $ rwork( jimag ) )
400 50 CONTINUE
401 60 CONTINUE
402*
403* Since B and BX are complex, the following call to DGEMM
404* is performed in two steps (real and imaginary parts).
405*
406* CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
407* $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
408*
409 j = nr*nrhs*2
410 DO 80 jcol = 1, nrhs
411 DO 70 jrow = nrf, nrf + nr - 1
412 j = j + 1
413 rwork( j ) = dble( b( jrow, jcol ) )
414 70 CONTINUE
415 80 CONTINUE
416 CALL dgemm( 'T', 'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
417 $ rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1 ), nr )
418 j = nr*nrhs*2
419 DO 100 jcol = 1, nrhs
420 DO 90 jrow = nrf, nrf + nr - 1
421 j = j + 1
422 rwork( j ) = dimag( b( jrow, jcol ) )
423 90 CONTINUE
424 100 CONTINUE
425 CALL dgemm( 'T', 'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
426 $ rwork( 1+nr*nrhs*2 ), nr, zero, rwork( 1+nr*nrhs ),
427 $ nr )
428 jreal = 0
429 jimag = nr*nrhs
430 DO 120 jcol = 1, nrhs
431 DO 110 jrow = nrf, nrf + nr - 1
432 jreal = jreal + 1
433 jimag = jimag + 1
434 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
435 $ rwork( jimag ) )
436 110 CONTINUE
437 120 CONTINUE
438*
439 130 CONTINUE
440*
441* Next copy the rows of B that correspond to unchanged rows
442* in the bidiagonal matrix to BX.
443*
444 DO 140 i = 1, nd
445 ic = iwork( inode+i-1 )
446 CALL zcopy( nrhs, b( ic, 1 ), ldb, bx( ic, 1 ), ldbx )
447 140 CONTINUE
448*
449* Finally go through the left singular vector matrices of all
450* the other subproblems bottom-up on the tree.
451*
452 j = 2**nlvl
453 sqre = 0
454*
455 DO 160 lvl = nlvl, 1, -1
456 lvl2 = 2*lvl - 1
457*
458* find the first node LF and last node LL on
459* the current level LVL
460*
461 IF( lvl.EQ.1 ) THEN
462 lf = 1
463 ll = 1
464 ELSE
465 lf = 2**( lvl-1 )
466 ll = 2*lf - 1
467 END IF
468 DO 150 i = lf, ll
469 im1 = i - 1
470 ic = iwork( inode+im1 )
471 nl = iwork( ndiml+im1 )
472 nr = iwork( ndimr+im1 )
473 nlf = ic - nl
474 nrf = ic + 1
475 j = j - 1
476 CALL zlals0( icompq, nl, nr, sqre, nrhs, bx( nlf, 1 ), ldbx,
477 $ b( nlf, 1 ), ldb, perm( nlf, lvl ),
478 $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
479 $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
480 $ difl( nlf, lvl ), difr( nlf, lvl2 ),
481 $ z( nlf, lvl ), k( j ), c( j ), s( j ), rwork,
482 $ info )
483 150 CONTINUE
484 160 CONTINUE
485 GO TO 330
486*
487* ICOMPQ = 1: applying back the right singular vector factors.
488*
489 170 CONTINUE
490*
491* First now go through the right singular vector matrices of all
492* the tree nodes top-down.
493*
494 j = 0
495 DO 190 lvl = 1, nlvl
496 lvl2 = 2*lvl - 1
497*
498* Find the first node LF and last node LL on
499* the current level LVL.
500*
501 IF( lvl.EQ.1 ) THEN
502 lf = 1
503 ll = 1
504 ELSE
505 lf = 2**( lvl-1 )
506 ll = 2*lf - 1
507 END IF
508 DO 180 i = ll, lf, -1
509 im1 = i - 1
510 ic = iwork( inode+im1 )
511 nl = iwork( ndiml+im1 )
512 nr = iwork( ndimr+im1 )
513 nlf = ic - nl
514 nrf = ic + 1
515 IF( i.EQ.ll ) THEN
516 sqre = 0
517 ELSE
518 sqre = 1
519 END IF
520 j = j + 1
521 CALL zlals0( icompq, nl, nr, sqre, nrhs, b( nlf, 1 ), ldb,
522 $ bx( nlf, 1 ), ldbx, perm( nlf, lvl ),
523 $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
524 $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
525 $ difl( nlf, lvl ), difr( nlf, lvl2 ),
526 $ z( nlf, lvl ), k( j ), c( j ), s( j ), rwork,
527 $ info )
528 180 CONTINUE
529 190 CONTINUE
530*
531* The nodes on the bottom level of the tree were solved
532* by DLASDQ. The corresponding right singular vector
533* matrices are in explicit form. Apply them back.
534*
535 ndb1 = ( nd+1 ) / 2
536 DO 320 i = ndb1, nd
537 i1 = i - 1
538 ic = iwork( inode+i1 )
539 nl = iwork( ndiml+i1 )
540 nr = iwork( ndimr+i1 )
541 nlp1 = nl + 1
542 IF( i.EQ.nd ) THEN
543 nrp1 = nr
544 ELSE
545 nrp1 = nr + 1
546 END IF
547 nlf = ic - nl
548 nrf = ic + 1
549*
550* Since B and BX are complex, the following call to DGEMM is
551* performed in two steps (real and imaginary parts).
552*
553* CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
554* $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
555*
556 j = nlp1*nrhs*2
557 DO 210 jcol = 1, nrhs
558 DO 200 jrow = nlf, nlf + nlp1 - 1
559 j = j + 1
560 rwork( j ) = dble( b( jrow, jcol ) )
561 200 CONTINUE
562 210 CONTINUE
563 CALL dgemm( 'T', 'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ), ldu,
564 $ rwork( 1+nlp1*nrhs*2 ), nlp1, zero, rwork( 1 ),
565 $ nlp1 )
566 j = nlp1*nrhs*2
567 DO 230 jcol = 1, nrhs
568 DO 220 jrow = nlf, nlf + nlp1 - 1
569 j = j + 1
570 rwork( j ) = dimag( b( jrow, jcol ) )
571 220 CONTINUE
572 230 CONTINUE
573 CALL dgemm( 'T', 'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ), ldu,
574 $ rwork( 1+nlp1*nrhs*2 ), nlp1, zero,
575 $ rwork( 1+nlp1*nrhs ), nlp1 )
576 jreal = 0
577 jimag = nlp1*nrhs
578 DO 250 jcol = 1, nrhs
579 DO 240 jrow = nlf, nlf + nlp1 - 1
580 jreal = jreal + 1
581 jimag = jimag + 1
582 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
583 $ rwork( jimag ) )
584 240 CONTINUE
585 250 CONTINUE
586*
587* Since B and BX are complex, the following call to DGEMM is
588* performed in two steps (real and imaginary parts).
589*
590* CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
591* $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
592*
593 j = nrp1*nrhs*2
594 DO 270 jcol = 1, nrhs
595 DO 260 jrow = nrf, nrf + nrp1 - 1
596 j = j + 1
597 rwork( j ) = dble( b( jrow, jcol ) )
598 260 CONTINUE
599 270 CONTINUE
600 CALL dgemm( 'T', 'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ), ldu,
601 $ rwork( 1+nrp1*nrhs*2 ), nrp1, zero, rwork( 1 ),
602 $ nrp1 )
603 j = nrp1*nrhs*2
604 DO 290 jcol = 1, nrhs
605 DO 280 jrow = nrf, nrf + nrp1 - 1
606 j = j + 1
607 rwork( j ) = dimag( b( jrow, jcol ) )
608 280 CONTINUE
609 290 CONTINUE
610 CALL dgemm( 'T', 'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ), ldu,
611 $ rwork( 1+nrp1*nrhs*2 ), nrp1, zero,
612 $ rwork( 1+nrp1*nrhs ), nrp1 )
613 jreal = 0
614 jimag = nrp1*nrhs
615 DO 310 jcol = 1, nrhs
616 DO 300 jrow = nrf, nrf + nrp1 - 1
617 jreal = jreal + 1
618 jimag = jimag + 1
619 bx( jrow, jcol ) = dcmplx( rwork( jreal ),
620 $ rwork( jimag ) )
621 300 CONTINUE
622 310 CONTINUE
623*
624 320 CONTINUE
625*
626 330 CONTINUE
627*
628 RETURN
629*
630* End of ZLALSA
631*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
zlals0
subroutine zlals0(icompq, nl, nr, sqre, nrhs, b, ldb, bx, ldbx, perm, givptr, givcol, ldgcol, givnum, ldgnum, poles, difl, difr, z, k, c, s, rwork, info)
ZLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer...
Definition zlals0.f:270
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:105
Here is the call graph for this function:
Here is the caller graph for this function: