From seager@zeus.llnl.gov Fri Apr  8 15:46:51 1988
Date: Fri, 8 Apr 88 13:45:57 PDT

sge.shar   C routines to do the  LINPACK functions   SGECO/SGEFA/SGESL
(dense  matrix conditon  number estimate,  factorization  and solution
routines) and some of  the BLAS in  C.  There is a driver  which shows
how to set up column oriented matrices in C for these routines.
Submitted by Mark K. Seager (seager@lll-crg.llnl.gov) 4/8/88.

========================== SGE.SHAR ====================================
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
#	READ.ME
#	makefile
#	ge.h
#	sample.c
#	sgeco.c
#	sgefa.c
#	sgesl.c
#	blas.c
#	sample.out
# This archive created: Fri Apr  8 12:48:29 1988
export PATH; PATH=/bin:/usr/bin:$PATH
if test -f 'READ.ME'
then
	echo shar: "will not over-write existing file 'READ.ME'"
else
cat << \SHAR_EOF > 'READ.ME'
**********************************************************************
**********************************************************************
**********************************************************************
********************	<<< DISCLAIMER >>>	**********************
**********************************************************************
**********************************************************************
**********************************************************************
This  document was prepared   as an account  of  work sponsored  by an
agency of the  United States  Government.   Neither  the United States
Government  nor the University     of  California nor   any   of their
employees, makes any warranty,  express  or implied,  or assumes   any
legal liability or responsibility for  the accuracy, completeness,  or
usefulness  of   any  information,  apparatus,    product,  or process
disclosed, or represents  that  its use would not   infringe privately
owned  rights.  Reference herein to  any specific commercial products,
process, or    service  by trade   name, trademark,  manufacturer,  or
otherwise, does not  necessarily constitute or  imply its endorsement,
recommendation, or  favoring by  the United States  Government or  the
University of California.  The views and opinions of authors expressed
herein do not necessarily state or reflect those of  the United States
Government  thereof, and shall not be  used for advertising or product
endorsement purposes.
**********************************************************************
**********************************************************************
**********************************************************************

	   This is the READ.ME for the SGEFA/SGESL package.

	SGECO/FA   is agaussian  elimination   (with partial pivoting)
package  written  in C  (single  precision).  The  SGECO routine  also
estimates the  condition number of the matrix.   It  is  a translation
from the FORTRAN LINPACK routines of the same  names.  The matrix data
structure was modified to account for the conventions  of C.  That is,
all elements of the matrix are referenced from 0 instead  of 1 and the
matrix is stored as a series of vectors  (which are the columns).  The
matrix data  structure is defined in  ge.h.   The user must supply the
storage for  the columns of  the matrix  by calls to the standard UNIX
memory allocation routine MALLOC.  See  the matgen routine in sample.c
for an example.  The necessary  BLAS routines are included in  blas.c.
Easy access to the matrix  elements can be  had by looking at the elem
and pelem macros defined in ge.h.

	A sample  output (from a  Sun-3/160   workstation with FPA) is
included in the file sample.out

	Send comments, suggestions and bug-reports to:
		Dr. Mark K. Seager
		Lawrence Livermore Nat. Lab.
		PO Box 808, L-316
		Livermore, CA 94550
		(415) 423-3141
		seager@lll-crg.llnl.gov
SHAR_EOF
if test 2935 -ne "`wc -c < 'READ.ME'`"
then
	echo shar: "error transmitting 'READ.ME'" '(should have been 2935 characters)'
fi
fi
if test -f 'makefile'
then
	echo shar: "will not over-write existing file 'makefile'"
else
cat << \SHAR_EOF > 'makefile'
#
#	This make file peices together the sgeco/fa/sl test routine(s)
#
#	@(#)makefile	1.3 4/8/88
#
#	To test this package compile with "make sample".  If you get
#	an unstatisfied external for "copysign" then change the compile
#	options for the blas.c file, see below, and recompile.
#
CFLAGS=-g -ffpa -fsingle
LDFLAGS=-lm
OBJS = sgeco.o sgefa.o sgesl.o blas.o

sample : sample.o $(OBJS)
	$(CC) $(CFLAGS) -o sample sample.o $(OBJS) $(LDFLAGS)
sgefat: driver.o $(OBJS)
	$(CC) $(CFLAGS) -o sgefat driver.o $(OBJS) $(LDFLAGS)

sample.o : sample.c ge.h
driver.o : driver.c ge.h

sgeco.o : sgeco.c ge.h 
sgefa.o: sgefa.c ge.h
sgesl.o: sgesl.c ge.h

blas.o:	blas.c
#
#	Use the following if you don't have copysign in your 3M library 
#	(i.e., uncomment the next two lines) and comment the line above out.
#blas.o : blas.c
#	$(CC) $(CFLAGS) -DCOPYSIGN -c blas.c

clean:
	rm -f *.o core gmon.out sgefat sample
SHAR_EOF
if test 901 -ne "`wc -c < 'makefile'`"
then
	echo shar: "error transmitting 'makefile'" '(should have been 901 characters)'
fi
fi
if test -f 'ge.h'
then
	echo shar: "will not over-write existing file 'ge.h'"
else
cat << \SHAR_EOF > 'ge.h'
/* SCCS ID @(#)ge.h	1.1		2/4/86 */
  /***************************************************************
 ******************************************************************
****	 Matrix data structure(s) for Gaussian Elimination	****
 ******************************************************************
  ****************************************************************/
/* 
   This file contains the definitions of the structures used in
   various algorithms for doing Gaussian Elimination.

   The following gives an array (of length 10) of pointers to floats.
     float *a[10]; 
   Now assume that each a[i] points to space for an array of floats (gotten
   by a call to malloc, say).
   Then the following is true:
        a[i]       can be thought of as a pointer to the i-th array of floats,
        *(a[i]+j)  is the j-th element of the i-th array.

   The following shows how to reference things for the definition of the FULL
   structure given below.
      a->cd is the value of (as apposed to a pointer to) the column dimension.
      a->rd is the value of (as apposed to a pointer to) the row dimension.
      a->pd[j] is a pointer to the j-th column (an array of floats).
      *(a->pd[j]+i) is the i-th element of the j-th column, viz., a(i,j).
   Here we think, as is natural in C, of all matrices and vectors indexed 
   from 0 insead of 1.
*/

#define MAXCOL 1000	/* Maximum number of Columns. */

struct FULL {		/* Struct definition for the FULL matrix structure. */
  int	cd;		/* Column dimension of the matrix. */
  int	rd;		/* Row Dimension of the matrix. */
  float	*pd[MAXCOL]; 	/* Array of pointers to the columns of a matrix. */
};

/* The following macro will get a(r,c) from a matrix in the FULL structure. */
#define elem(a,r,c)	(*(a.pd[(c)]+(r)))

/* The following macro will get a(r,c) from a pointer to a matrix 
   in the FULL structure. */
#define pelem(a,r,c)	(*(a->pd[(c)]+(r)))
SHAR_EOF
if test 1912 -ne "`wc -c < 'ge.h'`"
then
	echo shar: "error transmitting 'ge.h'" '(should have been 1912 characters)'
fi
fi
if test -f 'sample.c'
then
	echo shar: "will not over-write existing file 'sample.c'"
else
cat << \SHAR_EOF > 'sample.c'
  /***************************************************************
 ******************************************************************
****		    Sample driver for SGEFA/SL and CO		****
 ******************************************************************
  ****************************************************************/
/*
 * SYSTEM DEPENENT ROUTINES:
 *	double second();	Returns cpu time used in seconds
 */

#include <stdio.h>
#include <math.h>
#include "ge.h"
#define  MAXN	300			/* Maximum problem size. */

static char SAMPLESid[] = "@(#)sample.c	1.1  4/8/88";

main()
{
    /* Storage for the linear system and associates. */
    struct FULL a;
    float x[MAXCOL], b[MAXCOL], z[MAXCOL], rcond;
    int ipvt[MAXCOL], sgefa();
    void sgesl();

    /* Local vars. */
    double opsfa, opsco, opssl, xerrnrm, resnrm, xnrm;
    double tfai, tfao, tsli, tslo, tcoi, tcoo, second();
    int n, i, j, retval;
    void matgen();

#ifdef sun
  /* The following is a kludge for correct Sun FPA core dumps. */
  unsigned newmode = 62464;
  fpmode_(&newmode);
#endif

    for( n=50; n<=MAXN; n+=50 ) {

	/* Generate the linear system. */
	matgen( n, &a, b );
	opsfa = 1.0E-6*((2.0e0*(double)n*(double)n*(double)n)/3.0e0 + 
	  2.0e0*(double)n*(double)n);
	opsco = opsfa + 1.0E-6*(6.0e0*(double)n*(double)n);
	opssl = 1.0e-6*(2.0e0*(double)n*(double)n);

	/* Factor. */
	tfai = second();
	retval = sgefa( &a, ipvt );
	tfao = second();

	if( retval ) 
	  printf("Zero Column %d found.\n", retval );
	else {

	    /* Solve the system. */
	    tsli = second();
	    sgesl( &a, ipvt, b, 0 );
	    tslo = second();

	    /* Compute a residual to verify results. */
	    for( i=0; i<n; i++ ) x[i] = b[i];
	    matgen( n, &a, b );
	    for( j=0; j<n; j++ ) {
		register float *col = a.pd[j];
		register float xj = x[j];

		for( i=0; i<n; i++ ) {
		    b[i] -= col[i]*xj;
		}
	    }
	    xerrnrm = 0.0;
	    resnrm = 0.0;
	    for( i=0; i<n; i++ ) {
		xerrnrm += ((double)x[i]-1.0)*((double)x[i]-1.0);
		resnrm += (double)b[i]*(double)b[i];
	    }
	    xnrm = sqrt( (double)n );
	    xerrnrm = sqrt( xerrnrm )/xnrm;
	    resnrm = sqrt( resnrm )/xnrm;

	    /* Now factor and est condition number. */
	    tcoi = second();
	    retval = sgeco( &a, ipvt, &rcond, z );
	    tcoo = second();

	    printf("\n\n\tLinpack Benchmark in C of size %d\n",n);
	    printf("SGEFA time = %15e (Sec) = %15e (Mflops)\n",tfao-tfai,opsfa/(tfao-tfai) );
	    printf("SGECO time = %15e (Sec) = %15e (Mflops)\n",tcoo-tcoi,opsco/(tcoo-tcoi) );
	    if( tslo-tsli != 0.0 )
	      printf("SGESL time = %15e (Sec) = %15e (Mflops)\n",tslo-tsli,opssl/(tslo-tsli) );
	    else
	      printf("SGESL time = %15e (Sec)                \n",tslo-tsli );
	    printf(" 1/COND(A) (Condition number of A) = %15e\n",rcond);
	    printf("   ||x-X||/||X||  (Solution Error) = %15e\n",xerrnrm);
	    printf("   ||b-Ax||/||X|| (Residual Error) = %15e\n",resnrm);
	    printf("\tWhere X = True solution and and x = computed solution\n");
	}
    }
}
void matgen( n, a, b )
int n;
struct FULL *a;
float *b;
{

    /* This routine generates a struct FULL matrix */
    int i, j, init = 1325;
    char *malloc();
    int free();

    /* Allocate/Deallocate space for the columns. */
    for( i=0; i<n; i++ ) {
	if( a->pd[i] != NULL ) (void)free( (char *)a->pd[i] );
	if( (a->pd[i] = (float *)malloc( (unsigned)n*sizeof(float) )) == NULL ) {
	    fprintf(stderr, "MATGEN: Error allocating matrix for n=%d\n",n);
	    exit( 1 );
	}	   
    }

    /* Set the matrix elements and right hand side. */
    a->cd = n;
    a->rd = n;
    for( j=0; j<n; j++ ) {
	register float *col = a->pd[j];

	for( i=0; i<n; i++ ) {
	    init = 3125*init % 65536;
	    col[i] = ((float)init-32768.0)/16384.0;
	}
    }
    for( j=0; j<n; j++ ) b[j] = 0.0;
    for( j=0; j<n; j++ ) {
	register float *col = a->pd[j];

	for( i=0; i<n; i++ ) {
	    b[i] += col[i];
	}
    }
}

/****************************************************************************
 *		Routines for getting elapsed CPU usage.
 ****************************************************************************/
