Annotation of rpl/lapack/lapack/dlaed3.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b DLAED3
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLAED3 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed3.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed3.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed3.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
! 22: * CTOT, W, S, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * INTEGER INFO, K, LDQ, N, N1
! 26: * DOUBLE PRECISION RHO
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER CTOT( * ), INDX( * )
! 30: * DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
! 31: * $ S( * ), W( * )
! 32: * ..
! 33: *
! 34: *
! 35: *> \par Purpose:
! 36: * =============
! 37: *>
! 38: *> \verbatim
! 39: *>
! 40: *> DLAED3 finds the roots of the secular equation, as defined by the
! 41: *> values in D, W, and RHO, between 1 and K. It makes the
! 42: *> appropriate calls to DLAED4 and then updates the eigenvectors by
! 43: *> multiplying the matrix of eigenvectors of the pair of eigensystems
! 44: *> being combined by the matrix of eigenvectors of the K-by-K system
! 45: *> which is solved here.
! 46: *>
! 47: *> This code makes very mild assumptions about floating point
! 48: *> arithmetic. It will work on machines with a guard digit in
! 49: *> add/subtract, or on those binary machines without guard digits
! 50: *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
! 51: *> It could conceivably fail on hexadecimal or decimal machines
! 52: *> without guard digits, but we know of none.
! 53: *> \endverbatim
! 54: *
! 55: * Arguments:
! 56: * ==========
! 57: *
! 58: *> \param[in] K
! 59: *> \verbatim
! 60: *> K is INTEGER
! 61: *> The number of terms in the rational function to be solved by
! 62: *> DLAED4. K >= 0.
! 63: *> \endverbatim
! 64: *>
! 65: *> \param[in] N
! 66: *> \verbatim
! 67: *> N is INTEGER
! 68: *> The number of rows and columns in the Q matrix.
! 69: *> N >= K (deflation may result in N>K).
! 70: *> \endverbatim
! 71: *>
! 72: *> \param[in] N1
! 73: *> \verbatim
! 74: *> N1 is INTEGER
! 75: *> The location of the last eigenvalue in the leading submatrix.
! 76: *> min(1,N) <= N1 <= N/2.
! 77: *> \endverbatim
! 78: *>
! 79: *> \param[out] D
! 80: *> \verbatim
! 81: *> D is DOUBLE PRECISION array, dimension (N)
! 82: *> D(I) contains the updated eigenvalues for
! 83: *> 1 <= I <= K.
! 84: *> \endverbatim
! 85: *>
! 86: *> \param[out] Q
! 87: *> \verbatim
! 88: *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
! 89: *> Initially the first K columns are used as workspace.
! 90: *> On output the columns 1 to K contain
! 91: *> the updated eigenvectors.
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in] LDQ
! 95: *> \verbatim
! 96: *> LDQ is INTEGER
! 97: *> The leading dimension of the array Q. LDQ >= max(1,N).
! 98: *> \endverbatim
! 99: *>
! 100: *> \param[in] RHO
! 101: *> \verbatim
! 102: *> RHO is DOUBLE PRECISION
! 103: *> The value of the parameter in the rank one update equation.
! 104: *> RHO >= 0 required.
! 105: *> \endverbatim
! 106: *>
! 107: *> \param[in,out] DLAMDA
! 108: *> \verbatim
! 109: *> DLAMDA is DOUBLE PRECISION array, dimension (K)
! 110: *> The first K elements of this array contain the old roots
! 111: *> of the deflated updating problem. These are the poles
! 112: *> of the secular equation. May be changed on output by
! 113: *> having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
! 114: *> Cray-2, or Cray C-90, as described above.
! 115: *> \endverbatim
! 116: *>
! 117: *> \param[in] Q2
! 118: *> \verbatim
! 119: *> Q2 is DOUBLE PRECISION array, dimension (LDQ2, N)
! 120: *> The first K columns of this matrix contain the non-deflated
! 121: *> eigenvectors for the split problem.
! 122: *> \endverbatim
! 123: *>
! 124: *> \param[in] INDX
! 125: *> \verbatim
! 126: *> INDX is INTEGER array, dimension (N)
! 127: *> The permutation used to arrange the columns of the deflated
! 128: *> Q matrix into three groups (see DLAED2).
! 129: *> The rows of the eigenvectors found by DLAED4 must be likewise
! 130: *> permuted before the matrix multiply can take place.
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[in] CTOT
! 134: *> \verbatim
! 135: *> CTOT is INTEGER array, dimension (4)
! 136: *> A count of the total number of the various types of columns
! 137: *> in Q, as described in INDX. The fourth column type is any
! 138: *> column which has been deflated.
! 139: *> \endverbatim
! 140: *>
! 141: *> \param[in,out] W
! 142: *> \verbatim
! 143: *> W is DOUBLE PRECISION array, dimension (K)
! 144: *> The first K elements of this array contain the components
! 145: *> of the deflation-adjusted updating vector. Destroyed on
! 146: *> output.
! 147: *> \endverbatim
! 148: *>
! 149: *> \param[out] S
! 150: *> \verbatim
! 151: *> S is DOUBLE PRECISION array, dimension (N1 + 1)*K
! 152: *> Will contain the eigenvectors of the repaired matrix which
! 153: *> will be multiplied by the previously accumulated eigenvectors
! 154: *> to update the system.
! 155: *> \endverbatim
! 156: *>
! 157: *> \param[out] INFO
! 158: *> \verbatim
! 159: *> INFO is INTEGER
! 160: *> = 0: successful exit.
! 161: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 162: *> > 0: if INFO = 1, an eigenvalue did not converge
! 163: *> \endverbatim
! 164: *
! 165: * Authors:
! 166: * ========
! 167: *
! 168: *> \author Univ. of Tennessee
! 169: *> \author Univ. of California Berkeley
! 170: *> \author Univ. of Colorado Denver
! 171: *> \author NAG Ltd.
! 172: *
! 173: *> \date November 2011
! 174: *
! 175: *> \ingroup auxOTHERcomputational
! 176: *
! 177: *> \par Contributors:
! 178: * ==================
! 179: *>
! 180: *> Jeff Rutter, Computer Science Division, University of California
! 181: *> at Berkeley, USA \n
! 182: *> Modified by Francoise Tisseur, University of Tennessee
! 183: *>
! 184: * =====================================================================
1.1 bertrand 185: SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
186: $ CTOT, W, S, INFO )
187: *
1.8 ! bertrand 188: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 189: * -- LAPACK is a software package provided by Univ. of Tennessee, --
190: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 191: * November 2011
1.1 bertrand 192: *
193: * .. Scalar Arguments ..
194: INTEGER INFO, K, LDQ, N, N1
195: DOUBLE PRECISION RHO
196: * ..
197: * .. Array Arguments ..
198: INTEGER CTOT( * ), INDX( * )
199: DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
200: $ S( * ), W( * )
201: * ..
202: *
203: * =====================================================================
204: *
205: * .. Parameters ..
206: DOUBLE PRECISION ONE, ZERO
207: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
208: * ..
209: * .. Local Scalars ..
210: INTEGER I, II, IQ2, J, N12, N2, N23
211: DOUBLE PRECISION TEMP
212: * ..
213: * .. External Functions ..
214: DOUBLE PRECISION DLAMC3, DNRM2
215: EXTERNAL DLAMC3, DNRM2
216: * ..
217: * .. External Subroutines ..
218: EXTERNAL DCOPY, DGEMM, DLACPY, DLAED4, DLASET, XERBLA
219: * ..
220: * .. Intrinsic Functions ..
221: INTRINSIC MAX, SIGN, SQRT
222: * ..
223: * .. Executable Statements ..
224: *
225: * Test the input parameters.
226: *
227: INFO = 0
228: *
229: IF( K.LT.0 ) THEN
230: INFO = -1
231: ELSE IF( N.LT.K ) THEN
232: INFO = -2
233: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
234: INFO = -6
235: END IF
236: IF( INFO.NE.0 ) THEN
237: CALL XERBLA( 'DLAED3', -INFO )
238: RETURN
239: END IF
240: *
241: * Quick return if possible
242: *
243: IF( K.EQ.0 )
244: $ RETURN
245: *
246: * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
247: * be computed with high relative accuracy (barring over/underflow).
248: * This is a problem on machines without a guard digit in
249: * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
250: * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
251: * which on any of these machines zeros out the bottommost
252: * bit of DLAMDA(I) if it is 1; this makes the subsequent
253: * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
254: * occurs. On binary machines with a guard digit (almost all
255: * machines) it does not change DLAMDA(I) at all. On hexadecimal
256: * and decimal machines with a guard digit, it slightly
257: * changes the bottommost bits of DLAMDA(I). It does not account
258: * for hexadecimal or decimal machines without guard digits
259: * (we know of none). We use a subroutine call to compute
260: * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
261: * this code.
262: *
263: DO 10 I = 1, K
264: DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
265: 10 CONTINUE
266: *
267: DO 20 J = 1, K
268: CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
269: *
270: * If the zero finder fails, the computation is terminated.
271: *
272: IF( INFO.NE.0 )
273: $ GO TO 120
274: 20 CONTINUE
275: *
276: IF( K.EQ.1 )
277: $ GO TO 110
278: IF( K.EQ.2 ) THEN
279: DO 30 J = 1, K
280: W( 1 ) = Q( 1, J )
281: W( 2 ) = Q( 2, J )
282: II = INDX( 1 )
283: Q( 1, J ) = W( II )
284: II = INDX( 2 )
285: Q( 2, J ) = W( II )
286: 30 CONTINUE
287: GO TO 110
288: END IF
289: *
290: * Compute updated W.
291: *
292: CALL DCOPY( K, W, 1, S, 1 )
293: *
294: * Initialize W(I) = Q(I,I)
295: *
296: CALL DCOPY( K, Q, LDQ+1, W, 1 )
297: DO 60 J = 1, K
298: DO 40 I = 1, J - 1
299: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
300: 40 CONTINUE
301: DO 50 I = J + 1, K
302: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
303: 50 CONTINUE
304: 60 CONTINUE
305: DO 70 I = 1, K
306: W( I ) = SIGN( SQRT( -W( I ) ), S( I ) )
307: 70 CONTINUE
308: *
309: * Compute eigenvectors of the modified rank-1 modification.
310: *
311: DO 100 J = 1, K
312: DO 80 I = 1, K
313: S( I ) = W( I ) / Q( I, J )
314: 80 CONTINUE
315: TEMP = DNRM2( K, S, 1 )
316: DO 90 I = 1, K
317: II = INDX( I )
318: Q( I, J ) = S( II ) / TEMP
319: 90 CONTINUE
320: 100 CONTINUE
321: *
322: * Compute the updated eigenvectors.
323: *
324: 110 CONTINUE
325: *
326: N2 = N - N1
327: N12 = CTOT( 1 ) + CTOT( 2 )
328: N23 = CTOT( 2 ) + CTOT( 3 )
329: *
330: CALL DLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 )
331: IQ2 = N1*N12 + 1
332: IF( N23.NE.0 ) THEN
333: CALL DGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23,
334: $ ZERO, Q( N1+1, 1 ), LDQ )
335: ELSE
336: CALL DLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ )
337: END IF
338: *
339: CALL DLACPY( 'A', N12, K, Q, LDQ, S, N12 )
340: IF( N12.NE.0 ) THEN
341: CALL DGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q,
342: $ LDQ )
343: ELSE
344: CALL DLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ )
345: END IF
346: *
347: *
348: 120 CONTINUE
349: RETURN
350: *
351: * End of DLAED3
352: *
353: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack