LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlatdf.f
Go to the documentation of this file.
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
subroutine dlassq(n, x, incx, scl, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f90:126
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition: dgesc2.f:114
subroutine dgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
DGECON
Definition: dgecon.f:124
subroutine dlatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition: dlatdf.f:171
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
DLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: dlaswp.f:115