/***************************************************************
 ******************************************************************
****		      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. */
#ifndef SCALE
#define SCALE	 1		/* System sizes are scaled by this amount. */
#endif
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 ) **bt = 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");
}