#include <sys/time.h>
#include <sys/resource.h>
double second()
/****************************************************************************
 *		Returns the total cpu time used in seconds.
 *	Emulates the Cray fortran library function of the same name.
 ****************************************************************************/
{
    struct rusage buf;
    double temp;

    getrusage( RUSAGE_SELF, &buf );
    
    /* Get system time and user time in micro-seconds.*/
    temp = (double)buf.ru_utime.tv_sec*1.0e6 + (double)buf.ru_utime.tv_usec +
           (double)buf.ru_stime.tv_sec*1.0e6 + (double)buf.ru_stime.tv_usec;

    /* Return the sum of system and user time in SECONDS.*/
    return( temp*1.0e-6 );
}
SHAR_EOF
if test 4812 -ne "`wc -c < 'sample.c'`"
then
	echo shar: "error transmitting 'sample.c'" '(should have been 4812 characters)'
fi
fi
if test -f 'sgeco.c'
then
	echo shar: "will not over-write existing file 'sgeco.c'"
else
cat << \SHAR_EOF > 'sgeco.c'
  /***************************************************************
 ******************************************************************
****	    Gaussian Elimination with partial pivoting		****
**** 		 and condition number estimation.		****
****	  This file contains the factorization driver and 	****
****	      contion number estimation routine SGECO		****
 ******************************************************************
  ****************************************************************/

#include <math.h>
#include "ge.h"
static char SGECOSid[] = "@(#)sgeco.c	1.1  4/8/88";

int sgeco( a, ipvt, rcond, z )
struct FULL *a;
int	    *ipvt;
float       *rcond, *z;
/*
  PURPOSE
      SGECO factors a real matrix by gaussian elimination
      and estimates the condition of the matrix.
 
  REMARKS
      If  rcond  is not needed, SGEFA is slightly faster.
      to solve  A*x = b , follow SGECO by SGESL.
      To compute  inverse(A)*c , follow SGECO by SGESL.
      To compute  determinant(A) , follow SGECO by SGEDI.
      To compute  inverse(A) , follow SGECO by SGEDI.
 
  INPUT
 	a	A pointer to the FULL matrix structure.  
		See the definition in ge.h.
 
  OUTPUT
	a	A pointer to the FULL matrix structure containing 
		an upper triangular matrix and the multipliers
		which were used to obtain it.
		The factorization can be written  a = l*u  where
		l  is a product of permutation and unit lower
		triangular matrices and  u  is upper triangular.
	ipvt	An integer vector (of length a->cd) of pivot indices.
	rcond   A float estimate of the reciprocal condition of  A .
                for the system  A*x = b , relative perturbations
                in  A  and  b  of size  epsilon  may cause
                relative perturbations in  x  of size  epsilon/rcond .
                If  rcond  is so small that the logical expression
                            1.0 + rcond .eq. 1.0
                is true, then  a  may be singular to working
                precision.  In particular,  rcond  is zero  if
                exact singularity is detected or the estimate
                underflows.
         z      A float vector (of length a->cd) for a work vector 
	        whose contents are usually unimportant.  If  A  is 
		close to a singular matrix, then  z  is an approx-
		imate null vector in the sense that:
		         norm(a*z) = rcond*norm(a)*norm(z) .

    RETURNS
		= -1  Matrix is not square.
		=  0  Normal return value.
		=  k  if  u(k,k) .eq. 0.0 .  This is not an error
		      condition for this subroutine, but it does
		      indicate that sgesl or sgedi will divide by zero
		      if called.  Use  rcond  in sgeco for a reliable
		      indication of singularity.

    ROUTINES
	SGEFA(), blas sasum() and sdot(), copysign(), fabs();

    WARNINGS
        This routine uses the UN*X math library routines 
	copysign() and fabs().
*/
{
    void saxpy(), sscal();
    register int j;
    int k, n, info, sgefa();
    register float s;
    double sasum(), sdot();
    float ek, anorm, ynorm;

    n = a->cd;

    /* Compute 1-norm of A. */
    for( j=0, anorm=0.0; j<n; j++ ) {
	register double sum;

	sum = sasum( n, a->pd[j], 1 );
	anorm = (float)sum > anorm ? (float)sum : anorm;
    }

    /* Factor A. */
    info = sgefa( a, ipvt );
 
    /*
     * rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
     * estimate = norm(z)/norm(y) where  a*z = y  and  trans(a)*y = e .
     * trans(a)  is the transpose of a .  The components of  e  are
     * chosen to cause maximum local growth in the elements of w  where
     * trans(u)*w = e .  The vectors are frequently rescaled to avoid
     * overflow.
     */

    /* solve trans(u)*w = e */
    ek = 1.0;
    for( j=0; j<n; j++ ) {
	z[j] = 0.0e0;
    }
    for( k=0; k<n; k++ ) {
	register float zk = z[k];
	float wk, wkm, sm;
	int kp1 = k+1;

	if( zk != 0.0 ) ek = (float)copysign( (double)ek, (double)-zk );
	if( fabs( (double)(ek-zk) ) > fabs( (double)pelem(a,k,k) ) ) {
            s = (float)fabs( (double)pelem(a,k,k) )/(float)fabs( (double)(ek-zk) );
            sscal( n, (double)s, z, 1 );
	    zk = z[k];
            ek *= s;
	}
	wk = ek - zk;
	wkm = -ek - zk;
	s = (float)fabs( (double)wk );
	sm = (float)fabs( (double)wkm );
	if( pelem(a,k,k) == 0.0 ) {
            wk = 1.0;
            wkm = 1.0;
	}
	else {
            wk /= pelem(a,k,k);
            wkm /= pelem(a,k,k);
	}
	if( kp1<n ) {
	    for( j=kp1; j<n; j++ ) {
		sm = sm + (float)fabs( (double)(z[j]+wkm*pelem(a,k,j)) );
		z[j] += ((float)wk)*pelem(a,k,j);
		s += (float)fabs( (double)z[j] );
	    }
            if( s < sm ) {
		register float t = wkm-wk;

		wk = wkm;
		for( j=kp1; j<n; j++ ) {
		    z[j] += t*pelem(a,k,j);
		}
	    }
	}
	z[k] = wk;
    }
    s = 1.0/(float)sasum( n, z, 1 );
    sscal( n, (double)s, z, 1 );

    /* Solve trans(L)*y = w. */
    for( k=n-1; k>=0; k-- ) {
	register int   l;
	register float t;

	if( k < n-1 ) z[k] += (float)sdot( n-k-1, a->pd[k]+k+1, 1, (z+k+1), 1 );
	if( fabs( (double)z[k] ) > 1.0 ) {
            s = 1.0/(float)fabs( (double)z[k] );
            sscal( n, (double)s, z, 1 );
	}
	l = ipvt[k];
	t = z[l];
	z[l] = z[k];
	z[k] = t;
    }
    s = 1.0/(float)sasum( n, z, 1 );
    sscal( n, (double)s, z, 1);

    ynorm = 1.0;

    /* solve L*v = y. */
    for( k=0; k<n; k++ ) {
	register int   l;
	register float t;

	l = ipvt[k];
	t = z[l];
	z[l] = z[k];
	z[k] = t;
	if( k < n-1 ) saxpy( n-k-1, (double)t, (a->pd[k]+k+1), 1, (z+k+1), 1 );
	if( fabs( (double)z[k] ) > 1.0)  {
            s = 1.0/(float)fabs( (double)z[k] );
            sscal( n, (float)s, z, 1 );
            ynorm *= s;
	}
    }
    s = 1.0/(float)sasum( n, z, 1 );
    sscal( n, (double)s, z, 1 );
    ynorm *= s;

    /* Solve  U*z = v. */
    for( k=n-1; k>=0; k-- ) {
	register float t;

	if( fabs( (double)z[k] ) > fabs( (double)pelem(a,k,k) ) ) {
            s = (float)fabs( (double)pelem(a,k,k) )/(float)fabs( (double)z[k] );
            sscal( n, (float)s, z, 1 );
            ynorm *= s;
	}
	if( pelem(a,k,k) == 0.0 ) z[k] = 1.0;
	else z[k] /= pelem(a,k,k);
	t = -z[k];
	saxpy( k, (double)t, a->pd[k], 1, z, 1 );
    }

    /* Make znorm = 1.0. */
    s = 1.0/(float)sasum( n, z, 1 );
    sscal( n, (double)s, z, 1 );
    ynorm *= s;

    if( anorm == 0.0e0)  *rcond = 0.0;
    else *rcond = ynorm/anorm;

    return( info );
}
SHAR_EOF
if test 6315 -ne "`wc -c < 'sgeco.c'`"
then
	echo shar: "error transmitting 'sgeco.c'" '(should have been 6315 characters)'
