1 /* ../netlib/slatdf.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static integer c__1 = 1;
5 static integer c_n1 = -1;
6 static real c_b23 = 1.f;
7 static real c_b37 = -1.f;
8 /* > \brief \b SLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contrib ution to the reciprocal Dif-estimate. */
9 /* =========== DOCUMENTATION =========== */
10 /* Online html documentation available at */
11 /* http://www.netlib.org/lapack/explore-html/ */
12 /* > \htmlonly */
13 /* > Download SLATDF + dependencies */
14 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slatdf. f"> */
15 /* > [TGZ]</a> */
16 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slatdf. f"> */
17 /* > [ZIP]</a> */
18 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slatdf. f"> */
19 /* > [TXT]</a> */
20 /* > \endhtmlonly */
21 /* Definition: */
22 /* =========== */
23 /* SUBROUTINE SLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, */
24 /* JPIV ) */
25 /* .. Scalar Arguments .. */
26 /* INTEGER IJOB, LDZ, N */
27 /* REAL RDSCAL, RDSUM */
28 /* .. */
29 /* .. Array Arguments .. */
30 /* INTEGER IPIV( * ), JPIV( * ) */
31 /* REAL RHS( * ), Z( LDZ, * ) */
32 /* .. */
33 /* > \par Purpose: */
34 /* ============= */
35 /* > */
36 /* > \verbatim */
37 /* > */
38 /* > SLATDF uses the LU factorization of the n-by-n matrix Z computed by */
39 /* > SGETC2 and computes a contribution to the reciprocal Dif-estimate */
40 /* > by solving Z * x = b for x, and choosing the r.h.s. b such that */
41 /* > the norm of x is as large as possible. On entry RHS = b holds the */
42 /* > contribution from earlier solved sub-systems, and on return RHS = x. */
43 /* > */
44 /* > The factorization of Z returned by SGETC2 has the form Z = P*L*U*Q, */
45 /* > where P and Q are permutation matrices. L is lower triangular with */
46 /* > unit diagonal elements and U is upper triangular. */
47 /* > \endverbatim */
48 /* Arguments: */
49 /* ========== */
50 /* > \param[in] IJOB */
51 /* > \verbatim */
52 /* > IJOB is INTEGER */
53 /* > IJOB = 2: First compute an approximative null-vector e */
54 /* > of Z using SGECON, e is normalized and solve for */
55 /* > Zx = +-e - f with the sign giving the greater value */
56 /* > of 2-norm(x). About 5 times as expensive as Default. */
57 /* > IJOB .ne. 2: Local look ahead strategy where all entries of */
58 /* > the r.h.s. b is choosen as either +1 or -1 (Default). */
59 /* > \endverbatim */
60 /* > */
61 /* > \param[in] N */
62 /* > \verbatim */
63 /* > N is INTEGER */
64 /* > The number of columns of the matrix Z. */
65 /* > \endverbatim */
66 /* > */
67 /* > \param[in] Z */
68 /* > \verbatim */
69 /* > Z is REAL array, dimension (LDZ, N) */
70 /* > On entry, the LU part of the factorization of the n-by-n */
71 /* > matrix Z computed by SGETC2: Z = P * L * U * Q */
72 /* > \endverbatim */
73 /* > */
74 /* > \param[in] LDZ */
75 /* > \verbatim */
76 /* > LDZ is INTEGER */
77 /* > The leading dimension of the array Z. LDA >= max(1, N). */
78 /* > \endverbatim */
79 /* > */
80 /* > \param[in,out] RHS */
81 /* > \verbatim */
82 /* > RHS is REAL array, dimension N. */
83 /* > On entry, RHS contains contributions from other subsystems. */
84 /* > On exit, RHS contains the solution of the subsystem with */
85 /* > entries acoording to the value of IJOB (see above). */
86 /* > \endverbatim */
87 /* > */
88 /* > \param[in,out] RDSUM */
89 /* > \verbatim */
90 /* > RDSUM is REAL */
91 /* > On entry, the sum of squares of computed contributions to */
92 /* > the Dif-estimate under computation by STGSYL, where the */
93 /* > scaling factor RDSCAL (see below) has been factored out. */
94 /* > On exit, the corresponding sum of squares updated with the */
95 /* > contributions from the current sub-system. */
96 /* > If TRANS = 'T' RDSUM is not touched. */
97 /* > NOTE: RDSUM only makes sense when STGSY2 is called by STGSYL. */
98 /* > \endverbatim */
99 /* > */
100 /* > \param[in,out] RDSCAL */
101 /* > \verbatim */
102 /* > RDSCAL is REAL */
103 /* > On entry, scaling factor used to prevent overflow in RDSUM. */
104 /* > On exit, RDSCAL is updated w.r.t. the current contributions */
105 /* > in RDSUM. */
106 /* > If TRANS = 'T', RDSCAL is not touched. */
107 /* > NOTE: RDSCAL only makes sense when STGSY2 is called by */
108 /* > STGSYL. */
109 /* > \endverbatim */
110 /* > */
111 /* > \param[in] IPIV */
112 /* > \verbatim */
113 /* > IPIV is INTEGER array, dimension (N). */
114 /* > The pivot indices;
115 for 1 <= i <= N, row i of the */
116 /* > matrix has been interchanged with row IPIV(i). */
117 /* > \endverbatim */
118 /* > */
119 /* > \param[in] JPIV */
120 /* > \verbatim */
121 /* > JPIV is INTEGER array, dimension (N). */
122 /* > The pivot indices;
123 for 1 <= j <= N, column j of the */
124 /* > matrix has been interchanged with column JPIV(j). */
125 /* > \endverbatim */
126 /* Authors: */
127 /* ======== */
128 /* > \author Univ. of Tennessee */
129 /* > \author Univ. of California Berkeley */
130 /* > \author Univ. of Colorado Denver */
131 /* > \author NAG Ltd. */
132 /* > \date September 2012 */
133 /* > \ingroup realOTHERauxiliary */
134 /* > \par Further Details: */
135 /* ===================== */
136 /* > */
137 /* > This routine is a further developed implementation of algorithm */
138 /* > BSOLVE in [1] using complete pivoting in the LU factorization. */
139 /* > \par Contributors: */
140 /* ================== */
141 /* > */
142 /* > Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
143 /* > Umea University, S-901 87 Umea, Sweden. */
144 /* > \par References: */
145 /* ================ */
146 /* > */
147 /* > \verbatim */
148 /* > */
149 /* > */
150 /* > [1] Bo Kagstrom and Lars Westin, */
151 /* > Generalized Schur Methods with Condition Estimators for */
152 /* > Solving the Generalized Sylvester Equation, IEEE Transactions */
153 /* > on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751. */
154 /* > */
155 /* > [2] Peter Poromaa, */
156 /* > On Efficient and Robust Estimators for the Separation */
157 /* > between two Regular Matrix Pairs with Applications in */
158 /* > Condition Estimation. Report IMINF-95.05, Departement of */
159 /* > Computing Science, Umea University, S-901 87 Umea, Sweden, 1995. */
160 /* > \endverbatim */
161 /* > */
162 /* ===================================================================== */
163 /* Subroutine */
slatdf_(integer * ijob,integer * n,real * z__,integer * ldz,real * rhs,real * rdsum,real * rdscal,integer * ipiv,integer * jpiv)164 int slatdf_(integer *ijob, integer *n, real *z__, integer * ldz, real *rhs, real *rdsum, real *rdscal, integer *ipiv, integer * jpiv)
165 {
166 /* System generated locals */
167 integer z_dim1, z_offset, i__1, i__2;
168 real r__1;
169 /* Builtin functions */
170 double sqrt(doublereal);
171 /* Local variables */
172 integer i__, j, k;
173 real bm, bp, xm[8], xp[8];
174 integer info;
175 real temp;
176 extern real sdot_(integer *, real *, integer *, real *, integer *);
177 real work[32];
178 extern /* Subroutine */
179 int sscal_(integer *, real *, real *, integer *);
180 real pmone;
181 extern real sasum_(integer *, real *, integer *);
182 real sminu;
183 integer iwork[8];
184 extern /* Subroutine */
185 int scopy_(integer *, real *, integer *, real *, integer *), saxpy_(integer *, real *, real *, integer *, real *, integer *);
186 real splus;
187 extern /* Subroutine */
188 int sgesc2_(integer *, real *, integer *, real *, integer *, integer *, real *), sgecon_(char *, integer *, real *, integer *, real *, real *, real *, integer *, integer *), slassq_(integer *, real *, integer *, real *, real *), slaswp_( integer *, real *, integer *, integer *, integer *, integer *, integer *);
189 /* -- LAPACK auxiliary routine (version 3.4.2) -- */
190 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
191 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
192 /* September 2012 */
193 /* .. Scalar Arguments .. */
194 /* .. */
195 /* .. Array Arguments .. */
196 /* .. */
197 /* ===================================================================== */
198 /* .. Parameters .. */
199 /* .. */
200 /* .. Local Scalars .. */
201 /* .. */
202 /* .. Local Arrays .. */
203 /* .. */
204 /* .. External Subroutines .. */
205 /* .. */
206 /* .. External Functions .. */
207 /* .. */
208 /* .. Intrinsic Functions .. */
209 /* .. */
210 /* .. Executable Statements .. */
211 /* Parameter adjustments */
212 z_dim1 = *ldz;
213 z_offset = 1 + z_dim1;
214 z__ -= z_offset;
215 --rhs;
216 --ipiv;
217 --jpiv;
218 /* Function Body */
219 if (*ijob != 2)
220 {
221 /* Apply permutations IPIV to RHS */
222 i__1 = *n - 1;
223 slaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &ipiv[1], &c__1);
224 /* Solve for L-part choosing RHS either to +1 or -1. */
225 pmone = -1.f;
226 i__1 = *n - 1;
227 for (j = 1;
228 j <= i__1;
229 ++j)
230 {
231 bp = rhs[j] + 1.f;
232 bm = rhs[j] - 1.f;
233 splus = 1.f;
234 /* Look-ahead for L-part RHS(1:N-1) = + or -1, SPLUS and */
235 /* SMIN computed more efficiently than in BSOLVE [1]. */
236 i__2 = *n - j;
237 splus += sdot_(&i__2, &z__[j + 1 + j * z_dim1], &c__1, &z__[j + 1 + j * z_dim1], &c__1);
238 i__2 = *n - j;
239 sminu = sdot_(&i__2, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1], &c__1);
240 splus *= rhs[j];
241 if (splus > sminu)
242 {
243 rhs[j] = bp;
244 }
245 else if (sminu > splus)
246 {
247 rhs[j] = bm;
248 }
249 else
250 {
251 /* In this case the updating sums are equal and we can */
252 /* choose RHS(J) +1 or -1. The first time this happens */
253 /* we choose -1, thereafter +1. This is a simple way to */
254 /* get good estimates of matrices like Byers well-known */
255 /* example (see [1]). (Not done in BSOLVE.) */
256 rhs[j] += pmone;
257 pmone = 1.f;
258 }
259 /* Compute the remaining r.h.s. */
260 temp = -rhs[j];
261 i__2 = *n - j;
262 saxpy_(&i__2, &temp, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1], &c__1);
263 /* L10: */
264 }
265 /* Solve for U-part, look-ahead for RHS(N) = +-1. This is not done */
266 /* in BSOLVE and will hopefully give us a better estimate because */
267 /* any ill-conditioning of the original matrix is transfered to U */
268 /* and not to L. U(N, N) is an approximation to sigma_min(LU). */
269 i__1 = *n - 1;
270 scopy_(&i__1, &rhs[1], &c__1, xp, &c__1);
271 xp[*n - 1] = rhs[*n] + 1.f;
272 rhs[*n] += -1.f;
273 splus = 0.f;
274 sminu = 0.f;
275 for (i__ = *n;
276 i__ >= 1;
277 --i__)
278 {
279 temp = 1.f / z__[i__ + i__ * z_dim1];
280 xp[i__ - 1] *= temp;
281 rhs[i__] *= temp;
282 i__1 = *n;
283 for (k = i__ + 1;
284 k <= i__1;
285 ++k)
286 {
287 xp[i__ - 1] -= xp[k - 1] * (z__[i__ + k * z_dim1] * temp);
288 rhs[i__] -= rhs[k] * (z__[i__ + k * z_dim1] * temp);
289 /* L20: */
290 }
291 splus += (r__1 = xp[i__ - 1], f2c_abs(r__1));
292 sminu += (r__1 = rhs[i__], f2c_abs(r__1));
293 /* L30: */
294 }
295 if (splus > sminu)
296 {
297 scopy_(n, xp, &c__1, &rhs[1], &c__1);
298 }
299 /* Apply the permutations JPIV to the computed solution (RHS) */
300 i__1 = *n - 1;
301 slaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &jpiv[1], &c_n1);
302 /* Compute the sum of squares */
303 slassq_(n, &rhs[1], &c__1, rdscal, rdsum);
304 }
305 else
306 {
307 /* IJOB = 2, Compute approximate nullvector XM of Z */
308 sgecon_("I", n, &z__[z_offset], ldz, &c_b23, &temp, work, iwork, & info);
309 scopy_(n, &work[*n], &c__1, xm, &c__1);
310 /* Compute RHS */
311 i__1 = *n - 1;
312 slaswp_(&c__1, xm, ldz, &c__1, &i__1, &ipiv[1], &c_n1);
313 temp = 1.f / sqrt(sdot_(n, xm, &c__1, xm, &c__1));
314 sscal_(n, &temp, xm, &c__1);
315 scopy_(n, xm, &c__1, xp, &c__1);
316 saxpy_(n, &c_b23, &rhs[1], &c__1, xp, &c__1);
317 saxpy_(n, &c_b37, xm, &c__1, &rhs[1], &c__1);
318 sgesc2_(n, &z__[z_offset], ldz, &rhs[1], &ipiv[1], &jpiv[1], &temp);
319 sgesc2_(n, &z__[z_offset], ldz, xp, &ipiv[1], &jpiv[1], &temp);
320 if (sasum_(n, xp, &c__1) > sasum_(n, &rhs[1], &c__1))
321 {
322 scopy_(n, xp, &c__1, &rhs[1], &c__1);
323 }
324 /* Compute the sum of squares */
325 slassq_(n, &rhs[1], &c__1, rdscal, rdsum);
326 }
327 return 0;
328 /* End of SLATDF */
329 }
330 /* slatdf_ */
331