![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLAEDA 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLAEDA + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaeda.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaeda.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaeda.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLAEDA( 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 * REAL GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SLAEDA 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 REAL 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 REAL 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 REAL 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 REAL 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 SLAEDA( 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 REAL GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * ) 00181 * .. 00182 * 00183 * ===================================================================== 00184 * 00185 * .. Parameters .. 00186 REAL ZERO, HALF, ONE 00187 PARAMETER ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0 ) 00188 * .. 00189 * .. Local Scalars .. 00190 INTEGER BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2, 00191 $ PTR, ZPTR1 00192 * .. 00193 * .. External Subroutines .. 00194 EXTERNAL SCOPY, SGEMV, SROT, XERBLA 00195 * .. 00196 * .. Intrinsic Functions .. 00197 INTRINSIC INT, REAL, 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( 'SLAEDA', -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( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) 00236 BSIZ2 = INT( HALF+SQRT( REAL( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) ) 00237 DO 10 K = 1, MID - BSIZ1 - 1 00238 Z( K ) = ZERO 00239 10 CONTINUE 00240 CALL SCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1, 00241 $ Z( MID-BSIZ1 ), 1 ) 00242 CALL SCOPY( 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 SROT( 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 SROT( 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( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) 00286 BSIZ2 = INT( HALF+SQRT( REAL( QPTR( CURR+2 )-QPTR( CURR+ $ 1 ) ) ) ) 00287 IF( BSIZ1.GT.0 ) THEN 00288 CALL SGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ), 00289 $ BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 ) 00290 END IF 00291 CALL SCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ), 00292 $ 1 ) 00293 IF( BSIZ2.GT.0 ) THEN 00294 CALL SGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ), 00295 $ BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 ) 00296 END IF 00297 CALL SCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1, 00298 $ Z( MID+BSIZ2 ), 1 ) 00299 * 00300 PTR = PTR + 2**( TLVLS-K ) 00301 70 CONTINUE 00302 * 00303 RETURN 00304 * 00305 * End of SLAEDA 00306 * 00307 END 00308