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