kusano 7d535a
kusano 7d535a
/*! @file dlacon.c
kusano 7d535a
 * \brief Estimates the 1-norm
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 * -- SuperLU routine (version 2.0) --
kusano 7d535a
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
kusano 7d535a
 * and Lawrence Berkeley National Lab.
kusano 7d535a
 * November 15, 1997
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
#include <math.h></math.h>
kusano 7d535a
#include "slu_Cnames.h"
kusano 7d535a
kusano 7d535a
/*! \brief
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 *   Purpose   
kusano 7d535a
 *   =======   
kusano 7d535a
 *
kusano 7d535a
 *   DLACON estimates the 1-norm of a square matrix A.   
kusano 7d535a
 *   Reverse communication is used for evaluating matrix-vector products. 
kusano 7d535a
 * 
kusano 7d535a
 *
kusano 7d535a
 *   Arguments   
kusano 7d535a
 *   =========   
kusano 7d535a
 *
kusano 7d535a
 *   N      (input) INT
kusano 7d535a
 *          The order of the matrix.  N >= 1.   
kusano 7d535a
 *
kusano 7d535a
 *   V      (workspace) DOUBLE PRECISION array, dimension (N)   
kusano 7d535a
 *          On the final return, V = A*W,  where  EST = norm(V)/norm(W)   
kusano 7d535a
 *          (W is not returned).   
kusano 7d535a
 *
kusano 7d535a
 *   X      (input/output) DOUBLE PRECISION array, dimension (N)   
kusano 7d535a
 *          On an intermediate return, X should be overwritten by   
kusano 7d535a
 *                A * X,   if KASE=1,   
kusano 7d535a
 *                A' * X,  if KASE=2,
kusano 7d535a
 *         and DLACON must be re-called with all the other parameters   
kusano 7d535a
 *          unchanged.   
kusano 7d535a
 *
kusano 7d535a
 *   ISGN   (workspace) INT array, dimension (N)
kusano 7d535a
 *
kusano 7d535a
 *   EST    (output) DOUBLE PRECISION   
kusano 7d535a
 *          An estimate (a lower bound) for norm(A).   
kusano 7d535a
 *
kusano 7d535a
 *   KASE   (input/output) INT
kusano 7d535a
 *          On the initial call to DLACON, KASE should be 0.   
kusano 7d535a
 *          On an intermediate return, KASE will be 1 or 2, indicating   
kusano 7d535a
 *          whether X should be overwritten by A * X  or A' * X.   
kusano 7d535a
 *          On the final return from DLACON, KASE will again be 0.   
kusano 7d535a
 *
kusano 7d535a
 *   Further Details   
kusano 7d535a
 *   ======= =======   
kusano 7d535a
 *
kusano 7d535a
 *   Contributed by Nick Higham, University of Manchester.   
kusano 7d535a
 *   Originally named CONEST, dated March 16, 1988.   
kusano 7d535a
 *
kusano 7d535a
 *   Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of 
kusano 7d535a
 *   a real or complex matrix, with applications to condition estimation", 
kusano 7d535a
 *   ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.   
kusano 7d535a
 *   ===================================================================== 
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
int
kusano 7d535a
dlacon_(int *n, double *v, double *x, int *isgn, double *est, int *kase)
kusano 7d535a
kusano 7d535a
{
kusano 7d535a
kusano 7d535a
kusano 7d535a
    /* Table of constant values */
kusano 7d535a
    int c__1 = 1;
kusano 7d535a
    double      zero = 0.0;
kusano 7d535a
    double      one = 1.0;
kusano 7d535a
    
kusano 7d535a
    /* Local variables */
kusano 7d535a
    static int iter;
kusano 7d535a
    static int jump, jlast;
kusano 7d535a
    static double altsgn, estold;
kusano 7d535a
    static int i, j;
kusano 7d535a
    double temp;
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
    extern int ISAMAX(int *, double *, int *);
kusano 7d535a
    extern double SASUM(int *, double *, int *);
kusano 7d535a
    extern int SCOPY(int *, double *, int *, double *, int *);
kusano 7d535a
#else
kusano 7d535a
    extern int idamax_(int *, double *, int *);
kusano 7d535a
    extern double dasum_(int *, double *, int *);
kusano 7d535a
    extern int dcopy_(int *, double *, int *, double *, int *);
kusano 7d535a
#endif
kusano 7d535a
#define d_sign(a, b) (b >= 0 ? fabs(a) : -fabs(a))    /* Copy sign */
kusano 7d535a
#define i_dnnt(a) \
kusano 7d535a
	( a>=0 ? floor(a+.5) : -floor(.5-a) ) /* Round to nearest integer */
kusano 7d535a
kusano 7d535a
    if ( *kase == 0 ) {
kusano 7d535a
	for (i = 0; i < *n; ++i) {
kusano 7d535a
	    x[i] = 1. / (double) (*n);
kusano 7d535a
	}
kusano 7d535a
	*kase = 1;
kusano 7d535a
	jump = 1;
kusano 7d535a
	return 0;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    switch (jump) {
kusano 7d535a
	case 1:  goto L20;
kusano 7d535a
	case 2:  goto L40;
kusano 7d535a
	case 3:  goto L70;
kusano 7d535a
	case 4:  goto L110;
kusano 7d535a
	case 5:  goto L140;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /*     ................ ENTRY   (JUMP = 1)   
kusano 7d535a
	   FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X. */
kusano 7d535a
  L20:
kusano 7d535a
    if (*n == 1) {
kusano 7d535a
	v[0] = x[0];
kusano 7d535a
	*est = fabs(v[0]);
kusano 7d535a
	/*        ... QUIT */
kusano 7d535a
	goto L150;
kusano 7d535a
    }
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
    *est = SASUM(n, x, &c__1);
kusano 7d535a
#else
kusano 7d535a
    *est = dasum_(n, x, &c__1);
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
    for (i = 0; i < *n; ++i) {
kusano 7d535a
	x[i] = d_sign(one, x[i]);
kusano 7d535a
	isgn[i] = i_dnnt(x[i]);
kusano 7d535a
    }
kusano 7d535a
    *kase = 2;
kusano 7d535a
    jump = 2;
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
    /*     ................ ENTRY   (JUMP = 2)   
kusano 7d535a
	   FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. */
kusano 7d535a
L40:
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
    j = ISAMAX(n, &x[0], &c__1);
kusano 7d535a
#else
kusano 7d535a
    j = idamax_(n, &x[0], &c__1);
kusano 7d535a
#endif
kusano 7d535a
    --j;
kusano 7d535a
    iter = 2;
kusano 7d535a
kusano 7d535a
    /*     MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
kusano 7d535a
L50:
kusano 7d535a
    for (i = 0; i < *n; ++i) x[i] = zero;
kusano 7d535a
    x[j] = one;
kusano 7d535a
    *kase = 1;
kusano 7d535a
    jump = 3;
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
    /*     ................ ENTRY   (JUMP = 3)   
kusano 7d535a
	   X HAS BEEN OVERWRITTEN BY A*X. */
kusano 7d535a
L70:
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
    SCOPY(n, x, &c__1, v, &c__1);
kusano 7d535a
#else
kusano 7d535a
    dcopy_(n, x, &c__1, v, &c__1);
kusano 7d535a
#endif
kusano 7d535a
    estold = *est;
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
    *est = SASUM(n, v, &c__1);
kusano 7d535a
#else
kusano 7d535a
    *est = dasum_(n, v, &c__1);
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
    for (i = 0; i < *n; ++i)
kusano 7d535a
	if (i_dnnt(d_sign(one, x[i])) != isgn[i])
kusano 7d535a
	    goto L90;
kusano 7d535a
kusano 7d535a
    /*     REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. */
kusano 7d535a
    goto L120;
kusano 7d535a
kusano 7d535a
L90:
kusano 7d535a
    /*     TEST FOR CYCLING. */
kusano 7d535a
    if (*est <= estold) goto L120;
kusano 7d535a
kusano 7d535a
    for (i = 0; i < *n; ++i) {
kusano 7d535a
	x[i] = d_sign(one, x[i]);
kusano 7d535a
	isgn[i] = i_dnnt(x[i]);
kusano 7d535a
    }
kusano 7d535a
    *kase = 2;
kusano 7d535a
    jump = 4;
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
    /*     ................ ENTRY   (JUMP = 4)   
kusano 7d535a
	   X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
kusano 7d535a
L110:
kusano 7d535a
    jlast = j;
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
    j = ISAMAX(n, &x[0], &c__1);
kusano 7d535a
#else
kusano 7d535a
    j = idamax_(n, &x[0], &c__1);
kusano 7d535a
#endif
kusano 7d535a
    --j;
kusano 7d535a
    if (x[jlast] != fabs(x[j]) && iter < 5) {
kusano 7d535a
	++iter;
kusano 7d535a
	goto L50;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /*     ITERATION COMPLETE.  FINAL STAGE. */
kusano 7d535a
L120:
kusano 7d535a
    altsgn = 1.;
kusano 7d535a
    for (i = 1; i <= *n; ++i) {
kusano 7d535a
	x[i-1] = altsgn * ((double)(i - 1) / (double)(*n - 1) + 1.);
kusano 7d535a
	altsgn = -altsgn;
kusano 7d535a
    }
kusano 7d535a
    *kase = 1;
kusano 7d535a
    jump = 5;
kusano 7d535a
    return 0;
kusano 7d535a
    
kusano 7d535a
    /*     ................ ENTRY   (JUMP = 5)   
kusano 7d535a
	   X HAS BEEN OVERWRITTEN BY A*X. */
kusano 7d535a
L140:
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
    temp = SASUM(n, x, &c__1) / (double)(*n * 3) * 2.;
kusano 7d535a
#else
kusano 7d535a
    temp = dasum_(n, x, &c__1) / (double)(*n * 3) * 2.;
kusano 7d535a
#endif
kusano 7d535a
    if (temp > *est) {
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
	SCOPY(n, &x[0], &c__1, &v[0], &c__1);
kusano 7d535a
#else
kusano 7d535a
	dcopy_(n, &x[0], &c__1, &v[0], &c__1);
kusano 7d535a
#endif
kusano 7d535a
	*est = temp;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
L150:
kusano 7d535a
    *kase = 0;
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
} /* dlacon_ */