# This is a shell archive.  Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by lll-zeus!seager on Tue Sep 30 12:35:10 PDT 1986
# Contents:  READ.ME makefile driver.c sgefa.c sgesl.c ge.h blas.c sgefat.out
 
echo x - READ.ME
sed 's/^@//' > "READ.ME" <<'@//E*O*F 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.

	SGEFA is a   gaussian  elimination (with  partial    pivoting)
package written in C (single precision).  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 driver.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 sgefat.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.Arpa
@//E*O*F READ.ME//
chmod u=rw,g=r,o=r READ.ME
 
echo x - makefile
sed 's/^@//' > "makefile" <<'@//E*O*F makefile//'

#	This make file peices together the SGEFA.C test routine(s)

CFLAGS=-g
LFLAGS=-lm
OBJS = driver.o sgefa.o sgesl.o blas.o

sgefat: $(OBJS)
	cc $(CFLAGS) -o sgefat $(OBJS) $(LFLAGS)

driver.o: driver.c ge.h

sgefa.o: sgefa.c ge.h

sgesl.o: sgesl.c ge.h

blas.o:	blas.c
@//E*O*F makefile//
chmod u=rw,g=r,o=r makefile
 
echo x - driver.c
sed 's/^@//' > "driver.c" <<'@//E*O*F driver.c//'
  /***************************************************************
 ******************************************************************
****		      Test routine for SGEFA.C			****
 ******************************************************************
  ****************************************************************/

#include <stdio.h>
#include <math.h>
#include "ge.h"

#define BIG	 1.0e38		/* Largest representable float. */
#define SMALL	 1.0e-45	/* Smallest (in absolute value) representable float. */
#define MATPRINT 7		/* Matricies of size <= MATPRINT are printed out. */
#define VECPRINT 7	    	/* Vectors of size <= VECPRINT are printed out. */
#define SCALE	 1		/* System sizes are scaled by this amount. */
main()
{
  register int i, j, k;
  struct FULL a;

  /* Storage for rhs, rhs of trans(A)x and solution. */
  float *b, *bt, *x, anorm, *col, t;
  double err, snrm2();
  int 	*ipvt, retval, test_case = 0;
  void  matdump(), ivecdump(), fvecdump();
  
  /* Loop over the test cases. */
  /* Initilize matrix. */
  while( !matgen( &a, &x, &b, &bt, &ipvt, ++test_case, SCALE ) ) {
  
      /* Check that system size is reasonable. */
      if( a.rd > MAXCOL || a.cd > MAXCOL ) {
	  fprintf(stderr,"Matrix row dim (%d) or column dim (%d) too big.\n",a.rd,a.cd);
	  exit( 1 );
      }
      /* Calculate the 1 norm of A. */
      for( j=0, anorm=0.0; j<a.cd; j++ ) {
	  for( i=0, col=a.pd[j], t=0.0; i<a.rd; i++, col++ ) 
	    t += (*col<0.0 ? -*col : *col );
	  anorm = ( anorm > t ? anorm : t );
      }
      printf("One-Norm(A) ---------- %e.\n", anorm );


    /* Test SGEFA. */
    retval = sgefa( &a, ipvt );

    /* For a successful return from SGEFA test SGESL. */
    if( retval ) 
      printf("Zero Column %d found \n", retval );
    else {
      /* Solve system. */
      (void)sgesl( &a, ipvt, b, 0 );
      if( a.rd <= MATPRINT )
	  (void) matdump( &a, "FACTORED MATRIX FOLLOWS" );
      if( a.rd <= VECPRINT ) {
	  (void)fvecdump( x, a.rd, "True Solution");
	  (void)fvecdump( b, a.rd, "Solution");
      }
      (void)vexopy( a.rd, b, x, b, 2 );
      err = snrm2( a.rd, b, 1 );
      printf(" For Ax=b.    Absolute error = %e.  Relative error = %e.\n",
	     err, err/snrm2( a.rd, x, 1 ) );

      /* Solve transpose system. */
      (void)sgesl( &a, ipvt, bt, 1 );
      if( a.rd <= VECPRINT ) {
	  (void)fvecdump( bt, a.rd, "Solution to transposed system");
      }
      (void)vexopy( a.rd, bt, x, bt, 2 );
      err = snrm2( a.rd, bt, 1 );
      printf(" For A^Tx=b.  Absolute error = %e.  Relative error = %e.\n",
	     err, err/snrm2( a.rd, x, 1 ) );
    }
  }				/* End of while loop over test cases. */
}				/* End of MAIN */

