00001 /* dlatdf.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__1 = 1; 00019 static integer c_n1 = -1; 00020 static doublereal c_b23 = 1.; 00021 static doublereal c_b37 = -1.; 00022 00023 /* Subroutine */ int dlatdf_(integer *ijob, integer *n, doublereal *z__, 00024 integer *ldz, doublereal *rhs, doublereal *rdsum, doublereal *rdscal, 00025 integer *ipiv, integer *jpiv) 00026 { 00027 /* System generated locals */ 00028 integer z_dim1, z_offset, i__1, i__2; 00029 doublereal d__1; 00030 00031 /* Builtin functions */ 00032 double sqrt(doublereal); 00033 00034 /* Local variables */ 00035 integer i__, j, k; 00036 doublereal bm, bp, xm[8], xp[8]; 00037 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 00038 integer *); 00039 integer info; 00040 doublereal temp, work[32]; 00041 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 00042 integer *); 00043 extern doublereal dasum_(integer *, doublereal *, integer *); 00044 doublereal pmone; 00045 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 00046 doublereal *, integer *), daxpy_(integer *, doublereal *, 00047 doublereal *, integer *, doublereal *, integer *); 00048 doublereal sminu; 00049 integer iwork[8]; 00050 doublereal splus; 00051 extern /* Subroutine */ int dgesc2_(integer *, doublereal *, integer *, 00052 doublereal *, integer *, integer *, doublereal *), dgecon_(char *, 00053 integer *, doublereal *, integer *, doublereal *, doublereal *, 00054 doublereal *, integer *, integer *), dlassq_(integer *, 00055 doublereal *, integer *, doublereal *, doublereal *), dlaswp_( 00056 integer *, doublereal *, integer *, integer *, integer *, integer 00057 *, integer *); 00058 00059 00060 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00061 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00062 /* November 2006 */ 00063 00064 /* .. Scalar Arguments .. */ 00065 /* .. */ 00066 /* .. Array Arguments .. */ 00067 /* .. */ 00068 00069 /* Purpose */ 00070 /* ======= */ 00071 00072 /* DLATDF uses the LU factorization of the n-by-n matrix Z computed by */ 00073 /* DGETC2 and computes a contribution to the reciprocal Dif-estimate */ 00074 /* by solving Z * x = b for x, and choosing the r.h.s. b such that */ 00075 /* the norm of x is as large as possible. On entry RHS = b holds the */ 00076 /* contribution from earlier solved sub-systems, and on return RHS = x. */ 00077 00078 /* The factorization of Z returned by DGETC2 has the form Z = P*L*U*Q, */ 00079 /* where P and Q are permutation matrices. L is lower triangular with */ 00080 /* unit diagonal elements and U is upper triangular. */ 00081 00082 /* Arguments */ 00083 /* ========= */ 00084 00085 /* IJOB (input) INTEGER */ 00086 /* IJOB = 2: First compute an approximative null-vector e */ 00087 /* of Z using DGECON, e is normalized and solve for */ 00088 /* Zx = +-e - f with the sign giving the greater value */ 00089 /* of 2-norm(x). About 5 times as expensive as Default. */ 00090 /* IJOB .ne. 2: Local look ahead strategy where all entries of */ 00091 /* the r.h.s. b is choosen as either +1 or -1 (Default). */ 00092 00093 /* N (input) INTEGER */ 00094 /* The number of columns of the matrix Z. */ 00095 00096 /* Z (input) DOUBLE PRECISION array, dimension (LDZ, N) */ 00097 /* On entry, the LU part of the factorization of the n-by-n */ 00098 /* matrix Z computed by DGETC2: Z = P * L * U * Q */ 00099 00100 /* LDZ (input) INTEGER */ 00101 /* The leading dimension of the array Z. LDA >= max(1, N). */ 00102 00103 /* RHS (input/output) DOUBLE PRECISION array, dimension N. */ 00104 /* On entry, RHS contains contributions from other subsystems. */ 00105 /* On exit, RHS contains the solution of the subsystem with */ 00106 /* entries acoording to the value of IJOB (see above). */ 00107 00108 /* RDSUM (input/output) DOUBLE PRECISION */ 00109 /* On entry, the sum of squares of computed contributions to */ 00110 /* the Dif-estimate under computation by DTGSYL, where the */ 00111 /* scaling factor RDSCAL (see below) has been factored out. */ 00112 /* On exit, the corresponding sum of squares updated with the */ 00113 /* contributions from the current sub-system. */ 00114 /* If TRANS = 'T' RDSUM is not touched. */ 00115 /* NOTE: RDSUM only makes sense when DTGSY2 is called by STGSYL. */ 00116 00117 /* RDSCAL (input/output) DOUBLE PRECISION */ 00118 /* On entry, scaling factor used to prevent overflow in RDSUM. */ 00119 /* On exit, RDSCAL is updated w.r.t. the current contributions */ 00120 /* in RDSUM. */ 00121 /* If TRANS = 'T', RDSCAL is not touched. */ 00122 /* NOTE: RDSCAL only makes sense when DTGSY2 is called by */ 00123 /* DTGSYL. */ 00124 00125 /* IPIV (input) INTEGER array, dimension (N). */ 00126 /* The pivot indices; for 1 <= i <= N, row i of the */ 00127 /* matrix has been interchanged with row IPIV(i). */ 00128 00129 /* JPIV (input) INTEGER array, dimension (N). */ 00130 /* The pivot indices; for 1 <= j <= N, column j of the */ 00131 /* matrix has been interchanged with column JPIV(j). */ 00132 00133 /* Further Details */ 00134 /* =============== */ 00135 00136 /* Based on contributions by */ 00137 /* Bo Kagstrom and Peter Poromaa, Department of Computing Science, */ 00138 /* Umea University, S-901 87 Umea, Sweden. */ 00139 00140 /* This routine is a further developed implementation of algorithm */ 00141 /* BSOLVE in [1] using complete pivoting in the LU factorization. */ 00142 00143 /* [1] Bo Kagstrom and Lars Westin, */ 00144 /* Generalized Schur Methods with Condition Estimators for */ 00145 /* Solving the Generalized Sylvester Equation, IEEE Transactions */ 00146 /* on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751. */ 00147 00148 /* [2] Peter Poromaa, */ 00149 /* On Efficient and Robust Estimators for the Separation */ 00150 /* between two Regular Matrix Pairs with Applications in */ 00151 /* Condition Estimation. Report IMINF-95.05, Departement of */ 00152 /* Computing Science, Umea University, S-901 87 Umea, Sweden, 1995. */ 00153 00154 /* ===================================================================== */ 00155 00156 /* .. Parameters .. */ 00157 /* .. */ 00158 /* .. Local Scalars .. */ 00159 /* .. */ 00160 /* .. Local Arrays .. */ 00161 /* .. */ 00162 /* .. External Subroutines .. */ 00163 /* .. */ 00164 /* .. External Functions .. */ 00165 /* .. */ 00166 /* .. Intrinsic Functions .. */ 00167 /* .. */ 00168 /* .. Executable Statements .. */ 00169 00170 /* Parameter adjustments */ 00171 z_dim1 = *ldz; 00172 z_offset = 1 + z_dim1; 00173 z__ -= z_offset; 00174 --rhs; 00175 --ipiv; 00176 --jpiv; 00177 00178 /* Function Body */ 00179 if (*ijob != 2) { 00180 00181 /* Apply permutations IPIV to RHS */ 00182 00183 i__1 = *n - 1; 00184 dlaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &ipiv[1], &c__1); 00185 00186 /* Solve for L-part choosing RHS either to +1 or -1. */ 00187 00188 pmone = -1.; 00189 00190 i__1 = *n - 1; 00191 for (j = 1; j <= i__1; ++j) { 00192 bp = rhs[j] + 1.; 00193 bm = rhs[j] - 1.; 00194 splus = 1.; 00195 00196 /* Look-ahead for L-part RHS(1:N-1) = + or -1, SPLUS and */ 00197 /* SMIN computed more efficiently than in BSOLVE [1]. */ 00198 00199 i__2 = *n - j; 00200 splus += ddot_(&i__2, &z__[j + 1 + j * z_dim1], &c__1, &z__[j + 1 00201 + j * z_dim1], &c__1); 00202 i__2 = *n - j; 00203 sminu = ddot_(&i__2, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1], 00204 &c__1); 00205 splus *= rhs[j]; 00206 if (splus > sminu) { 00207 rhs[j] = bp; 00208 } else if (sminu > splus) { 00209 rhs[j] = bm; 00210 } else { 00211 00212 /* In this case the updating sums are equal and we can */ 00213 /* choose RHS(J) +1 or -1. The first time this happens */ 00214 /* we choose -1, thereafter +1. This is a simple way to */ 00215 /* get good estimates of matrices like Byers well-known */ 00216 /* example (see [1]). (Not done in BSOLVE.) */ 00217 00218 rhs[j] += pmone; 00219 pmone = 1.; 00220 } 00221 00222 /* Compute the remaining r.h.s. */ 00223 00224 temp = -rhs[j]; 00225 i__2 = *n - j; 00226 daxpy_(&i__2, &temp, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1], 00227 &c__1); 00228 00229 /* L10: */ 00230 } 00231 00232 /* Solve for U-part, look-ahead for RHS(N) = +-1. This is not done */ 00233 /* in BSOLVE and will hopefully give us a better estimate because */ 00234 /* any ill-conditioning of the original matrix is transfered to U */ 00235 /* and not to L. U(N, N) is an approximation to sigma_min(LU). */ 00236 00237 i__1 = *n - 1; 00238 dcopy_(&i__1, &rhs[1], &c__1, xp, &c__1); 00239 xp[*n - 1] = rhs[*n] + 1.; 00240 rhs[*n] += -1.; 00241 splus = 0.; 00242 sminu = 0.; 00243 for (i__ = *n; i__ >= 1; --i__) { 00244 temp = 1. / z__[i__ + i__ * z_dim1]; 00245 xp[i__ - 1] *= temp; 00246 rhs[i__] *= temp; 00247 i__1 = *n; 00248 for (k = i__ + 1; k <= i__1; ++k) { 00249 xp[i__ - 1] -= xp[k - 1] * (z__[i__ + k * z_dim1] * temp); 00250 rhs[i__] -= rhs[k] * (z__[i__ + k * z_dim1] * temp); 00251 /* L20: */ 00252 } 00253 splus += (d__1 = xp[i__ - 1], abs(d__1)); 00254 sminu += (d__1 = rhs[i__], abs(d__1)); 00255 /* L30: */ 00256 } 00257 if (splus > sminu) { 00258 dcopy_(n, xp, &c__1, &rhs[1], &c__1); 00259 } 00260 00261 /* Apply the permutations JPIV to the computed solution (RHS) */ 00262 00263 i__1 = *n - 1; 00264 dlaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &jpiv[1], &c_n1); 00265 00266 /* Compute the sum of squares */ 00267 00268 dlassq_(n, &rhs[1], &c__1, rdscal, rdsum); 00269 00270 } else { 00271 00272 /* IJOB = 2, Compute approximate nullvector XM of Z */ 00273 00274 dgecon_("I", n, &z__[z_offset], ldz, &c_b23, &temp, work, iwork, & 00275 info); 00276 dcopy_(n, &work[*n], &c__1, xm, &c__1); 00277 00278 /* Compute RHS */ 00279 00280 i__1 = *n - 1; 00281 dlaswp_(&c__1, xm, ldz, &c__1, &i__1, &ipiv[1], &c_n1); 00282 temp = 1. / sqrt(ddot_(n, xm, &c__1, xm, &c__1)); 00283 dscal_(n, &temp, xm, &c__1); 00284 dcopy_(n, xm, &c__1, xp, &c__1); 00285 daxpy_(n, &c_b23, &rhs[1], &c__1, xp, &c__1); 00286 daxpy_(n, &c_b37, xm, &c__1, &rhs[1], &c__1); 00287 dgesc2_(n, &z__[z_offset], ldz, &rhs[1], &ipiv[1], &jpiv[1], &temp); 00288 dgesc2_(n, &z__[z_offset], ldz, xp, &ipiv[1], &jpiv[1], &temp); 00289 if (dasum_(n, xp, &c__1) > dasum_(n, &rhs[1], &c__1)) { 00290 dcopy_(n, xp, &c__1, &rhs[1], &c__1); 00291 } 00292 00293 /* Compute the sum of squares */ 00294 00295 dlassq_(n, &rhs[1], &c__1, rdscal, rdsum); 00296 00297 } 00298 00299 return 0; 00300 00301 /* End of DLATDF */ 00302 00303 } /* dlatdf_ */