fi
fi
if test -f 'sgefa.c'
then
	echo shar: "will not over-write existing file 'sgefa.c'"
else
cat << \SHAR_EOF > 'sgefa.c'
  /***************************************************************
 ******************************************************************
****	    Gaussian Elimination with partial pivoting.		****
****	This file contains the factorization routine SGEFA	****
 ******************************************************************
  ****************************************************************/

#include "ge.h"
static char SGEFASid[] = "@(#)sgefa.c	1.3  4/8/88";

int sgefa( a, ipvt )
struct FULL *a;
int	    *ipvt;
/*
    PURPOSE
	SGEFA factors a real matrix by gaussian elimination.

    REMARKS
	SGEFA is usually called by SGECO, but it can be called
	directly with a saving in time if  rcond  is not needed.
	(time for SGECO) = (1 + 9/n)*(time for SGEFA) .

    INPUT
	a	A pointer to the FULL matrix structure.  
		See the definition in ge.h.

    OUTPUT
	a	A pointer to the FULL matrix structure containing 
		an upper triangular matrix and the multipliers
		which were used to obtain it.
		The factorization can be written  a = l*u  where
		l  is a product of permutation and unit lower
		triangular matrices and  u  is upper triangular.
	ipvt	An integer vector (of length a->cd) of pivot indices.

    RETURNS
		= -1  Matrix is not square.
		=  0  Normal return value.
		=  k  if  u(k,k) .eq. 0.0 .  This is not an error
		      condition for this subroutine, but it does
		      indicate that sgesl or sgedi will divide by zero
		      if called.  Use  rcond  in sgeco for a reliable
		      indication of singularity.

    ROUTINES
	blas ISAMAX
*/
{
    register int i, j;
    int		isamax(), k, l, nm1, info, n;
    float	*akk, *alk;
    register float t, *mik;
    
    /* Gaussian elimination with partial pivoting. */
    if( a->cd != a->rd ) return( -1 );
    n    = a->cd;
    nm1  = n - 1;
    akk  = a->pd[0];
    info = 0;				/* Assume nothing will go wrong! */
    if( n < 2 ) goto CLEAN_UP;
    
    /*  Loop over Diagional */
    for( k=0; k<nm1; k++, ipvt++ ) {
	
	/* Find index of max elem in col below the diagional (l = pivot index). */
	akk   = a->pd[k] + k;
	l     = isamax( n-k, akk, 1 ) + k;
	*ipvt = l;
	
	/* Zero pivot implies this column already triangularized. */
	alk = a->pd[k] + l;
	if( *alk == 0.0e0) {
	    info = k;
	    continue;
	}
	
	/* Interchange a(k,k) and a(l,k) if necessary. */
	if( l != k ) {
	    t    = *alk;
	    *alk = *akk;
	    *akk = t;
	}
	
	/* Compute multipliers for this column. */
	t = -1.0e0 / (*akk);
	for( i=k+1, mik=a->pd[k]; i<n; i++ )
	  mik[i] *= t;
	
	/* Column elimination with row indexing. */
	
	if( l != k ) {
	    /* Interchange a(k,j) and a(l,j) if necessary. */
	    for( j=k+1; j<n; j++ ) {
		t = pelem(a,k,j);
		pelem(a,k,j) = pelem(a,l,j);
		pelem(a,l,j) = t;
	    }
	}
	for( j=k+1; j<n; j++ ) {
	    register float *aij = a->pd[j];

	    t = pelem(a,k,j);
	    for( i=k+1, mik=a->pd[k]; i<n; i++ ) 
	      aij[i] += t*mik[i];
	}
    }			/* End of for k loop */
    
  CLEAN_UP: 
    *ipvt = nm1;
    if( *akk == 0.0e0 ) info = n;
    return( info );
}
SHAR_EOF
if test 3028 -ne "`wc -c < 'sgefa.c'`"
then
	echo shar: "error transmitting 'sgefa.c'" '(should have been 3028 characters)'
fi
fi
if test -f 'sgesl.c'
then
	echo shar: "will not over-write existing file 'sgesl.c'"
else
cat << \SHAR_EOF > 'sgesl.c'
  /***************************************************************
 ******************************************************************
****	    Gaussian Elimination with partial pivoting.		****
****	This file contains the solution routine SGESL		****
 ******************************************************************
  ****************************************************************/

