zlatdf.c
Go to the documentation of this file.
00001 /* zlatdf.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 doublecomplex c_b1 = {1.,0.};
00019 static integer c__1 = 1;
00020 static integer c_n1 = -1;
00021 static doublereal c_b24 = 1.;
00022 
00023 /* Subroutine */ int zlatdf_(integer *ijob, integer *n, doublecomplex *z__, 
00024         integer *ldz, doublecomplex *rhs, doublereal *rdsum, doublereal *
00025         rdscal, integer *ipiv, integer *jpiv)
00026 {
00027     /* System generated locals */
00028     integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;
00029     doublecomplex z__1, z__2, z__3;
00030 
00031     /* Builtin functions */
00032     void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
00033     double z_abs(doublecomplex *);
00034     void z_sqrt(doublecomplex *, doublecomplex *);
00035 
00036     /* Local variables */
00037     integer i__, j, k;
00038     doublecomplex bm, bp, xm[2], xp[2];
00039     integer info;
00040     doublecomplex temp, work[8];
00041     doublereal scale;
00042     extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 
00043             doublecomplex *, integer *);
00044     doublecomplex pmone;
00045     extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 
00046             doublecomplex *, integer *, doublecomplex *, integer *);
00047     doublereal rtemp, sminu, rwork[2];
00048     extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 
00049             doublecomplex *, integer *);
00050     doublereal splus;
00051     extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *, 
00052             doublecomplex *, integer *, doublecomplex *, integer *), zgesc2_(
00053             integer *, doublecomplex *, integer *, doublecomplex *, integer *, 
00054              integer *, doublereal *), zgecon_(char *, integer *, 
00055             doublecomplex *, integer *, doublereal *, doublereal *, 
00056             doublecomplex *, doublereal *, integer *);
00057     extern doublereal dzasum_(integer *, doublecomplex *, integer *);
00058     extern /* Subroutine */ int zlassq_(integer *, doublecomplex *, integer *, 
00059              doublereal *, doublereal *), zlaswp_(integer *, doublecomplex *, 
00060             integer *, integer *, integer *, integer *, integer *);
00061 
00062 
00063 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00064 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00065 /*     November 2006 */
00066 
00067 /*     .. Scalar Arguments .. */
00068 /*     .. */
00069 /*     .. Array Arguments .. */
00070 /*     .. */
00071 
00072 /*  Purpose */
00073 /*  ======= */
00074 
00075 /*  ZLATDF computes the contribution to the reciprocal Dif-estimate */
00076 /*  by solving for x in Z * x = b, where b is chosen such that the norm */
00077 /*  of x is as large as possible. It is assumed that LU decomposition */
00078 /*  of Z has been computed by ZGETC2. On entry RHS = f holds the */
00079 /*  contribution from earlier solved sub-systems, and on return RHS = x. */
00080 
00081 /*  The factorization of Z returned by ZGETC2 has the form */
00082 /*  Z = P * L * U * Q, where P and Q are permutation matrices. L is lower */
00083 /*  triangular with unit diagonal elements and U is upper triangular. */
00084 
00085 /*  Arguments */
00086 /*  ========= */
00087 
00088 /*  IJOB    (input) INTEGER */
00089 /*          IJOB = 2: First compute an approximative null-vector e */
00090 /*              of Z using ZGECON, e is normalized and solve for */
00091 /*              Zx = +-e - f with the sign giving the greater value of */
00092 /*              2-norm(x).  About 5 times as expensive as Default. */
00093 /*          IJOB .ne. 2: Local look ahead strategy where */
00094 /*              all entries of the r.h.s. b is choosen as either +1 or */
00095 /*              -1.  Default. */
00096 
00097 /*  N       (input) INTEGER */
00098 /*          The number of columns of the matrix Z. */
00099 
00100 /*  Z       (input) DOUBLE PRECISION array, dimension (LDZ, N) */
00101 /*          On entry, the LU part of the factorization of the n-by-n */
00102 /*          matrix Z computed by ZGETC2:  Z = P * L * U * Q */
00103 
00104 /*  LDZ     (input) INTEGER */
00105 /*          The leading dimension of the array Z.  LDA >= max(1, N). */
00106 
00107 /*  RHS     (input/output) DOUBLE PRECISION array, dimension (N). */
00108 /*          On entry, RHS contains contributions from other subsystems. */
00109 /*          On exit, RHS contains the solution of the subsystem with */
00110 /*          entries according to the value of IJOB (see above). */
00111 
00112 /*  RDSUM   (input/output) DOUBLE PRECISION */
00113 /*          On entry, the sum of squares of computed contributions to */
00114 /*          the Dif-estimate under computation by ZTGSYL, where the */
00115 /*          scaling factor RDSCAL (see below) has been factored out. */
00116 /*          On exit, the corresponding sum of squares updated with the */
00117 /*          contributions from the current sub-system. */
00118 /*          If TRANS = 'T' RDSUM is not touched. */
00119 /*          NOTE: RDSUM only makes sense when ZTGSY2 is called by CTGSYL. */
00120 
00121 /*  RDSCAL  (input/output) DOUBLE PRECISION */
00122 /*          On entry, scaling factor used to prevent overflow in RDSUM. */
00123 /*          On exit, RDSCAL is updated w.r.t. the current contributions */
00124 /*          in RDSUM. */
00125 /*          If TRANS = 'T', RDSCAL is not touched. */
00126 /*          NOTE: RDSCAL only makes sense when ZTGSY2 is called by */
00127 /*          ZTGSYL. */
00128 
00129 /*  IPIV    (input) INTEGER array, dimension (N). */
00130 /*          The pivot indices; for 1 <= i <= N, row i of the */
00131 /*          matrix has been interchanged with row IPIV(i). */
00132 
00133 /*  JPIV    (input) INTEGER array, dimension (N). */
00134 /*          The pivot indices; for 1 <= j <= N, column j of the */
00135 /*          matrix has been interchanged with column JPIV(j). */
00136 
00137 /*  Further Details */
00138 /*  =============== */
00139 
00140 /*  Based on contributions by */
00141 /*     Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
00142 /*     Umea University, S-901 87 Umea, Sweden. */
00143 
00144 /*  This routine is a further developed implementation of algorithm */
00145 /*  BSOLVE in [1] using complete pivoting in the LU factorization. */
00146 
00147 /*   [1]   Bo Kagstrom and Lars Westin, */
00148 /*         Generalized Schur Methods with Condition Estimators for */
00149 /*         Solving the Generalized Sylvester Equation, IEEE Transactions */
00150 /*         on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751. */
00151 
00152 /*   [2]   Peter Poromaa, */
00153 /*         On Efficient and Robust Estimators for the Separation */
00154 /*         between two Regular Matrix Pairs with Applications in */
00155 /*         Condition Estimation. Report UMINF-95.05, Department of */
00156 /*         Computing Science, Umea University, S-901 87 Umea, Sweden, */
00157 /*         1995. */
00158 
00159 /*  ===================================================================== */
00160 
00161 /*     .. Parameters .. */
00162 /*     .. */
00163 /*     .. Local Scalars .. */
00164 /*     .. */
00165 /*     .. Local Arrays .. */
00166 /*     .. */
00167 /*     .. External Subroutines .. */
00168 /*     .. */
00169 /*     .. External Functions .. */
00170 /*     .. */
00171 /*     .. Intrinsic Functions .. */
00172 /*     .. */
00173 /*     .. Executable Statements .. */
00174 
00175     /* Parameter adjustments */
00176     z_dim1 = *ldz;
00177     z_offset = 1 + z_dim1;
00178     z__ -= z_offset;
00179     --rhs;
00180     --ipiv;
00181     --jpiv;
00182 
00183     /* Function Body */
00184     if (*ijob != 2) {
00185 
00186 /*        Apply permutations IPIV to RHS */
00187 
00188         i__1 = *n - 1;
00189         zlaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &ipiv[1], &c__1);
00190 
00191 /*        Solve for L-part choosing RHS either to +1 or -1. */
00192 
00193         z__1.r = -1., z__1.i = -0.;
00194         pmone.r = z__1.r, pmone.i = z__1.i;
00195         i__1 = *n - 1;
00196         for (j = 1; j <= i__1; ++j) {
00197             i__2 = j;
00198             z__1.r = rhs[i__2].r + 1., z__1.i = rhs[i__2].i + 0.;
00199             bp.r = z__1.r, bp.i = z__1.i;
00200             i__2 = j;
00201             z__1.r = rhs[i__2].r - 1., z__1.i = rhs[i__2].i - 0.;
00202             bm.r = z__1.r, bm.i = z__1.i;
00203             splus = 1.;
00204 
00205 /*           Lockahead for L- part RHS(1:N-1) = +-1 */
00206 /*           SPLUS and SMIN computed more efficiently than in BSOLVE[1]. */
00207 
00208             i__2 = *n - j;
00209             zdotc_(&z__1, &i__2, &z__[j + 1 + j * z_dim1], &c__1, &z__[j + 1 
00210                     + j * z_dim1], &c__1);
00211             splus += z__1.r;
00212             i__2 = *n - j;
00213             zdotc_(&z__1, &i__2, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1], 
00214                      &c__1);
00215             sminu = z__1.r;
00216             i__2 = j;
00217             splus *= rhs[i__2].r;
00218             if (splus > sminu) {
00219                 i__2 = j;
00220                 rhs[i__2].r = bp.r, rhs[i__2].i = bp.i;
00221             } else if (sminu > splus) {
00222                 i__2 = j;
00223                 rhs[i__2].r = bm.r, rhs[i__2].i = bm.i;
00224             } else {
00225 
00226 /*              In this case the updating sums are equal and we can */
00227 /*              choose RHS(J) +1 or -1. The first time this happens we */
00228 /*              choose -1, thereafter +1. This is a simple way to get */
00229 /*              good estimates of matrices like Byers well-known example */
00230 /*              (see [1]). (Not done in BSOLVE.) */
00231 
00232                 i__2 = j;
00233                 i__3 = j;
00234                 z__1.r = rhs[i__3].r + pmone.r, z__1.i = rhs[i__3].i + 
00235                         pmone.i;
00236                 rhs[i__2].r = z__1.r, rhs[i__2].i = z__1.i;
00237                 pmone.r = 1., pmone.i = 0.;
00238             }
00239 
00240 /*           Compute the remaining r.h.s. */
00241 
00242             i__2 = j;
00243             z__1.r = -rhs[i__2].r, z__1.i = -rhs[i__2].i;
00244             temp.r = z__1.r, temp.i = z__1.i;
00245             i__2 = *n - j;
00246             zaxpy_(&i__2, &temp, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1], 
00247                      &c__1);
00248 /* L10: */
00249         }
00250 
00251 /*        Solve for U- part, lockahead for RHS(N) = +-1. This is not done */
00252 /*        In BSOLVE and will hopefully give us a better estimate because */
00253 /*        any ill-conditioning of the original matrix is transfered to U */
00254 /*        and not to L. U(N, N) is an approximation to sigma_min(LU). */
00255 
00256         i__1 = *n - 1;
00257         zcopy_(&i__1, &rhs[1], &c__1, work, &c__1);
00258         i__1 = *n - 1;
00259         i__2 = *n;
00260         z__1.r = rhs[i__2].r + 1., z__1.i = rhs[i__2].i + 0.;
00261         work[i__1].r = z__1.r, work[i__1].i = z__1.i;
00262         i__1 = *n;
00263         i__2 = *n;
00264         z__1.r = rhs[i__2].r - 1., z__1.i = rhs[i__2].i - 0.;
00265         rhs[i__1].r = z__1.r, rhs[i__1].i = z__1.i;
00266         splus = 0.;
00267         sminu = 0.;
00268         for (i__ = *n; i__ >= 1; --i__) {
00269             z_div(&z__1, &c_b1, &z__[i__ + i__ * z_dim1]);
00270             temp.r = z__1.r, temp.i = z__1.i;
00271             i__1 = i__ - 1;
00272             i__2 = i__ - 1;
00273             z__1.r = work[i__2].r * temp.r - work[i__2].i * temp.i, z__1.i = 
00274                     work[i__2].r * temp.i + work[i__2].i * temp.r;
00275             work[i__1].r = z__1.r, work[i__1].i = z__1.i;
00276             i__1 = i__;
00277             i__2 = i__;
00278             z__1.r = rhs[i__2].r * temp.r - rhs[i__2].i * temp.i, z__1.i = 
00279                     rhs[i__2].r * temp.i + rhs[i__2].i * temp.r;
00280             rhs[i__1].r = z__1.r, rhs[i__1].i = z__1.i;
00281             i__1 = *n;
00282             for (k = i__ + 1; k <= i__1; ++k) {
00283                 i__2 = i__ - 1;
00284                 i__3 = i__ - 1;
00285                 i__4 = k - 1;
00286                 i__5 = i__ + k * z_dim1;
00287                 z__3.r = z__[i__5].r * temp.r - z__[i__5].i * temp.i, z__3.i =
00288                          z__[i__5].r * temp.i + z__[i__5].i * temp.r;
00289                 z__2.r = work[i__4].r * z__3.r - work[i__4].i * z__3.i, 
00290                         z__2.i = work[i__4].r * z__3.i + work[i__4].i * 
00291                         z__3.r;
00292                 z__1.r = work[i__3].r - z__2.r, z__1.i = work[i__3].i - 
00293                         z__2.i;
00294                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00295                 i__2 = i__;
00296                 i__3 = i__;
00297                 i__4 = k;
00298                 i__5 = i__ + k * z_dim1;
00299                 z__3.r = z__[i__5].r * temp.r - z__[i__5].i * temp.i, z__3.i =
00300                          z__[i__5].r * temp.i + z__[i__5].i * temp.r;
00301                 z__2.r = rhs[i__4].r * z__3.r - rhs[i__4].i * z__3.i, z__2.i =
00302                          rhs[i__4].r * z__3.i + rhs[i__4].i * z__3.r;
00303                 z__1.r = rhs[i__3].r - z__2.r, z__1.i = rhs[i__3].i - z__2.i;
00304                 rhs[i__2].r = z__1.r, rhs[i__2].i = z__1.i;
00305 /* L20: */
00306             }
00307             splus += z_abs(&work[i__ - 1]);
00308             sminu += z_abs(&rhs[i__]);
00309 /* L30: */
00310         }
00311         if (splus > sminu) {
00312             zcopy_(n, work, &c__1, &rhs[1], &c__1);
00313         }
00314 
00315 /*        Apply the permutations JPIV to the computed solution (RHS) */
00316 
00317         i__1 = *n - 1;
00318         zlaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &jpiv[1], &c_n1);
00319 
00320 /*        Compute the sum of squares */
00321 
00322         zlassq_(n, &rhs[1], &c__1, rdscal, rdsum);
00323         return 0;
00324     }
00325 
00326 /*     ENTRY IJOB = 2 */
00327 
00328 /*     Compute approximate nullvector XM of Z */
00329 
00330     zgecon_("I", n, &z__[z_offset], ldz, &c_b24, &rtemp, work, rwork, &info);
00331     zcopy_(n, &work[*n], &c__1, xm, &c__1);
00332 
00333 /*     Compute RHS */
00334 
00335     i__1 = *n - 1;
00336     zlaswp_(&c__1, xm, ldz, &c__1, &i__1, &ipiv[1], &c_n1);
00337     zdotc_(&z__3, n, xm, &c__1, xm, &c__1);
00338     z_sqrt(&z__2, &z__3);
00339     z_div(&z__1, &c_b1, &z__2);
00340     temp.r = z__1.r, temp.i = z__1.i;
00341     zscal_(n, &temp, xm, &c__1);
00342     zcopy_(n, xm, &c__1, xp, &c__1);
00343     zaxpy_(n, &c_b1, &rhs[1], &c__1, xp, &c__1);
00344     z__1.r = -1., z__1.i = -0.;
00345     zaxpy_(n, &z__1, xm, &c__1, &rhs[1], &c__1);
00346     zgesc2_(n, &z__[z_offset], ldz, &rhs[1], &ipiv[1], &jpiv[1], &scale);
00347     zgesc2_(n, &z__[z_offset], ldz, xp, &ipiv[1], &jpiv[1], &scale);
00348     if (dzasum_(n, xp, &c__1) > dzasum_(n, &rhs[1], &c__1)) {
00349         zcopy_(n, xp, &c__1, &rhs[1], &c__1);
00350     }
00351 
00352 /*     Compute the sum of squares */
00353 
00354     zlassq_(n, &rhs[1], &c__1, rdscal, rdsum);
00355     return 0;
00356 
00357 /*     End of ZLATDF */
00358 
00359 } /* zlatdf_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:42