1: *> \brief \b DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLATDF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlatdf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlatdf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlatdf.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
22: * JPIV )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER IJOB, LDZ, N
26: * DOUBLE PRECISION RDSCAL, RDSUM
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IPIV( * ), JPIV( * )
30: * DOUBLE PRECISION RHS( * ), Z( LDZ, * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DLATDF uses the LU factorization of the n-by-n matrix Z computed by
40: *> DGETC2 and computes a contribution to the reciprocal Dif-estimate
41: *> by solving Z * x = b for x, and choosing the r.h.s. b such that
42: *> the norm of x is as large as possible. On entry RHS = b holds the
43: *> contribution from earlier solved sub-systems, and on return RHS = x.
44: *>
45: *> The factorization of Z returned by DGETC2 has the form Z = P*L*U*Q,
46: *> where P and Q are permutation matrices. L is lower triangular with
47: *> unit diagonal elements and U is upper triangular.
48: *> \endverbatim
49: *
50: * Arguments:
51: * ==========
52: *
53: *> \param[in] IJOB
54: *> \verbatim
55: *> IJOB is INTEGER
56: *> IJOB = 2: First compute an approximative null-vector e
57: *> of Z using DGECON, e is normalized and solve for
58: *> Zx = +-e - f with the sign giving the greater value
59: *> of 2-norm(x). About 5 times as expensive as Default.
60: *> IJOB .ne. 2: Local look ahead strategy where all entries of
61: *> the r.h.s. b is chosen as either +1 or -1 (Default).
62: *> \endverbatim
63: *>
64: *> \param[in] N
65: *> \verbatim
66: *> N is INTEGER
67: *> The number of columns of the matrix Z.
68: *> \endverbatim
69: *>
70: *> \param[in] Z
71: *> \verbatim
72: *> Z is DOUBLE PRECISION array, dimension (LDZ, N)
73: *> On entry, the LU part of the factorization of the n-by-n
74: *> matrix Z computed by DGETC2: Z = P * L * U * Q
75: *> \endverbatim
76: *>
77: *> \param[in] LDZ
78: *> \verbatim
79: *> LDZ is INTEGER
80: *> The leading dimension of the array Z. LDA >= max(1, N).
81: *> \endverbatim
82: *>
83: *> \param[in,out] RHS
84: *> \verbatim
85: *> RHS is DOUBLE PRECISION array, dimension (N)
86: *> On entry, RHS contains contributions from other subsystems.
87: *> On exit, RHS contains the solution of the subsystem with
88: *> entries according to the value of IJOB (see above).
89: *> \endverbatim
90: *>
91: *> \param[in,out] RDSUM
92: *> \verbatim
93: *> RDSUM is DOUBLE PRECISION
94: *> On entry, the sum of squares of computed contributions to
95: *> the Dif-estimate under computation by DTGSYL, where the
96: *> scaling factor RDSCAL (see below) has been factored out.
97: *> On exit, the corresponding sum of squares updated with the
98: *> contributions from the current sub-system.
99: *> If TRANS = 'T' RDSUM is not touched.
100: *> NOTE: RDSUM only makes sense when DTGSY2 is called by STGSYL.
101: *> \endverbatim
102: *>
103: *> \param[in,out] RDSCAL
104: *> \verbatim
105: *> RDSCAL is DOUBLE PRECISION
106: *> On entry, scaling factor used to prevent overflow in RDSUM.
107: *> On exit, RDSCAL is updated w.r.t. the current contributions
108: *> in RDSUM.
109: *> If TRANS = 'T', RDSCAL is not touched.
110: *> NOTE: RDSCAL only makes sense when DTGSY2 is called by
111: *> DTGSYL.
112: *> \endverbatim
113: *>
114: *> \param[in] IPIV
115: *> \verbatim
116: *> IPIV is INTEGER array, dimension (N).
117: *> The pivot indices; for 1 <= i <= N, row i of the
118: *> matrix has been interchanged with row IPIV(i).
119: *> \endverbatim
120: *>
121: *> \param[in] JPIV
122: *> \verbatim
123: *> JPIV is INTEGER array, dimension (N).
124: *> The pivot indices; for 1 <= j <= N, column j of the
125: *> matrix has been interchanged with column JPIV(j).
126: *> \endverbatim
127: *
128: * Authors:
129: * ========
130: *
131: *> \author Univ. of Tennessee
132: *> \author Univ. of California Berkeley
133: *> \author Univ. of Colorado Denver
134: *> \author NAG Ltd.
135: *
136: *> \ingroup doubleOTHERauxiliary
137: *
138: *> \par Further Details:
139: * =====================
140: *>
141: *> This routine is a further developed implementation of algorithm
142: *> BSOLVE in [1] using complete pivoting in the LU factorization.
143: *
144: *> \par Contributors:
145: * ==================
146: *>
147: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
148: *> Umea University, S-901 87 Umea, Sweden.
149: *
150: *> \par References:
151: * ================
152: *>
153: *> \verbatim
154: *>
155: *>
156: *> [1] Bo Kagstrom and Lars Westin,
157: *> Generalized Schur Methods with Condition Estimators for
158: *> Solving the Generalized Sylvester Equation, IEEE Transactions
159: *> on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.
160: *>
161: *> [2] Peter Poromaa,
162: *> On Efficient and Robust Estimators for the Separation
163: *> between two Regular Matrix Pairs with Applications in
164: *> Condition Estimation. Report IMINF-95.05, Departement of
165: *> Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.
166: *> \endverbatim
167: *>
168: * =====================================================================
169: SUBROUTINE DLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
170: $ JPIV )
171: *
172: * -- LAPACK auxiliary routine --
173: * -- LAPACK is a software package provided by Univ. of Tennessee, --
174: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
175: *
176: * .. Scalar Arguments ..
177: INTEGER IJOB, LDZ, N
178: DOUBLE PRECISION RDSCAL, RDSUM
179: * ..
180: * .. Array Arguments ..
181: INTEGER IPIV( * ), JPIV( * )
182: DOUBLE PRECISION RHS( * ), Z( LDZ, * )
183: * ..
184: *
185: * =====================================================================
186: *
187: * .. Parameters ..
188: INTEGER MAXDIM
189: PARAMETER ( MAXDIM = 8 )
190: DOUBLE PRECISION ZERO, ONE
191: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
192: * ..
193: * .. Local Scalars ..
194: INTEGER I, INFO, J, K
195: DOUBLE PRECISION BM, BP, PMONE, SMINU, SPLUS, TEMP
196: * ..
197: * .. Local Arrays ..
198: INTEGER IWORK( MAXDIM )
199: DOUBLE PRECISION WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM )
200: * ..
201: * .. External Subroutines ..
202: EXTERNAL DAXPY, DCOPY, DGECON, DGESC2, DLASSQ, DLASWP,
203: $ DSCAL
204: * ..
205: * .. External Functions ..
206: DOUBLE PRECISION DASUM, DDOT
207: EXTERNAL DASUM, DDOT
208: * ..
209: * .. Intrinsic Functions ..
210: INTRINSIC ABS, SQRT
211: * ..
212: * .. Executable Statements ..
213: *
214: IF( IJOB.NE.2 ) THEN
215: *
216: * Apply permutations IPIV to RHS
217: *
218: CALL DLASWP( 1, RHS, LDZ, 1, N-1, IPIV, 1 )
219: *
220: * Solve for L-part choosing RHS either to +1 or -1.
221: *
222: PMONE = -ONE
223: *
224: DO 10 J = 1, N - 1
225: BP = RHS( J ) + ONE
226: BM = RHS( J ) - ONE
227: SPLUS = ONE
228: *
229: * Look-ahead for L-part RHS(1:N-1) = + or -1, SPLUS and
230: * SMIN computed more efficiently than in BSOLVE [1].
231: *
232: SPLUS = SPLUS + DDOT( N-J, Z( J+1, J ), 1, Z( J+1, J ), 1 )
233: SMINU = DDOT( N-J, Z( J+1, J ), 1, RHS( J+1 ), 1 )
234: SPLUS = SPLUS*RHS( J )
235: IF( SPLUS.GT.SMINU ) THEN
236: RHS( J ) = BP
237: ELSE IF( SMINU.GT.SPLUS ) THEN
238: RHS( J ) = BM
239: ELSE
240: *
241: * In this case the updating sums are equal and we can
242: * choose RHS(J) +1 or -1. The first time this happens
243: * we choose -1, thereafter +1. This is a simple way to
244: * get good estimates of matrices like Byers well-known
245: * example (see [1]). (Not done in BSOLVE.)
246: *
247: RHS( J ) = RHS( J ) + PMONE
248: PMONE = ONE
249: END IF
250: *
251: * Compute the remaining r.h.s.
252: *
253: TEMP = -RHS( J )
254: CALL DAXPY( N-J, TEMP, Z( J+1, J ), 1, RHS( J+1 ), 1 )
255: *
256: 10 CONTINUE
257: *
258: * Solve for U-part, look-ahead for RHS(N) = +-1. This is not done
259: * in BSOLVE and will hopefully give us a better estimate because
260: * any ill-conditioning of the original matrix is transferred to U
261: * and not to L. U(N, N) is an approximation to sigma_min(LU).
262: *
263: CALL DCOPY( N-1, RHS, 1, XP, 1 )
264: XP( N ) = RHS( N ) + ONE
265: RHS( N ) = RHS( N ) - ONE
266: SPLUS = ZERO
267: SMINU = ZERO
268: DO 30 I = N, 1, -1
269: TEMP = ONE / Z( I, I )
270: XP( I ) = XP( I )*TEMP
271: RHS( I ) = RHS( I )*TEMP
272: DO 20 K = I + 1, N
273: XP( I ) = XP( I ) - XP( K )*( Z( I, K )*TEMP )
274: RHS( I ) = RHS( I ) - RHS( K )*( Z( I, K )*TEMP )
275: 20 CONTINUE
276: SPLUS = SPLUS + ABS( XP( I ) )
277: SMINU = SMINU + ABS( RHS( I ) )
278: 30 CONTINUE
279: IF( SPLUS.GT.SMINU )
280: $ CALL DCOPY( N, XP, 1, RHS, 1 )
281: *
282: * Apply the permutations JPIV to the computed solution (RHS)
283: *
284: CALL DLASWP( 1, RHS, LDZ, 1, N-1, JPIV, -1 )
285: *
286: * Compute the sum of squares
287: *
288: CALL DLASSQ( N, RHS, 1, RDSCAL, RDSUM )
289: *
290: ELSE
291: *
292: * IJOB = 2, Compute approximate nullvector XM of Z
293: *
294: CALL DGECON( 'I', N, Z, LDZ, ONE, TEMP, WORK, IWORK, INFO )
295: CALL DCOPY( N, WORK( N+1 ), 1, XM, 1 )
296: *
297: * Compute RHS
298: *
299: CALL DLASWP( 1, XM, LDZ, 1, N-1, IPIV, -1 )
300: TEMP = ONE / SQRT( DDOT( N, XM, 1, XM, 1 ) )
301: CALL DSCAL( N, TEMP, XM, 1 )
302: CALL DCOPY( N, XM, 1, XP, 1 )
303: CALL DAXPY( N, ONE, RHS, 1, XP, 1 )
304: CALL DAXPY( N, -ONE, XM, 1, RHS, 1 )
305: CALL DGESC2( N, Z, LDZ, RHS, IPIV, JPIV, TEMP )
306: CALL DGESC2( N, Z, LDZ, XP, IPIV, JPIV, TEMP )
307: IF( DASUM( N, XP, 1 ).GT.DASUM( N, RHS, 1 ) )
308: $ CALL DCOPY( N, XP, 1, RHS, 1 )
309: *
310: * Compute the sum of squares
311: *
312: CALL DLASSQ( N, RHS, 1, RDSCAL, RDSUM )
313: *
314: END IF
315: *
316: RETURN
317: *
318: * End of DLATDF
319: *
320: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlatdf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack