Annotation of rpl/lapack/lapack/dlaeda.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLAEDA
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLAEDA + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaeda.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaeda.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaeda.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
! 22: * GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * INTEGER CURLVL, CURPBM, INFO, N, TLVLS
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * INTEGER GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
! 29: * $ PRMPTR( * ), QPTR( * )
! 30: * DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DLAEDA computes the Z vector corresponding to the merge step in the
! 40: *> CURLVLth step of the merge process with TLVLS steps for the CURPBMth
! 41: *> problem.
! 42: *> \endverbatim
! 43: *
! 44: * Arguments:
! 45: * ==========
! 46: *
! 47: *> \param[in] N
! 48: *> \verbatim
! 49: *> N is INTEGER
! 50: *> The dimension of the symmetric tridiagonal matrix. N >= 0.
! 51: *> \endverbatim
! 52: *>
! 53: *> \param[in] TLVLS
! 54: *> \verbatim
! 55: *> TLVLS is INTEGER
! 56: *> The total number of merging levels in the overall divide and
! 57: *> conquer tree.
! 58: *> \endverbatim
! 59: *>
! 60: *> \param[in] CURLVL
! 61: *> \verbatim
! 62: *> CURLVL is INTEGER
! 63: *> The current level in the overall merge routine,
! 64: *> 0 <= curlvl <= tlvls.
! 65: *> \endverbatim
! 66: *>
! 67: *> \param[in] CURPBM
! 68: *> \verbatim
! 69: *> CURPBM is INTEGER
! 70: *> The current problem in the current level in the overall
! 71: *> merge routine (counting from upper left to lower right).
! 72: *> \endverbatim
! 73: *>
! 74: *> \param[in] PRMPTR
! 75: *> \verbatim
! 76: *> PRMPTR is INTEGER array, dimension (N lg N)
! 77: *> Contains a list of pointers which indicate where in PERM a
! 78: *> level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
! 79: *> indicates the size of the permutation and incidentally the
! 80: *> size of the full, non-deflated problem.
! 81: *> \endverbatim
! 82: *>
! 83: *> \param[in] PERM
! 84: *> \verbatim
! 85: *> PERM is INTEGER array, dimension (N lg N)
! 86: *> Contains the permutations (from deflation and sorting) to be
! 87: *> applied to each eigenblock.
! 88: *> \endverbatim
! 89: *>
! 90: *> \param[in] GIVPTR
! 91: *> \verbatim
! 92: *> GIVPTR is INTEGER array, dimension (N lg N)
! 93: *> Contains a list of pointers which indicate where in GIVCOL a
! 94: *> level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
! 95: *> indicates the number of Givens rotations.
! 96: *> \endverbatim
! 97: *>
! 98: *> \param[in] GIVCOL
! 99: *> \verbatim
! 100: *> GIVCOL is INTEGER array, dimension (2, N lg N)
! 101: *> Each pair of numbers indicates a pair of columns to take place
! 102: *> in a Givens rotation.
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[in] GIVNUM
! 106: *> \verbatim
! 107: *> GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)
! 108: *> Each number indicates the S value to be used in the
! 109: *> corresponding Givens rotation.
! 110: *> \endverbatim
! 111: *>
! 112: *> \param[in] Q
! 113: *> \verbatim
! 114: *> Q is DOUBLE PRECISION array, dimension (N**2)
! 115: *> Contains the square eigenblocks from previous levels, the
! 116: *> starting positions for blocks are given by QPTR.
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[in] QPTR
! 120: *> \verbatim
! 121: *> QPTR is INTEGER array, dimension (N+2)
! 122: *> Contains a list of pointers which indicate where in Q an
! 123: *> eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates
! 124: *> the size of the block.
! 125: *> \endverbatim
! 126: *>
! 127: *> \param[out] Z
! 128: *> \verbatim
! 129: *> Z is DOUBLE PRECISION array, dimension (N)
! 130: *> On output this vector contains the updating vector (the last
! 131: *> row of the first sub-eigenvector matrix and the first row of
! 132: *> the second sub-eigenvector matrix).
! 133: *> \endverbatim
! 134: *>
! 135: *> \param[out] ZTEMP
! 136: *> \verbatim
! 137: *> ZTEMP is DOUBLE PRECISION array, dimension (N)
! 138: *> \endverbatim
! 139: *>
! 140: *> \param[out] INFO
! 141: *> \verbatim
! 142: *> INFO is INTEGER
! 143: *> = 0: successful exit.
! 144: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 145: *> \endverbatim
! 146: *
! 147: * Authors:
! 148: * ========
! 149: *
! 150: *> \author Univ. of Tennessee
! 151: *> \author Univ. of California Berkeley
! 152: *> \author Univ. of Colorado Denver
! 153: *> \author NAG Ltd.
! 154: *
! 155: *> \date November 2011
! 156: *
! 157: *> \ingroup auxOTHERcomputational
! 158: *
! 159: *> \par Contributors:
! 160: * ==================
! 161: *>
! 162: *> Jeff Rutter, Computer Science Division, University of California
! 163: *> at Berkeley, USA
! 164: *
! 165: * =====================================================================
1.1 bertrand 166: SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
167: $ GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
168: *
1.9 ! bertrand 169: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 170: * -- LAPACK is a software package provided by Univ. of Tennessee, --
171: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 172: * November 2011
1.1 bertrand 173: *
174: * .. Scalar Arguments ..
175: INTEGER CURLVL, CURPBM, INFO, N, TLVLS
176: * ..
177: * .. Array Arguments ..
178: INTEGER GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
179: $ PRMPTR( * ), QPTR( * )
180: DOUBLE PRECISION GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
181: * ..
182: *
183: * =====================================================================
184: *
185: * .. Parameters ..
186: DOUBLE PRECISION ZERO, HALF, ONE
187: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 )
188: * ..
189: * .. Local Scalars ..
190: INTEGER BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
191: $ PTR, ZPTR1
192: * ..
193: * .. External Subroutines ..
194: EXTERNAL DCOPY, DGEMV, DROT, XERBLA
195: * ..
196: * .. Intrinsic Functions ..
197: INTRINSIC DBLE, INT, SQRT
198: * ..
199: * .. Executable Statements ..
200: *
201: * Test the input parameters.
202: *
203: INFO = 0
204: *
205: IF( N.LT.0 ) THEN
206: INFO = -1
207: END IF
208: IF( INFO.NE.0 ) THEN
209: CALL XERBLA( 'DLAEDA', -INFO )
210: RETURN
211: END IF
212: *
213: * Quick return if possible
214: *
215: IF( N.EQ.0 )
216: $ RETURN
217: *
218: * Determine location of first number in second half.
219: *
220: MID = N / 2 + 1
221: *
222: * Gather last/first rows of appropriate eigenblocks into center of Z
223: *
224: PTR = 1
225: *
226: * Determine location of lowest level subproblem in the full storage
227: * scheme
228: *
229: CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1
230: *
231: * Determine size of these matrices. We add HALF to the value of
232: * the SQRT in case the machine underestimates one of these square
233: * roots.
234: *
235: BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
236: BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
237: DO 10 K = 1, MID - BSIZ1 - 1
238: Z( K ) = ZERO
239: 10 CONTINUE
240: CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
241: $ Z( MID-BSIZ1 ), 1 )
242: CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
243: DO 20 K = MID + BSIZ2, N
244: Z( K ) = ZERO
245: 20 CONTINUE
246: *
1.5 bertrand 247: * Loop through remaining levels 1 -> CURLVL applying the Givens
1.1 bertrand 248: * rotations and permutation and then multiplying the center matrices
249: * against the current Z.
250: *
251: PTR = 2**TLVLS + 1
252: DO 70 K = 1, CURLVL - 1
253: CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1
254: PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
255: PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
256: ZPTR1 = MID - PSIZ1
257: *
258: * Apply Givens at CURR and CURR+1
259: *
260: DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1
261: CALL DROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1,
262: $ Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ),
263: $ GIVNUM( 2, I ) )
264: 30 CONTINUE
265: DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1
266: CALL DROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1,
267: $ Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ),
268: $ GIVNUM( 2, I ) )
269: 40 CONTINUE
270: PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
271: PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
272: DO 50 I = 0, PSIZ1 - 1
273: ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 )
274: 50 CONTINUE
275: DO 60 I = 0, PSIZ2 - 1
276: ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 )
277: 60 CONTINUE
278: *
279: * Multiply Blocks at CURR and CURR+1
280: *
281: * Determine size of these matrices. We add HALF to the value of
282: * the SQRT in case the machine underestimates one of these
283: * square roots.
284: *
285: BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
286: BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+
287: $ 1 ) ) ) )
288: IF( BSIZ1.GT.0 ) THEN
289: CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
290: $ BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
291: END IF
292: CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
293: $ 1 )
294: IF( BSIZ2.GT.0 ) THEN
295: CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
296: $ BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
297: END IF
298: CALL DCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1,
299: $ Z( MID+BSIZ2 ), 1 )
300: *
301: PTR = PTR + 2**( TLVLS-K )
302: 70 CONTINUE
303: *
304: RETURN
305: *
306: * End of DLAEDA
307: *
308: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaeda.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack