 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLAEDA( 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 REAL GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * SLAEDA 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) REAL 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) REAL 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) REAL 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) REAL 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 REAL ZERO, HALF, ONE 00097 PARAMETER ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0 ) 00098 * .. 00099 * .. Local Scalars .. 00100 INTEGER BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2, 00101 $ PTR, ZPTR1 00102 * .. 00103 * .. External Subroutines .. 00104 EXTERNAL SCOPY, SGEMV, SROT, XERBLA 00105 * .. 00106 * .. Intrinsic Functions .. 00107 INTRINSIC INT, REAL, 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( 'SLAEDA', -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( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) 00146 BSIZ2 = INT( HALF+SQRT( REAL( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) ) 00147 DO 10 K = 1, MID - BSIZ1 - 1 00148 Z( K ) = ZERO 00149 10 CONTINUE 00150 CALL SCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1, 00151 $ Z( MID-BSIZ1 ), 1 ) 00152 CALL SCOPY( 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 SROT( 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 SROT( 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( REAL( QPTR( CURR+1 )-QPTR( CURR ) ) ) ) 00196 BSIZ2 = INT( HALF+SQRT( REAL( QPTR( CURR+2 )-QPTR( CURR+ $ 1 ) ) ) ) 00197 IF( BSIZ1.GT.0 ) THEN 00198 CALL SGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ), 00199 $ BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 ) 00200 END IF 00201 CALL SCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ), 00202 $ 1 ) 00203 IF( BSIZ2.GT.0 ) THEN 00204 CALL SGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ), 00205 $ BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 ) 00206 END IF 00207 CALL SCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1, 00208 $ Z( MID+BSIZ2 ), 1 ) 00209 * 00210 PTR = PTR + 2**( TLVLS-K ) 00211 70 CONTINUE 00212 * 00213 RETURN 00214 * 00215 * End of SLAEDA 00216 * 00217 END 00218
1.7.3 
