swiftnav: claic1.c Source File

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