LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
slatdf.f
Go to the documentation of this file.
1 *> \brief \b SLATDF 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 SLATDF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slatdf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slatdf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slatdf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
22 * JPIV )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER IJOB, LDZ, N
26 * REAL RDSCAL, RDSUM
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IPIV( * ), JPIV( * )
30 * REAL RHS( * ), Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SLATDF uses the LU factorization of the n-by-n matrix Z computed by
40 *> SGETC2 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 SGETC2 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 SGECON, 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 REAL array, dimension (LDZ, N)
73 *> On entry, the LU part of the factorization of the n-by-n
74 *> matrix Z computed by SGETC2: 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 REAL 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 REAL
94 *> On entry, the sum of squares of computed contributions to
95 *> the Dif-estimate under computation by STGSYL, 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 STGSY2 is called by STGSYL.
101 *> \endverbatim
102 *>
103 *> \param[in,out] RDSCAL
104 *> \verbatim
105 *> RDSCAL is REAL
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 STGSY2 is called by
111 *> STGSYL.
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 realOTHERauxiliary
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 slatdf( 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  REAL RDSCAL, RDSUM
179 * ..
180 * .. Array Arguments ..
181  INTEGER IPIV( * ), JPIV( * )
182  REAL RHS( * ), Z( LDZ, * )
183 * ..
184 *
185 * =====================================================================
186 *
187 * .. Parameters ..
188  INTEGER MAXDIM
189  parameter( maxdim = 8 )
190  REAL ZERO, ONE
191  parameter( zero = 0.0e+0, one = 1.0e+0 )
192 * ..
193 * .. Local Scalars ..
194  INTEGER I, INFO, J, K
195  REAL BM, BP, PMONE, SMINU, SPLUS, TEMP
196 * ..
197 * .. Local Arrays ..
198  INTEGER IWORK( MAXDIM )
199  REAL WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM )
200 * ..
201 * .. External Subroutines ..
202  EXTERNAL saxpy, scopy, sgecon, sgesc2, slassq, slaswp,
203  $ sscal
204 * ..
205 * .. External Functions ..
206  REAL SASUM, SDOT
207  EXTERNAL sasum, sdot
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 slaswp( 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 + sdot( n-j, z( j+1, j ), 1, z( j+1, j ), 1 )
233  sminu = sdot( 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 saxpy( 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 scopy( 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 scopy( n, xp, 1, rhs, 1 )
281 *
282 * Apply the permutations JPIV to the computed solution (RHS)
283 *
284  CALL slaswp( 1, rhs, ldz, 1, n-1, jpiv, -1 )
285 *
286 * Compute the sum of squares
287 *
288  CALL slassq( n, rhs, 1, rdscal, rdsum )
289 *
290  ELSE
291 *
292 * IJOB = 2, Compute approximate nullvector XM of Z
293 *
294  CALL sgecon( 'I', n, z, ldz, one, temp, work, iwork, info )
295  CALL scopy( n, work( n+1 ), 1, xm, 1 )
296 *
297 * Compute RHS
298 *
299  CALL slaswp( 1, xm, ldz, 1, n-1, ipiv, -1 )
300  temp = one / sqrt( sdot( n, xm, 1, xm, 1 ) )
301  CALL sscal( n, temp, xm, 1 )
302  CALL scopy( n, xm, 1, xp, 1 )
303  CALL saxpy( n, one, rhs, 1, xp, 1 )
304  CALL saxpy( n, -one, xm, 1, rhs, 1 )
305  CALL sgesc2( n, z, ldz, rhs, ipiv, jpiv, temp )
306  CALL sgesc2( n, z, ldz, xp, ipiv, jpiv, temp )
307  IF( sasum( n, xp, 1 ).GT.sasum( n, rhs, 1 ) )
308  $ CALL scopy( n, xp, 1, rhs, 1 )
309 *
310 * Compute the sum of squares
311 *
312  CALL slassq( n, rhs, 1, rdscal, rdsum )
313 *
314  END IF
315 *
316  RETURN
317 *
318 * End of SLATDF
319 *
320  END
subroutine slassq(n, x, incx, scl, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition: slassq.f90:126
subroutine sgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
SGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition: sgesc2.f:114
subroutine sgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
SGECON
Definition: sgecon.f:124
subroutine slaswp(N, A, LDA, K1, K2, IPIV, INCX)
SLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: slaswp.f:115
subroutine slatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
SLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition: slatdf.f:171
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89