Mercurial > octave-nkf
diff libcruft/lapack/slasda.f @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Sun, 27 Apr 2008 22:34:17 +0200 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/slasda.f Sun Apr 27 22:34:17 2008 +0200 @@ -0,0 +1,389 @@ + SUBROUTINE SLASDA( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, + $ DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, + $ PERM, GIVNUM, C, S, WORK, IWORK, INFO ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE +* .. +* .. Array Arguments .. + INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ), + $ K( * ), PERM( LDGCOL, * ) + REAL C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ), + $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ), + $ S( * ), U( LDU, * ), VT( LDU, * ), WORK( * ), + $ Z( LDU, * ) +* .. +* +* 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. +* +* Arguments +* ========= +* +* ICOMPQ (input) 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. +* +* SMLSIZ (input) INTEGER +* The maximum size of the subproblems at the bottom of the +* computation tree. +* +* N (input) INTEGER +* The row dimension of the upper bidiagonal matrix. This is +* also the dimension of the main diagonal array D. +* +* SQRE (input) 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. +* +* D (input/output) 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. +* +* E (input) REAL array, dimension ( M-1 ) +* Contains the subdiagonal entries of the bidiagonal matrix. +* On exit, E has been destroyed. +* +* U (output) 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. +* +* LDU (input) INTEGER, LDU = > N. +* The leading dimension of arrays U, VT, DIFL, DIFR, POLES, +* GIVNUM, and Z. +* +* VT (output) REAL array, +* dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced +* if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT' contains the right +* singular vector matrices of all subproblems at the bottom +* level. +* +* K (output) 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. +* +* DIFL (output) REAL array, dimension ( LDU, NLVL ), +* where NLVL = floor(log_2 (N/SMLSIZ))). +* +* DIFR (output) 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. +* +* Z (output) 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. +* +* POLES (output) 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. +* +* GIVPTR (output) 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. +* +* GIVCOL (output) 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. +* +* LDGCOL (input) INTEGER, LDGCOL = > N. +* The leading dimension of arrays GIVCOL and PERM. +* +* PERM (output) 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. +* +* GIVNUM (output) 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. +* +* C (output) 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. +* +* S (output) 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. +* +* WORK (workspace) REAL array, dimension +* (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)). +* +* IWORK (workspace) INTEGER array, dimension (7*N). +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* > 0: if INFO = 1, an singular value did not converge +* +* Further Details +* =============== +* +* Based on contributions by +* Ming Gu and Huan Ren, Computer Science Division, University of +* California at Berkeley, USA +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK, + $ J, LF, LL, LVL, LVL2, M, NCC, ND, NDB1, NDIML, + $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, NRU, + $ NWORK1, NWORK2, SMLSZP, SQREI, VF, VFI, VL, VLI + REAL ALPHA, BETA +* .. +* .. External Subroutines .. + EXTERNAL SCOPY, SLASD6, SLASDQ, SLASDT, SLASET, XERBLA +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 +* + IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN + INFO = -1 + ELSE IF( SMLSIZ.LT.3 ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN + INFO = -4 + ELSE IF( LDU.LT.( N+SQRE ) ) THEN + INFO = -8 + ELSE IF( LDGCOL.LT.N ) THEN + INFO = -17 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SLASDA', -INFO ) + RETURN + END IF +* + M = N + SQRE +* +* If the input matrix is too small, call SLASDQ to find the SVD. +* + IF( N.LE.SMLSIZ ) THEN + IF( ICOMPQ.EQ.0 ) THEN + CALL SLASDQ( 'U', SQRE, N, 0, 0, 0, D, E, VT, LDU, U, LDU, + $ U, LDU, WORK, INFO ) + ELSE + CALL SLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDU, U, LDU, + $ U, LDU, WORK, INFO ) + END IF + RETURN + END IF +* +* Book-keeping and set up the computation tree. +* + INODE = 1 + NDIML = INODE + N + NDIMR = NDIML + N + IDXQ = NDIMR + N + IWK = IDXQ + N +* + NCC = 0 + NRU = 0 +* + SMLSZP = SMLSIZ + 1 + VF = 1 + VL = VF + M + NWORK1 = VL + M + NWORK2 = NWORK1 + SMLSZP*SMLSZP +* + CALL SLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), + $ IWORK( NDIMR ), SMLSIZ ) +* +* for the nodes on bottom level of the tree, solve +* their subproblems by SLASDQ. +* + NDB1 = ( ND+1 ) / 2 + DO 30 I = NDB1, ND +* +* IC : center row of each node +* NL : number of rows of left subproblem +* NR : number of rows of right subproblem +* NLF: starting row of the left subproblem +* NRF: starting row of the right subproblem +* + I1 = I - 1 + IC = IWORK( INODE+I1 ) + NL = IWORK( NDIML+I1 ) + NLP1 = NL + 1 + NR = IWORK( NDIMR+I1 ) + NLF = IC - NL + NRF = IC + 1 + IDXQI = IDXQ + NLF - 2 + VFI = VF + NLF - 1 + VLI = VL + NLF - 1 + SQREI = 1 + IF( ICOMPQ.EQ.0 ) THEN + CALL SLASET( 'A', NLP1, NLP1, ZERO, ONE, WORK( NWORK1 ), + $ SMLSZP ) + CALL SLASDQ( 'U', SQREI, NL, NLP1, NRU, NCC, D( NLF ), + $ E( NLF ), WORK( NWORK1 ), SMLSZP, + $ WORK( NWORK2 ), NL, WORK( NWORK2 ), NL, + $ WORK( NWORK2 ), INFO ) + ITEMP = NWORK1 + NL*SMLSZP + CALL SCOPY( NLP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 ) + CALL SCOPY( NLP1, WORK( ITEMP ), 1, WORK( VLI ), 1 ) + ELSE + CALL SLASET( 'A', NL, NL, ZERO, ONE, U( NLF, 1 ), LDU ) + CALL SLASET( 'A', NLP1, NLP1, ZERO, ONE, VT( NLF, 1 ), LDU ) + CALL SLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), + $ E( NLF ), VT( NLF, 1 ), LDU, U( NLF, 1 ), LDU, + $ U( NLF, 1 ), LDU, WORK( NWORK1 ), INFO ) + CALL SCOPY( NLP1, VT( NLF, 1 ), 1, WORK( VFI ), 1 ) + CALL SCOPY( NLP1, VT( NLF, NLP1 ), 1, WORK( VLI ), 1 ) + END IF + IF( INFO.NE.0 ) THEN + RETURN + END IF + DO 10 J = 1, NL + IWORK( IDXQI+J ) = J + 10 CONTINUE + IF( ( I.EQ.ND ) .AND. ( SQRE.EQ.0 ) ) THEN + SQREI = 0 + ELSE + SQREI = 1 + END IF + IDXQI = IDXQI + NLP1 + VFI = VFI + NLP1 + VLI = VLI + NLP1 + NRP1 = NR + SQREI + IF( ICOMPQ.EQ.0 ) THEN + CALL SLASET( 'A', NRP1, NRP1, ZERO, ONE, WORK( NWORK1 ), + $ SMLSZP ) + CALL SLASDQ( 'U', SQREI, NR, NRP1, NRU, NCC, D( NRF ), + $ E( NRF ), WORK( NWORK1 ), SMLSZP, + $ WORK( NWORK2 ), NR, WORK( NWORK2 ), NR, + $ WORK( NWORK2 ), INFO ) + ITEMP = NWORK1 + ( NRP1-1 )*SMLSZP + CALL SCOPY( NRP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 ) + CALL SCOPY( NRP1, WORK( ITEMP ), 1, WORK( VLI ), 1 ) + ELSE + CALL SLASET( 'A', NR, NR, ZERO, ONE, U( NRF, 1 ), LDU ) + CALL SLASET( 'A', NRP1, NRP1, ZERO, ONE, VT( NRF, 1 ), LDU ) + CALL SLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), + $ E( NRF ), VT( NRF, 1 ), LDU, U( NRF, 1 ), LDU, + $ U( NRF, 1 ), LDU, WORK( NWORK1 ), INFO ) + CALL SCOPY( NRP1, VT( NRF, 1 ), 1, WORK( VFI ), 1 ) + CALL SCOPY( NRP1, VT( NRF, NRP1 ), 1, WORK( VLI ), 1 ) + END IF + IF( INFO.NE.0 ) THEN + RETURN + END IF + DO 20 J = 1, NR + IWORK( IDXQI+J ) = J + 20 CONTINUE + 30 CONTINUE +* +* Now conquer each subproblem bottom-up. +* + J = 2**NLVL + DO 50 LVL = NLVL, 1, -1 + LVL2 = LVL*2 - 1 +* +* Find the first node LF and last node LL on +* the current level LVL. +* + IF( LVL.EQ.1 ) THEN + LF = 1 + LL = 1 + ELSE + LF = 2**( LVL-1 ) + LL = 2*LF - 1 + END IF + DO 40 I = LF, LL + IM1 = I - 1 + IC = IWORK( INODE+IM1 ) + NL = IWORK( NDIML+IM1 ) + NR = IWORK( NDIMR+IM1 ) + NLF = IC - NL + NRF = IC + 1 + IF( I.EQ.LL ) THEN + SQREI = SQRE + ELSE + SQREI = 1 + END IF + VFI = VF + NLF - 1 + VLI = VL + NLF - 1 + IDXQI = IDXQ + NLF - 1 + ALPHA = D( IC ) + BETA = E( IC ) + IF( ICOMPQ.EQ.0 ) THEN + CALL SLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ), + $ WORK( VFI ), WORK( VLI ), ALPHA, BETA, + $ IWORK( IDXQI ), PERM, GIVPTR( 1 ), GIVCOL, + $ LDGCOL, GIVNUM, LDU, POLES, DIFL, DIFR, Z, + $ K( 1 ), C( 1 ), S( 1 ), WORK( NWORK1 ), + $ IWORK( IWK ), INFO ) + ELSE + J = J - 1 + CALL SLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ), + $ WORK( VFI ), WORK( VLI ), ALPHA, BETA, + $ IWORK( IDXQI ), PERM( NLF, LVL ), + $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, + $ GIVNUM( NLF, LVL2 ), LDU, + $ POLES( NLF, LVL2 ), DIFL( NLF, LVL ), + $ DIFR( NLF, LVL2 ), Z( NLF, LVL ), K( J ), + $ C( J ), S( J ), WORK( NWORK1 ), + $ IWORK( IWK ), INFO ) + END IF + IF( INFO.NE.0 ) THEN + RETURN + END IF + 40 CONTINUE + 50 CONTINUE +* + RETURN +* +* End of SLASDA +* + END