LAPACK  3.6.1 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
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 acoording 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 *> \date June 2016
137 *
138 *> \ingroup doubleOTHERauxiliary
139 *
140 *> \par Further Details:
141 * =====================
142 *>
143 *> This routine is a further developed implementation of algorithm
144 *> BSOLVE in [1] using complete pivoting in the LU factorization.
145 *
146 *> \par Contributors:
147 * ==================
148 *>
149 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
150 *> Umea University, S-901 87 Umea, Sweden.
151 *
152 *> \par References:
153 * ================
154 *>
155 *> \verbatim
156 *>
157 *>
158 *> [1] Bo Kagstrom and Lars Westin,
159 *> Generalized Schur Methods with Condition Estimators for
160 *> Solving the Generalized Sylvester Equation, IEEE Transactions
161 *> on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.
162 *>
163 *> [2] Peter Poromaa,
164 *> On Efficient and Robust Estimators for the Separation
165 *> between two Regular Matrix Pairs with Applications in
166 *> Condition Estimation. Report IMINF-95.05, Departement of
167 *> Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.
168 *> \endverbatim
169 *>
170 * =====================================================================
171  SUBROUTINE dlatdf( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
172  \$ jpiv )
173 *
174 * -- LAPACK auxiliary routine (version 3.6.1) --
175 * -- LAPACK is a software package provided by Univ. of Tennessee, --
176 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
177 * June 2016
178 *
179 * .. Scalar Arguments ..
180  INTEGER IJOB, LDZ, N
181  DOUBLE PRECISION RDSCAL, RDSUM
182 * ..
183 * .. Array Arguments ..
184  INTEGER IPIV( * ), JPIV( * )
185  DOUBLE PRECISION RHS( * ), Z( ldz, * )
186 * ..
187 *
188 * =====================================================================
189 *
190 * .. Parameters ..
191  INTEGER MAXDIM
192  parameter ( maxdim = 8 )
193  DOUBLE PRECISION ZERO, ONE
194  parameter ( zero = 0.0d+0, one = 1.0d+0 )
195 * ..
196 * .. Local Scalars ..
197  INTEGER I, INFO, J, K
198  DOUBLE PRECISION BM, BP, PMONE, SMINU, SPLUS, TEMP
199 * ..
200 * .. Local Arrays ..
201  INTEGER IWORK( maxdim )
202  DOUBLE PRECISION WORK( 4*maxdim ), XM( maxdim ), XP( maxdim )
203 * ..
204 * .. External Subroutines ..
205  EXTERNAL daxpy, dcopy, dgecon, dgesc2, dlassq, dlaswp,
206  \$ dscal
207 * ..
208 * .. External Functions ..
209  DOUBLE PRECISION DASUM, DDOT
210  EXTERNAL dasum, ddot
211 * ..
212 * .. Intrinsic Functions ..
213  INTRINSIC abs, sqrt
214 * ..
215 * .. Executable Statements ..
216 *
217  IF( ijob.NE.2 ) THEN
218 *
219 * Apply permutations IPIV to RHS
220 *
221  CALL dlaswp( 1, rhs, ldz, 1, n-1, ipiv, 1 )
222 *
223 * Solve for L-part choosing RHS either to +1 or -1.
224 *
225  pmone = -one
226 *
227  DO 10 j = 1, n - 1
228  bp = rhs( j ) + one
229  bm = rhs( j ) - one
230  splus = one
231 *
232 * Look-ahead for L-part RHS(1:N-1) = + or -1, SPLUS and
233 * SMIN computed more efficiently than in BSOLVE [1].
234 *
235  splus = splus + ddot( n-j, z( j+1, j ), 1, z( j+1, j ), 1 )
236  sminu = ddot( n-j, z( j+1, j ), 1, rhs( j+1 ), 1 )
237  splus = splus*rhs( j )
238  IF( splus.GT.sminu ) THEN
239  rhs( j ) = bp
240  ELSE IF( sminu.GT.splus ) THEN
241  rhs( j ) = bm
242  ELSE
243 *
244 * In this case the updating sums are equal and we can
245 * choose RHS(J) +1 or -1. The first time this happens
246 * we choose -1, thereafter +1. This is a simple way to
247 * get good estimates of matrices like Byers well-known
248 * example (see [1]). (Not done in BSOLVE.)
249 *
250  rhs( j ) = rhs( j ) + pmone
251  pmone = one
252  END IF
253 *
254 * Compute the remaining r.h.s.
255 *
256  temp = -rhs( j )
257  CALL daxpy( n-j, temp, z( j+1, j ), 1, rhs( j+1 ), 1 )
258 *
259  10 CONTINUE
260 *
261 * Solve for U-part, look-ahead for RHS(N) = +-1. This is not done
262 * in BSOLVE and will hopefully give us a better estimate because
263 * any ill-conditioning of the original matrix is transfered to U
264 * and not to L. U(N, N) is an approximation to sigma_min(LU).
265 *
266  CALL dcopy( n-1, rhs, 1, xp, 1 )
267  xp( n ) = rhs( n ) + one
268  rhs( n ) = rhs( n ) - one
269  splus = zero
270  sminu = zero
271  DO 30 i = n, 1, -1
272  temp = one / z( i, i )
273  xp( i ) = xp( i )*temp
274  rhs( i ) = rhs( i )*temp
275  DO 20 k = i + 1, n
276  xp( i ) = xp( i ) - xp( k )*( z( i, k )*temp )
277  rhs( i ) = rhs( i ) - rhs( k )*( z( i, k )*temp )
278  20 CONTINUE
279  splus = splus + abs( xp( i ) )
280  sminu = sminu + abs( rhs( i ) )
281  30 CONTINUE
282  IF( splus.GT.sminu )
283  \$ CALL dcopy( n, xp, 1, rhs, 1 )
284 *
285 * Apply the permutations JPIV to the computed solution (RHS)
286 *
287  CALL dlaswp( 1, rhs, ldz, 1, n-1, jpiv, -1 )
288 *
289 * Compute the sum of squares
290 *
291  CALL dlassq( n, rhs, 1, rdscal, rdsum )
292 *
293  ELSE
294 *
295 * IJOB = 2, Compute approximate nullvector XM of Z
296 *
297  CALL dgecon( 'I', n, z, ldz, one, temp, work, iwork, info )
298  CALL dcopy( n, work( n+1 ), 1, xm, 1 )
299 *
300 * Compute RHS
301 *
302  CALL dlaswp( 1, xm, ldz, 1, n-1, ipiv, -1 )
303  temp = one / sqrt( ddot( n, xm, 1, xm, 1 ) )
304  CALL dscal( n, temp, xm, 1 )
305  CALL dcopy( n, xm, 1, xp, 1 )
306  CALL daxpy( n, one, rhs, 1, xp, 1 )
307  CALL daxpy( n, -one, xm, 1, rhs, 1 )
308  CALL dgesc2( n, z, ldz, rhs, ipiv, jpiv, temp )
309  CALL dgesc2( n, z, ldz, xp, ipiv, jpiv, temp )
310  IF( dasum( n, xp, 1 ).GT.dasum( n, rhs, 1 ) )
311  \$ CALL dcopy( n, xp, 1, rhs, 1 )
312 *
313 * Compute the sum of squares
314 *
315  CALL dlassq( n, rhs, 1, rdscal, rdsum )
316 *
317  END IF
318 *
319  RETURN
320 *
321 * End of DLATDF
322 *
323  END
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
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:116
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:173
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
DLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: dlaswp.f:116
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105
subroutine dgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
DGECON
Definition: dgecon.f:126