218 lines
7.3 KiB
FortranFixed
218 lines
7.3 KiB
FortranFixed
SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
|
|
$ GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.1) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|
* November 2006
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER CURLVL, CURPBM, INFO, N, TLVLS
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
|
|
$ PRMPTR( * ), QPTR( * )
|
|
DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DLAEDA computes the Z vector corresponding to the merge step in the
|
|
* CURLVLth step of the merge process with TLVLS steps for the CURPBMth
|
|
* problem.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* N (input) INTEGER
|
|
* The dimension of the symmetric tridiagonal matrix. N >= 0.
|
|
*
|
|
* TLVLS (input) INTEGER
|
|
* The total number of merging levels in the overall divide and
|
|
* conquer tree.
|
|
*
|
|
* CURLVL (input) INTEGER
|
|
* The current level in the overall merge routine,
|
|
* 0 <= curlvl <= tlvls.
|
|
*
|
|
* CURPBM (input) INTEGER
|
|
* The current problem in the current level in the overall
|
|
* merge routine (counting from upper left to lower right).
|
|
*
|
|
* PRMPTR (input) INTEGER array, dimension (N lg N)
|
|
* Contains a list of pointers which indicate where in PERM a
|
|
* level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
|
|
* indicates the size of the permutation and incidentally the
|
|
* size of the full, non-deflated problem.
|
|
*
|
|
* PERM (input) INTEGER array, dimension (N lg N)
|
|
* Contains the permutations (from deflation and sorting) to be
|
|
* applied to each eigenblock.
|
|
*
|
|
* GIVPTR (input) INTEGER array, dimension (N lg N)
|
|
* Contains a list of pointers which indicate where in GIVCOL a
|
|
* level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
|
|
* indicates the number of Givens rotations.
|
|
*
|
|
* GIVCOL (input) INTEGER array, dimension (2, N lg N)
|
|
* Each pair of numbers indicates a pair of columns to take place
|
|
* in a Givens rotation.
|
|
*
|
|
* GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
|
|
* Each number indicates the S value to be used in the
|
|
* corresponding Givens rotation.
|
|
*
|
|
* Q (input) DOUBLE PRECISION array, dimension (N**2)
|
|
* Contains the square eigenblocks from previous levels, the
|
|
* starting positions for blocks are given by QPTR.
|
|
*
|
|
* QPTR (input) INTEGER array, dimension (N+2)
|
|
* Contains a list of pointers which indicate where in Q an
|
|
* eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates
|
|
* the size of the block.
|
|
*
|
|
* Z (output) DOUBLE PRECISION array, dimension (N)
|
|
* On output this vector contains the updating vector (the last
|
|
* row of the first sub-eigenvector matrix and the first row of
|
|
* the second sub-eigenvector matrix).
|
|
*
|
|
* ZTEMP (workspace) DOUBLE PRECISION array, dimension (N)
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit.
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value.
|
|
*
|
|
* Further Details
|
|
* ===============
|
|
*
|
|
* Based on contributions by
|
|
* Jeff Rutter, Computer Science Division, University of California
|
|
* at Berkeley, USA
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ZERO, HALF, ONE
|
|
PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
|
|
$ PTR, ZPTR1
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DCOPY, DGEMV, DROT, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC DBLE, INT, SQRT
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
*
|
|
IF( N.LT.0 ) THEN
|
|
INFO = -1
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DLAEDA', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* Determine location of first number in second half.
|
|
*
|
|
MID = N / 2 + 1
|
|
*
|
|
* Gather last/first rows of appropriate eigenblocks into center of Z
|
|
*
|
|
PTR = 1
|
|
*
|
|
* Determine location of lowest level subproblem in the full storage
|
|
* scheme
|
|
*
|
|
CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1
|
|
*
|
|
* Determine size of these matrices. We add HALF to the value of
|
|
* the SQRT in case the machine underestimates one of these square
|
|
* roots.
|
|
*
|
|
BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
|
|
BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
|
|
DO 10 K = 1, MID - BSIZ1 - 1
|
|
Z( K ) = ZERO
|
|
10 CONTINUE
|
|
CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
|
|
$ Z( MID-BSIZ1 ), 1 )
|
|
CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
|
|
DO 20 K = MID + BSIZ2, N
|
|
Z( K ) = ZERO
|
|
20 CONTINUE
|
|
*
|
|
* Loop thru remaining levels 1 -> CURLVL applying the Givens
|
|
* rotations and permutation and then multiplying the center matrices
|
|
* against the current Z.
|
|
*
|
|
PTR = 2**TLVLS + 1
|
|
DO 70 K = 1, CURLVL - 1
|
|
CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1
|
|
PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
|
|
PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
|
|
ZPTR1 = MID - PSIZ1
|
|
*
|
|
* Apply Givens at CURR and CURR+1
|
|
*
|
|
DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1
|
|
CALL DROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1,
|
|
$ Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ),
|
|
$ GIVNUM( 2, I ) )
|
|
30 CONTINUE
|
|
DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1
|
|
CALL DROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1,
|
|
$ Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ),
|
|
$ GIVNUM( 2, I ) )
|
|
40 CONTINUE
|
|
PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
|
|
PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
|
|
DO 50 I = 0, PSIZ1 - 1
|
|
ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 )
|
|
50 CONTINUE
|
|
DO 60 I = 0, PSIZ2 - 1
|
|
ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 )
|
|
60 CONTINUE
|
|
*
|
|
* Multiply Blocks at CURR and CURR+1
|
|
*
|
|
* Determine size of these matrices. We add HALF to the value of
|
|
* the SQRT in case the machine underestimates one of these
|
|
* square roots.
|
|
*
|
|
BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
|
|
BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+
|
|
$ 1 ) ) ) )
|
|
IF( BSIZ1.GT.0 ) THEN
|
|
CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
|
|
$ BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
|
|
END IF
|
|
CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
|
|
$ 1 )
|
|
IF( BSIZ2.GT.0 ) THEN
|
|
CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
|
|
$ BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
|
|
END IF
|
|
CALL DCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1,
|
|
$ Z( MID+BSIZ2 ), 1 )
|
|
*
|
|
PTR = PTR + 2**( TLVLS-K )
|
|
70 CONTINUE
|
|
*
|
|
RETURN
|
|
*
|
|
* End of DLAEDA
|
|
*
|
|
END
|