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


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