zlaic1.c
Go to the documentation of this file.
00001 /* zlaic1.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 
00020 /* Subroutine */ int zlaic1_(integer *job, integer *j, doublecomplex *x, 
00021         doublereal *sest, doublecomplex *w, doublecomplex *gamma, doublereal *
00022         sestpr, doublecomplex *s, doublecomplex *c__)
00023 {
00024     /* System generated locals */
00025     doublereal d__1, d__2;
00026     doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;
00027 
00028     /* Builtin functions */
00029     double z_abs(doublecomplex *);
00030     void d_cnjg(doublecomplex *, doublecomplex *), z_sqrt(doublecomplex *, 
00031             doublecomplex *);
00032     double sqrt(doublereal);
00033     void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
00034 
00035     /* Local variables */
00036     doublereal b, t, s1, s2, scl, eps, tmp;
00037     doublecomplex sine;
00038     doublereal test, zeta1, zeta2;
00039     doublecomplex alpha;
00040     doublereal norma;
00041     extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 
00042             doublecomplex *, integer *, doublecomplex *, integer *);
00043     extern doublereal dlamch_(char *);
00044     doublereal absgam, absalp;
00045     doublecomplex cosine;
00046     doublereal absest;
00047 
00048 
00049 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00050 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00051 /*     November 2006 */
00052 
00053 /*     .. Scalar Arguments .. */
00054 /*     .. */
00055 /*     .. Array Arguments .. */
00056 /*     .. */
00057 
00058 /*  Purpose */
00059 /*  ======= */
00060 
00061 /*  ZLAIC1 applies one step of incremental condition estimation in */
00062 /*  its simplest version: */
00063 
00064 /*  Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j */
00065 /*  lower triangular matrix L, such that */
00066 /*           twonorm(L*x) = sest */
00067 /*  Then ZLAIC1 computes sestpr, s, c such that */
00068 /*  the vector */
00069 /*                  [ s*x ] */
00070 /*           xhat = [  c  ] */
00071 /*  is an approximate singular vector of */
00072 /*                  [ L     0  ] */
00073 /*           Lhat = [ w' gamma ] */
00074 /*  in the sense that */
00075 /*           twonorm(Lhat*xhat) = sestpr. */
00076 
00077 /*  Depending on JOB, an estimate for the largest or smallest singular */
00078 /*  value is computed. */
00079 
00080 /*  Note that [s c]' and sestpr**2 is an eigenpair of the system */
00081 
00082 /*      diag(sest*sest, 0) + [alpha  gamma] * [ conjg(alpha) ] */
00083 /*                                            [ conjg(gamma) ] */
00084 
00085 /*  where  alpha =  conjg(x)'*w. */
00086 
00087 /*  Arguments */
00088 /*  ========= */
00089 
00090 /*  JOB     (input) INTEGER */
00091 /*          = 1: an estimate for the largest singular value is computed. */
00092 /*          = 2: an estimate for the smallest singular value is computed. */
00093 
00094 /*  J       (input) INTEGER */
00095 /*          Length of X and W */
00096 
00097 /*  X       (input) COMPLEX*16 array, dimension (J) */
00098 /*          The j-vector x. */
00099 
00100 /*  SEST    (input) DOUBLE PRECISION */
00101 /*          Estimated singular value of j by j matrix L */
00102 
00103 /*  W       (input) COMPLEX*16 array, dimension (J) */
00104 /*          The j-vector w. */
00105 
00106 /*  GAMMA   (input) COMPLEX*16 */
00107 /*          The diagonal element gamma. */
00108 
00109 /*  SESTPR  (output) DOUBLE PRECISION */
00110 /*          Estimated singular value of (j+1) by (j+1) matrix Lhat. */
00111 
00112 /*  S       (output) COMPLEX*16 */
00113 /*          Sine needed in forming xhat. */
00114 
00115 /*  C       (output) COMPLEX*16 */
00116 /*          Cosine needed in forming xhat. */
00117 
00118 /*  ===================================================================== */
00119 
00120 /*     .. Parameters .. */
00121 /*     .. */
00122 /*     .. Local Scalars .. */
00123 /*     .. */
00124 /*     .. Intrinsic Functions .. */
00125 /*     .. */
00126 /*     .. External Functions .. */
00127 /*     .. */
00128 /*     .. Executable Statements .. */
00129 
00130     /* Parameter adjustments */
00131     --w;
00132     --x;
00133 
00134     /* Function Body */
00135     eps = dlamch_("Epsilon");
00136     zdotc_(&z__1, j, &x[1], &c__1, &w[1], &c__1);
00137     alpha.r = z__1.r, alpha.i = z__1.i;
00138 
00139     absalp = z_abs(&alpha);
00140     absgam = z_abs(gamma);
00141     absest = abs(*sest);
00142 
00143     if (*job == 1) {
00144 
00145 /*        Estimating largest singular value */
00146 
00147 /*        special cases */
00148 
00149         if (*sest == 0.) {
00150             s1 = max(absgam,absalp);
00151             if (s1 == 0.) {
00152                 s->r = 0., s->i = 0.;
00153                 c__->r = 1., c__->i = 0.;
00154                 *sestpr = 0.;
00155             } else {
00156                 z__1.r = alpha.r / s1, z__1.i = alpha.i / s1;
00157                 s->r = z__1.r, s->i = z__1.i;
00158                 z__1.r = gamma->r / s1, z__1.i = gamma->i / s1;
00159                 c__->r = z__1.r, c__->i = z__1.i;
00160                 d_cnjg(&z__4, s);
00161                 z__3.r = s->r * z__4.r - s->i * z__4.i, z__3.i = s->r * 
00162                         z__4.i + s->i * z__4.r;
00163                 d_cnjg(&z__6, c__);
00164                 z__5.r = c__->r * z__6.r - c__->i * z__6.i, z__5.i = c__->r * 
00165                         z__6.i + c__->i * z__6.r;
00166                 z__2.r = z__3.r + z__5.r, z__2.i = z__3.i + z__5.i;
00167                 z_sqrt(&z__1, &z__2);
00168                 tmp = z__1.r;
00169                 z__1.r = s->r / tmp, z__1.i = s->i / tmp;
00170                 s->r = z__1.r, s->i = z__1.i;
00171                 z__1.r = c__->r / tmp, z__1.i = c__->i / tmp;
00172                 c__->r = z__1.r, c__->i = z__1.i;
00173                 *sestpr = s1 * tmp;
00174             }
00175             return 0;
00176         } else if (absgam <= eps * absest) {
00177             s->r = 1., s->i = 0.;
00178             c__->r = 0., c__->i = 0.;
00179             tmp = max(absest,absalp);
00180             s1 = absest / tmp;
00181             s2 = absalp / tmp;
00182             *sestpr = tmp * sqrt(s1 * s1 + s2 * s2);
00183             return 0;
00184         } else if (absalp <= eps * absest) {
00185             s1 = absgam;
00186             s2 = absest;
00187             if (s1 <= s2) {
00188                 s->r = 1., s->i = 0.;
00189                 c__->r = 0., c__->i = 0.;
00190                 *sestpr = s2;
00191             } else {
00192                 s->r = 0., s->i = 0.;
00193                 c__->r = 1., c__->i = 0.;
00194                 *sestpr = s1;
00195             }
00196             return 0;
00197         } else if (absest <= eps * absalp || absest <= eps * absgam) {
00198             s1 = absgam;
00199             s2 = absalp;
00200             if (s1 <= s2) {
00201                 tmp = s1 / s2;
00202                 scl = sqrt(tmp * tmp + 1.);
00203                 *sestpr = s2 * scl;
00204                 z__2.r = alpha.r / s2, z__2.i = alpha.i / s2;
00205                 z__1.r = z__2.r / scl, z__1.i = z__2.i / scl;
00206                 s->r = z__1.r, s->i = z__1.i;
00207                 z__2.r = gamma->r / s2, z__2.i = gamma->i / s2;
00208                 z__1.r = z__2.r / scl, z__1.i = z__2.i / scl;
00209                 c__->r = z__1.r, c__->i = z__1.i;
00210             } else {
00211                 tmp = s2 / s1;
00212                 scl = sqrt(tmp * tmp + 1.);
00213                 *sestpr = s1 * scl;
00214                 z__2.r = alpha.r / s1, z__2.i = alpha.i / s1;
00215                 z__1.r = z__2.r / scl, z__1.i = z__2.i / scl;
00216                 s->r = z__1.r, s->i = z__1.i;
00217                 z__2.r = gamma->r / s1, z__2.i = gamma->i / s1;
00218                 z__1.r = z__2.r / scl, z__1.i = z__2.i / scl;
00219                 c__->r = z__1.r, c__->i = z__1.i;
00220             }
00221             return 0;
00222         } else {
00223 
00224 /*           normal case */
00225 
00226             zeta1 = absalp / absest;
00227             zeta2 = absgam / absest;
00228 
00229             b = (1. - zeta1 * zeta1 - zeta2 * zeta2) * .5;
00230             d__1 = zeta1 * zeta1;
00231             c__->r = d__1, c__->i = 0.;
00232             if (b > 0.) {
00233                 d__1 = b * b;
00234                 z__4.r = d__1 + c__->r, z__4.i = c__->i;
00235                 z_sqrt(&z__3, &z__4);
00236                 z__2.r = b + z__3.r, z__2.i = z__3.i;
00237                 z_div(&z__1, c__, &z__2);
00238                 t = z__1.r;
00239             } else {
00240                 d__1 = b * b;
00241                 z__3.r = d__1 + c__->r, z__3.i = c__->i;
00242                 z_sqrt(&z__2, &z__3);
00243                 z__1.r = z__2.r - b, z__1.i = z__2.i;
00244                 t = z__1.r;
00245             }
00246 
00247             z__3.r = alpha.r / absest, z__3.i = alpha.i / absest;
00248             z__2.r = -z__3.r, z__2.i = -z__3.i;
00249             z__1.r = z__2.r / t, z__1.i = z__2.i / t;
00250             sine.r = z__1.r, sine.i = z__1.i;
00251             z__3.r = gamma->r / absest, z__3.i = gamma->i / absest;
00252             z__2.r = -z__3.r, z__2.i = -z__3.i;
00253             d__1 = t + 1.;
00254             z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;
00255             cosine.r = z__1.r, cosine.i = z__1.i;
00256             d_cnjg(&z__4, &sine);
00257             z__3.r = sine.r * z__4.r - sine.i * z__4.i, z__3.i = sine.r * 
00258                     z__4.i + sine.i * z__4.r;
00259             d_cnjg(&z__6, &cosine);
00260             z__5.r = cosine.r * z__6.r - cosine.i * z__6.i, z__5.i = cosine.r 
00261                     * z__6.i + cosine.i * z__6.r;
00262             z__2.r = z__3.r + z__5.r, z__2.i = z__3.i + z__5.i;
00263             z_sqrt(&z__1, &z__2);
00264             tmp = z__1.r;
00265             z__1.r = sine.r / tmp, z__1.i = sine.i / tmp;
00266             s->r = z__1.r, s->i = z__1.i;
00267             z__1.r = cosine.r / tmp, z__1.i = cosine.i / tmp;
00268             c__->r = z__1.r, c__->i = z__1.i;
00269             *sestpr = sqrt(t + 1.) * absest;
00270             return 0;
00271         }
00272 
00273     } else if (*job == 2) {
00274 
00275 /*        Estimating smallest singular value */
00276 
00277 /*        special cases */
00278 
00279         if (*sest == 0.) {
00280             *sestpr = 0.;
00281             if (max(absgam,absalp) == 0.) {
00282                 sine.r = 1., sine.i = 0.;
00283                 cosine.r = 0., cosine.i = 0.;
00284             } else {
00285                 d_cnjg(&z__2, gamma);
00286                 z__1.r = -z__2.r, z__1.i = -z__2.i;
00287                 sine.r = z__1.r, sine.i = z__1.i;
00288                 d_cnjg(&z__1, &alpha);
00289                 cosine.r = z__1.r, cosine.i = z__1.i;
00290             }
00291 /* Computing MAX */
00292             d__1 = z_abs(&sine), d__2 = z_abs(&cosine);
00293             s1 = max(d__1,d__2);
00294             z__1.r = sine.r / s1, z__1.i = sine.i / s1;
00295             s->r = z__1.r, s->i = z__1.i;
00296             z__1.r = cosine.r / s1, z__1.i = cosine.i / s1;
00297             c__->r = z__1.r, c__->i = z__1.i;
00298             d_cnjg(&z__4, s);
00299             z__3.r = s->r * z__4.r - s->i * z__4.i, z__3.i = s->r * z__4.i + 
00300                     s->i * z__4.r;
00301             d_cnjg(&z__6, c__);
00302             z__5.r = c__->r * z__6.r - c__->i * z__6.i, z__5.i = c__->r * 
00303                     z__6.i + c__->i * z__6.r;
00304             z__2.r = z__3.r + z__5.r, z__2.i = z__3.i + z__5.i;
00305             z_sqrt(&z__1, &z__2);
00306             tmp = z__1.r;
00307             z__1.r = s->r / tmp, z__1.i = s->i / tmp;
00308             s->r = z__1.r, s->i = z__1.i;
00309             z__1.r = c__->r / tmp, z__1.i = c__->i / tmp;
00310             c__->r = z__1.r, c__->i = z__1.i;
00311             return 0;
00312         } else if (absgam <= eps * absest) {
00313             s->r = 0., s->i = 0.;
00314             c__->r = 1., c__->i = 0.;
00315             *sestpr = absgam;
00316             return 0;
00317         } else if (absalp <= eps * absest) {
00318             s1 = absgam;
00319             s2 = absest;
00320             if (s1 <= s2) {
00321                 s->r = 0., s->i = 0.;
00322                 c__->r = 1., c__->i = 0.;
00323                 *sestpr = s1;
00324             } else {
00325                 s->r = 1., s->i = 0.;
00326                 c__->r = 0., c__->i = 0.;
00327                 *sestpr = s2;
00328             }
00329             return 0;
00330         } else if (absest <= eps * absalp || absest <= eps * absgam) {
00331             s1 = absgam;
00332             s2 = absalp;
00333             if (s1 <= s2) {
00334                 tmp = s1 / s2;
00335                 scl = sqrt(tmp * tmp + 1.);
00336                 *sestpr = absest * (tmp / scl);
00337                 d_cnjg(&z__4, gamma);
00338                 z__3.r = z__4.r / s2, z__3.i = z__4.i / s2;
00339                 z__2.r = -z__3.r, z__2.i = -z__3.i;
00340                 z__1.r = z__2.r / scl, z__1.i = z__2.i / scl;
00341                 s->r = z__1.r, s->i = z__1.i;
00342                 d_cnjg(&z__3, &alpha);
00343                 z__2.r = z__3.r / s2, z__2.i = z__3.i / s2;
00344                 z__1.r = z__2.r / scl, z__1.i = z__2.i / scl;
00345                 c__->r = z__1.r, c__->i = z__1.i;
00346             } else {
00347                 tmp = s2 / s1;
00348                 scl = sqrt(tmp * tmp + 1.);
00349                 *sestpr = absest / scl;
00350                 d_cnjg(&z__4, gamma);
00351                 z__3.r = z__4.r / s1, z__3.i = z__4.i / s1;
00352                 z__2.r = -z__3.r, z__2.i = -z__3.i;
00353                 z__1.r = z__2.r / scl, z__1.i = z__2.i / scl;
00354                 s->r = z__1.r, s->i = z__1.i;
00355                 d_cnjg(&z__3, &alpha);
00356                 z__2.r = z__3.r / s1, z__2.i = z__3.i / s1;
00357                 z__1.r = z__2.r / scl, z__1.i = z__2.i / scl;
00358                 c__->r = z__1.r, c__->i = z__1.i;
00359             }
00360             return 0;
00361         } else {
00362 
00363 /*           normal case */
00364 
00365             zeta1 = absalp / absest;
00366             zeta2 = absgam / absest;
00367 
00368 /* Computing MAX */
00369             d__1 = zeta1 * zeta1 + 1. + zeta1 * zeta2, d__2 = zeta1 * zeta2 + 
00370                     zeta2 * zeta2;
00371             norma = max(d__1,d__2);
00372 
00373 /*           See if root is closer to zero or to ONE */
00374 
00375             test = (zeta1 - zeta2) * 2. * (zeta1 + zeta2) + 1.;
00376             if (test >= 0.) {
00377 
00378 /*              root is close to zero, compute directly */
00379 
00380                 b = (zeta1 * zeta1 + zeta2 * zeta2 + 1.) * .5;
00381                 d__1 = zeta2 * zeta2;
00382                 c__->r = d__1, c__->i = 0.;
00383                 d__2 = b * b;
00384                 z__2.r = d__2 - c__->r, z__2.i = -c__->i;
00385                 d__1 = b + sqrt(z_abs(&z__2));
00386                 z__1.r = c__->r / d__1, z__1.i = c__->i / d__1;
00387                 t = z__1.r;
00388                 z__2.r = alpha.r / absest, z__2.i = alpha.i / absest;
00389                 d__1 = 1. - t;
00390                 z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;
00391                 sine.r = z__1.r, sine.i = z__1.i;
00392                 z__3.r = gamma->r / absest, z__3.i = gamma->i / absest;
00393                 z__2.r = -z__3.r, z__2.i = -z__3.i;
00394                 z__1.r = z__2.r / t, z__1.i = z__2.i / t;
00395                 cosine.r = z__1.r, cosine.i = z__1.i;
00396                 *sestpr = sqrt(t + eps * 4. * eps * norma) * absest;
00397             } else {
00398 
00399 /*              root is closer to ONE, shift by that amount */
00400 
00401                 b = (zeta2 * zeta2 + zeta1 * zeta1 - 1.) * .5;
00402                 d__1 = zeta1 * zeta1;
00403                 c__->r = d__1, c__->i = 0.;
00404                 if (b >= 0.) {
00405                     z__2.r = -c__->r, z__2.i = -c__->i;
00406                     d__1 = b * b;
00407                     z__5.r = d__1 + c__->r, z__5.i = c__->i;
00408                     z_sqrt(&z__4, &z__5);
00409                     z__3.r = b + z__4.r, z__3.i = z__4.i;
00410                     z_div(&z__1, &z__2, &z__3);
00411                     t = z__1.r;
00412                 } else {
00413                     d__1 = b * b;
00414                     z__3.r = d__1 + c__->r, z__3.i = c__->i;
00415                     z_sqrt(&z__2, &z__3);
00416                     z__1.r = b - z__2.r, z__1.i = -z__2.i;
00417                     t = z__1.r;
00418                 }
00419                 z__3.r = alpha.r / absest, z__3.i = alpha.i / absest;
00420                 z__2.r = -z__3.r, z__2.i = -z__3.i;
00421                 z__1.r = z__2.r / t, z__1.i = z__2.i / t;
00422                 sine.r = z__1.r, sine.i = z__1.i;
00423                 z__3.r = gamma->r / absest, z__3.i = gamma->i / absest;
00424                 z__2.r = -z__3.r, z__2.i = -z__3.i;
00425                 d__1 = t + 1.;
00426                 z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;
00427                 cosine.r = z__1.r, cosine.i = z__1.i;
00428                 *sestpr = sqrt(t + 1. + eps * 4. * eps * norma) * absest;
00429             }
00430             d_cnjg(&z__4, &sine);
00431             z__3.r = sine.r * z__4.r - sine.i * z__4.i, z__3.i = sine.r * 
00432                     z__4.i + sine.i * z__4.r;
00433             d_cnjg(&z__6, &cosine);
00434             z__5.r = cosine.r * z__6.r - cosine.i * z__6.i, z__5.i = cosine.r 
00435                     * z__6.i + cosine.i * z__6.r;
00436             z__2.r = z__3.r + z__5.r, z__2.i = z__3.i + z__5.i;
00437             z_sqrt(&z__1, &z__2);
00438             tmp = z__1.r;
00439             z__1.r = sine.r / tmp, z__1.i = sine.i / tmp;
00440             s->r = z__1.r, s->i = z__1.i;
00441             z__1.r = cosine.r / tmp, z__1.i = cosine.i / tmp;
00442             c__->r = z__1.r, c__->i = z__1.i;
00443             return 0;
00444 
00445         }
00446     }
00447     return 0;
00448 
00449 /*     End of ZLAIC1 */
00450 
00451 } /* zlaic1_ */

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