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