LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ clalsa()

 subroutine clalsa ( integer ICOMPQ, integer SMLSIZ, integer N, integer NRHS, complex, dimension( ldb, * ) B, integer LDB, complex, dimension( ldbx, * ) BX, integer LDBX, 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( * ) RWORK, integer, dimension( * ) IWORK, integer INFO )

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

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

Purpose:
``` CLALSA is an itermediate 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 orthorgonal
matrices.).

If ICOMPQ = 0, CLALSA applies the inverse of the left singular vector
matrix of an upper bidiagonal matrix to the right hand side; and if
ICOMPQ = 1, CLALSA applies the right singular vector matrix to the
right hand side. The singular vector matrices were generated in
compact form by CLALSA.```
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 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 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 REAL 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 REAL 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 REAL array, dimension ( LDU, NLVL ). where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.``` [in] DIFR ``` DIFR is REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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.```
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 clalsa.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 REAL C( * ), DIFL( LDU, * ), DIFR( LDU, * ),
280 \$ GIVNUM( LDU, * ), POLES( LDU, * ), RWORK( * ),
281 \$ S( * ), U( LDU, * ), VT( LDU, * ), Z( LDU, * )
282 COMPLEX B( LDB, * ), BX( LDBX, * )
283* ..
284*
285* =====================================================================
286*
287* .. Parameters ..
288 REAL ZERO, ONE
289 parameter( zero = 0.0e0, one = 1.0e0 )
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 ccopy, clals0, sgemm, slasdt, xerbla
298* ..
299* .. Intrinsic Functions ..
300 INTRINSIC aimag, cmplx, real
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( 'CLALSA', -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 slasdt( 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 SLASDQ. 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 SGEMM
368* is performed in two steps (real and imaginary parts).
369*
370* CALL SGEMM( '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 ) = real( b( jrow, jcol ) )
378 10 CONTINUE
379 20 CONTINUE
380 CALL sgemm( '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 ) = aimag( b( jrow, jcol ) )
387 30 CONTINUE
388 40 CONTINUE
389 CALL sgemm( '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 ) = cmplx( rwork( jreal ),
399 \$ rwork( jimag ) )
400 50 CONTINUE
401 60 CONTINUE
402*
403* Since B and BX are complex, the following call to SGEMM
404* is performed in two steps (real and imaginary parts).
405*
406* CALL SGEMM( '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 ) = real( b( jrow, jcol ) )
414 70 CONTINUE
415 80 CONTINUE
416 CALL sgemm( '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 ) = aimag( b( jrow, jcol ) )
423 90 CONTINUE
424 100 CONTINUE
425 CALL sgemm( '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 ) = cmplx( 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 ccopy( 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 clals0( 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 clals0( 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 SLASDQ. 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 SGEMM is
551* performed in two steps (real and imaginary parts).
552*
553* CALL SGEMM( '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 ) = real( b( jrow, jcol ) )
561 200 CONTINUE
562 210 CONTINUE
563 CALL sgemm( '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 ) = aimag( b( jrow, jcol ) )
571 220 CONTINUE
572 230 CONTINUE
573 CALL sgemm( '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 ) = cmplx( rwork( jreal ),
583 \$ rwork( jimag ) )
584 240 CONTINUE
585 250 CONTINUE
586*
587* Since B and BX are complex, the following call to SGEMM is
588* performed in two steps (real and imaginary parts).
589*
590* CALL SGEMM( '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 ) = real( b( jrow, jcol ) )
598 260 CONTINUE
599 270 CONTINUE
600 CALL sgemm( '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 ) = aimag( b( jrow, jcol ) )
608 280 CONTINUE
609 290 CONTINUE
610 CALL sgemm( '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 ) = cmplx( 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 CLALSA
631*
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine clals0(ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO)
CLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer...
Definition: clals0.f:270
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
Here is the call graph for this function:
Here is the caller graph for this function: