dlatdf.c 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304
  1. /* dlatdf.f -- translated by f2c (version 20061008).
  2. You must link the resulting object file with libf2c:
  3. on Microsoft Windows system, link with libf2c.lib;
  4. on Linux or Unix systems, link with .../path/to/libf2c.a -lm
  5. or, if you install libf2c.a in a standard place, with -lf2c -lm
  6. -- in that order, at the end of the command line, as in
  7. cc *.o -lf2c -lm
  8. Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
  9. http://www.netlib.org/f2c/libf2c.zip
  10. */
  11. #include "f2c.h"
  12. #include "blaswrap.h"
  13. /* Table of constant values */
  14. static integer c__1 = 1;
  15. static integer c_n1 = -1;
  16. static doublereal c_b23 = 1.;
  17. static doublereal c_b37 = -1.;
  18. /* Subroutine */ int dlatdf_(integer *ijob, integer *n, doublereal *z__,
  19. integer *ldz, doublereal *rhs, doublereal *rdsum, doublereal *rdscal,
  20. integer *ipiv, integer *jpiv)
  21. {
  22. /* System generated locals */
  23. integer z_dim1, z_offset, i__1, i__2;
  24. doublereal d__1;
  25. /* Builtin functions */
  26. double sqrt(doublereal);
  27. /* Local variables */
  28. integer i__, j, k;
  29. doublereal bm, bp, xm[8], xp[8];
  30. extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
  31. integer *);
  32. integer info;
  33. doublereal temp, work[32];
  34. extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
  35. integer *);
  36. extern doublereal dasum_(integer *, doublereal *, integer *);
  37. doublereal pmone;
  38. extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  39. doublereal *, integer *), daxpy_(integer *, doublereal *,
  40. doublereal *, integer *, doublereal *, integer *);
  41. doublereal sminu;
  42. integer iwork[8];
  43. doublereal splus;
  44. extern /* Subroutine */ int dgesc2_(integer *, doublereal *, integer *,
  45. doublereal *, integer *, integer *, doublereal *), dgecon_(char *,
  46. integer *, doublereal *, integer *, doublereal *, doublereal *,
  47. doublereal *, integer *, integer *), dlassq_(integer *,
  48. doublereal *, integer *, doublereal *, doublereal *), dlaswp_(
  49. integer *, doublereal *, integer *, integer *, integer *, integer
  50. *, integer *);
  51. /* -- LAPACK auxiliary routine (version 3.2) -- */
  52. /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  53. /* November 2006 */
  54. /* .. Scalar Arguments .. */
  55. /* .. */
  56. /* .. Array Arguments .. */
  57. /* .. */
  58. /* Purpose */
  59. /* ======= */
  60. /* DLATDF uses the LU factorization of the n-by-n matrix Z computed by */
  61. /* DGETC2 and computes a contribution to the reciprocal Dif-estimate */
  62. /* by solving Z * x = b for x, and choosing the r.h.s. b such that */
  63. /* the norm of x is as large as possible. On entry RHS = b holds the */
  64. /* contribution from earlier solved sub-systems, and on return RHS = x. */
  65. /* The factorization of Z returned by DGETC2 has the form Z = P*L*U*Q, */
  66. /* where P and Q are permutation matrices. L is lower triangular with */
  67. /* unit diagonal elements and U is upper triangular. */
  68. /* Arguments */
  69. /* ========= */
  70. /* IJOB (input) INTEGER */
  71. /* IJOB = 2: First compute an approximative null-vector e */
  72. /* of Z using DGECON, e is normalized and solve for */
  73. /* Zx = +-e - f with the sign giving the greater value */
  74. /* of 2-norm(x). About 5 times as expensive as Default. */
  75. /* IJOB .ne. 2: Local look ahead strategy where all entries of */
  76. /* the r.h.s. b is choosen as either +1 or -1 (Default). */
  77. /* N (input) INTEGER */
  78. /* The number of columns of the matrix Z. */
  79. /* Z (input) DOUBLE PRECISION array, dimension (LDZ, N) */
  80. /* On entry, the LU part of the factorization of the n-by-n */
  81. /* matrix Z computed by DGETC2: Z = P * L * U * Q */
  82. /* LDZ (input) INTEGER */
  83. /* The leading dimension of the array Z. LDA >= max(1, N). */
  84. /* RHS (input/output) DOUBLE PRECISION array, dimension N. */
  85. /* On entry, RHS contains contributions from other subsystems. */
  86. /* On exit, RHS contains the solution of the subsystem with */
  87. /* entries acoording to the value of IJOB (see above). */
  88. /* RDSUM (input/output) DOUBLE PRECISION */
  89. /* On entry, the sum of squares of computed contributions to */
  90. /* the Dif-estimate under computation by DTGSYL, where the */
  91. /* scaling factor RDSCAL (see below) has been factored out. */
  92. /* On exit, the corresponding sum of squares updated with the */
  93. /* contributions from the current sub-system. */
  94. /* If TRANS = 'T' RDSUM is not touched. */
  95. /* NOTE: RDSUM only makes sense when DTGSY2 is called by STGSYL. */
  96. /* RDSCAL (input/output) DOUBLE PRECISION */
  97. /* On entry, scaling factor used to prevent overflow in RDSUM. */
  98. /* On exit, RDSCAL is updated w.r.t. the current contributions */
  99. /* in RDSUM. */
  100. /* If TRANS = 'T', RDSCAL is not touched. */
  101. /* NOTE: RDSCAL only makes sense when DTGSY2 is called by */
  102. /* DTGSYL. */
  103. /* IPIV (input) INTEGER array, dimension (N). */
  104. /* The pivot indices; for 1 <= i <= N, row i of the */
  105. /* matrix has been interchanged with row IPIV(i). */
  106. /* JPIV (input) INTEGER array, dimension (N). */
  107. /* The pivot indices; for 1 <= j <= N, column j of the */
  108. /* matrix has been interchanged with column JPIV(j). */
  109. /* Further Details */
  110. /* =============== */
  111. /* Based on contributions by */
  112. /* Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
  113. /* Umea University, S-901 87 Umea, Sweden. */
  114. /* This routine is a further developed implementation of algorithm */
  115. /* BSOLVE in [1] using complete pivoting in the LU factorization. */
  116. /* [1] Bo Kagstrom and Lars Westin, */
  117. /* Generalized Schur Methods with Condition Estimators for */
  118. /* Solving the Generalized Sylvester Equation, IEEE Transactions */
  119. /* on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751. */
  120. /* [2] Peter Poromaa, */
  121. /* On Efficient and Robust Estimators for the Separation */
  122. /* between two Regular Matrix Pairs with Applications in */
  123. /* Condition Estimation. Report IMINF-95.05, Departement of */
  124. /* Computing Science, Umea University, S-901 87 Umea, Sweden, 1995. */
  125. /* ===================================================================== */
  126. /* .. Parameters .. */
  127. /* .. */
  128. /* .. Local Scalars .. */
  129. /* .. */
  130. /* .. Local Arrays .. */
  131. /* .. */
  132. /* .. External Subroutines .. */
  133. /* .. */
  134. /* .. External Functions .. */
  135. /* .. */
  136. /* .. Intrinsic Functions .. */
  137. /* .. */
  138. /* .. Executable Statements .. */
  139. /* Parameter adjustments */
  140. z_dim1 = *ldz;
  141. z_offset = 1 + z_dim1;
  142. z__ -= z_offset;
  143. --rhs;
  144. --ipiv;
  145. --jpiv;
  146. /* Function Body */
  147. if (*ijob != 2) {
  148. /* Apply permutations IPIV to RHS */
  149. i__1 = *n - 1;
  150. dlaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &ipiv[1], &c__1);
  151. /* Solve for L-part choosing RHS either to +1 or -1. */
  152. pmone = -1.;
  153. i__1 = *n - 1;
  154. for (j = 1; j <= i__1; ++j) {
  155. bp = rhs[j] + 1.;
  156. bm = rhs[j] - 1.;
  157. splus = 1.;
  158. /* Look-ahead for L-part RHS(1:N-1) = + or -1, SPLUS and */
  159. /* SMIN computed more efficiently than in BSOLVE [1]. */
  160. i__2 = *n - j;
  161. splus += ddot_(&i__2, &z__[j + 1 + j * z_dim1], &c__1, &z__[j + 1
  162. + j * z_dim1], &c__1);
  163. i__2 = *n - j;
  164. sminu = ddot_(&i__2, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1],
  165. &c__1);
  166. splus *= rhs[j];
  167. if (splus > sminu) {
  168. rhs[j] = bp;
  169. } else if (sminu > splus) {
  170. rhs[j] = bm;
  171. } else {
  172. /* In this case the updating sums are equal and we can */
  173. /* choose RHS(J) +1 or -1. The first time this happens */
  174. /* we choose -1, thereafter +1. This is a simple way to */
  175. /* get good estimates of matrices like Byers well-known */
  176. /* example (see [1]). (Not done in BSOLVE.) */
  177. rhs[j] += pmone;
  178. pmone = 1.;
  179. }
  180. /* Compute the remaining r.h.s. */
  181. temp = -rhs[j];
  182. i__2 = *n - j;
  183. daxpy_(&i__2, &temp, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1],
  184. &c__1);
  185. /* L10: */
  186. }
  187. /* Solve for U-part, look-ahead for RHS(N) = +-1. This is not done */
  188. /* in BSOLVE and will hopefully give us a better estimate because */
  189. /* any ill-conditioning of the original matrix is transfered to U */
  190. /* and not to L. U(N, N) is an approximation to sigma_min(LU). */
  191. i__1 = *n - 1;
  192. dcopy_(&i__1, &rhs[1], &c__1, xp, &c__1);
  193. xp[*n - 1] = rhs[*n] + 1.;
  194. rhs[*n] += -1.;
  195. splus = 0.;
  196. sminu = 0.;
  197. for (i__ = *n; i__ >= 1; --i__) {
  198. temp = 1. / z__[i__ + i__ * z_dim1];
  199. xp[i__ - 1] *= temp;
  200. rhs[i__] *= temp;
  201. i__1 = *n;
  202. for (k = i__ + 1; k <= i__1; ++k) {
  203. xp[i__ - 1] -= xp[k - 1] * (z__[i__ + k * z_dim1] * temp);
  204. rhs[i__] -= rhs[k] * (z__[i__ + k * z_dim1] * temp);
  205. /* L20: */
  206. }
  207. splus += (d__1 = xp[i__ - 1], abs(d__1));
  208. sminu += (d__1 = rhs[i__], abs(d__1));
  209. /* L30: */
  210. }
  211. if (splus > sminu) {
  212. dcopy_(n, xp, &c__1, &rhs[1], &c__1);
  213. }
  214. /* Apply the permutations JPIV to the computed solution (RHS) */
  215. i__1 = *n - 1;
  216. dlaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &jpiv[1], &c_n1);
  217. /* Compute the sum of squares */
  218. dlassq_(n, &rhs[1], &c__1, rdscal, rdsum);
  219. } else {
  220. /* IJOB = 2, Compute approximate nullvector XM of Z */
  221. dgecon_("I", n, &z__[z_offset], ldz, &c_b23, &temp, work, iwork, &
  222. info);
  223. dcopy_(n, &work[*n], &c__1, xm, &c__1);
  224. /* Compute RHS */
  225. i__1 = *n - 1;
  226. dlaswp_(&c__1, xm, ldz, &c__1, &i__1, &ipiv[1], &c_n1);
  227. temp = 1. / sqrt(ddot_(n, xm, &c__1, xm, &c__1));
  228. dscal_(n, &temp, xm, &c__1);
  229. dcopy_(n, xm, &c__1, xp, &c__1);
  230. daxpy_(n, &c_b23, &rhs[1], &c__1, xp, &c__1);
  231. daxpy_(n, &c_b37, xm, &c__1, &rhs[1], &c__1);
  232. dgesc2_(n, &z__[z_offset], ldz, &rhs[1], &ipiv[1], &jpiv[1], &temp);
  233. dgesc2_(n, &z__[z_offset], ldz, xp, &ipiv[1], &jpiv[1], &temp);
  234. if (dasum_(n, xp, &c__1) > dasum_(n, &rhs[1], &c__1)) {
  235. dcopy_(n, xp, &c__1, &rhs[1], &c__1);
  236. }
  237. /* Compute the sum of squares */
  238. dlassq_(n, &rhs[1], &c__1, rdscal, rdsum);
  239. }
  240. return 0;
  241. /* End of DLATDF */
  242. } /* dlatdf_ */