int matgen( a, x, b, bt, ipvt, test_case, scale )
struct FULL *a;
float	    **x, **b, **bt;
int	    **ipvt, test_case, scale;
/*
  This routine generates test matrices for the SGE routines.
  
  INPUT
 test_case	Switch to type of matrix to generate.
 		1, 2, 3 => Hilbert slices of various sizes.  Note
		due to the memory allocalion local to this routine
		test_case=1 must be run before any of the others.
		4, 5    => monoelemental test.
  
 scale          Sizes of systems are scaled according to scale.

  OUTPUT
  a		The matrix stored in the FULL structure.
  x		Generated solution.
  b		Generated right hand side.
  bt		Generated rhs of trans(A)x.
  ipvt		A pointer to an array of ints.

  RETURN VALUE
  0 => everything went O.K.
*/
{
  register int i, j;
  int  n;
  float *col, tl, tu;
  char *malloc();
  void free(), matdump(), ivecdump(), fvecdump();

  if( test_case>1 ) {			/* Free up memory used in the last test. */
    printf("\n\n**********************************************************************\n");
    for( j=0; j<a->rd; j++ )
      free( a->pd[j] );
    free( *x );
    free( *b );
    free( *bt );
  }

  /* Determine the test to generate. */
  switch( test_case ) {
  case 1:				/* Hilbert slice of various sizes. */
  case 2:
  case 3:
    n = 3*test_case*scale;		/* Set system size. */
    a->cd = n;				/* Set column and row dimensions. */
    a->rd = n;
  
    if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */
      return( 1 );

    printf("Hilbert Slice.  Test case %d of size %d.\n", test_case, n );
    for( j=0; j<n; j++ ) {
      for( i=0, col=a->pd[j]; i<n; i++, col++ ) {
	*col = 0.0;
	if( i>=(j-3) && i<=(j+2) ) *col = 1.0/(float)(i+j+1);
      }
    }
    break;

  case 4:				/* Monoemenental test (NOT SCALED). */
  case 5:
    n=1;
    a->cd = n;				/* Set column and row dimensions. */
    a->rd = n;
  
    if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */
      return( 1 );

    printf("Monoelemental.  Test case %d of size %d.\n", test_case, n );
    *a->pd[0] = ( test_case == 4 ? 3.0 : 0.0 ); 
    break;

  case 6:				/* Tridiagional of various types. */
  case 7:
  case 8:
    n = 15*scale;
    a->cd = n;				/* Set column and row dimensions. */
    a->rd = n;
  
    if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */
      return( 1 );

    printf("Tridiagional.  Test case %d of size %d.\n", test_case, n );
    tu = 1.0;
    tl = 1.0;
    if( test_case == 7 ) tl = 100.0;
    if( test_case == 8 ) tu = 100.0;
    for( j=0; j<n; j++ ) {
      for( i=0, col=a->pd[j]; i<n; i++, col++ ) {
	*col = 0.0;
	if( i==j ) *col = 4.0;
	else if( i==j-1 ) *col = tl;
	else if( i==j+1 ) *col = tu;
      }
    }
    break;

  case 9:				/* Rank one. */
    n = 5*scale;
    a->cd = n;				/* Set column and row dimensions. */
    a->rd = n;
  
    if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */
      return( 1 );

    printf("Rank One.  Test case %d of size %d.\n", test_case, n );
    for( j=0; j<n; j++ ) {
      for( i=0, col=a->pd[j]; i<n; i++, col++ ) {
	*col = (float) pow( 10.0, (double)(i-j) );
      }
    }
    break;

  case 10:				/* Zero column. */
    n = 4*scale;
    a->cd = n;				/* Set column and row dimensions. */
    a->rd = n;
  
    if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */
      return( 1 );

    printf("Zero Column.  Test case %d of size %d.\n", test_case, n );
    for( j=0; j<n; j++ ) {
      for( i=0, col=a->pd[j]; i<n; i++, col++ ) {
	tu   = (float)(j-2);
	tl   = (float)(i+1);
	*col = tu/tl;
      }
    }
    break;

  case 11:				/* Upper Triangular. */
    n = 6*scale;
    a->cd = n;				/* Set column and row dimensions. */
    a->rd = n;
  
    if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */
      return( 1 );

    printf("Upper Triangular.  Test case %d of size %d.\n", test_case, n );
    for( j=0; j<n; j++ ) 
      for( i=0, col=a->pd[j]; i<n; i++, col++ ) 
	*col = ( i>j ? 0.0 : (float)(i-j+1) );
    break;

  case 12:				/* Lower Triangular. */
    n = 6*scale;
    a->cd = n;				/* Set column and row dimensions. */
    a->rd = n;
  
    if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */
      return( 1 );

    printf("Lower Triangular.  Test case %d of size %d.\n", test_case, n );
    for( j=0; j<n; j++ ) 
      for( i=0, col=a->pd[j]; i<n; i++, col++ ) 
	*col = ( i<j ? 0.0 : (float)(i-j+1) );
    break;

  case 13:				/* Near Overflow. */
    n = 5*scale;
    a->cd = n;				/* Set column and row dimensions. */
    a->rd = n;
  
    if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */
      return( 1 );

    printf("Near Overflow.  Test case %d of size %d. BIG = %e\n", test_case, n, BIG );
    tl = (float)(n*n);
    for( j=0; j<n; j++ ) 
      for( i=0, col=a->pd[j]; i<n; i++, col++ ) {
	tu = (float)(j+1)/(float)( i>j ? i+1 : j+1 );		/* A number <= 1.0 */
	*col = BIG*tu/tl;
      }
    break;

  case 14:				/* Near Underflow. */
    n = 5*scale;
    a->cd = n;				/* Set column and row dimensions. */
    a->rd = n;
  
    if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */
      return( 1 );

    /* Assumes that BIG < 1/SMALL */
    printf("Near Underflow.  Test case %d of size %d. SMALL = %e\n", test_case, n, 1.0/BIG );
    tl = (float)(n*n);
    for( j=0; j<n; j++ ) 
      for( i=0, col=a->pd[j]; i<n; i++, col++ ) {
	tu = (float)( i>j ? i+1 : j+1 )/(float)(j+1);	/* A number >= 1.0 */
	*col = tu*tl/BIG;
      }
    break;

  default:
    printf("MATGEN: All tests complete.\n");
    return( 1 );
    break;

  }

  /* Generate solution. */
  **x = 1.0;
  if( n>1 ) **(x+1) = 0.0;
  if( n>2 ) {
    for( i=2, col=(*x)+2; i<n; i++, col++ )
      *col = - *(col-2);
  }

  /* Generate rhs. */
  if( matvec( a, *x, *b, 0 ) ) {
    printf("MATGEN: Error in matvec.\n");
    return( 1 );
  }

  /* Generate rhs of transposed system. */
  if( matvec( a, *x, *bt, 1 ) ) {
    printf("MATGEN: Error in matvec.\n");
    return( 1 );
  }
  if( n<=MATPRINT ) 
    (void) matdump( a, "MATRIX FOLLOWS" );
  if( n<=VECPRINT ) {
    (void) fvecdump( *x,  n, "SOLUTION" );
    (void) fvecdump( *b,  n, "RIGHT HAND SIDE" );
    (void) fvecdump( *bt, n, "TRANSPOSE RIGHT HAND SIDE" );
  }
  return( 0 );
}				/* End of MATGEN */

