1
SUBROUTINE ZLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
4
* -- LAPACK auxiliary routine (version 3.0) --
5
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6
* Courant Institute, Argonne National Lab, and Rice University
9
* .. Scalar Arguments ..
11
DOUBLE PRECISION RDSCAL, RDSUM
13
* .. Array Arguments ..
14
INTEGER IPIV( * ), JPIV( * )
15
COMPLEX*16 RHS( * ), Z( LDZ, * )
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.
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.
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
44
* The number of columns of the matrix Z.
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
51
* The leading dimension of the array Z. LDA >= max(1, N).
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).
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.
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
71
* If TRANS = 'T', RDSCAL is not touched.
72
* NOTE: RDSCAL only makes sense when ZTGSY2 is called by
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).
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).
86
* Based on contributions by
87
* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
88
* Umea University, S-901 87 Umea, Sweden.
90
* This routine is a further developed implementation of algorithm
91
* BSOLVE in [1] using complete pivoting in the LU factorization.
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.
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,
105
* =====================================================================
109
PARAMETER ( MAXDIM = 2 )
110
DOUBLE PRECISION ZERO, ONE
111
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
113
PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
115
* .. Local Scalars ..
116
INTEGER I, INFO, J, K
117
DOUBLE PRECISION RTEMP, SCALE, SMINU, SPLUS
118
COMPLEX*16 BM, BP, PMONE, TEMP
121
DOUBLE PRECISION RWORK( MAXDIM )
122
COMPLEX*16 WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM )
124
* .. External Subroutines ..
125
EXTERNAL ZAXPY, ZCOPY, ZGECON, ZGESC2, ZLASSQ, ZLASWP,
128
* .. External Functions ..
129
DOUBLE PRECISION DZASUM
131
EXTERNAL DZASUM, ZDOTC
133
* .. Intrinsic Functions ..
134
INTRINSIC ABS, DBLE, SQRT
136
* .. Executable Statements ..
140
* Apply permutations IPIV to RHS
142
CALL ZLASWP( 1, RHS, LDZ, 1, N-1, IPIV, 1 )
144
* Solve for L-part choosing RHS either to +1 or -1.
152
* Lockahead for L- part RHS(1:N-1) = +-1
153
* SPLUS and SMIN computed more efficiently than in BSOLVE[1].
155
SPLUS = SPLUS + DBLE( ZDOTC( N-J, Z( J+1, J ), 1, Z( 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
161
ELSE IF( SMINU.GT.SPLUS ) THEN
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.)
171
RHS( J ) = RHS( J ) + PMONE
175
* Compute the remaining r.h.s.
178
CALL ZAXPY( N-J, TEMP, Z( J+1, J ), 1, RHS( J+1 ), 1 )
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).
186
CALL ZCOPY( N-1, RHS, 1, WORK, 1 )
187
WORK( N ) = RHS( N ) + CONE
188
RHS( N ) = RHS( N ) - CONE
192
TEMP = CONE / Z( I, I )
193
WORK( I ) = WORK( I )*TEMP
194
RHS( I ) = RHS( I )*TEMP
196
WORK( I ) = WORK( I ) - WORK( K )*( Z( I, K )*TEMP )
197
RHS( I ) = RHS( I ) - RHS( K )*( Z( I, K )*TEMP )
199
SPLUS = SPLUS + ABS( WORK( I ) )
200
SMINU = SMINU + ABS( RHS( I ) )
203
$ CALL ZCOPY( N, WORK, 1, RHS, 1 )
205
* Apply the permutations JPIV to the computed solution (RHS)
207
CALL ZLASWP( 1, RHS, LDZ, 1, N-1, JPIV, -1 )
209
* Compute the sum of squares
211
CALL ZLASSQ( N, RHS, 1, RDSCAL, RDSUM )
217
* Compute approximate nullvector XM of Z
219
CALL ZGECON( 'I', N, Z, LDZ, ONE, RTEMP, WORK, RWORK, INFO )
220
CALL ZCOPY( N, WORK( N+1 ), 1, XM, 1 )
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 )
235
* Compute the sum of squares
237
CALL ZLASSQ( N, RHS, 1, RDSCAL, RDSUM )