A C++ String Class

#include "stdafx.h"

/*
 * -- SuperLU routine (version 2.0) --
 * Lawrence Berkeley National Lab, Univ. of California Berkeley,
 * and Xerox Palo Alto Research Center.
 * September 10, 2007
 *
 */
#include <math.h>
#include "hnum_slu_mt_Cnames.h"
#include "hnum_slu_mt_util.h"
#include "hnum_slu_dcomplex.h"

#include "hnum_slucmn.h"

namespace harlinn
{
    namespace numerics
    {
        namespace SuperLU
        {
            int zlacon_(int *n, doublecomplex *v, doublecomplex *x, double *est, int *kase)

            {
            /*
                Purpose   
                =======   

                ZLACON estimates the 1-norm of a square matrix A.   
                Reverse communication is used for evaluating matrix-vector products. 
  

                Arguments   
                =========   

                N      (input) INT
                        The order of the matrix.  N >= 1.   

                V      (workspace) DOUBLE COMPLEX PRECISION array, dimension (N)   
                        On the final return, V = A*W,  where  EST = norm(V)/norm(W)   
                        (W is not returned).   

                X      (input/output) DOUBLE COMPLEX PRECISION array, dimension (N)   
                        On an intermediate return, X should be overwritten by   
                                A * X,   if KASE=1,   
                                A' * X,  if KASE=2,
                        where A' is the conjugate transpose of A,
                        and ZLACON must be re-called with all the other parameters   
                        unchanged.   


                EST    (output) DOUBLE PRECISION   
                        An estimate (a lower bound) for norm(A).   

                KASE   (input/output) INT
                        On the initial call to ZLACON, KASE should be 0.   
                        On an intermediate return, KASE will be 1 or 2, indicating   
                        whether X should be overwritten by A * X  or A' * X.   
                        On the final return from ZLACON, KASE will again be 0.   

                Further Details   
                ======= =======   

                Contributed by Nick Higham, University of Manchester.   
                Originally named CONEST, dated March 16, 1988.   

                Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of 
                a real or complex matrix, with applications to condition estimation", 
                ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.   
                ===================================================================== 
            */

                /* Table of constant values */
                int c__1 = 1;
                doublecomplex      zero = {0.0, 0.0};
                doublecomplex      one = {1.0, 0.0};

                /* System generated locals */
                double d__1;
    
                /* Local variables */
                static int iter;
                static int jump, jlast;
                static double altsgn, estold;
                static int i, j;
                double temp;
                double safmin;
                    

                safmin = dlamch_("Safe minimum");
                if ( *kase == 0 ) {
	            for (i = 0; i < *n; ++i) {
	                x[i].r = 1. / (double) (*n);
	                x[i].i = 0.;
	            }
	            *kase = 1;
	            jump = 1;
	            return 0;
                }

                switch (jump) {
	            case 1:  goto L20;
	            case 2:  goto L40;
	            case 3:  goto L70;
	            case 4:  goto L110;
	            case 5:  goto L140;
                }

                /*     ................ ENTRY   (JUMP = 1)   
	                FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X. */
                L20:
                if (*n == 1) {
	            v[0] = x[0];
	            *est = z_abs(&v[0]);
	            /*        ... QUIT */
	            goto L150;
                }
                *est = dzsum1_(n, x, &c__1);

                for (i = 0; i < *n; ++i) {
	            d__1 = z_abs(&x[i]);
	            if (d__1 > safmin) {
	                d__1 = 1 / d__1;
	                x[i].r *= d__1;
	                x[i].i *= d__1;
	            } else {
	                x[i] = one;
	            }
                }
                *kase = 2;
                jump = 2;
                return 0;

                /*     ................ ENTRY   (JUMP = 2)   
	                FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. */
            L40:
                j = izmax1_(n, &x[0], &c__1);
                --j;
                iter = 2;

                /*     MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
            L50:
                for (i = 0; i < *n; ++i) x[i] = zero;
                x[j] = one;
                *kase = 1;
                jump = 3;
                return 0;

                /*     ................ ENTRY   (JUMP = 3)   
	                X HAS BEEN OVERWRITTEN BY A*X. */
            L70:
#ifdef USE_VENDOR_BLAS
                zcopy_(n, &x[0], &c__1, &v[0], &c__1);
#else
                zcopy_(*n, &x[0],&v[0]);
#endif
                estold = *est;
                *est = dzsum1_(n, v, &c__1);


            L90:
                /*     TEST FOR CYCLING. */
                if (*est <= estold) goto L120;

                for (i = 0; i < *n; ++i) {
	            d__1 = z_abs(&x[i]);
	            if (d__1 > safmin) {
	                d__1 = 1 / d__1;
	                x[i].r *= d__1;
	                x[i].i *= d__1;
	            } else {
	                x[i] = one;
	            }
                }
                *kase = 2;
                jump = 4;
                return 0;

                /*     ................ ENTRY   (JUMP = 4)   
	                X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
            L110:
                jlast = j;
                j = izmax1_(n, &x[0], &c__1);
                --j;
                if (x[jlast].r != (d__1 = x[j].r, fabs(d__1)) && iter < 5) {
	            ++iter;
	            goto L50;
                }

                /*     ITERATION COMPLETE.  FINAL STAGE. */
            L120:
                altsgn = 1.;
                for (i = 1; i <= *n; ++i) {
	            x[i-1].r = altsgn * ((double) (i - 1) / (double) (*n - 1) + 1.);
	            x[i-1].i = 0.;
	            altsgn = -altsgn;
                }
                *kase = 1;
                jump = 5;
                return 0;
    
                /*     ................ ENTRY   (JUMP = 5)   
	                X HAS BEEN OVERWRITTEN BY A*X. */
            L140:
                temp = dzsum1_(n, x, &c__1) / (double) (*n * 3) * 2.;
                if (temp > *est) {
#ifdef USE_VENDOR_BLAS
	            zcopy_(n, &x[0], &c__1, &v[0], &c__1);
#else
                zcopy_(*n, &x[0], &v[0]);
#endif
	            *est = temp;
                }

            L150:
                *kase = 0;
                return 0;

            } /* zlacon_ */

        };
    };
};