int get_space( a, x, b, bt, ipvt )
struct FULL *a;
float       **x, **b, **bt;
int         **ipvt;
/*
  This routine gets space for the above vectors.
*/
{
  char *malloc();
  register int i,j;

  /* Get space for the columns. */
  for( j=0; j<a->rd; j++ ) {
    a->pd[j] = (float *)malloc( a->cd*sizeof( float ) );
    if( a->pd[j] == NULL ) {
      printf("GET_SPACE: Can't get enouph space for matricies...\n");
      return( 1 );
    }
  }
  *x    = (float *)malloc( a->cd*sizeof( float ) );
  *b    = (float *)malloc( a->cd*sizeof( float ) );
  *bt   = (float *)malloc( a->cd*sizeof( float ) );
  *ipvt = (int *)malloc( a->cd*sizeof( int ) );
  if( *x == NULL || *b == NULL || *bt == NULL || *ipvt == NULL) {
      printf("GET_SPACE: Can't get enouph space for vectors...\n");
      return( 1 );
    }
  return( 0 );
}

int matvec( a, x, b, job)
struct FULL *a;
float	    *x, *b;
int	    job;
/*
  This routine calculates b = a*x (if job=0 and b=trans(a)x else)
  via a column oriented approach.  It is most efficient if the matrix 
  is stored by columns.   a is a matrix (not its transpose, even when 
  job is nonzero) and x, b are vectors of the appropriate size.

  RETURNS
  Nonzero if something goes wrong.
*/
{
  register int i, j;
  float *px, *pb, *col, *row;
  
  /* Check input. */
  if( (a->cd < 1) || (a->rd < 1) ) return( 1 );

  /* Job non-zero => do b = trans(A)x. */
  if( job ) {
    /* Loop over (transposed columns) rows. */
    for( i=0, pb=b; i<a->rd; i++, pb++ ) {
      for( j=0, row=a->pd[i], px=x, *pb=0.0; j<a->cd; j++, px++, row++ )
	*pb += (*row)*(*px);
    }
    return( 0 );
  }

  /* Job zero => do b = Ax. */
  /* Loop over columns. */
  for( i=0, px=x, pb=b, col=a->pd[0]; i<a->rd; i++, pb++, col++)
    *pb = (*col)*(*px);
  for( j=1; j<a->cd; j++ ) {
    for( i=0, px=x+j, pb=b, col=a->pd[j]; i<a->rd; i++, pb++, col++)
      *pb += (*col)*(*px);
  }
  return( 0 );
}				/* End of MATVEC */

void matdump( a, head )
struct FULL *a;
char	    *head;
/*
  This routine prints out a FULL matrix with headding head.
*/
{
  register int i, j;
  int k, jj, ncolmod, ncolrem;

  ncolmod = (a->cd)/6;
  ncolrem = a->cd - ncolmod*6;

  printf("%s\n", head );
  for( i=0;  i<a->rd;  i++) {
    printf("%3d|", i );
    j = 0;
    for( k=0; k<ncolmod; k++ ) {
      if( k ) printf("    ");
      for( jj=0;  jj<6;  jj++, j++) 
	printf("%12.4e", pelem(a,i,j));
      printf("\n");
    }
    
    /* Clean up remainder of this row. */
    for( jj=0;  jj<ncolrem;  jj++, j++) 
      printf("%12.4e", pelem(a,i,j));
    printf("\n");
    }
  printf("\n");
}				/* End of MATDUMP */

void fvecdump( vec, len, head )
float *vec;
int   len;
char  *head;
/*
  This routine prints out a float vector `vec' of length `len' and
  a headding of `head'.
*/
{
  register int i, j, count;
  int lenmod, lenrem;

  lenmod = len/6;
  lenrem = len - lenmod*6;

  printf("%s\n", head );
  count = 0;
  for( j=0; j<lenmod; j++ ) {
    printf("%3d|", count );
    for( i=0; i<6; i++, vec++, count++ )
      printf("%12.4e", *vec );
    printf("\n");
  }
  /* Clean up loop. */
  if( lenrem ) {
    printf("%3d|", count );
    for( i=0; i<lenrem; i++, vec++ )
      printf("%12.4e", *vec );
    printf("\n");
  }
  printf("\n");
}				/* End of FVECDUMP */

void ivecdump( vec, len, head )
int  *vec;
int  len;
char *head;
/*
  This routine prints out a float vector `vec' of length `len' and
  a headding of `head'.
*/
{
  register int i, j, count;
  int lenmod, lenrem;

  lenmod = len/9;
  lenrem = len - lenmod*6;

  printf("%s\n", head );
  count = 0;
  for( j=0; j<lenmod; j++ ) {
    printf("%3d|", count );
    for( i=0; i<9; i++, vec++, count++ )
      printf("%8d", *vec );
    printf("\n");
  }
  /* Clean up loop. */
  if( lenrem ) {
    printf("%3d|", count );
    for( i=0; i<lenrem; i++, vec++ )
      printf("%8d", *vec );
    printf("\n");
  }
  printf("\n");
}
@//E*O*F driver.c//
chmod u=rw,g=r,o=r driver.c
 