#include "ge.h"
static char SGESLSid[] = "@(#)sgesl.c	1.3  4/8/88";

int sgesl( a, ipvt, b, job )
struct FULL *a;
int 	    *ipvt, job;
float	    b[];
/*
    PURPOSE
	SGESL solves the real system
	a * x = b  or  trans(a) * x = b
	using the factors computed by SGECO or SGEFA.

    INPUT
	a	A pointer to the FULL matrix structure containing the factored
		matrix.  See the definition of FULL in ge.h.
	ipvt    The pivot vector (of length a->cd) from SGECO or SGEFA.
	b       The right hand side vector (of length a->cd).
	job     = 0         to solve  a*x = b ,
		= nonzero   to solve  trans(a)*x = b  where
			    trans(a)  is the transpose.

    OUTPUT
	b       The solution vector x. 

    REMARKS
	Error condition:
	A division by zero will occur if the input factor contains a
	zero on the diagonal.  Technically this indicates singularity
	but it is often caused by improper arguments or improper
	setting of lda .  It will not occur if the subroutines are
	called correctly and if sgeco has set rcond .gt. 0.0
	or sgefa has set info .eq. 0 .
*/
{
    register float t, *mik;
    float	*akk;
    register int i, k;
    int	l, n, nm1;
    
    n   = a->cd;
    nm1 = n - 1;
    
    /* job = 0 , solve  A * x = b.  */
    if( job == 0 ) {

	/* Forward elimination. Solve L*y = b. */
	for( k=0; k<nm1; k++ ) {
	    akk = a->pd[k] + k;		/* akk points to a(k,k). */
	    
	    /* Interchange b[k] and b[l] if necessary. */
	    l = ipvt[k];
	    t = b[l];
	    if( l != k ) {
		b[l] = b[k];
		b[k] = t;
	    }
	    for( i=k+1, mik=a->pd[k]; i<n; i++ )
	      b[i] += t*mik[i];
	}
	
	/* Back substitution.  Solve  U*x = y. */
	for( k=nm1; k>=0; k-- ) {
	    register float *uik = a->pd[k];
	    akk = uik + k;
	    b[k] /= (*akk);
	    for( i=0; i<k; i++ )
	      b[i] -= uik[i]*b[k];
	}
	return;
    }
    
    /* job = nonzero.  Solve  trans(A) * x = b. */
    /* First solve trans(U)*y = b. */
    for( k=0; k<n; k++ ) {
	register float *uik = a->pd[k];
	akk = uik + k;
	t = 0.0;
	for( i=0; i<k; i++ )
	  t += uik[i]*b[i];
	b[k] = (b[k] - t) / (*akk);
    }
    
    /* b now contains y. */
    /* Solve trans(L)*x = y. */
    for( k=n-2; k>=0; k-- ) {
	mik = a->pd[k];
	t = 0.0;
	for( i=k+1; i<n; i++ )
	  t += mik[i]*b[i];
	b[k] += t;
	
	/* Interchange b(k) and b(ipvt(k)) if necessary. */
	l    = ipvt[k];
	if( l == k ) continue; 
	t    = b[l];
	b[l] = b[k];
	b[k] = t;
    }
    return;
}
SHAR_EOF
if test 2782 -ne "`wc -c < 'sgesl.c'`"
then
	echo shar: "error transmitting 'sgesl.c'" '(should have been 2782 characters)'
fi
fi
if test -f 'blas.c'
then
	echo shar: "will not over-write existing file 'blas.c'"
else
cat << \SHAR_EOF > 'blas.c'
  /***************************************************************
  *****************************************************************
 *******************************************************************
*****								*****
*****				BLAS				*****
*****		Basic Linear Algebra Subroutines		*****
*****	      Written in the C Programming Language.		*****
*****								*****
*****	Functions include:					*****
*****	isamax, sasum, saxpy, scopy, sdot, snrm2,		*****
*****								*****
*****   In addition a few other routines are included:		*****
*****	vexopy, vfill						*****
*****								*****
*****	If your 3M library does not have the copysign function	*****
*****   then compile this file with -DCOPYSIGN and one will be	*****
*****   be supplied.						*****
 *******************************************************************
  *****************************************************************
   ***************************************************************/
#include <math.h>
#ifdef vax
#define HUGE	1.701411733192644270e38
#endif
#ifdef mc68000
#ifndef HUGE
#define HUGE	1.79769313486231+308	/* IEEE infinity */
#endif
#define HUGEsp	3.40282346638529E+38	/* Largest binary float. */
#define SMALLsp	1.17549435082229E-38	/* Smallest (positive) binary float. */
#endif
#ifndef HUGEsp
#define HUGEsp 1.0e+38
#endif
#ifndef SMALLsp
#define SMALLsp 1.0e-45
#endif

static char BLASSid[] = "@(#)blas.c	1.3 4/8/88";

int isamax( n, sx, incx )
float	*sx;
int	n, incx;
/*
    PURPOSE
        Finds the index of element having max. absolute value.

    INPUT
	n	Number of elements to check.
	sx	Vector to be checked.
	incx	Every incx-th element is checked.

*/
{
    register float smax = 0.0e0;
    register int i, istmp = 0;

#ifndef abs
#define abs(x) ((x)<0.0?-(x):(x)) 
#endif
#ifdef alliant
  int isamax_();

  istmp = isamax_( &n, sx, &incx )-1;	/* Remember 0 is first. */
  return( istmp );
#else
  if( n <= 1 ) return( istmp );
  if( incx != 1 ) {
    /* Code for increment not equal to 1. */
    if( incx < 0 ) sx = sx + ((-n+1)*incx + 1);
    istmp = 0;
    smax  = abs( *sx );
    sx += incx;
    for( i=1; i<n; i++, sx+=incx ) 
      if( abs( *sx ) > smax ) {
	istmp = i;
	smax  = abs( *sx );
      }
    return( istmp );
  }
  /* Code for increment equal to 1. */
  istmp = 0;
  smax  = abs(*sx);
  sx++;
  for( i=1; i<n; i++, sx++ ) 
    if( abs( *sx ) > smax ) { 
      istmp = i;
      smax  = abs( *sx );
    }
  return( istmp );
#endif
}
 
