LAPACK: SRC/dlaeda.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
00002      $                   GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     June 2010
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            CURLVL, CURPBM, INFO, N, TLVLS
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
00014      $                   PRMPTR( * ), QPTR( * )
00015       DOUBLE PRECISION   GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DLAEDA computes the Z vector corresponding to the merge step in the
00022 *  CURLVLth step of the merge process with TLVLS steps for the CURPBMth
00023 *  problem.
00024 *
00025 *  Arguments
00026 *  =========
00027 *
00028 *  N      (input) INTEGER
00029 *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
00030 *
00031 *  TLVLS  (input) INTEGER
00032 *         The total number of merging levels in the overall divide and
00033 *         conquer tree.
00034 *
00035 *  CURLVL (input) INTEGER
00036 *         The current level in the overall merge routine,
00037 *         0 <= curlvl <= tlvls.
00038 *
00039 *  CURPBM (input) INTEGER
00040 *         The current problem in the current level in the overall
00041 *         merge routine (counting from upper left to lower right).
00042 *
00043 *  PRMPTR (input) INTEGER array, dimension (N lg N)
00044 *         Contains a list of pointers which indicate where in PERM a
00045 *         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
00046 *         indicates the size of the permutation and incidentally the
00047 *         size of the full, non-deflated problem.
00048 *
00049 *  PERM   (input) INTEGER array, dimension (N lg N)
00050 *         Contains the permutations (from deflation and sorting) to be
00051 *         applied to each eigenblock.
00052 *
00053 *  GIVPTR (input) INTEGER array, dimension (N lg N)
00054 *         Contains a list of pointers which indicate where in GIVCOL a
00055 *         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
00056 *         indicates the number of Givens rotations.
00057 *
00058 *  GIVCOL (input) INTEGER array, dimension (2, N lg N)
00059 *         Each pair of numbers indicates a pair of columns to take place
00060 *         in a Givens rotation.
00061 *
00062 *  GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
00063 *         Each number indicates the S value to be used in the
00064 *         corresponding Givens rotation.
00065 *
00066 *  Q      (input) DOUBLE PRECISION array, dimension (N**2)
00067 *         Contains the square eigenblocks from previous levels, the
00068 *         starting positions for blocks are given by QPTR.
00069 *
00070 *  QPTR   (input) INTEGER array, dimension (N+2)
00071 *         Contains a list of pointers which indicate where in Q an
00072 *         eigenblock is stored.  SQRT( QPTR(i+1) - QPTR(i) ) indicates
00073 *         the size of the block.
00074 *
00075 *  Z      (output) DOUBLE PRECISION array, dimension (N)
00076 *         On output this vector contains the updating vector (the last
00077 *         row of the first sub-eigenvector matrix and the first row of
00078 *         the second sub-eigenvector matrix).
00079 *
00080 *  ZTEMP  (workspace) DOUBLE PRECISION array, dimension (N)
00081 *
00082 *  INFO   (output) INTEGER
00083 *          = 0:  successful exit.
00084 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00085 *
00086 *  Further Details
00087 *  ===============
00088 *
00089 *  Based on contributions by
00090 *     Jeff Rutter, Computer Science Division, University of California
00091 *     at Berkeley, USA
00092 *
00093 *  =====================================================================
00094 *
00095 *     .. Parameters ..
00096       DOUBLE PRECISION   ZERO, HALF, ONE
00097       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 )
00098 *     ..
00099 *     .. Local Scalars ..
00100       INTEGER            BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
00101      $                   PTR, ZPTR1
00102 *     ..
00103 *     .. External Subroutines ..
00104       EXTERNAL           DCOPY, DGEMV, DROT, XERBLA
00105 *     ..
00106 *     .. Intrinsic Functions ..
00107       INTRINSIC          DBLE, INT, SQRT
00108 *     ..
00109 *     .. Executable Statements ..
00110 *
00111 *     Test the input parameters.
00112 *
00113       INFO = 0
00114 *
00115       IF( N.LT.0 ) THEN
00116          INFO = -1
00117       END IF
00118       IF( INFO.NE.0 ) THEN
00119          CALL XERBLA( 'DLAEDA', -INFO )
00120          RETURN
00121       END IF
00122 *
00123 *     Quick return if possible
00124 *
00125       IF( N.EQ.0 )
00126      $   RETURN
00127 *
00128 *     Determine location of first number in second half.
00129 *
00130       MID = N / 2 + 1
00131 *
00132 *     Gather last/first rows of appropriate eigenblocks into center of Z
00133 *
00134       PTR = 1
00135 *
00136 *     Determine location of lowest level subproblem in the full storage
00137 *     scheme
00138 *
00139       CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1
00140 *
00141 *     Determine size of these matrices.  We add HALF to the value of
00142 *     the SQRT in case the machine underestimates one of these square
00143 *     roots.
00144 *
00145       BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
00146       BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
00147       DO 10 K = 1, MID - BSIZ1 - 1
00148          Z( K ) = ZERO
00149    10 CONTINUE
00150       CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
00151      $            Z( MID-BSIZ1 ), 1 )
00152       CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
00153       DO 20 K = MID + BSIZ2, N
00154          Z( K ) = ZERO
00155    20 CONTINUE
00156 *
00157 *     Loop through remaining levels 1 -> CURLVL applying the Givens
00158 *     rotations and permutation and then multiplying the center matrices
00159 *     against the current Z.
00160 *
00161       PTR = 2**TLVLS + 1
00162       DO 70 K = 1, CURLVL - 1
00163          CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1
00164          PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
00165          PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
00166          ZPTR1 = MID - PSIZ1
00167 *
00168 *       Apply Givens at CURR and CURR+1
00169 *
00170          DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1
00171             CALL DROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1,
00172      $                 Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ),
00173      $                 GIVNUM( 2, I ) )
00174    30    CONTINUE
00175          DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1
00176             CALL DROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1,
00177      $                 Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ),
00178      $                 GIVNUM( 2, I ) )
00179    40    CONTINUE
00180          PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
00181          PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
00182          DO 50 I = 0, PSIZ1 - 1
00183             ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 )
00184    50    CONTINUE
00185          DO 60 I = 0, PSIZ2 - 1
00186             ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 )
00187    60    CONTINUE
00188 *
00189 *        Multiply Blocks at CURR and CURR+1
00190 *
00191 *        Determine size of these matrices.  We add HALF to the value of
00192 *        the SQRT in case the machine underestimates one of these
00193 *        square roots.
00194 *
00195          BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
00196          BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+
00197      $           1 ) ) ) )
00198          IF( BSIZ1.GT.0 ) THEN
00199             CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
00200      $                  BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
00201          END IF
00202          CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
00203      $               1 )
00204          IF( BSIZ2.GT.0 ) THEN
00205             CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
00206      $                  BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
00207          END IF
00208          CALL DCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1,
00209      $               Z( MID+BSIZ2 ), 1 )
00210 *
00211          PTR = PTR + 2**( TLVLS-K )
00212    70 CONTINUE
00213 *
00214       RETURN
00215 *
00216 *     End of DLAEDA
00217 *
00218       END