clacon.c
Go to the documentation of this file.
00001 /* clacon.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 clacon_(integer *n, complex *v, complex *x, real *est, 
00021         integer *kase)
00022 {
00023     /* System generated locals */
00024     integer i__1, i__2, i__3;
00025     real r__1, r__2;
00026     complex q__1;
00027 
00028     /* Builtin functions */
00029     double c_abs(complex *), r_imag(complex *);
00030 
00031     /* Local variables */
00032     static integer i__, j, iter;
00033     static real temp;
00034     static integer jump;
00035     static real absxi;
00036     static integer jlast;
00037     extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 
00038             complex *, integer *);
00039     extern integer icmax1_(integer *, complex *, integer *);
00040     extern doublereal scsum1_(integer *, complex *, integer *), slamch_(char *
00041 );
00042     static real safmin, altsgn, estold;
00043 
00044 
00045 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00046 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00047 /*     November 2006 */
00048 
00049 /*     .. Scalar Arguments .. */
00050 /*     .. */
00051 /*     .. Array Arguments .. */
00052 /*     .. */
00053 
00054 /*  Purpose */
00055 /*  ======= */
00056 
00057 /*  CLACON estimates the 1-norm of a square, complex matrix A. */
00058 /*  Reverse communication is used for evaluating matrix-vector products. */
00059 
00060 /*  Arguments */
00061 /*  ========= */
00062 
00063 /*  N      (input) INTEGER */
00064 /*         The order of the matrix.  N >= 1. */
00065 
00066 /*  V      (workspace) COMPLEX array, dimension (N) */
00067 /*         On the final return, V = A*W,  where  EST = norm(V)/norm(W) */
00068 /*         (W is not returned). */
00069 
00070 /*  X      (input/output) COMPLEX array, dimension (N) */
00071 /*         On an intermediate return, X should be overwritten by */
00072 /*               A * X,   if KASE=1, */
00073 /*               A' * X,  if KASE=2, */
00074 /*         where A' is the conjugate transpose of A, and CLACON must be */
00075 /*         re-called with all the other parameters unchanged. */
00076 
00077 /*  EST    (input/output) REAL */
00078 /*         On entry with KASE = 1 or 2 and JUMP = 3, EST should be */
00079 /*         unchanged from the previous call to CLACON. */
00080 /*         On exit, EST is an estimate (a lower bound) for norm(A). */
00081 
00082 /*  KASE   (input/output) INTEGER */
00083 /*         On the initial call to CLACON, KASE should be 0. */
00084 /*         On an intermediate return, KASE will be 1 or 2, indicating */
00085 /*         whether X should be overwritten by A * X  or A' * X. */
00086 /*         On the final return from CLACON, KASE will again be 0. */
00087 
00088 /*  Further Details */
00089 /*  ======= ======= */
00090 
00091 /*  Contributed by Nick Higham, University of Manchester. */
00092 /*  Originally named CONEST, dated March 16, 1988. */
00093 
00094 /*  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of */
00095 /*  a real or complex matrix, with applications to condition estimation", */
00096 /*  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. */
00097 
00098 /*  Last modified:  April, 1999 */
00099 
00100 /*  ===================================================================== */
00101 
00102 /*     .. Parameters .. */
00103 /*     .. */
00104 /*     .. Local Scalars .. */
00105 /*     .. */
00106 /*     .. External Functions .. */
00107 /*     .. */
00108 /*     .. External Subroutines .. */
00109 /*     .. */
00110 /*     .. Intrinsic Functions .. */
00111 /*     .. */
00112 /*     .. Save statement .. */
00113 /*     .. */
00114 /*     .. Executable Statements .. */
00115 
00116     /* Parameter adjustments */
00117     --x;
00118     --v;
00119 
00120     /* Function Body */
00121     safmin = slamch_("Safe minimum");
00122     if (*kase == 0) {
00123         i__1 = *n;
00124         for (i__ = 1; i__ <= i__1; ++i__) {
00125             i__2 = i__;
00126             r__1 = 1.f / (real) (*n);
00127             q__1.r = r__1, q__1.i = 0.f;
00128             x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00129 /* L10: */
00130         }
00131         *kase = 1;
00132         jump = 1;
00133         return 0;
00134     }
00135 
00136     switch (jump) {
00137         case 1:  goto L20;
00138         case 2:  goto L40;
00139         case 3:  goto L70;
00140         case 4:  goto L90;
00141         case 5:  goto L120;
00142     }
00143 
00144 /*     ................ ENTRY   (JUMP = 1) */
00145 /*     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X. */
00146 
00147 L20:
00148     if (*n == 1) {
00149         v[1].r = x[1].r, v[1].i = x[1].i;
00150         *est = c_abs(&v[1]);
00151 /*        ... QUIT */
00152         goto L130;
00153     }
00154     *est = scsum1_(n, &x[1], &c__1);
00155 
00156     i__1 = *n;
00157     for (i__ = 1; i__ <= i__1; ++i__) {
00158         absxi = c_abs(&x[i__]);
00159         if (absxi > safmin) {
00160             i__2 = i__;
00161             i__3 = i__;
00162             r__1 = x[i__3].r / absxi;
00163             r__2 = r_imag(&x[i__]) / absxi;
00164             q__1.r = r__1, q__1.i = r__2;
00165             x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00166         } else {
00167             i__2 = i__;
00168             x[i__2].r = 1.f, x[i__2].i = 0.f;
00169         }
00170 /* L30: */
00171     }
00172     *kase = 2;
00173     jump = 2;
00174     return 0;
00175 
00176 /*     ................ ENTRY   (JUMP = 2) */
00177 /*     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY CTRANS(A)*X. */
00178 
00179 L40:
00180     j = icmax1_(n, &x[1], &c__1);
00181     iter = 2;
00182 
00183 /*     MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
00184 
00185 L50:
00186     i__1 = *n;
00187     for (i__ = 1; i__ <= i__1; ++i__) {
00188         i__2 = i__;
00189         x[i__2].r = 0.f, x[i__2].i = 0.f;
00190 /* L60: */
00191     }
00192     i__1 = j;
00193     x[i__1].r = 1.f, x[i__1].i = 0.f;
00194     *kase = 1;
00195     jump = 3;
00196     return 0;
00197 
00198 /*     ................ ENTRY   (JUMP = 3) */
00199 /*     X HAS BEEN OVERWRITTEN BY A*X. */
00200 
00201 L70:
00202     ccopy_(n, &x[1], &c__1, &v[1], &c__1);
00203     estold = *est;
00204     *est = scsum1_(n, &v[1], &c__1);
00205 
00206 /*     TEST FOR CYCLING. */
00207     if (*est <= estold) {
00208         goto L100;
00209     }
00210 
00211     i__1 = *n;
00212     for (i__ = 1; i__ <= i__1; ++i__) {
00213         absxi = c_abs(&x[i__]);
00214         if (absxi > safmin) {
00215             i__2 = i__;
00216             i__3 = i__;
00217             r__1 = x[i__3].r / absxi;
00218             r__2 = r_imag(&x[i__]) / absxi;
00219             q__1.r = r__1, q__1.i = r__2;
00220             x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00221         } else {
00222             i__2 = i__;
00223             x[i__2].r = 1.f, x[i__2].i = 0.f;
00224         }
00225 /* L80: */
00226     }
00227     *kase = 2;
00228     jump = 4;
00229     return 0;
00230 
00231 /*     ................ ENTRY   (JUMP = 4) */
00232 /*     X HAS BEEN OVERWRITTEN BY CTRANS(A)*X. */
00233 
00234 L90:
00235     jlast = j;
00236     j = icmax1_(n, &x[1], &c__1);
00237     if (c_abs(&x[jlast]) != c_abs(&x[j]) && iter < 5) {
00238         ++iter;
00239         goto L50;
00240     }
00241 
00242 /*     ITERATION COMPLETE.  FINAL STAGE. */
00243 
00244 L100:
00245     altsgn = 1.f;
00246     i__1 = *n;
00247     for (i__ = 1; i__ <= i__1; ++i__) {
00248         i__2 = i__;
00249         r__1 = altsgn * ((real) (i__ - 1) / (real) (*n - 1) + 1.f);
00250         q__1.r = r__1, q__1.i = 0.f;
00251         x[i__2].r = q__1.r, x[i__2].i = q__1.i;
00252         altsgn = -altsgn;
00253 /* L110: */
00254     }
00255     *kase = 1;
00256     jump = 5;
00257     return 0;
00258 
00259 /*     ................ ENTRY   (JUMP = 5) */
00260 /*     X HAS BEEN OVERWRITTEN BY A*X. */
00261 
00262 L120:
00263     temp = scsum1_(n, &x[1], &c__1) / (real) (*n * 3) * 2.f;
00264     if (temp > *est) {
00265         ccopy_(n, &x[1], &c__1, &v[1], &c__1);
00266         *est = temp;
00267     }
00268 
00269 L130:
00270     *kase = 0;
00271     return 0;
00272 
00273 /*     End of CLACON */
00274 
00275 } /* clacon_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:29