double sasum( n, sx, incx )
float *sx;
int   n, incx;
/*
  PURPOSE
      Returns sum of magnitudes of single precision SX.
      sasum = sum from 0 to n-1 of  ABS(SX(1+I*INCX))

  INPUT
    n		Number of elements to multiply.
    sx		Pointer to float vector to take abs sum of.
    incx	Storage incrament for sx.

  RETURNS
    sasum       Double variable with the result.

  WARNINGS
    This routine uses the UN*X math library function fabs().

*/
{
#ifndef abs
#define abs(x) ((x)<0.0?-(x):(x)) 
#endif
    register float sum = 0.0;
    
    if( n<= 0 ) return( sum );
    if( incx == 1 ) {
	register int i, m;

        /* Code for increments equal to 1. */

        /* Clean-up loop so remaining vector length is a multiple of 6. */
	m = n % 6;
	if( m != 0 ) {
	    for( i=0; i<m; i++ )
	      sum += abs( sx[i] );
	if( n < 6 ) return( (double)sum );
	}
	for( i=m; i<n; i+=6 ) {
	    sum += abs( sx[i] ) + abs( sx[i+1]) + abs( sx[i+2] ) + 
	      abs( sx[i+3] ) + abs( sx[i+4] ) + abs( sx[i+5] );
	}
	return( (double)sum );
    } 
    else {
	register int i, ns;

        /* Code for increments not equal to 1. */
	ns = n*incx;
	for( i=0; i<ns; i+=incx ) {
	    sum += abs( sx[i] );
	}
	return( (double)sum );
    }
}

void saxpy( n, sa, sx, incx, sy, incy )
float *sx, *sy;
double sa;
int   n, incx, incy;
/*
  PURPOSE
    Vector times a scalar plus a vector.  sy = sy + sa*sx.

  INPUT
    n		Number of elements to multiply.
    sa		Scalar to multiply by (note that this is a double).
    sx		Pointer to float vector to scale.
    incx	Storage incrament for sx.
    sy		Pointer to float vector to add.
    incy	Storage incrament for sy.

  OUTPUT
    sy		sy = sy + sa*sx
*/
{
    register int i;
    register float ssa = (float)sa;

#ifdef alliant
    void saxpy_();
    float ssaa = (float)sa;
    
    saxpy_( &n, &ssaa, sx, &incx, sy, &incy );
    return;
#else
    
    if( n<=0 || ssa==0.0 ) return;
    if( incx == incy ) {
	if( incx == 1 ) {

	    /* Both increments = 1 */
	    for( i=0; i<n; i++ ) 
	      sy[i] += ssa*sx[i];
	    return;
	}
	if( incx>0 ) {

	    /* Equal, positive, non-unit increments. */
	    for( i=0; i<n; i++,sx+=incx,sy+=incx )
	      *sy += ssa*(*sx);
	    return;
	}
    }

    /* Unequal or negative increments. */
    if( incx < 0 ) sx += ((-n+1)*incx + 1);
    if( incy < 0 ) sy += ((-n+1)*incy + 1);
    for( i=0; i<n; i++,sx+=incx,sy+=incy ) 
      *sy += ssa*(*sx);
#endif
}
 
void saxpyx( n, sa, sx, incx, sy, incy )
float *sx, *sy;
double sa;
int   n, incx, incy;
/*
  PURPOSE
    Vector times a scalar plus a vector.  sx = sy + sa*sx.

  INPUT
    n		Number of elements to multiply.
    sa		Scalar to multiply by (note that this is a double).
    sx		Pointer to float vector to scale.
    incx	Storage incrament for sx.
    sy		Pointer to float vector to add.
    incy	Storage incrament for sy.

  OUTPUT
    sx		sx = sy + sa*sx
*/
{
    register i;
    register float ssa = (float)sa;

    if( n<=0 || ssa==0.0 ) return;
    if( incx == incy ) {
	if( incx == 1 ) {

	    /* Both increments = 1 */
	    for( i=0; i<n; i++ ) 
	      sx[i] = sy[i] + ssa*sx[i];
	    return;
	}
	if( incx>0 ) {

	    /* Equal, positive, non-unit increments. */
	    for( i=0; i<n; i++, sx+=incx, sy+=incx)
	      *sx = *sy + ssa*(*sx);
	    return;
	}
    }

    /* Unequal or negative increments. */
    if( incx < 0 ) sx += ((-n+1)*incx + 1);
    if( incy < 0 ) sy += ((-n+1)*incy + 1);
    for( i=0; i<n; i++,sx+=incx,sy+=incy ) 
      *sx = *sy + ssa*(*sx);
}

void scopy( n, sx, incx, sy, incy )
float  *sx, *sy;
int     n, incx, incy;
/*
    PURPOSE
        Copies vector sx into vector sy.
 
    INPUT
        n    Number of components to copy.
	sx   Source vector
	incx Index increment for sx.
        incy Index increment for sy.
 
    OUTPUT
        sy   Destination vector.
*/
{
    register int i;

#ifdef alliant
    void scopy_();

    scopy_( &n, sx, &incx, sy, &incy );
    return;
#else

    if( n<1  ) return;
    if( incx == incy ) {
	if( incx == 1 ) {

	    /* Both increments = 1 */
	    for( i=0; i<n; i++ )
	      sy[i] = sx[i];
	    return;
	}
	if( incx > 0 ) {

	    /* Equal, positive, non-unit increments. */
	    for( i=0; i<n; i++, sx+=incx, sy+=incx)
	      *sy = *sx;
	    return;
	}
    }

    /* Non-equal or negative increments. */
    if( incx < 0 ) sx += ((-n+1)*incx + 1);
    if( incy < 0 ) sy += ((-n+1)*incy + 1);
    for( i=0; i<n; i++,sx+=incx,sy+=incy ) 
      (*sx) = (*sy);
    return;
#endif
}

double sdot( n, sx, incx, sy, incy )
float	*sx, *sy;
int	n, incx, incy;
/*
    PURPOSE
        Forms the dot product of a vector.

    INPUT
        n       Number of elements to sum.
        sx      Address of first element of x vector.
        incx    Incrament for the x vector.
        sy      Address of first element of y vector.
        incy    incrament for the y vector.

    OUPUT
        sdot    Dot product x and y.  Double returned
		due to `C' language features.
*/
{
    register int i;
    register float stemp = 0.0e0;

    if( n<1 ) return( (double)stemp );
    if( incx == incy ) {
	if( incx == 1 ) {

	    /* Both increments = 1 */
	    for( i=0; i<n; i++ )
	      stemp += sx[i]*sy[i];
	    return( (double)stemp );
	}
	if( incx>0 ) {

	    /* Equal, positive, non-unit increments. */
	    for( i=0; i<n; i++, sx+=incx, sy+=incx)
	      stemp += (*sx)*(*sy);
	    return( (double)stemp );
	}
    }

    /* Unequal or negative increments. */
    if( incx < 0 ) sx += ((-n+1)*incx + 1);
    if( incy < 0 ) sy += ((-n+1)*incy + 1);
    for( i=0; i<n; i++,sx+=incx,sy+=incy ) 
      stemp += (*sx)*(*sy);
    return( (double)stemp );
}				/* End of ---SDOT--- */

