[] add metering mode to CLI

1 files changed, 440 insertions, 0 deletions

diff --git a/contrib/libs/clapack/slar1v.c b/contrib/libs/clapack/slar1v.c
new file mode 100644
index 0000000000..59de90199a
--- /dev/null
+++ b/contrib/libs/clapack/slar1v.c
@@ -0,0 +1,440 @@
+/* slar1v.f -- translated by f2c (version 20061008).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "f2c.h"
+#include "blaswrap.h"
+
+/* Subroutine */ int slar1v_(integer *n, integer *b1, integer *bn, real *
+	lambda, real *d__, real *l, real *ld, real *lld, real *pivmin, real *
+	gaptol, real *z__, logical *wantnc, integer *negcnt, real *ztz, real *
+	mingma, integer *r__, integer *isuppz, real *nrminv, real *resid, 
+	real *rqcorr, real *work)
+{
+    /* System generated locals */
+    integer i__1;
+    real r__1, r__2, r__3;
+
+    /* Builtin functions */
+    double sqrt(doublereal);
+
+    /* Local variables */
+    integer i__;
+    real s;
+    integer r1, r2;
+    real eps, tmp;
+    integer neg1, neg2, indp, inds;
+    real dplus;
+    extern doublereal slamch_(char *);
+    integer indlpl, indumn;
+    extern logical sisnan_(real *);
+    real dminus;
+    logical sawnan1, sawnan2;
+
+
+/*  -- LAPACK auxiliary routine (version 3.2) -- */
+/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/*     November 2006 */
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  SLAR1V computes the (scaled) r-th column of the inverse of */
+/*  the sumbmatrix in rows B1 through BN of the tridiagonal matrix */
+/*  L D L^T - sigma I. When sigma is close to an eigenvalue, the */
+/*  computed vector is an accurate eigenvector. Usually, r corresponds */
+/*  to the index where the eigenvector is largest in magnitude. */
+/*  The following steps accomplish this computation : */
+/*  (a) Stationary qd transform,  L D L^T - sigma I = L(+) D(+) L(+)^T, */
+/*  (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T, */
+/*  (c) Computation of the diagonal elements of the inverse of */
+/*      L D L^T - sigma I by combining the above transforms, and choosing */
+/*      r as the index where the diagonal of the inverse is (one of the) */
+/*      largest in magnitude. */
+/*  (d) Computation of the (scaled) r-th column of the inverse using the */
+/*      twisted factorization obtained by combining the top part of the */
+/*      the stationary and the bottom part of the progressive transform. */
+
+/*  Arguments */
+/*  ========= */
+
+/*  N        (input) INTEGER */
+/*           The order of the matrix L D L^T. */
+
+/*  B1       (input) INTEGER */
+/*           First index of the submatrix of L D L^T. */
+
+/*  BN       (input) INTEGER */
+/*           Last index of the submatrix of L D L^T. */
+
+/*  LAMBDA    (input) REAL */
+/*           The shift. In order to compute an accurate eigenvector, */
+/*           LAMBDA should be a good approximation to an eigenvalue */
+/*           of L D L^T. */
+
+/*  L        (input) REAL             array, dimension (N-1) */
+/*           The (n-1) subdiagonal elements of the unit bidiagonal matrix */
+/*           L, in elements 1 to N-1. */
+
+/*  D        (input) REAL             array, dimension (N) */
+/*           The n diagonal elements of the diagonal matrix D. */
+
+/*  LD       (input) REAL             array, dimension (N-1) */
+/*           The n-1 elements L(i)*D(i). */
+
+/*  LLD      (input) REAL             array, dimension (N-1) */
+/*           The n-1 elements L(i)*L(i)*D(i). */
+
+/*  PIVMIN   (input) REAL */
+/*           The minimum pivot in the Sturm sequence. */
+
+/*  GAPTOL   (input) REAL */
+/*           Tolerance that indicates when eigenvector entries are negligible */
+/*           w.r.t. their contribution to the residual. */
+
+/*  Z        (input/output) REAL             array, dimension (N) */
+/*           On input, all entries of Z must be set to 0. */
+/*           On output, Z contains the (scaled) r-th column of the */
+/*           inverse. The scaling is such that Z(R) equals 1. */
+
+/*  WANTNC   (input) LOGICAL */
+/*           Specifies whether NEGCNT has to be computed. */
+
+/*  NEGCNT   (output) INTEGER */
+/*           If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin */
+/*           in the  matrix factorization L D L^T, and NEGCNT = -1 otherwise. */
+
+/*  ZTZ      (output) REAL */
+/*           The square of the 2-norm of Z. */
+
+/*  MINGMA   (output) REAL */
+/*           The reciprocal of the largest (in magnitude) diagonal */
+/*           element of the inverse of L D L^T - sigma I. */
+
+/*  R        (input/output) INTEGER */
+/*           The twist index for the twisted factorization used to */
+/*           compute Z. */
+/*           On input, 0 <= R <= N. If R is input as 0, R is set to */
+/*           the index where (L D L^T - sigma I)^{-1} is largest */
+/*           in magnitude. If 1 <= R <= N, R is unchanged. */
+/*           On output, R contains the twist index used to compute Z. */
+/*           Ideally, R designates the position of the maximum entry in the */
+/*           eigenvector. */
+
+/*  ISUPPZ   (output) INTEGER array, dimension (2) */
+/*           The support of the vector in Z, i.e., the vector Z is */
+/*           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ). */
+
+/*  NRMINV   (output) REAL */
+/*           NRMINV = 1/SQRT( ZTZ ) */
+
+/*  RESID    (output) REAL */
+/*           The residual of the FP vector. */
+/*           RESID = ABS( MINGMA )/SQRT( ZTZ ) */
+
+/*  RQCORR   (output) REAL */
+/*           The Rayleigh Quotient correction to LAMBDA. */
+/*           RQCORR = MINGMA*TMP */
+
+/*  WORK     (workspace) REAL             array, dimension (4*N) */
+
+/*  Further Details */
+/*  =============== */
+
+/*  Based on contributions by */
+/*     Beresford Parlett, University of California, Berkeley, USA */
+/*     Jim Demmel, University of California, Berkeley, USA */
+/*     Inderjit Dhillon, University of Texas, Austin, USA */
+/*     Osni Marques, LBNL/NERSC, USA */
+/*     Christof Voemel, University of California, Berkeley, USA */
+
+/*  ===================================================================== */
+
+/*     .. Parameters .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+    /* Parameter adjustments */
+    --work;
+    --isuppz;
+    --z__;
+    --lld;
+    --ld;
+    --l;
+    --d__;
+
+    /* Function Body */
+    eps = slamch_("Precision");
+    if (*r__ == 0) {
+	r1 = *b1;
+	r2 = *bn;
+    } else {
+	r1 = *r__;
+	r2 = *r__;
+    }
+/*     Storage for LPLUS */
+    indlpl = 0;
+/*     Storage for UMINUS */
+    indumn = *n;
+    inds = (*n << 1) + 1;
+    indp = *n * 3 + 1;
+    if (*b1 == 1) {
+	work[inds] = 0.f;
+    } else {
+	work[inds + *b1 - 1] = lld[*b1 - 1];
+    }
+
+/*     Compute the stationary transform (using the differential form) */
+/*     until the index R2. */
+
+    sawnan1 = FALSE_;
+    neg1 = 0;
+    s = work[inds + *b1 - 1] - *lambda;
+    i__1 = r1 - 1;
+    for (i__ = *b1; i__ <= i__1; ++i__) {
+	dplus = d__[i__] + s;
+	work[indlpl + i__] = ld[i__] / dplus;
+	if (dplus < 0.f) {
+	    ++neg1;
+	}
+	work[inds + i__] = s * work[indlpl + i__] * l[i__];
+	s = work[inds + i__] - *lambda;
+/* L50: */
+    }
+    sawnan1 = sisnan_(&s);
+    if (sawnan1) {
+	goto L60;
+    }
+    i__1 = r2 - 1;
+    for (i__ = r1; i__ <= i__1; ++i__) {
+	dplus = d__[i__] + s;
+	work[indlpl + i__] = ld[i__] / dplus;
+	work[inds + i__] = s * work[indlpl + i__] * l[i__];
+	s = work[inds + i__] - *lambda;
+/* L51: */
+    }
+    sawnan1 = sisnan_(&s);
+
+L60:
+    if (sawnan1) {
+/*        Runs a slower version of the above loop if a NaN is detected */
+	neg1 = 0;
+	s = work[inds + *b1 - 1] - *lambda;
+	i__1 = r1 - 1;
+	for (i__ = *b1; i__ <= i__1; ++i__) {
+	    dplus = d__[i__] + s;
+	    if (dabs(dplus) < *pivmin) {
+		dplus = -(*pivmin);
+	    }
+	    work[indlpl + i__] = ld[i__] / dplus;
+	    if (dplus < 0.f) {
+		++neg1;
+	    }
+	    work[inds + i__] = s * work[indlpl + i__] * l[i__];
+	    if (work[indlpl + i__] == 0.f) {
+		work[inds + i__] = lld[i__];
+	    }
+	    s = work[inds + i__] - *lambda;
+/* L70: */
+	}
+	i__1 = r2 - 1;
+	for (i__ = r1; i__ <= i__1; ++i__) {
+	    dplus = d__[i__] + s;
+	    if (dabs(dplus) < *pivmin) {
+		dplus = -(*pivmin);
+	    }
+	    work[indlpl + i__] = ld[i__] / dplus;
+	    work[inds + i__] = s * work[indlpl + i__] * l[i__];
+	    if (work[indlpl + i__] == 0.f) {
+		work[inds + i__] = lld[i__];
+	    }
+	    s = work[inds + i__] - *lambda;
+/* L71: */
+	}
+    }
+
+/*     Compute the progressive transform (using the differential form) */
+/*     until the index R1 */
+
+    sawnan2 = FALSE_;
+    neg2 = 0;
+    work[indp + *bn - 1] = d__[*bn] - *lambda;
+    i__1 = r1;
+    for (i__ = *bn - 1; i__ >= i__1; --i__) {
+	dminus = lld[i__] + work[indp + i__];
+	tmp = d__[i__] / dminus;
+	if (dminus < 0.f) {
+	    ++neg2;
+	}
+	work[indumn + i__] = l[i__] * tmp;
+	work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
+/* L80: */
+    }
+    tmp = work[indp + r1 - 1];
+    sawnan2 = sisnan_(&tmp);
+    if (sawnan2) {
+/*        Runs a slower version of the above loop if a NaN is detected */
+	neg2 = 0;
+	i__1 = r1;
+	for (i__ = *bn - 1; i__ >= i__1; --i__) {
+	    dminus = lld[i__] + work[indp + i__];
+	    if (dabs(dminus) < *pivmin) {
+		dminus = -(*pivmin);
+	    }
+	    tmp = d__[i__] / dminus;
+	    if (dminus < 0.f) {
+		++neg2;
+	    }
+	    work[indumn + i__] = l[i__] * tmp;
+	    work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
+	    if (tmp == 0.f) {
+		work[indp + i__ - 1] = d__[i__] - *lambda;
+	    }
+/* L100: */
+	}
+    }
+
+/*     Find the index (from R1 to R2) of the largest (in magnitude) */
+/*     diagonal element of the inverse */
+
+    *mingma = work[inds + r1 - 1] + work[indp + r1 - 1];
+    if (*mingma < 0.f) {
+	++neg1;
+    }
+    if (*wantnc) {
+	*negcnt = neg1 + neg2;
+    } else {
+	*negcnt = -1;
+    }
+    if (dabs(*mingma) == 0.f) {
+	*mingma = eps * work[inds + r1 - 1];
+    }
+    *r__ = r1;
+    i__1 = r2 - 1;
+    for (i__ = r1; i__ <= i__1; ++i__) {
+	tmp = work[inds + i__] + work[indp + i__];
+	if (tmp == 0.f) {
+	    tmp = eps * work[inds + i__];
+	}
+	if (dabs(tmp) <= dabs(*mingma)) {
+	    *mingma = tmp;
+	    *r__ = i__ + 1;
+	}
+/* L110: */
+    }
+
+/*     Compute the FP vector: solve N^T v = e_r */
+
+    isuppz[1] = *b1;
+    isuppz[2] = *bn;
+    z__[*r__] = 1.f;
+    *ztz = 1.f;
+
+/*     Compute the FP vector upwards from R */
+
+    if (! sawnan1 && ! sawnan2) {
+	i__1 = *b1;
+	for (i__ = *r__ - 1; i__ >= i__1; --i__) {
+	    z__[i__] = -(work[indlpl + i__] * z__[i__ + 1]);
+	    if (((r__1 = z__[i__], dabs(r__1)) + (r__2 = z__[i__ + 1], dabs(
+		    r__2))) * (r__3 = ld[i__], dabs(r__3)) < *gaptol) {
+		z__[i__] = 0.f;
+		isuppz[1] = i__ + 1;
+		goto L220;
+	    }
+	    *ztz += z__[i__] * z__[i__];
+/* L210: */
+	}
+L220:
+	;
+    } else {
+/*        Run slower loop if NaN occurred. */
+	i__1 = *b1;
+	for (i__ = *r__ - 1; i__ >= i__1; --i__) {
+	    if (z__[i__ + 1] == 0.f) {
+		z__[i__] = -(ld[i__ + 1] / ld[i__]) * z__[i__ + 2];
+	    } else {
+		z__[i__] = -(work[indlpl + i__] * z__[i__ + 1]);
+	    }
+	    if (((r__1 = z__[i__], dabs(r__1)) + (r__2 = z__[i__ + 1], dabs(
+		    r__2))) * (r__3 = ld[i__], dabs(r__3)) < *gaptol) {
+		z__[i__] = 0.f;
+		isuppz[1] = i__ + 1;
+		goto L240;
+	    }
+	    *ztz += z__[i__] * z__[i__];
+/* L230: */
+	}
+L240:
+	;
+    }
+/*     Compute the FP vector downwards from R in blocks of size BLKSIZ */
+    if (! sawnan1 && ! sawnan2) {
+	i__1 = *bn - 1;
+	for (i__ = *r__; i__ <= i__1; ++i__) {
+	    z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
+	    if (((r__1 = z__[i__], dabs(r__1)) + (r__2 = z__[i__ + 1], dabs(
+		    r__2))) * (r__3 = ld[i__], dabs(r__3)) < *gaptol) {
+		z__[i__ + 1] = 0.f;
+		isuppz[2] = i__;
+		goto L260;
+	    }
+	    *ztz += z__[i__ + 1] * z__[i__ + 1];
+/* L250: */
+	}
+L260:
+	;
+    } else {
+/*        Run slower loop if NaN occurred. */
+	i__1 = *bn - 1;
+	for (i__ = *r__; i__ <= i__1; ++i__) {
+	    if (z__[i__] == 0.f) {
+		z__[i__ + 1] = -(ld[i__ - 1] / ld[i__]) * z__[i__ - 1];
+	    } else {
+		z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
+	    }
+	    if (((r__1 = z__[i__], dabs(r__1)) + (r__2 = z__[i__ + 1], dabs(
+		    r__2))) * (r__3 = ld[i__], dabs(r__3)) < *gaptol) {
+		z__[i__ + 1] = 0.f;
+		isuppz[2] = i__;
+		goto L280;
+	    }
+	    *ztz += z__[i__ + 1] * z__[i__ + 1];
+/* L270: */
+	}
+L280:
+	;
+    }
+
+/*     Compute quantities for convergence test */
+
+    tmp = 1.f / *ztz;
+    *nrminv = sqrt(tmp);
+    *resid = dabs(*mingma) * *nrminv;
+    *rqcorr = *mingma * tmp;
+
+
+    return 0;
+
+/*     End of SLAR1V */
+
+} /* slar1v_ */