echo x - sgefa.c
sed 's/^@//' > "sgefa.c" <<'@//E*O*F sgefa.c//'
  /***************************************************************
 ******************************************************************
****	    Gaussian Elimination with partial pivoting.		****
****	This file contains the factorization routine SGEFA	****
 ******************************************************************
  ****************************************************************/

#include "ge.h"
static char SGEFASid[] = "@(#)sgefa.c	1.1  2/4/86";

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		t, *akk, *alk, *aij, *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 = akk+1; i<n; i++, mik++ )
      *mik *= t;

    /* Column elimination with row indexing. */
    for( j=k+1; j<n; j++ ) {

      /* Interchange a(k,j) and a(l,j) if necessary. */
      t = pelem(a,k,j);
      if( l != k ) {
	pelem(a,k,j) = pelem(a,l,j);
	pelem(a,l,j) = t;
	t = pelem(a,k,j);
      }
      for( i=k+1, aij=a->pd[j]+k+1, mik=akk+1; i<n; i++, aij++, mik++ ) 
	*aij += t*(*mik);
    }
  }				/* End of for k loop */
  
 CLEAN_UP: 
  *ipvt = nm1;
  if( *akk == 0.0e0 ) info = n;
  return( info );
}
@//E*O*F sgefa.c//
chmod u=r,g=r,o=r sgefa.c
 
echo x - sgesl.c
sed 's/^@//' > "sgesl.c" <<'@//E*O*F sgesl.c//'
  /***************************************************************
 ******************************************************************
****	    Gaussian Elimination with partial pivoting.		****
****	This file contains the solution routine SGESL		****
 ******************************************************************
  ****************************************************************/

#include "ge.h"
static char SGESLSid[] = "@(#)sgesl.c	1.1  2/4/86";

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 .
*/
{
  float	t;
  float	*akk, *mik, *uik, *bi;
  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++, ipvt++ ) {
      akk = a->pd[k] + k;		/* akk points to a(k,k). */

      /* Interchange b[k] and b[l] if necessary. */
      l = *ipvt;
      t = b[l];
      if( l != k ) {
	b[l] = b[k];
	b[k] = t;
      }
      for( i=k+1, mik=akk+1; i<n; i++, mik++ )
	b[i] += (*mik)*t;
    }

    /* Back substitution.  Solve  u*x = y. */
    for( k=nm1; k>=0; k-- ) {
      akk = a->pd[k] +k;
      b[k] /= (*akk);
      for( i=0, uik=a->pd[k]; i<k; i++, uik++ )
	b[i] -= (*uik)*b[k];
    }
    return;
  }

  /* job = nonzero.  Solve  trans(a) * x = b. */
  /* First solve trans(u)*y = b. */
  for( k=0; k<n; k++ ) {
    akk = a->pd[k] + k;
    for( i=0, t=0.0, uik=a->pd[k], bi=b; i<k; i++, uik++, bi++ )
      t += (*uik)*(*bi);
    b[k] = (b[k] - t) / (*akk);
  }

  /* b now contains y. */
  /* Solve trans(l)*x = y. */
  ipvt += n-2;
  for( k=n-2; k>=0; k--, ipvt-- ) {
    for( i=k+1, t=0.0, mik=a->pd[k]+k+1, bi=b+k+1; i<n; i++, mik++, bi++ )
      t += (*mik)*(*bi);
    b[k] += t;
    
    /* Interchange b(k) and b(ipvt(k)) if necessary. */
    l    = *ipvt;
    if( l == k ) continue; 
    t    = b[l];
    b[l] = b[k];
    b[k] = t;
  }
  return;
}
@//E*O*F sgesl.c//
chmod u=r,g=r,o=r sgesl.c
 
echo x - ge.h
sed 's/^@//' > "ge.h" <<'@//E*O*F 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)))
@//E*O*F ge.h//
chmod u=r,g=r,o=r ge.h
 
echo x - blas.c
sed 's/^@//' > "blas.c" <<'@//E*O*F blas.c//'
  /***************************************************************
  *****************************************************************
 *******************************************************************
*****								*****
*****				BLAS				*****
*****		Basic Linear Algebra Subroutines		*****
*****	      Written in the C Programming Language.		*****
*****								*****
*****	Functions include:					*****
*****	isamax, saxpy, saxpyx, scopy, sdot, snrm2,		****
*****	vexopy, vfill						*****
*****								*****
 *******************************************************************
  *****************************************************************
   ***************************************************************/
#include <math.h>
#ifndef SMALLsp
#define SMALLsp		1.0e-45		/* Smallest (pos) binary float */
#endif
#ifndef HUGEsp
#define HUGEsp		1.0e+38		/* Largest binary float */
#endif
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.

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

#ifndef abs
#define abs(x) ((x)>0.0?(x):-(x)) 
#endif
  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 );
}
 