double snrm2( n, sx, incx )
float	*sx;
int	n, incx;
/*
    PURPOSE
        Computes the Euclidean norm of sx while being
	very careful of distructive underflow and overflow.

    INPUT
        n       Number of elements to use.
        sx      Address of first element of x vector.
        incx    Incrament for the x vector (>0).

    OUPUT
        snrm2   Euclidean norm of sx.  Returns double
		due to `C' language features.
    REMARKS
        This algorithm proceeds in four steps.
	1) scan zero components.
	2) do phase 2 when component is near underflow.
*/
{
  register int i;
  static double cutlo, cuthi;
  double sum = 0.0e0, hitst, r1mach();
  float xmax;

  if( n<1 || incx<1 ) return( sum );

  /* Calculate near underflow */
  if( cutlo == 0.0 ) cutlo = sqrt( SMALLsp/r1mach() );
  /* Calculate near  overflow */
  if( cuthi == 0.0 ) cuthi = sqrt( HUGEsp );
  hitst = cuthi/(double) n;
  i     = 0;

  /* Zero Sum. */
  while( *sx == 0.0 && i<n ) {
    i++;
    sx += incx;
  }
  if( i>=n ) return( sum );

START:
  if( abs( *sx ) > cutlo ) {
    for( ; i<n; i++, sx+=incx ) {		/* Loop over elements. */
      if( abs( *sx ) > hitst ) goto GOT_LARGE;  
      sum += (*sx) * (*sx);
    }
    sum = sqrt( sum );
    return( sum );				/* Sum completed normaly. */
  }
  else {					/* Small sum prepare for phase 2. */
    xmax  = abs( *sx );
    sx += incx;
    i++;
    sum += 1.0;
    for( ; i<n; i++, sx+=incx ) {
      if( abs( *sx ) > cutlo ) {		/* Got normal elem.  Rescale and process. */
	sum = (sum*xmax)*xmax;
	goto START;
      }
      if( abs( *sx ) > xmax ) {
	sum  = 1.0 + sum*(xmax /(*sx))*(xmax /(*sx));
	xmax = abs( *sx );
	continue;
      }
      sum += ((*sx)/xmax)*((*sx)/xmax);
    }
    return( (double)xmax*sqrt( sum ) );
  }						/* End of small sum. */

 GOT_LARGE:
  sum  = 1.0 + (sum/(*sx))/(*sx);		/* Rescale and process. */
  xmax = abs( *sx );
  sx   += incx;
  i++;
  for( ; i<n; i++, sx+=incx ) {
    if( abs( *sx ) > xmax ) {
      sum  = 1.0 + sum*(xmax /(*sx))*(xmax /(*sx));
      xmax = abs( *sx );
      continue;
    }
    sum += ((*sx)/xmax)*((*sx)/xmax);
  }
  return( (double)xmax*sqrt( sum ) );		/* End of small sum. */
}						/* End of ---SDOT--- */

double r1mach()
/***********************************************************************
 ****   This routine computes the unit roundoff for single precision 
 ****	of the machine.  This is defined as the smallest positive 
 ****	machine number u such that  1.0 + u .ne. 1.0
 ****	Returns a double due to `C' language features.
 **********************************************************************/
{
    float u = 1.0e0, comp;
 
    do {
        u *= 0.5e0;
        comp = 1.0e0 + u;
    }
    while( comp != 1.0e0 );
    return( (double)u*2.0e0 );
} /*-------------------- end of function r1mach ------------------------*/
 
