LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaeda.f
Go to the documentation of this file.
00001 *> \brief \b DLAEDA
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLAEDA + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaeda.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaeda.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaeda.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
00022 *                          GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            CURLVL, CURPBM, INFO, N, TLVLS
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
00029 *      $                   PRMPTR( * ), QPTR( * )
00030 *       DOUBLE PRECISION   GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DLAEDA computes the Z vector corresponding to the merge step in the
00040 *> CURLVLth step of the merge process with TLVLS steps for the CURPBMth
00041 *> problem.
00042 *> \endverbatim
00043 *
00044 *  Arguments:
00045 *  ==========
00046 *
00047 *> \param[in] N
00048 *> \verbatim
00049 *>          N is INTEGER
00050 *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
00051 *> \endverbatim
00052 *>
00053 *> \param[in] TLVLS
00054 *> \verbatim
00055 *>          TLVLS is INTEGER
00056 *>         The total number of merging levels in the overall divide and
00057 *>         conquer tree.
00058 *> \endverbatim
00059 *>
00060 *> \param[in] CURLVL
00061 *> \verbatim
00062 *>          CURLVL is INTEGER
00063 *>         The current level in the overall merge routine,
00064 *>         0 <= curlvl <= tlvls.
00065 *> \endverbatim
00066 *>
00067 *> \param[in] CURPBM
00068 *> \verbatim
00069 *>          CURPBM is INTEGER
00070 *>         The current problem in the current level in the overall
00071 *>         merge routine (counting from upper left to lower right).
00072 *> \endverbatim
00073 *>
00074 *> \param[in] PRMPTR
00075 *> \verbatim
00076 *>          PRMPTR is INTEGER array, dimension (N lg N)
00077 *>         Contains a list of pointers which indicate where in PERM a
00078 *>         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
00079 *>         indicates the size of the permutation and incidentally the
00080 *>         size of the full, non-deflated problem.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] PERM
00084 *> \verbatim
00085 *>          PERM is INTEGER array, dimension (N lg N)
00086 *>         Contains the permutations (from deflation and sorting) to be
00087 *>         applied to each eigenblock.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] GIVPTR
00091 *> \verbatim
00092 *>          GIVPTR is INTEGER array, dimension (N lg N)
00093 *>         Contains a list of pointers which indicate where in GIVCOL a
00094 *>         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
00095 *>         indicates the number of Givens rotations.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] GIVCOL
00099 *> \verbatim
00100 *>          GIVCOL is INTEGER array, dimension (2, N lg N)
00101 *>         Each pair of numbers indicates a pair of columns to take place
00102 *>         in a Givens rotation.
00103 *> \endverbatim
00104 *>
00105 *> \param[in] GIVNUM
00106 *> \verbatim
00107 *>          GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)
00108 *>         Each number indicates the S value to be used in the
00109 *>         corresponding Givens rotation.
00110 *> \endverbatim
00111 *>
00112 *> \param[in] Q
00113 *> \verbatim
00114 *>          Q is DOUBLE PRECISION array, dimension (N**2)
00115 *>         Contains the square eigenblocks from previous levels, the
00116 *>         starting positions for blocks are given by QPTR.
00117 *> \endverbatim
00118 *>
00119 *> \param[in] QPTR
00120 *> \verbatim
00121 *>          QPTR is INTEGER array, dimension (N+2)
00122 *>         Contains a list of pointers which indicate where in Q an
00123 *>         eigenblock is stored.  SQRT( QPTR(i+1) - QPTR(i) ) indicates
00124 *>         the size of the block.
00125 *> \endverbatim
00126 *>
00127 *> \param[out] Z
00128 *> \verbatim
00129 *>          Z is DOUBLE PRECISION array, dimension (N)
00130 *>         On output this vector contains the updating vector (the last
00131 *>         row of the first sub-eigenvector matrix and the first row of
00132 *>         the second sub-eigenvector matrix).
00133 *> \endverbatim
00134 *>
00135 *> \param[out] ZTEMP
00136 *> \verbatim
00137 *>          ZTEMP is DOUBLE PRECISION array, dimension (N)
00138 *> \endverbatim
00139 *>
00140 *> \param[out] INFO
00141 *> \verbatim
00142 *>          INFO is INTEGER
00143 *>          = 0:  successful exit.
00144 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00145 *> \endverbatim
00146 *
00147 *  Authors:
00148 *  ========
00149 *
00150 *> \author Univ. of Tennessee 
00151 *> \author Univ. of California Berkeley 
00152 *> \author Univ. of Colorado Denver 
00153 *> \author NAG Ltd. 
00154 *
00155 *> \date November 2011
00156 *
00157 *> \ingroup auxOTHERcomputational
00158 *
00159 *> \par Contributors:
00160 *  ==================
00161 *>
00162 *> Jeff Rutter, Computer Science Division, University of California
00163 *> at Berkeley, USA
00164 *
00165 *  =====================================================================
00166       SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
00167      $                   GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
00168 *
00169 *  -- LAPACK computational routine (version 3.4.0) --
00170 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00171 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00172 *     November 2011
00173 *
00174 *     .. Scalar Arguments ..
00175       INTEGER            CURLVL, CURPBM, INFO, N, TLVLS
00176 *     ..
00177 *     .. Array Arguments ..
00178       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
00179      $                   PRMPTR( * ), QPTR( * )
00180       DOUBLE PRECISION   GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
00181 *     ..
00182 *
00183 *  =====================================================================
00184 *
00185 *     .. Parameters ..
00186       DOUBLE PRECISION   ZERO, HALF, ONE
00187       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 )
00188 *     ..
00189 *     .. Local Scalars ..
00190       INTEGER            BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
00191      $                   PTR, ZPTR1
00192 *     ..
00193 *     .. External Subroutines ..
00194       EXTERNAL           DCOPY, DGEMV, DROT, XERBLA
00195 *     ..
00196 *     .. Intrinsic Functions ..
00197       INTRINSIC          DBLE, INT, SQRT
00198 *     ..
00199 *     .. Executable Statements ..
00200 *
00201 *     Test the input parameters.
00202 *
00203       INFO = 0
00204 *
00205       IF( N.LT.0 ) THEN
00206          INFO = -1
00207       END IF
00208       IF( INFO.NE.0 ) THEN
00209          CALL XERBLA( 'DLAEDA', -INFO )
00210          RETURN
00211       END IF
00212 *
00213 *     Quick return if possible
00214 *
00215       IF( N.EQ.0 )
00216      $   RETURN
00217 *
00218 *     Determine location of first number in second half.
00219 *
00220       MID = N / 2 + 1
00221 *
00222 *     Gather last/first rows of appropriate eigenblocks into center of Z
00223 *
00224       PTR = 1
00225 *
00226 *     Determine location of lowest level subproblem in the full storage
00227 *     scheme
00228 *
00229       CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1
00230 *
00231 *     Determine size of these matrices.  We add HALF to the value of
00232 *     the SQRT in case the machine underestimates one of these square
00233 *     roots.
00234 *
00235       BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
00236       BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
00237       DO 10 K = 1, MID - BSIZ1 - 1
00238          Z( K ) = ZERO
00239    10 CONTINUE
00240       CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
00241      $            Z( MID-BSIZ1 ), 1 )
00242       CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
00243       DO 20 K = MID + BSIZ2, N
00244          Z( K ) = ZERO
00245    20 CONTINUE
00246 *
00247 *     Loop through remaining levels 1 -> CURLVL applying the Givens
00248 *     rotations and permutation and then multiplying the center matrices
00249 *     against the current Z.
00250 *
00251       PTR = 2**TLVLS + 1
00252       DO 70 K = 1, CURLVL - 1
00253          CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1
00254          PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
00255          PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
00256          ZPTR1 = MID - PSIZ1
00257 *
00258 *       Apply Givens at CURR and CURR+1
00259 *
00260          DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1
00261             CALL DROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1,
00262      $                 Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ),
00263      $                 GIVNUM( 2, I ) )
00264    30    CONTINUE
00265          DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1
00266             CALL DROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1,
00267      $                 Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ),
00268      $                 GIVNUM( 2, I ) )
00269    40    CONTINUE
00270          PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
00271          PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
00272          DO 50 I = 0, PSIZ1 - 1
00273             ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 )
00274    50    CONTINUE
00275          DO 60 I = 0, PSIZ2 - 1
00276             ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 )
00277    60    CONTINUE
00278 *
00279 *        Multiply Blocks at CURR and CURR+1
00280 *
00281 *        Determine size of these matrices.  We add HALF to the value of
00282 *        the SQRT in case the machine underestimates one of these
00283 *        square roots.
00284 *
00285          BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
00286          BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+
00287      $           1 ) ) ) )
00288          IF( BSIZ1.GT.0 ) THEN
00289             CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
00290      $                  BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
00291          END IF
00292          CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
00293      $               1 )
00294          IF( BSIZ2.GT.0 ) THEN
00295             CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
00296      $                  BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
00297          END IF
00298          CALL DCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1,
00299      $               Z( MID+BSIZ2 ), 1 )
00300 *
00301          PTR = PTR + 2**( TLVLS-K )
00302    70 CONTINUE
00303 *
00304       RETURN
00305 *
00306 *     End of DLAEDA
00307 *
00308       END
 All Files Functions