233 lines
7.6 KiB
FortranFixed
233 lines
7.6 KiB
FortranFixed
SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
|
|
$ IDXQ, IWORK, WORK, INFO )
|
|
*
|
|
* -- LAPACK auxiliary routine (version 3.1) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|
* November 2006
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDU, LDVT, NL, NR, SQRE
|
|
DOUBLE PRECISION ALPHA, BETA
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IDXQ( * ), IWORK( * )
|
|
DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
|
|
* where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
|
|
*
|
|
* A related subroutine DLASD7 handles the case in which the singular
|
|
* values (and the singular vectors in factored form) are desired.
|
|
*
|
|
* DLASD1 computes the SVD as follows:
|
|
*
|
|
* ( D1(in) 0 0 0 )
|
|
* B = U(in) * ( Z1' a Z2' b ) * VT(in)
|
|
* ( 0 0 D2(in) 0 )
|
|
*
|
|
* = U(out) * ( D(out) 0) * VT(out)
|
|
*
|
|
* where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
|
|
* with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
|
|
* elsewhere; and the entry b is empty if SQRE = 0.
|
|
*
|
|
* The left singular vectors of the original matrix are stored in U, and
|
|
* the transpose of the right singular vectors are stored in VT, and the
|
|
* singular values are in D. The algorithm consists of three stages:
|
|
*
|
|
* The first stage consists of deflating the size of the problem
|
|
* when there are multiple singular values or when there are zeros in
|
|
* the Z vector. For each such occurence the dimension of the
|
|
* secular equation problem is reduced by one. This stage is
|
|
* performed by the routine DLASD2.
|
|
*
|
|
* The second stage consists of calculating the updated
|
|
* singular values. This is done by finding the square roots of the
|
|
* roots of the secular equation via the routine DLASD4 (as called
|
|
* by DLASD3). This routine also calculates the singular vectors of
|
|
* the current problem.
|
|
*
|
|
* The final stage consists of computing the updated singular vectors
|
|
* directly using the updated singular values. The singular vectors
|
|
* for the current problem are multiplied with the singular vectors
|
|
* from the overall problem.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* NL (input) INTEGER
|
|
* The row dimension of the upper block. NL >= 1.
|
|
*
|
|
* NR (input) INTEGER
|
|
* The row dimension of the lower block. NR >= 1.
|
|
*
|
|
* SQRE (input) INTEGER
|
|
* = 0: the lower block is an NR-by-NR square matrix.
|
|
* = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
|
|
*
|
|
* The bidiagonal matrix has row dimension N = NL + NR + 1,
|
|
* and column dimension M = N + SQRE.
|
|
*
|
|
* D (input/output) DOUBLE PRECISION array,
|
|
* dimension (N = NL+NR+1).
|
|
* On entry D(1:NL,1:NL) contains the singular values of the
|
|
* upper block; and D(NL+2:N) contains the singular values of
|
|
* the lower block. On exit D(1:N) contains the singular values
|
|
* of the modified matrix.
|
|
*
|
|
* ALPHA (input/output) DOUBLE PRECISION
|
|
* Contains the diagonal element associated with the added row.
|
|
*
|
|
* BETA (input/output) DOUBLE PRECISION
|
|
* Contains the off-diagonal element associated with the added
|
|
* row.
|
|
*
|
|
* U (input/output) DOUBLE PRECISION array, dimension(LDU,N)
|
|
* On entry U(1:NL, 1:NL) contains the left singular vectors of
|
|
* the upper block; U(NL+2:N, NL+2:N) contains the left singular
|
|
* vectors of the lower block. On exit U contains the left
|
|
* singular vectors of the bidiagonal matrix.
|
|
*
|
|
* LDU (input) INTEGER
|
|
* The leading dimension of the array U. LDU >= max( 1, N ).
|
|
*
|
|
* VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
|
|
* where M = N + SQRE.
|
|
* On entry VT(1:NL+1, 1:NL+1)' contains the right singular
|
|
* vectors of the upper block; VT(NL+2:M, NL+2:M)' contains
|
|
* the right singular vectors of the lower block. On exit
|
|
* VT' contains the right singular vectors of the
|
|
* bidiagonal matrix.
|
|
*
|
|
* LDVT (input) INTEGER
|
|
* The leading dimension of the array VT. LDVT >= max( 1, M ).
|
|
*
|
|
* IDXQ (output) INTEGER array, dimension(N)
|
|
* This contains the permutation which will reintegrate the
|
|
* subproblem just solved back into sorted order, i.e.
|
|
* D( IDXQ( I = 1, N ) ) will be in ascending order.
|
|
*
|
|
* IWORK (workspace) INTEGER array, dimension( 4 * N )
|
|
*
|
|
* WORK (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
|
|
*
|
|
* 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 ..
|
|
*
|
|
DOUBLE PRECISION ONE, ZERO
|
|
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
|
|
$ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
|
|
DOUBLE PRECISION ORGNRM
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC ABS, MAX
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
*
|
|
IF( NL.LT.1 ) THEN
|
|
INFO = -1
|
|
ELSE IF( NR.LT.1 ) THEN
|
|
INFO = -2
|
|
ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
|
|
INFO = -3
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DLASD1', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
N = NL + NR + 1
|
|
M = N + SQRE
|
|
*
|
|
* The following values are for bookkeeping purposes only. They are
|
|
* integer pointers which indicate the portion of the workspace
|
|
* used by a particular array in DLASD2 and DLASD3.
|
|
*
|
|
LDU2 = N
|
|
LDVT2 = M
|
|
*
|
|
IZ = 1
|
|
ISIGMA = IZ + M
|
|
IU2 = ISIGMA + N
|
|
IVT2 = IU2 + LDU2*N
|
|
IQ = IVT2 + LDVT2*M
|
|
*
|
|
IDX = 1
|
|
IDXC = IDX + N
|
|
COLTYP = IDXC + N
|
|
IDXP = COLTYP + N
|
|
*
|
|
* Scale.
|
|
*
|
|
ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
|
|
D( NL+1 ) = ZERO
|
|
DO 10 I = 1, N
|
|
IF( ABS( D( I ) ).GT.ORGNRM ) THEN
|
|
ORGNRM = ABS( D( I ) )
|
|
END IF
|
|
10 CONTINUE
|
|
CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
|
|
ALPHA = ALPHA / ORGNRM
|
|
BETA = BETA / ORGNRM
|
|
*
|
|
* Deflate singular values.
|
|
*
|
|
CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
|
|
$ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
|
|
$ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
|
|
$ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
|
|
*
|
|
* Solve Secular Equation and update singular vectors.
|
|
*
|
|
LDQ = K
|
|
CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
|
|
$ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
|
|
$ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
|
|
$ INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Unscale.
|
|
*
|
|
CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
|
|
*
|
|
* Prepare the IDXQ sorting permutation.
|
|
*
|
|
N1 = K
|
|
N2 = N - K
|
|
CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
|
|
*
|
|
RETURN
|
|
*
|
|
* End of DLASD1
|
|
*
|
|
END
|