int min0( n, a, b, c, d, e, f, g, h, i, j, k, l, m, o, p )
/*
    PURPOSE
        Determine the minimum of the arguments a-p.
 
    INPUT
        n       Number of arguments to check 1 <= n <= 15.
        a-p     Integer arguments of which the minumum is desired.
 
    RETURNS
        min0    Minimum of a thru p.
*/
int n, a, b, c, d, e, f, g, h, i, j, k, l, m, o, p;
{
    int mt;
 
    if( n < 1 || n > 15 ) return( -1 );
    mt = a;
    if( n == 1 ) return( mt );
 
    if( mt > b ) mt = b;
    if( n == 2 ) return( mt );
 
    if( mt > c ) mt = c;
    if( n == 3 ) return( mt );
 
    if( mt > d ) mt = d;
    if( n == 4 ) return( mt );
 
    if( mt > e ) mt = e;
    if( n == 5 ) return( mt );
 
    if( mt > f ) mt = f;
    if( n == 6 ) return( mt );
 
    if( mt > g ) mt = g;
    if( n == 7 ) return( mt );
 
    if( mt > h ) mt = h;
    if( n == 8 ) return( mt );
 
    if( mt > i ) mt = i;
    if( n == 9 ) return( mt );
 
    if( mt > j  ) mt = j;
    if( n == 10 ) return( mt );
 
    if( mt > k  ) mt = k;
    if( n == 11 ) return( mt );
 
    if( mt > l  ) mt = l;
    if( n == 12 ) return( mt );
 
    if( mt > m ) mt = m;
    if( n == 13 ) return( mt );
 
    if( mt > o  ) mt = o;
    if( n == 14 ) return( mt );
 
    if( mt > p  ) mt = p;
    return( mt );
}
void sscal( n, sa, sx, incx )
int     n, incx;
double  sa;
float   *sx;
/*
    PURPOSE
        Scales a vector by a constant.
 
    INPUT
        n    Number of components to scale.
        sa   Scale value (note that this is a double).
        sx   Vector to scale.
        incx Every incx-th element of sx will be scaled.
 
    OUTPUT
        sx   Scaled vector.
*/
{
    register float ssa = (float)sa;
    int i;

#ifdef alliant
    void sscal_();
    float ssaa = (float)sa;

    sscal_( &n, &ssaa, sx, &incx );
    return;
#else
 
    if( n < 1 ) return;

    /* Code for increment not equal to 1.*/
    if( incx != 1 ) {
	if( incx < 0 ) sx += (-n+1)*incx;
	for( i=0; i<n; i++, sx+=incx )
	  *sx *= ssa;
        return;
    }

    /*  Code for unit increment. */
    for( i=0; i<n; i++ ) 
      sx[i] *= ssa;
    return;
#endif
}
void vexopy( n, v, x, y, itype )
int	n, itype;
float	*v, *x, *y;
/*
  Purpose:
    To operate on the vectors x and y.

  Input:
    n		Number of elements to scale.
    x		First operand vector.
    y		Second operand vector.
    itype	Type of operation to perform:
		itype = 1 => '+'
		itype = 2 => '-'

  Output:
    v		Result vector of x op y.
*/
{
    register int i;

    if( n<1 ) return;

    if( itype == 1 )		/* ADDITION. */
      for( i=0; i<n; i++ )
	v[i] = x[i] + y[i];
    else			/* SUBTRACTION. */
      for( i=0; i<n; i++ )
	v[i] = x[i] - y[i];
}
void vfill( n, v, val )
int	n;
float	*v;
double val;
/*
  Purpose
    To fill the FLOAT vector v with the value val.
    Make sure that if you pass a value for val that you cast
    it to double, viz.,
		vfill( 100, vector, (double)0.0 );
*/
{
    register int i;
    register float vval = (float)val;

    if( n<1 ) return;
    for( i=0; i<n; i++ ) 
      v[i] = vval;
}
#ifdef COPYSIGN
double copysign( x, y )
double x, y;
/*
  PURPOSE
     This routine is recomended by the IEEE standard 754.  Most C 
     3M libraries contain this function.  If yours does not then
     compile this file with the -DCOPYSIGN option and the following
     code will be generated.  This routine copies the sign from y
     onto x.

  INPUT
     x		Thing to have it's sign changed.
     y		Variable to supply the sign.

  RETURNS
     copysign	Returns x with the sign of y (the fortran intrinsic
     		sign(x,y)).
*/
{
    /* Registers are used for reduced memory traffic. */
    register double xx = x;
    register double yy = y;

    if( xx*yy < 0.0 ) return( -xx );
    else return( xx );
}
#endif
SHAR_EOF
if test 14706 -ne "`wc -c < 'blas.c'`"
then
	echo shar: "error transmitting 'blas.c'" '(should have been 14706 characters)'
fi
fi
if test -f 'sample.out'
then
	echo shar: "will not over-write existing file 'sample.out'"
else
cat << \SHAR_EOF > 'sample.out'


	Linpack Benchmark in C of size 50
SGEFA time =    2.400000e-01 (Sec) =    3.680556e-01 (Mflops)
SGECO time =    3.600000e-01 (Sec) =    2.870370e-01 (Mflops)
SGESL time =    2.000000e-02 (Sec) =    2.500000e-01 (Mflops)
 1/COND(A) (Condition number of A) =    6.184206e-04
   ||x-X||/||X||  (Solution Error) =    6.839908e-06
   ||b-Ax||/||X|| (Residual Error) =    3.501223e-06
	Where X = True solution and and x = computed solution


	Linpack Benchmark in C of size 100
SGEFA time =    1.760000e+00 (Sec) =    3.901515e-01 (Mflops)
SGECO time =    2.200000e+00 (Sec) =    3.393939e-01 (Mflops)
SGESL time =    8.000000e-02 (Sec) =    2.500000e-01 (Mflops)
 1/COND(A) (Condition number of A) =    1.554807e-04
   ||x-X||/||X||  (Solution Error) =    1.603013e-05
   ||b-Ax||/||X|| (Residual Error) =    1.260402e-05
	Where X = True solution and and x = computed solution


	Linpack Benchmark in C of size 150
SGEFA time =    5.800000e+00 (Sec) =    3.956897e-01 (Mflops)
SGECO time =    6.680000e+00 (Sec) =    3.637725e-01 (Mflops)
SGESL time =    1.400000e-01 (Sec) =    3.214286e-01 (Mflops)
 1/COND(A) (Condition number of A) =    1.645726e-03
   ||x-X||/||X||  (Solution Error) =    1.207713e-05
   ||b-Ax||/||X|| (Residual Error) =    2.412636e-05
	Where X = True solution and and x = computed solution


	Linpack Benchmark in C of size 200
SGEFA time =    1.346000e+01 (Sec) =    4.021793e-01 (Mflops)
SGECO time =    1.522000e+01 (Sec) =    3.714411e-01 (Mflops)
SGESL time =    2.600000e-01 (Sec) =    3.076923e-01 (Mflops)
 1/COND(A) (Condition number of A) =    1.040173e-04
   ||x-X||/||X||  (Solution Error) =    1.950866e-05
   ||b-Ax||/||X|| (Residual Error) =    3.119719e-05
	Where X = True solution and and x = computed solution


	Linpack Benchmark in C of size 250
SGEFA time =    2.634000e+01 (Sec) =    4.002151e-01 (Mflops)
SGECO time =    2.850000e+01 (Sec) =    3.830409e-01 (Mflops)
SGESL time =    4.000000e-01 (Sec) =    3.125000e-01 (Mflops)
 1/COND(A) (Condition number of A) =    2.600729e-04
   ||x-X||/||X||  (Solution Error) =    1.682144e-05
   ||b-Ax||/||X|| (Residual Error) =    5.737117e-05
	Where X = True solution and and x = computed solution


	Linpack Benchmark in C of size 300
SGEFA time =    4.620000e+01 (Sec) =    3.935065e-01 (Mflops)
SGECO time =    5.242000e+01 (Sec) =    3.571156e-01 (Mflops)
SGESL time =    5.600000e-01 (Sec) =    3.214286e-01 (Mflops)
 1/COND(A) (Condition number of A) =    1.293807e-04
   ||x-X||/||X||  (Solution Error) =    3.347697e-05
   ||b-Ax||/||X|| (Residual Error) =    9.074241e-05
	Where X = True solution and and x = computed solution
SHAR_EOF
if test 2627 -ne "`wc -c < 'sample.out'`"
then
	echo shar: "error transmitting 'sample.out'" '(should have been 2627 characters)'
fi
fi
exit 0
#	End of shell archive