void saxpy( n, sa, sx, incx, sy, incy )
float *sx, *sy, 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.
    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;

  if( n<=0 || sa==0.0 ) return;
  if( incx == incy ) {
    if( incx == 1 ) {
      /* Both increments = 1 */
      for( i=0; i<n; i++,sy++,sx++ )
	*sy += sa*(*sx);
      return;
    }
    if( incx>0 ) {
      /* Equal, positive, non-unit increments. */
      for( i=0; i<n; i++,sx+=incx,sy+=incx )
	*sy += sa*(*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 += sa*(*sx);
}
 
void saxpyx( n, sa, sx, incx, sy, incy )
float *sx, *sy, 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.
    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;

  if( n<=0 || sa==0.0 ) return;
  if( incx == incy ) {
    if( incx == 1 ) {
      /* Both increments = 1 */
      for( i=0; i<n; i++, sx++, sy++ )
	*sx = *sy + sa*(*sx);
      return;
    }
    if( incx>0 ) {
      /* Equal, positive, non-unit increments. */
      for( i=0; i<n; i++, sx+=incx, sy+=incx)
	*sx = *sy + sa*(*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 + sa*(*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;

  if( n<1  ) return;
  if( incx == incy ) {
    if( incx == 1 ) {
      /* Both increments = 1 */
      for( i=0; i<n; i++ )
	*(sy++) = *(sx++);
      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;
}
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 i;
  float	stemp = 0.0e0;

  if( n<1 ) return( stemp );
  if( incx == incy ) {
    if( incx == 1 ) {
      /* Both increments = 1 */
      for( i=0; i<n; i++, sx++, sy++ )
	stemp += (*sx)*(*sy);
      return( stemp );
    }
    if( incx>0 ) {
      /* Equal, positive, non-unit increments. */
      for( i=0; i<n; i++, sx+=incx, sy+=incx)
	stemp += (*sx)*(*sy);
      return( 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( 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;
  int	 phase = 3;
  double sum = 0.0e0, cutlo, cuthi, hitst, r1mach();
  float xmax;

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

  cutlo = sqrt( SMALLsp/r1mach() );			/* Calculate near underflow */
  cuthi = sqrt( HUGEsp );				/* Calculate near  overflow */
  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 );
}
int sscal( n, sa, sx, incx )
int     n, incx;
float  sa, *sx;
/*
    PURPOSE
        Scales a vector by a constant.
 
    INPUT
        n    Number of components to scale.
        sa   Scale value.
        sx   Vector to scale.
        incx Every incx-th element of sx will be scaled.
 
    OUTPUT
        sx   Scaled vector.
*/
{
  int i;
 
  if( n < 1 ) return( 1 );

  /* 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 *= sa;
        return;
  }
  /*  Code for unit increment. */
  for( i=0; i<n; i++, sx++ )
    *sx *= sa;
  return;
}
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++, x++, y++, v++ )
      *v = *x + *y;
  else			/* SUBTRACTION. */
    for( i=0; i<n; i++, x++, y++, v++ )
      *v = *x - *y;
}
void vfill( n, v, val )
int	n;
float	*v, val;
/*
  Purpose
    To fill the FLOAT vector v with the FLOAT value val.
    Make sure that if you pass a value for val that you cast
    it to float, viz.,
		vfill( 100, vector, (float)0.0 );
*/
{
  register int i;

  if( n<1 ) return;
  for( i=0; i<n; i++, v++ )
    *v = val;
}
@//E*O*F blas.c//
chmod u=rw,g=r,o=r blas.c
 
echo x - sgefat.out
sed 's/^@//' > "sgefat.out" <<'@//E*O*F sgefat.out//'
Hilbert Slice.  Test case 1 of size 3.
MATRIX FOLLOWS
  0|  1.0000e+00  5.0000e-01  3.3333e-01
  1|  5.0000e-01  3.3333e-01  2.5000e-01
  2|  3.3333e-01  2.5000e-01  2.0000e-01

SOLUTION
  0|  1.0000e+00  0.0000e+00 -1.0000e+00

RIGHT HAND SIDE
  0|  6.6667e-01  2.5000e-01  1.3333e-01

TRANSPOSE RIGHT HAND SIDE
  0|  6.6667e-01  2.5000e-01  1.3333e-01

One-Norm(A) ---------- 1.833333e+00.
FACTORED MATRIX FOLLOWS
  0|  1.0000e+00  5.0000e-01  3.3333e-01
  1| -5.0000e-01  8.3333e-02  8.3333e-02
  2| -3.3333e-01 -1.0000e+00  5.5556e-03

True Solution
  0|  1.0000e+00  0.0000e+00 -1.0000e+00

Solution
  0|  1.0000e+00  1.7881e-07 -1.0000e+00

 For Ax=b.    Absolute error = 1.884864e-07.  Relative error = 1.332800e-07.
Solution to transposed system
  0|  1.0000e+00  1.5497e-06 -1.0000e+00

 For A^Tx=b.  Absolute error = 2.061321e-06.  Relative error = 1.457574e-06.


**********************************************************************
Hilbert Slice.  Test case 2 of size 6.
MATRIX FOLLOWS
  0|  1.0000e+00  5.0000e-01  3.3333e-01  2.5000e-01  0.0000e+00  0.0000e+00

  1|  5.0000e-01  3.3333e-01  2.5000e-01  2.0000e-01  1.6667e-01  0.0000e+00

  2|  3.3333e-01  2.5000e-01  2.0000e-01  1.6667e-01  1.4286e-01  1.2500e-01

  3|  0.0000e+00  2.0000e-01  1.6667e-01  1.4286e-01  1.2500e-01  1.1111e-01

  4|  0.0000e+00  0.0000e+00  1.4286e-01  1.2500e-01  1.1111e-01  1.0000e-01

  5|  0.0000e+00  0.0000e+00  0.0000e+00  1.1111e-01  1.0000e-01  9.0909e-02


SOLUTION
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

RIGHT HAND SIDE
  0|  6.6667e-01  4.1667e-01  2.7619e-01 -4.1667e-02 -3.1746e-02  1.0000e-01

TRANSPOSE RIGHT HAND SIDE
  0|  6.6667e-01  2.5000e-01  2.7619e-01  2.0833e-01 -3.1746e-02 -2.5000e-02

One-Norm(A) ---------- 1.833333e+00.
FACTORED MATRIX FOLLOWS
  0|  1.0000e+00  5.0000e-01  3.3333e-01  2.5000e-01  0.0000e+00  0.0000e+00

  1| -5.0000e-01  2.0000e-01  1.6667e-01  1.4286e-01  1.2500e-01  1.1111e-01

  2| -3.3333e-01 -4.1667e-01  1.4286e-01  1.2500e-01  1.1111e-01  1.0000e-01

  3|  0.0000e+00 -4.1667e-01 -9.7222e-02  1.1111e-01  1.0000e-01  9.0909e-02

  4|  0.0000e+00  0.0000e+00 -1.3611e-01 -6.1161e-02  1.0079e-01 -5.8738e-02

  5|  0.0000e+00  0.0000e+00  0.0000e+00 -2.9911e-02 -6.8989e-01  1.0006e-01


True Solution
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

Solution
  0|  1.0000e+00  0.0000e+00 -1.0000e+00 -6.7055e-08  1.0000e+00 -1.4893e-07

 For Ax=b.    Absolute error = 2.889978e-07.  Relative error = 1.668530e-07.
Solution to transposed system
  0|  1.0000e+00 -5.9605e-08 -1.0000e+00  5.9605e-08  1.0000e+00 -1.5274e-07

 For A^Tx=b.  Absolute error = 1.843548e-07.  Relative error = 1.064373e-07.


**********************************************************************
Hilbert Slice.  Test case 3 of size 9.
One-Norm(A) ---------- 1.833333e+00.
 For Ax=b.    Absolute error = 1.155101e-05.  Relative error = 5.165769e-06.
 For A^Tx=b.  Absolute error = 4.500883e-06.  Relative error = 2.012856e-06.


**********************************************************************
Monoelemental.  Test case 4 of size 1.
MATRIX FOLLOWS
  0|  3.0000e+00

SOLUTION
  0|  1.0000e+00

RIGHT HAND SIDE
  0|  3.0000e+00

TRANSPOSE RIGHT HAND SIDE
  0|  3.0000e+00

One-Norm(A) ---------- 3.000000e+00.
FACTORED MATRIX FOLLOWS
  0|  3.0000e+00

True Solution
  0|  1.0000e+00

Solution
  0|  1.0000e+00

 For Ax=b.    Absolute error = 0.000000e+00.  Relative error = 0.000000e+00.
Solution to transposed system
  0|  1.0000e+00

 For A^Tx=b.  Absolute error = 0.000000e+00.  Relative error = 0.000000e+00.


**********************************************************************
Monoelemental.  Test case 5 of size 1.
MATRIX FOLLOWS
  0|  0.0000e+00

SOLUTION
  0|  1.0000e+00

RIGHT HAND SIDE
  0|  0.0000e+00

TRANSPOSE RIGHT HAND SIDE
  0|  0.0000e+00

One-Norm(A) ---------- 0.000000e+00.
Zero Column 1 found 


**********************************************************************
Tridiagional.  Test case 6 of size 15.
One-Norm(A) ---------- 6.000000e+00.
 For Ax=b.    Absolute error = 0.000000e+00.  Relative error = 0.000000e+00.
 For A^Tx=b.  Absolute error = 0.000000e+00.  Relative error = 0.000000e+00.


**********************************************************************
Tridiagional.  Test case 7 of size 15.
One-Norm(A) ---------- 1.050000e+02.
 For Ax=b.    Absolute error = 0.000000e+00.  Relative error = 0.000000e+00.
 For A^Tx=b.  Absolute error = 0.000000e+00.  Relative error = 0.000000e+00.


**********************************************************************
Tridiagional.  Test case 8 of size 15.
One-Norm(A) ---------- 1.050000e+02.
 For Ax=b.    Absolute error = 2.625374e+07.  Relative error = 9.282099e+06.
 For A^Tx=b.  Absolute error = 1.000835e+00.  Relative error = 3.538487e-01.


**********************************************************************
Rank One.  Test case 9 of size 5.
MATRIX FOLLOWS
  0|  1.0000e+00  1.0000e-01  1.0000e-02  1.0000e-03  1.0000e-04
  1|  1.0000e+01  1.0000e+00  1.0000e-01  1.0000e-02  1.0000e-03
  2|  1.0000e+02  1.0000e+01  1.0000e+00  1.0000e-01  1.0000e-02
  3|  1.0000e+03  1.0000e+02  1.0000e+01  1.0000e+00  1.0000e-01
  4|  1.0000e+04  1.0000e+03  1.0000e+02  1.0000e+01  1.0000e+00

SOLUTION
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00

RIGHT HAND SIDE
  0|  9.9010e-01  9.9010e+00  9.9010e+01  9.9010e+02  9.9010e+03

TRANSPOSE RIGHT HAND SIDE
  0|  9.9010e+03  9.9010e+02  9.9010e+01  9.9010e+00  9.9010e-01

One-Norm(A) ---------- 1.111100e+04.
FACTORED MATRIX FOLLOWS
  0|  1.0000e+04  1.0000e+03  1.0000e+02  1.0000e+01  1.0000e+00
  1| -1.0000e-03  7.6294e-06  9.5367e-07  5.9605e-08  7.4506e-09
  2| -1.0000e-02  0.0000e+00 -9.3132e-10  5.8208e-11 -7.2760e-12
  3| -1.0000e-01 -7.8125e-03  0.0000e+00  7.4506e-09  0.0000e+00
  4| -1.0000e-04 -9.7656e-04  0.0000e+00  6.2500e-02  5.8208e-11

True Solution
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00

Solution
  0| -2.3230e-01 -8.0000e+00 -6.4000e+01  1.0240e+03  1.6384e+04

 For Ax=b.    Absolute error = 1.641509e+04.  Relative error = 9.477259e+03.
Solution to transposed system
  0| -8.1920e+03  0.0000e+00  1.2800e+02  0.0000e+00  5.2930e-01

 For A^Tx=b.  Absolute error = 8.194016e+03.  Relative error = 4.730817e+03.


**********************************************************************
Zero Column.  Test case 10 of size 4.
MATRIX FOLLOWS
  0| -2.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00
  1| -1.0000e+00 -5.0000e-01  0.0000e+00  5.0000e-01
  2| -6.6667e-01 -3.3333e-01  0.0000e+00  3.3333e-01
  3| -5.0000e-01 -2.5000e-01  0.0000e+00  2.5000e-01

SOLUTION
  0|  1.0000e+00  1.0000e+01 -1.0000e+00 -1.0000e+01

RIGHT HAND SIDE
  0| -2.2000e+01 -1.1000e+01 -7.3333e+00 -5.5000e+00

TRANSPOSE RIGHT HAND SIDE
  0| -6.3333e+00 -3.1667e+00  0.0000e+00  3.1667e+00

One-Norm(A) ---------- 4.166667e+00.
Zero Column 4 found 


**********************************************************************
Upper Triangular.  Test case 11 of size 6.
MATRIX FOLLOWS
  0|  1.0000e+00  0.0000e+00 -1.0000e+00 -2.0000e+00 -3.0000e+00 -4.0000e+00

  1|  0.0000e+00  1.0000e+00  0.0000e+00 -1.0000e+00 -2.0000e+00 -3.0000e+00

  2|  0.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00 -1.0000e+00 -2.0000e+00

  3|  0.0000e+00  0.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00 -1.0000e+00

  4|  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

  5|  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  1.0000e+00


SOLUTION
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

RIGHT HAND SIDE
  0| -1.0000e+00 -2.0000e+00 -2.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

TRANSPOSE RIGHT HAND SIDE
  0|  1.0000e+00  0.0000e+00 -2.0000e+00 -2.0000e+00 -1.0000e+00 -2.0000e+00

One-Norm(A) ---------- 1.100000e+01.
FACTORED MATRIX FOLLOWS
  0|  1.0000e+00  0.0000e+00 -1.0000e+00 -2.0000e+00 -3.0000e+00 -4.0000e+00

  1|  0.0000e+00  1.0000e+00  0.0000e+00 -1.0000e+00 -2.0000e+00 -3.0000e+00

  2|  0.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00 -1.0000e+00 -2.0000e+00

  3|  0.0000e+00  0.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00 -1.0000e+00

  4|  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

  5|  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  1.0000e+00


True Solution
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

Solution
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

 For Ax=b.    Absolute error = 0.000000e+00.  Relative error = 0.000000e+00.
Solution to transposed system
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

 For A^Tx=b.  Absolute error = 0.000000e+00.  Relative error = 0.000000e+00.


**********************************************************************
Lower Triangular.  Test case 12 of size 6.
MATRIX FOLLOWS
  0|  1.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00

  1|  2.0000e+00  1.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00

  2|  3.0000e+00  2.0000e+00  1.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00

  3|  4.0000e+00  3.0000e+00  2.0000e+00  1.0000e+00  0.0000e+00  0.0000e+00

  4|  5.0000e+00  4.0000e+00  3.0000e+00  2.0000e+00  1.0000e+00  0.0000e+00

  5|  6.0000e+00  5.0000e+00  4.0000e+00  3.0000e+00  2.0000e+00  1.0000e+00


SOLUTION
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

RIGHT HAND SIDE
  0|  1.0000e+00  2.0000e+00  2.0000e+00  2.0000e+00  3.0000e+00  4.0000e+00

TRANSPOSE RIGHT HAND SIDE
  0|  3.0000e+00  2.0000e+00  2.0000e+00  2.0000e+00  1.0000e+00  0.0000e+00

One-Norm(A) ---------- 2.100000e+01.
FACTORED MATRIX FOLLOWS
  0|  6.0000e+00  5.0000e+00  4.0000e+00  3.0000e+00  2.0000e+00  1.0000e+00

  1| -3.3333e-01 -8.3333e-01 -6.6667e-01 -5.0000e-01 -3.3333e-01 -1.6667e-01

  2| -5.0000e-01 -6.0000e-01 -8.0000e-01 -6.0000e-01 -4.0000e-01 -2.0000e-01

  3| -6.6667e-01 -4.0000e-01 -5.0000e-01 -7.5000e-01 -5.0000e-01 -2.5000e-01

  4| -8.3333e-01 -2.0000e-01 -2.5000e-01 -3.3333e-01 -6.6667e-01 -3.3333e-01

  5| -1.6667e-01 -8.0000e-01 -7.5000e-01 -6.6667e-01 -5.0000e-01 -5.0000e-01


True Solution
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00  0.0000e+00

Solution
  0|  1.0000e+00  0.0000e+00 -1.0000e+00 -3.9736e-08  1.0000e+00  3.5763e-07

 For Ax=b.    Absolute error = 3.647319e-07.  Relative error = 2.105781e-07.
Solution to transposed system
  0|  1.0000e+00  5.9605e-08 -1.0000e+00 -2.9802e-07  1.0000e+00  5.9605e-08

 For A^Tx=b.  Absolute error = 7.609812e-07.  Relative error = 4.393527e-07.


**********************************************************************
Near Overflow.  Test case 13 of size 5. BIG = 1.000000e+38
MATRIX FOLLOWS
  0|  4.0000e+36  4.0000e+36  4.0000e+36  4.0000e+36  4.0000e+36
  1|  2.0000e+36  4.0000e+36  4.0000e+36  4.0000e+36  4.0000e+36
  2|  1.3333e+36  2.6667e+36  4.0000e+36  4.0000e+36  4.0000e+36
  3|  1.0000e+36  2.0000e+36  3.0000e+36  4.0000e+36  4.0000e+36
  4|  8.0000e+35  1.6000e+36  2.4000e+36  3.2000e+36  4.0000e+36

SOLUTION
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00

RIGHT HAND SIDE
  0|  4.0000e+36  2.0000e+36  1.3333e+36  2.0000e+36  2.4000e+36

TRANSPOSE RIGHT HAND SIDE
  0|  3.4667e+36  2.9333e+36  2.4000e+36  3.2000e+36  4.0000e+36

One-Norm(A) ---------- 2.000000e+37.
FACTORED MATRIX FOLLOWS
  0|  4.0000e+36  4.0000e+36  4.0000e+36  4.0000e+36  4.0000e+36
  1| -5.0000e-01  2.0000e+36  2.0000e+36  2.0000e+36  2.0000e+36
  2| -3.3333e-01 -6.6667e-01  1.3333e+36  1.3333e+36  1.3333e+36
  3| -2.5000e-01 -5.0000e-01 -7.5000e-01  1.0000e+36  1.0000e+36
  4| -2.0000e-01 -4.0000e-01 -6.0000e-01 -8.0000e-01  8.0000e+35

True Solution
  0|  1.0000e+00  0.0000e+00 -1.0000e+00  0.0000e+00  1.0000e+00

Solution
  0|  1.0000e+00  7.9228e-08 -1.0000e+00  5.5460e-07  1.0000e+00

 For Ax=b.    Absolute error = 7.571020e-07.  Relative error = 4.371130e-07.
Solution to transposed system
  0|  1.0000e+00 -1.4901e-07 -1.0000e+00  0.0000e+00  1.0000e+00

 For A^Tx=b.  Absolute error = 2.250026e-07.  Relative error = 1.299053e-07.


**********************************************************************
Near Underflow.  Test case 14 of size 5. SMALL = 1.000000e-38
MATRIX FOLLOWS
  0|  2.5000e-37  2.5000e-37  2.5000e-37  2.5000e-37  2.5000e-37
  1|  5.0000e-37  2.5000e-37  2.5000e-37  2.5000e-37  2.5000e-37
  2|  7.5000e-37  3.7500e-37  2.5000e-37  2.5000e-37  2.5000e-37
  3|  1.0000e-36  5.0000e-37  3.3333e-37  2.5000e-37  2.5000e-37
  4|  1.2500e-36  6.2500e-37  4.1667e-37  3.1250e-37  2.5000e-37

SOLUTION
  0|  1.0000e+00  1.4901e-07 -1.0000e+00 -1.4901e-07  1.0000e+00

RIGHT HAND SIDE
  0|  2.5000e-37  5.0000e-37  7.5000e-37  9.1667e-37  1.0833e-36

TRANSPOSE RIGHT HAND SIDE
  0|  7.5000e-37  5.0000e-37  4.1667e-37  3.1250e-37  2.5000e-37

One-Norm(A) ---------- 3.750000e-36.
FACTORED MATRIX FOLLOWS
  0|  1.2500e-36  6.2500e-37  4.1667e-37  3.1250e-37  2.5000e-37
  1| -4.0000e-01  1.2500e-37  1.6667e-37  1.8750e-37  2.0000e-37
  2| -6.0000e-01 -1.7937e-07  8.3333e-38  1.2500e-37  1.5000e-37
  3| -8.0000e-01  0.0000e+00 -2.6905e-07  6.2500e-38  1.0000e-37
  4| -2.0000e-01  0.0000e+00  8.4078e-08  5.3810e-07  5.0000e-38

True Solution
  0|  1.0000e+00  1.4901e-07 -1.0000e+00 -1.4901e-07  1.0000e+00

Solution
  0|  1.0000e+00  1.7937e-07 -1.0000e+00 -3.0492e-06  1.0000e+00

 For Ax=b.    Absolute error = 4.261550e-06.  Relative error = 2.460407e-06.
Solution to transposed system
  0|  1.0000e+00  4.5402e-07 -1.0000e+00  4.4842e-07  1.0000e+00

 For A^Tx=b.  Absolute error = 1.016167e-06.  Relative error = 5.866841e-07.


**********************************************************************
MATGEN: All tests complete.
@//E*O*F sgefat.out//
chmod u=rw,g=r,o=r sgefat.out
 
echo Inspecting for damage in transit...
temp=/tmp/shar$$; dtemp=/tmp/.shar$$
trap "rm -f $temp $dtemp; exit" 0 1 2 3 15
cat > $temp <<\!!!
      51     337    2860 READ.ME
      17      37     269 makefile
     479    2106   12921 driver.c
     107     465    2988 sgefa.c
     103     458    2828 sgesl.c
      42     277    1912 ge.h
     487    1968   11320 blas.c
     410    1374   13778 sgefat.out
    1696    7022   48876 total
!!!
wc  READ.ME makefile driver.c sgefa.c sgesl.c ge.h blas.c sgefat.out | sed 's=[^ ]*/==' | diff -b $temp - >$dtemp
if [ -s $dtemp ]
then echo "Ouch [diff of wc output]:" ; cat $dtemp
else echo "No problems found."
fi
exit 0