Annotation of rpl/lapack/lapack/zlatdf.f, revision 1.3
1.1 bertrand 1: SUBROUTINE ZLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
2: $ JPIV )
3: *
4: * -- LAPACK auxiliary routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER IJOB, LDZ, N
11: DOUBLE PRECISION RDSCAL, RDSUM
12: * ..
13: * .. Array Arguments ..
14: INTEGER IPIV( * ), JPIV( * )
15: COMPLEX*16 RHS( * ), Z( LDZ, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * ZLATDF computes the contribution to the reciprocal Dif-estimate
22: * by solving for x in Z * x = b, where b is chosen such that the norm
23: * of x is as large as possible. It is assumed that LU decomposition
24: * of Z has been computed by ZGETC2. On entry RHS = f holds the
25: * contribution from earlier solved sub-systems, and on return RHS = x.
26: *
27: * The factorization of Z returned by ZGETC2 has the form
28: * Z = P * L * U * Q, where P and Q are permutation matrices. L is lower
29: * triangular with unit diagonal elements and U is upper triangular.
30: *
31: * Arguments
32: * =========
33: *
34: * IJOB (input) INTEGER
35: * IJOB = 2: First compute an approximative null-vector e
36: * of Z using ZGECON, e is normalized and solve for
37: * Zx = +-e - f with the sign giving the greater value of
38: * 2-norm(x). About 5 times as expensive as Default.
39: * IJOB .ne. 2: Local look ahead strategy where
40: * all entries of the r.h.s. b is choosen as either +1 or
41: * -1. Default.
42: *
43: * N (input) INTEGER
44: * The number of columns of the matrix Z.
45: *
46: * Z (input) DOUBLE PRECISION array, dimension (LDZ, N)
47: * On entry, the LU part of the factorization of the n-by-n
48: * matrix Z computed by ZGETC2: Z = P * L * U * Q
49: *
50: * LDZ (input) INTEGER
51: * The leading dimension of the array Z. LDA >= max(1, N).
52: *
53: * RHS (input/output) DOUBLE PRECISION array, dimension (N).
54: * On entry, RHS contains contributions from other subsystems.
55: * On exit, RHS contains the solution of the subsystem with
56: * entries according to the value of IJOB (see above).
57: *
58: * RDSUM (input/output) DOUBLE PRECISION
59: * On entry, the sum of squares of computed contributions to
60: * the Dif-estimate under computation by ZTGSYL, where the
61: * scaling factor RDSCAL (see below) has been factored out.
62: * On exit, the corresponding sum of squares updated with the
63: * contributions from the current sub-system.
64: * If TRANS = 'T' RDSUM is not touched.
65: * NOTE: RDSUM only makes sense when ZTGSY2 is called by CTGSYL.
66: *
67: * RDSCAL (input/output) DOUBLE PRECISION
68: * On entry, scaling factor used to prevent overflow in RDSUM.
69: * On exit, RDSCAL is updated w.r.t. the current contributions
70: * in RDSUM.
71: * If TRANS = 'T', RDSCAL is not touched.
72: * NOTE: RDSCAL only makes sense when ZTGSY2 is called by
73: * ZTGSYL.
74: *
75: * IPIV (input) INTEGER array, dimension (N).
76: * The pivot indices; for 1 <= i <= N, row i of the
77: * matrix has been interchanged with row IPIV(i).
78: *
79: * JPIV (input) INTEGER array, dimension (N).
80: * The pivot indices; for 1 <= j <= N, column j of the
81: * matrix has been interchanged with column JPIV(j).
82: *
83: * Further Details
84: * ===============
85: *
86: * Based on contributions by
87: * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
88: * Umea University, S-901 87 Umea, Sweden.
89: *
90: * This routine is a further developed implementation of algorithm
91: * BSOLVE in [1] using complete pivoting in the LU factorization.
92: *
93: * [1] Bo Kagstrom and Lars Westin,
94: * Generalized Schur Methods with Condition Estimators for
95: * Solving the Generalized Sylvester Equation, IEEE Transactions
96: * on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.
97: *
98: * [2] Peter Poromaa,
99: * On Efficient and Robust Estimators for the Separation
100: * between two Regular Matrix Pairs with Applications in
101: * Condition Estimation. Report UMINF-95.05, Department of
102: * Computing Science, Umea University, S-901 87 Umea, Sweden,
103: * 1995.
104: *
105: * =====================================================================
106: *
107: * .. Parameters ..
108: INTEGER MAXDIM
109: PARAMETER ( MAXDIM = 2 )
110: DOUBLE PRECISION ZERO, ONE
111: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
112: COMPLEX*16 CONE
113: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
114: * ..
115: * .. Local Scalars ..
116: INTEGER I, INFO, J, K
117: DOUBLE PRECISION RTEMP, SCALE, SMINU, SPLUS
118: COMPLEX*16 BM, BP, PMONE, TEMP
119: * ..
120: * .. Local Arrays ..
121: DOUBLE PRECISION RWORK( MAXDIM )
122: COMPLEX*16 WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM )
123: * ..
124: * .. External Subroutines ..
125: EXTERNAL ZAXPY, ZCOPY, ZGECON, ZGESC2, ZLASSQ, ZLASWP,
126: $ ZSCAL
127: * ..
128: * .. External Functions ..
129: DOUBLE PRECISION DZASUM
130: COMPLEX*16 ZDOTC
131: EXTERNAL DZASUM, ZDOTC
132: * ..
133: * .. Intrinsic Functions ..
134: INTRINSIC ABS, DBLE, SQRT
135: * ..
136: * .. Executable Statements ..
137: *
138: IF( IJOB.NE.2 ) THEN
139: *
140: * Apply permutations IPIV to RHS
141: *
142: CALL ZLASWP( 1, RHS, LDZ, 1, N-1, IPIV, 1 )
143: *
144: * Solve for L-part choosing RHS either to +1 or -1.
145: *
146: PMONE = -CONE
147: DO 10 J = 1, N - 1
148: BP = RHS( J ) + CONE
149: BM = RHS( J ) - CONE
150: SPLUS = ONE
151: *
152: * Lockahead for L- part RHS(1:N-1) = +-1
153: * SPLUS and SMIN computed more efficiently than in BSOLVE[1].
154: *
155: SPLUS = SPLUS + DBLE( ZDOTC( N-J, Z( J+1, J ), 1, Z( J+1,
156: $ J ), 1 ) )
157: SMINU = DBLE( ZDOTC( N-J, Z( J+1, J ), 1, RHS( J+1 ), 1 ) )
158: SPLUS = SPLUS*DBLE( RHS( J ) )
159: IF( SPLUS.GT.SMINU ) THEN
160: RHS( J ) = BP
161: ELSE IF( SMINU.GT.SPLUS ) THEN
162: RHS( J ) = BM
163: ELSE
164: *
165: * In this case the updating sums are equal and we can
166: * choose RHS(J) +1 or -1. The first time this happens we
167: * choose -1, thereafter +1. This is a simple way to get
168: * good estimates of matrices like Byers well-known example
169: * (see [1]). (Not done in BSOLVE.)
170: *
171: RHS( J ) = RHS( J ) + PMONE
172: PMONE = CONE
173: END IF
174: *
175: * Compute the remaining r.h.s.
176: *
177: TEMP = -RHS( J )
178: CALL ZAXPY( N-J, TEMP, Z( J+1, J ), 1, RHS( J+1 ), 1 )
179: 10 CONTINUE
180: *
181: * Solve for U- part, lockahead for RHS(N) = +-1. This is not done
182: * In BSOLVE and will hopefully give us a better estimate because
183: * any ill-conditioning of the original matrix is transfered to U
184: * and not to L. U(N, N) is an approximation to sigma_min(LU).
185: *
186: CALL ZCOPY( N-1, RHS, 1, WORK, 1 )
187: WORK( N ) = RHS( N ) + CONE
188: RHS( N ) = RHS( N ) - CONE
189: SPLUS = ZERO
190: SMINU = ZERO
191: DO 30 I = N, 1, -1
192: TEMP = CONE / Z( I, I )
193: WORK( I ) = WORK( I )*TEMP
194: RHS( I ) = RHS( I )*TEMP
195: DO 20 K = I + 1, N
196: WORK( I ) = WORK( I ) - WORK( K )*( Z( I, K )*TEMP )
197: RHS( I ) = RHS( I ) - RHS( K )*( Z( I, K )*TEMP )
198: 20 CONTINUE
199: SPLUS = SPLUS + ABS( WORK( I ) )
200: SMINU = SMINU + ABS( RHS( I ) )
201: 30 CONTINUE
202: IF( SPLUS.GT.SMINU )
203: $ CALL ZCOPY( N, WORK, 1, RHS, 1 )
204: *
205: * Apply the permutations JPIV to the computed solution (RHS)
206: *
207: CALL ZLASWP( 1, RHS, LDZ, 1, N-1, JPIV, -1 )
208: *
209: * Compute the sum of squares
210: *
211: CALL ZLASSQ( N, RHS, 1, RDSCAL, RDSUM )
212: RETURN
213: END IF
214: *
215: * ENTRY IJOB = 2
216: *
217: * Compute approximate nullvector XM of Z
218: *
219: CALL ZGECON( 'I', N, Z, LDZ, ONE, RTEMP, WORK, RWORK, INFO )
220: CALL ZCOPY( N, WORK( N+1 ), 1, XM, 1 )
221: *
222: * Compute RHS
223: *
224: CALL ZLASWP( 1, XM, LDZ, 1, N-1, IPIV, -1 )
225: TEMP = CONE / SQRT( ZDOTC( N, XM, 1, XM, 1 ) )
226: CALL ZSCAL( N, TEMP, XM, 1 )
227: CALL ZCOPY( N, XM, 1, XP, 1 )
228: CALL ZAXPY( N, CONE, RHS, 1, XP, 1 )
229: CALL ZAXPY( N, -CONE, XM, 1, RHS, 1 )
230: CALL ZGESC2( N, Z, LDZ, RHS, IPIV, JPIV, SCALE )
231: CALL ZGESC2( N, Z, LDZ, XP, IPIV, JPIV, SCALE )
232: IF( DZASUM( N, XP, 1 ).GT.DZASUM( N, RHS, 1 ) )
233: $ CALL ZCOPY( N, XP, 1, RHS, 1 )
234: *
235: * Compute the sum of squares
236: *
237: CALL ZLASSQ( N, RHS, 1, RDSCAL, RDSUM )
238: RETURN
239: *
240: * End of ZLATDF
241: *
242: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlatdf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack