# 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 11:14:49 PDT 1986
# Contents:  READ.ME makefile dcgdrv.c dcg.c dcg.h dbatimes.c dbmsolve.c
#	dcgutil.c dblas.c
 
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 DCG package.

	DCG is a Preconditioned Conjugate Gradient  package written in
C (double precision).  It uses the modular  approach  of LINPACK.  The
necessary BLAS routines are included in blas.c.  This package is  also
set up so that the user can define his/her  own  matrix data structure
in dcg.h and write matrix multiply and matrix  solve routines based on
that   matrix  structure.   Examples of   both   are found in   dcg.h,
dbatimes.c and dbmsolve.c for a banded matrix.  A driver (dcgdrv.c) is
also included and runs several simple minded examples.   Feel  free to
modify all of the above routines for your purposes.

	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//'
#	SID = "@(#)makefile	1.2  9/30/86"
#	This file makes the Preconditioned Conjugate Gradient routine.
#	Also, it makes the non-symmetric system solvers: ORTHOMIN(k)
#	Define DEBUG if it is necessary to get debug output.
#
CFLAGS=-g
LDFLAGS=-lm
SOBJ=scgdrv.o scg.o sbatimes.o sbmsolve.o scgutil.o
DOBJ=dcgdrv.o dcg.o dbatimes.o dbmsolve.o dcgutil.o dblas.o
DOBJN = dnsymdrv.o dorthomin.o dorthatimes.o

# Single precsision (float) version....
scg: $(SOBJ)
	cc $(CFLAGS) -o scg $(SOBJ) $(LDFLAGS) 

scg.o: scg.c scg.h

scgdrv.o: scgdrv.c scg.h

scgutil.o : scgutil.c scg.h

sbatimes.o : sbatimes.c scg.h

sbmsolve.o : sbmsolve.c scg.h

# Double precsision (double) version....
dcg: $(DOBJ)
	cc $(CFLAGS) -o dcg $(DOBJ) $(LDFLAGS) 

dcg.o: dcg.c dcg.h

dcgdrv.o: dcgdrv.c dcg.h

dcgutil.o : dcgutil.c dcg.h

dbatimes.o : dbatimes.c dcg.h

dbmsolve.o : dbmsolve.c dcg.h

# Double precision Non-symmetric solver ORTHOMIN.
dorthomin: $(DOBJN)
	cc $(CFLAGS) -o dorthomin $(DOBJN) $(LDFLAGS)

dnsymdrv.o : dnsymdrv.c dorthomin.h

dorthomin.o : dorthomin.c dorthomin.h

dorthatimes.o : dorthatimes.c dorthomin.h

# Clean up stuff.
clean:
	rm -f scg dcg $(SOBJ) $(DOBJ) $(DOBJN)

@//E*O*F makefile//
chmod u=r,g=r,o=r makefile
 
echo x - dcgdrv.c
sed 's/^@//' > "dcgdrv.c" <<'@//E*O*F dcgdrv.c//'
/*
  This routine generates a test problem for the preconditionded CG routine.
*/
#include <stdio.h>
#include <math.h>
#include "dcg.h"
static char DCGDRVsid[] = "@(#)dcgdrv.c	1.3  9/30/86";

#define MAXMAT  1000	/* Maximum Matrix size. */
#define MATPRINT 26	/* Matricies <= MATPRINT in size are printed. */
#define VECPRINT 26	/* Vectors <= VECPRINT in size are printed. */

main(argc, argv)
int	argc;
char	*argv[];
{
  extern int *(dbatimes)(), *(dbmsolve)();
  int 	 i, j, n, itmax, test_type = 0;
  double  *sol, *x, *b, err, err_est;
  struct DCGAMAT a;
  struct DCGMMAT m;
  double etime(), tary[2];				/* Total cpu time used in SECONDS. */
  double job_time;					/* Time to complete a task. */
  double dnrm2();

  /* Set up structures. */
  if( argc <= 1 ) test_type++;				/* Skip first test if no file is given. */
  while( !matgen( &a, &m, &sol, &x, &b, ++test_type, argv[1] ) ) {
      n = a.rd;
      if( n>MAXMAT || n<1 ) {
	  printf("Error return in matrix size.  n=%d.\n",n);
	  exit( 1 );
      }

      /* Solve the system with Diagional Scaling. */
      itmax   = n;
      m.itype = 1;
      job_time = etime( tary );
      j = dcg( n, dbatimes, &a, dbmsolve, &m, x, b, itmax, 1.0e-5, &err_est );
      job_time = etime( tary ) - job_time;
      printf("\nCG(DS)  Est err = %e.  Took %d it (%e in sec).\n", err_est, j, job_time );
      if( n<VECPRINT ) {
	  printf("CGDRV:     j          x[j],          sol[j],        error\n");
	  for(i=0; i<n; i++ ) {
	      printf("        %4d %15e %15e %15e\n",i,x[i],sol[i],x[i]-sol[i]);
	  }
      }
      
      /* Calculate error. */
      (void)dexopy( n, x, sol, x, 2 );
      err = dnrm2( n, x, 1 );
      printf("Error = %e relative err = %e.\n", err, err/dnrm2(n, sol, 1 ) );

      /* Solve the system with Parameterized Incomplete Inverse with 2 terms. */
      itmax   = n;
      m.itype = 2;
      dfill( n, x, 0.0 );
      
      /* Scale system. */
      dsbscale( n, &a, &m, sol, b, 1, 1 );
      job_time = etime( tary );
      j = dcg( n, dbatimes, &a, dbmsolve, &m, x, b, itmax, 1.0e-5, &err_est );
      job_time = etime( tary ) - job_time;
      printf("\nCG(PII)  Est err = %e.  Took %d it (%e in sec).\n", err_est, j, job_time );
      if( n<VECPRINT ) {
	  printf("CGDRV:     j          x[j],          sol[j],        error\n");
	  for(i=0; i<n; i++ ) {
	      printf("        %4d %15e %15e %15e\n",i,x[i],sol[i],x[i]-sol[i]);
	  }
      }

      /* Calculate error. */
      (void)dexopy( n, x, sol, x, 2 );
      err = dnrm2( n, x, 1 );
      printf("Error = %e relative err = %e.\n", err, err/dnrm2(n, sol, 1 ) );
      
      /* UNScale system. */
      dsbscale( n, &a, &m, sol, b, 0, 1 );
  }				/* End of test matrix loop. */
}				/* End of driver. */

int matgen( aptr, mptr, sol, x, b, how, filenam )
struct DCGAMAT *aptr;
struct DCGMMAT *mptr;
double	    **sol, **x, **b;
int	    how;
char	    *filenam;
{
  char 	 chtmp[81];
  int	 n;				/* System size. */
  int	 dosoln = 1;
  int	 intmp, i, j;
  double  *tmp1, *tmp2;
  double ranf();			/* Random number generator ala FORTRAN. */
  double etime(), tary[2];		/* Total cpu time used in SECONDS. */
  double job_time;			/* Time to complete a task. */
  FILE 	 *fptr, *fopen();
  char *malloc();

  switch( how )	{			/* Type of matrix to generate. */
  case 1:				/* Read in from file. */
      if( (fptr = fopen( filenam, "r" )) == NULL ) {
	  printf("MATGEN: Error opening %s\n",filenam );
	  return( 1 );
      }
      printf("\n\n**************************************************\n");
      printf("Test case 1.  Matrix read in from %s.\n", filenam );
      fscanf(fptr,"%d %d\n", &aptr->hbw, &n );
      mptr->hbw = aptr->hbw;
      aptr->rd = mptr->rd = n;
      printf("MATGEN: half-band-width = %d, System Size = %d\n",aptr->hbw, n);
      if( aptr->hbw >= MAXBND || n < 1 ) {
	  printf("MATGEN: Bad hbw (%d) or n (%d).\n",aptr->hbw, n);
	  return( 1 );
      }
      fgets( chtmp, 81, fptr );
      mptr->itype = 1;
      dosoln = 0;
      *sol = (double *)malloc( (unsigned)n*sizeof( double ) );
      *b   = (double *)malloc( (unsigned)n*sizeof( double ) );
      if( *sol == NULL || *b == NULL ) {
	  printf("MATGEN: Error in geting space for sol, b vectors.\n");
	  return( 1 );
      }
      for( j=0; j<=aptr->hbw; j++ ) {
	  aptr->diagptr[j] = (double *)malloc( (unsigned)n*sizeof( double ) );
	  mptr->diagptr[j] = (double *)malloc( (unsigned)n*sizeof( double ) );
	  if( aptr->diagptr[j] == NULL || mptr->diagptr[j] == NULL ) {
	      printf("MATGEN: Can't get enouph space for matricies...\n");
	      return( 1 );
	  }
      }
      
      for( i=0; i<n; i++ ) {				/* Read in diagionals line by line. */
	  fscanf(fptr,"%d,",&intmp);
	  for( j=0; j<=aptr->hbw; j++ ) 
	    fscanf(fptr,"%le",aptr->diagptr[j]+i);
	  fscanf( fptr, "%le",*sol+i);
	  fscanf( fptr, "%le",*b+i);
/*	  fgets( chtmp, 81, fptr );	*/		/* Clear out rest of line. */
	  *(mptr->diagptr[0]+i) = *(aptr->diagptr[0]+i);	/* Diag scaling pre-conditioning. */
      }
      break;
    case 2:				/* Small Tridiagional Matrix. */
      n = 10;
      printf("\n\n**************************************************\n");
      printf("Test case 2.  Tridiagional system of size %d.\n", n );
      aptr->hbw = 1;
      goto gen_space;
      
    case 3:				/* Large Matrix. */
      n = 400;
      aptr->hbw = 20;
      printf("\n\n**************************************************\n");
      printf("Test case 3.  System with %d upper diagionals and size %d.\n", aptr->hbw, n );
      
    gen_space:
      mptr->hbw   = aptr->hbw;
      mptr->itype = 1;
      aptr->rd    = mptr->rd = n;
      if( (n>MAXMAT || n<2) || (aptr->hbw<1 || aptr->hbw>MAXBND)) {
	  printf("MATGEN: Error in size of matrix.\n");
	  return( 1 );
      }
      for( i=0; i<=aptr->hbw; i++ ) {
	  aptr->diagptr[i] = (double *)malloc( (unsigned)n*sizeof( double ) );
	  mptr->diagptr[i] = (double *)malloc( (unsigned)n*sizeof( double ) );
	  if( aptr->diagptr[i] == NULL || mptr->diagptr[i] == NULL ) {
	      printf("MATGEN: Can't get enouph space for matricies...\n");
	      return( 1 );
	  }
      }
      
      /* Set up Matrix main diagional elements and M (diag scaling). */
      tmp1 = aptr->diagptr[0];
      tmp2 = mptr->diagptr[0];
      for( i=0; i<n; i++, tmp1++, tmp2++ ) {
	  *tmp1 = 1.0 + 100.0*ranf();
	  *tmp2 = *tmp1;
      }
      
      /* Off diagional elements. */
      for( j=1; j<=aptr->hbw; j++ ) {
	  tmp1 = aptr->diagptr[j];
	  for( i=0; i<n-j; i++, tmp1++ ) {
	      *tmp1 = -ranf();
	  }
      }
      break;
      
    default:				/* Unknown type to generate. */
      return( 1 );
  }					/* End of Switch */

  /* Set true solution, get rhs & initial guess. */
  if( (*x   = (double *)malloc( (unsigned)n*sizeof( double ) )) == NULL ) {
      printf("MATGEN: Error in geting space for x vector.\n");
      return( 1 );
  }
  dfill( n, *x, 0.0 );
  if( dosoln ) {
      *sol = (double *)malloc( (unsigned)n*sizeof( double ) );
      *b   = (double *)malloc( (unsigned)n*sizeof( double ) );
      if( *sol == NULL || *b == NULL ) {
	  printf("MATGEN: Error in geting space for vectors.\n");
	  return( 1 );
      }
      for( tmp1 = *sol, i=0; i<n; i++, tmp1++ )
	*tmp1 = (double)i;

      job_time = etime( tary );
      if( dbatimes( n, aptr, *sol, *b ) ) {
	  printf("CGDRV: Error return from dbatimes.\n");
	  exit( 1 );
      }
      job_time = etime( tary ) - job_time;
  }
  printf("Matrix generation complete.  DBATIMES call = %e (sec).\n",job_time);

  /* Echo out generated system if desired. */
  if( n <= MATPRINT ) 
    (void) dsbdumpmat( aptr, "Matrix diagionals follow...");
  if( n <= VECPRINT ) {
      (void) dumpdvec( n, *sol, "Solution follows...");
      (void) dumpdvec( n, *b,   "Rhs follows...");
  }
      
  return( 0 );
}
/**********************************************************************
****
****  This routine generates a uniform psudo-random number in [0,1].
****  It is designed to model the FORTRAN function of the same name.
****
 **********************************************************************/
double ranf()
{
    long	random();			/* random number generator. */
    double	rng = 2147483647.;		/* random num gen range = 2**31-1 */

    return( (double)random()/rng );
}

#include <sys/types.h>
#include <sys/times.h>

double etime( tary )
double tary[2];
/*********************************************************************** 
   This routine returns the elapsed total (system+user) time used by
   this process in SECONDS.

 INPUT
   NONE

 OUTPUT
   tary[2]	This array should be declaired by the user.
		The first element returned is the user time used.
		The second element returned is the system time used.
 RETURNS
 		The total of system and user time in SECONDS.
***********************************************************************/
{
    struct tms buf;			/* Buffer to hold system time. */

    (void)times( &buf );		/* Get times (in HZ=1/60 sec) from the system. */
    tary[0] = (double)buf.tms_utime / 60.0;	/* User time in SECONDS. */
    tary[1] = (double)buf.tms_stime / 60.0;	/* System time in SECONDS. */
    return( tary[0] + tary[1] );		/* Return Total time in SECONDS. */
}

/****************************************************************************
****	    This routine prints the (double) vector vec of lenght	  ****
****		    len to stdout with a heading head.		  	  ****
 ****************************************************************************/
int dumpdvec( len, vec, head )
int   len;
double *vec;
char  *head;
{
    register int i;

    printf("\n%s\n  0>",head );
    for( i=0; i<len; i++ ) {
	printf("%13.6e|",vec[i]);
	if( (i+1) % 5 == 0 ) {
	    printf("\n");
	    if( (i+1) < len ) printf("%3d>",i);
	}
    }
    if( i % 5 != 0 ) printf("\n");
}    

/****************************************************************************
****	    This routine prints the (float) vector vec of lenght	  ****
****		    len to stdout with a heading head.		  	  ****
 ****************************************************************************/
int dumpfvec( len, vec, head )
int   len;
float *vec;
char  *head;
{
    register int i;

    printf("\n%s\n  0>",head );
    for( i=0; i<len; i++ ) {
	printf("%13.6e|",vec[i]);
	if( (i+1) % 5 == 0 ) {
	    printf("\n");
	    if( (i+1) < len ) printf("%3d>",i);
	}
    }
    if( i % 5 != 0 ) printf("\n");
}    
@//E*O*F dcgdrv.c//
chmod u=r,g=r,o=r dcgdrv.c
 
echo x - dcg.c
sed 's/^@//' > "dcg.c" <<'@//E*O*F dcg.c//'
  /************************************************************************
 ***************************************************************************
*****************************************************************************
*****					 				*****
*****		    Preconditioned Conjugate Gradient			*****
*****		    Double Precision (double) version.			*****
*****					 				*****
*****************************************************************************
 ***************************************************************************
  ************************************************************************/
#include "dcg.h"
#include <math.h>
static char DCGsid[] = "@(#)dcg.c	1.2  5/19/86";

#define NULL 0

int dcg( n, atimes, aptr, msolve, mptr, x, b, itmax, eps, err )
struct DCGAMAT *aptr;
struct DCGMMAT *mptr;
int (*atimes)(), (*msolve)();
int	n, itmax;
double	*x, *b, eps, *err;
/*
  Purpose:
    This routine does preconditioned conjugate gradient iteration
    on the symmetric positive definte system Ax = b.

  Input:
    n		Size of system to be iterated (ie A is nxn).
    aptr	Pointer to the structure defining the matrix A.
		This structure need not be known by this routine.
		The structure is passed on to the atimes routine which does
		the matrix multiply.  See below for the defintion of atimes.
    mptr	Pointer to the structure defining the preconditioning matrix 
	        M.  M should be chosen so that:
		  Mz = r is easy to solve for z, given r (compared to Ax = b).
		  Inv(M) is close to Inv(A).  Inv = inverse.
		Obviously, the better choice of M the faster the convergence.
    x		Initial guess of the solution.  If no information on the 
	        solution is available then use a zero vector or b.
    b		Right hand side.
    eps		Requested error tollerence.  System is iterated until
		||b-Ax||/||b|| < eps.  Normal choice is 1.0e-6.
    itmax	Maximum number of iterations the user is willing to allow.

  Output:
    dcg		Returns number of iterations taken.  Error Returns follow:
    		dcg = itmax 	May have not converged.  Check err.
		dcg = -1	System size too small (n<2).
		dcg = -2	Could not get memory for tempory vectors.
		dcg = -3	Error tollerence requested too tight.
		dcg = -4	Divergence of iteration detected.
		dcg = -10	User declared error in atimes.
		dcg = -11	User declared error in msolve.
    err		L2 norm of the relative residual error estimate 
	        ||b-Ax||/||b||.  This estimate of the true error is not always
		very accurate.  Note that err is a return value so the user
		must call dcg with a pointer to err.

  User Defined Functions and Structures.
    aptr	The user defines an external structure DCGAMAT and hands this
		routine a pointer to that structure.  The structure should
		contain all the information the user needs to calculate
		A times a vector.  See DCG.H for an example.

    mptr	The user defines an external structure DCGMMAT and hands this
		routine a pointer to that structure.  The structure should
		contain all the information the user needs to solve the system
		Mz=r for z.  See DCG.H for an example.

    int atimes( n, aptr, x, y )  
	        The user should write a routine that calculates y = Ax, given
		an x vector.  Return non-zero if an error occures.
    int msolve( n, mptr, z, r, atimes )
	        The user should write a routine that solves Mz=r for z, given
		a vector r.  Return non-zero if an error occures.
*/
{
  char   *malloc();
  register int i;
  double  *r, *z, *p;
  extern double ddot();
  double bnorm, bknum, bkden, bk;
  double akden, ak;
  double toobig;	/* When err est gets this big we signal divergence! */

  /* Check input data and allocate temporary storage. */
  if( n<2 ) return( -1 );
  if( (r = (double *)malloc( n*sizeof(double) )) == NULL ) return( -2 );
  if( (z = (double *)malloc( n*sizeof(double) )) == NULL ) return( -2 );
  if( (p = (double *)malloc( n*sizeof(double) )) == NULL ) return( -2 );
  if( eps < 1.0e-12 ) return( -3 );

  /* Calculate norm of right hand size. */
  bnorm = sqrt( ddot( n, b, 1, b, 1 ) );

  /* Calculate r = b - Ax and initinal error. */
  if( atimes( n, aptr, x, r ) ) return( -10 );
  dexopy( n, r, b, r, 2 );
  *err = sqrt( ddot(n, r, 1, r, 1 ) )/bnorm;
  if( *err < eps ) return( 0 );
  toobig = *err * 1.0e2;

  /* Iterate!!! */
  for( i=0; i<itmax; i++ ) {

    /* Solve Mz = r. */
    if( msolve( n, mptr, aptr, z, r, atimes ) ) return( -11 );

    /* Calculate bknum = (z, Mz) and p = z (first iteration). */
    if( i == 0 ) {
      dcopy( n, z, 1, p, 1 );
      bknum = bkden = ddot( n, z, 1, r, 1 );
    } 
    else {
      /* Calculate bknum = (z, r), bkden and bk. */
      bknum = ddot( n, z, 1, r, 1 );
      bk    = bknum/bkden;
      bkden = bknum;

      /* Calculate p = z + bk*p, z = Ap, akden = (p, Ap) and ak. */
      daxpyx( n, bk, p, 1, z, 1 );
    }
    if( atimes( n, aptr, p, z ) ) return( -10 );
    akden = ddot( n, p, 1, z, 1 );
    ak    = bknum/akden;

    /* Update x and r. Calculate error. */
    daxpy( n,  ak, p, 1, x, 1 );
    daxpy( n, -ak, z, 1, r, 1 );
    *err = sqrt( ddot( n, r, 1, r, 1 ) )/bnorm;
    if( *err < eps ) return( i+1 );
    if( *err > toobig ) return( -4 );
  }				/* end for loop. */
  return( itmax );
}				/* end dcg. */
@//E*O*F dcg.c//
chmod u=r,g=r,o=r dcg.c
 
echo x - dcg.h
sed 's/^@//' > "dcg.h" <<'@//E*O*F dcg.h//'
/* SID  "@(#)dcg.h	1.1  2/4/86" */
 /************************************************************************
*****									*****
*****   Structure defintions for Preconditioned Conjugate Gradient	*****
*****    		Double precision (double).		 	*****
*****	This file uses the example of a Symmetric Banded Matrix.	*****
*****									*****
  ************************************************************************/
/* 
   The following gives an array (of length 10) of pointers to doubles.
     double *a[10]; 
   Now assume that each a[i] has been given space for an array of doubles.  
   Then the following is true:
       a[i]       can be thought of as a pointer to the i-th array of doubles,
       *(a[i]+j)  is the j-th element of the i-th array.
   In your user code you are at some point going to need to pass 
   the DCGAMAT and DCGMMAT structures to a routine to do this declare 
struct DCGAMAT a;
struct DCGMMAT m;
   Which sets aside storage for the matrix pointers.  Then pass &a and &m
   to routines (e.g., dcg).  Note that you must use the address operator here.
   Inside the called routines you declare:
struct DCGAMAT *a;
struct DCGMMAT *m;
   The following shows how to reference things for the definitions of 
   DCGAMAT and DCGMMAT structures given below.
     a->hbw is the (value of, as apposed to a pointer to) half band width.
     a->diagptr[i] is a pointer to the i-th diagional (an array of doubles).
     *(a->diagptr[i]+j) is the j-th element of the i-th diagional (array).
   Here we think of the 0-th diagional as the main diagional, the 1-st 
   diagional as the 1-st upper diagional and so on.
*/

#define MAXBND 100	/* Maximum number of bands(including diagional). */

struct DCGAMAT {	/* Structure definition for the A matrix structure. */
  int   rd;		/* Number of rows in the matrix. */
  int	hbw;		/* Half band width.  Number of bands = 2hbw + 1. */
			/* Pointers to the main and upper diagionals. */
  double *diagptr[MAXBND];
};

struct DCGMMAT {	/* Struct definition for preconditioner matrix M. */
  int	itype;		/* Type of msolve to do. The user may want code up 
	  		   several choices for M and then set this parameter
			   before calling the iteration routine to specify
			   which type to use. */
  int   rd;		/* Number of rows in the matrix. */
  int	hbw;		/* Half band width.  Number of bands = 2hbw + 1. */
			/* Pointers to diagionals. */
  double *diagptr[MAXBND];
};
@//E*O*F dcg.h//
chmod u=r,g=r,o=r dcg.h
 
echo x - dbatimes.c
sed 's/^@//' > "dbatimes.c" <<'@//E*O*F dbatimes.c//'
 /**************************************************************************
*****			Sample atimes routine				*****
*****		    Designed for banded Matrices.			*****
 **************************************************************************/

#include "dcg.h"
static char DBATIMESsid[] = "@(#)dbatimes.c	1.2  5/19/86";

int dbatimes( n, aptr, x, y )
int	n;
struct DCGAMAT *aptr;
double	*x, *y;
/*
  Purpose
    To multiply A times a vector, viz., y = Ax.
*/
{
  register int i, j;
  register double *diag, *ytmp, *xtmp;

  /* Check input data.*/
  if( n<1 ) return( 1 );
  if( aptr->hbw < 1 ) return( 1 );

  /* Do main diagional */
  diag = aptr->diagptr[0];
  xtmp = x;
  ytmp = y;
  for( i=0; i<n; i++, xtmp++, ytmp++, diag++ )
    *ytmp = *diag * (*xtmp);

  /* Loop over other diagionals. */
  for( j=1; j<=aptr->hbw; j++ ) {

    /* Upper diagional. */
    diag = aptr->diagptr[j];
    xtmp = x+j;
    ytmp = y;
    for( i=0; i<n-j; i++, ytmp++, xtmp++, diag++ )
      *ytmp += *diag * (*xtmp);

    /* Lower diagional. */
    diag = aptr->diagptr[j];
    xtmp = x;
    ytmp = y+j;
    for( i=0; i<n-j; i++, ytmp++, xtmp++, diag++ )
      *ytmp += *diag * (*xtmp);
  }					/* End of diagional loop. */
  return( 0 );
}					/* End of dbatimes. */
@//E*O*F dbatimes.c//
chmod u=r,g=r,o=r dbatimes.c
 
echo x - dbmsolve.c
sed 's/^@//' > "dbmsolve.c" <<'@//E*O*F dbmsolve.c//'
 /**************************************************************************
*****			Sample msolve routine				*****
*****		    Designed for BANDED matrices.			*****
 **************************************************************************/

#include "dcg.h"
static char DBMSOLVEsid[] = "@(#)dbmsolve.c	1.2  5/19/86";

int dbmsolve( n, mptr, aptr, z, r, atimes )
int	n; 
struct DCGMMAT *mptr;
struct DCGAMAT *aptr;
double	*z, *r;
int (*atimes)();
/*
  Purpose
    To calculate an approximate inverse of A times a vector.

  Types of inverses
    itype = 1	Diagional scaling.   M=DIAG(A)
          = 2   Parameterized Incomplete Inverse (one matrix multiply).
	        WARNING: The User must first diagionaly scale A before
		this approximat inverse is useful.

 User supplied routines
    int atimes( n, aptr, x, y )  The user should write a routine that calculates
    		y = Ax.  Return non-zero if an error occures.
*/
{
  register int i;
  register double *diag;
  double gam0, gam1, gam01;

  /* Check input data. */
  if( n<1 ) return( 1 );
  if( mptr->hbw < 1 ) return( 1 );
  
  /* Options are as follows: */
  switch( mptr->itype ) {
  case 1:				/* Diagional Scaling. */
    diag = mptr->diagptr[0];		/* z = r/DIAG(A). */
    for( i=0; i<n; i++, z++, r++, diag++ ) {
      *z = *r / (*diag);
    }
    break;

  case 2:				/* Prameterized Incomplete Inverse p=1 */
    gam0  = 7.0/6.0;			/* z = INV(m)r = (gam0*I+gam1G)r, G=I-A */
    gam1  = 5.0/6.0;
    gam01 = gam0+gam1;

    if( atimes( n, aptr, r, z ) ) return( 1 );
    diag = z;
    for( i=0; i<n; i++, diag++ )
      *diag *= gam1;
    for( i=0; i<n; i++, z++, r++ )
      *z = gam01 * (*r) - *z;
    break;
  default:				/* Wrong itype given. */
    return( 1 );
  }
  return( 0 );
}				/* End of sbmsolve. */
@//E*O*F dbmsolve.c//
chmod u=r,g=r,o=r dbmsolve.c
 
echo x - dcgutil.c
sed 's/^@//' > "dcgutil.c" <<'@//E*O*F dcgutil.c//'
  /***************************************************************
 ******************************************************************
****		  Utilitie routines for CG.C			****
 ******************************************************************
  ****************************************************************/

#include "dcg.h"
static char DCGUTILsid[] = "@(#)dcgutil.c	1.2  5/19/86";

void dsbdumpmat( a, head )
struct DCGAMAT *a;
char	    *head;
/*
  This routine prints out a Double Symmetric Banded MATRIX with headding head.
*/
{
  register int i, j;
  int k, jj, ncolmod, ncolrem;

  jj = a->hbw+1;				/* Must include diagional. */
  ncolmod = jj / 6;
  ncolrem = jj - 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", *(a->diagptr[j]+i));
      printf("\n");
    }
    
    /* Clean up remainder of this row. */
    for( jj=0;  jj<ncolrem;  jj++, j++) 
      printf("%12.4e", *(a->diagptr[j]+i));
    printf("\n");
    }
  printf("\n");
}				/* End of SYM_BND_MATDUMP */

int dsbscale( n, aptr, mptr, x, b, scale, dox )
int n, scale, dox;
struct DCGAMAT *aptr;
struct DCGMMAT *mptr;
double *x, *b;
/*
  Purpose
    To (diagionaly) scale or unscale the matrix problem Ax=b.
    The structure of a is assumed to be Double precision (double) 
    Symmetric Banded (SSB).

  Input
    n	  System size.
    aptr  Pointer to the A matrix structure (this routine assumes a
          DIAGIONAL structure).
    mptr  Pointer to the M matrix structure.  This is used to hold
	  the sqrt of the diagional diagional for unscaling. 
    x     Pointer to the solution (if known. See dox).
    b     Pointer to the right hand side.
    scale Boolian which when true (non-zero) => scale, false => unscale.
    dox   Boolian which when true (non-zero) => scale/unscale solution x.
*/
{
  register int i, j;
  double *t1, *t2;
  double sqrt();

  if( n<2 ) return( 1 );		/* Check input parmaters. */
  if( aptr->hbw<1 ) return( 1 );

  if( scale ) {				/* Set up for SCALING. */
    t1 = mptr->diagptr[0];
    t2 = aptr->diagptr[0];
    for( i=0; i<n; i++, t1++, t2++ ) {
      *t1 = 1.0/sqrt( *t2 );
      *t2 = 1.0;
    }
  }
  else {				/* Set up for UNSCALING. */
    t1 = mptr->diagptr[0];
    t2 = aptr->diagptr[0];
    for( i=0; i<n; i++, t1++, t2++ ) {
      *t1 = 1.0/(*t1);
      *t2 = (*t1) * (*t1);
    }
  }

  /* Do scaling/unscaling of A. */
  for( j=1; j<=aptr->hbw; j++ ) {
    t1 = mptr->diagptr[0];
    t2 = aptr->diagptr[j];
    for( i=0; i<n-j; i++, t1++, t2++ ) 
      *t2 *= (*t1)*(*(t1+j));
  }

  /* Do scaling/unscaling of RHS. */
  t1 = mptr->diagptr[0];
  t2 = b;
  for( i=0; i<n; i++, t1++, t2++ )
    *t2 *= *t1;

  /* Do scaling/unscaling of Solution (if requtested). */
  if( dox ) {
    t1 = mptr->diagptr[0];
    t2 = x;
    for( i=0; i<n; i++, t1++, t2++ )
      *t2 /= *t1;
  }

#ifdef DEBUG
  /* Echo out matrix. */
  printf("Matrix diagionals,rhs,solution follow...\n");
  for( i=0, t1=b, t2=x; i<n; i++, t1++, t2++ ) {
    printf("%3d",i);
    for( j=0; j<=aptr->hbw; j++ ) 
	printf("%15e",*(aptr->diagptr[j]+i));
    printf("|%15e|%15e\n", *t1, *t2);
  }
#endif

 return( 0 );
}
@//E*O*F dcgutil.c//
chmod u=r,g=r,o=r dcgutil.c
 
echo x - dblas.c
sed 's/^@//' > "dblas.c" <<'@//E*O*F dblas.c//'
  /***************************************************************
  *****************************************************************
 *******************************************************************
*****								*****
*****				BLAS				*****
*****		Basic Linear Algebra Subroutines		*****
*****	      Written in the C Programming Language.		*****
*****		 << DOUBLE PRECISION VERSION >>			*****
*****								*****
*****	Functions include:					*****
*****	idamax, daxpy, daxpyx, dcopy, ddot, dnrm2,		*****
*****	dexopy, dfill						*****
*****								*****
 *******************************************************************
  *****************************************************************
   ***************************************************************/
#include <math.h>
#ifndef SMALLdp
#define SMALLdp		2.22507385850720e-308		/* Smallest (pos) binary float */
#endif
#ifndef HUGEdp
#define HUGEdp		1.79769313486232e+308		/* Largest binary float */
#endif

int idamax( n, sx, incx )
double	*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.

*/
{
  double 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 daxpy( n, sa, sx, incx, sy, incy )
double *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 double vector to scale.
    incx	Storage incrament for sx.
    sy		Pointer to double 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 daxpyx( n, sa, sx, incx, sy, incy )
double *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 double vector to scale.
    incx	Storage incrament for sx.
    sy		Pointer to double 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 dcopy( n, sx, incx, sy, incy )
double  *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 ddot( n, sx, incx, sy, incy )
double	*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;
  double 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 ---DDOT--- */

double dnrm2( n, sx, incx )
double	*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, r2mach();
  double xmax;

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

  cutlo = sqrt( SMALLdp/r2mach() );			/* Calculate near underflow */
  cuthi = sqrt( HUGEdp );				/* 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 ---DNRM2--- */

int dscal( n, sa, sx, incx )
double  sa, *sx;
int     n, incx;
/*
    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.
*/
{
  register 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 dexopy( n, v, x, y, itype )
int	n, itype;
double	*v, *x, *y;
/*
  Purpose: v = x op y
    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++, x++, y++ )
      *v = *x + *y;
  else			/* SUBTRACTION. */
    for( i=0; i<n; i++, v++, x++, y++ )
      *v = *x - *y;
}

void dfill( n, v, val )
int	n;
double	*v, val;
/*
  Purpose
    To fill the DOUBLE vector v with the DOUBLE value val.
*/
{
  register int i;

  if( n<1 ) return;
  for( i=0; i<n; i++, v++ )
    *v = val;
}

double r2mach()
/* ---------------------------------------------------------------------
        This routine computes the unit roundoff for double 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.
----------------------------------------------------------------------*/
{
    double u = 1.0e0, comp;
 
    do {
        u *= 0.5e0;
        comp = 1.0e0 + u;
    }
    while( comp != 1.0e0 );
    return( (double)u*2.0e0 );
}
/*-------------------- end of function r2mach ------------------------*/
 
@//E*O*F dblas.c//
chmod u=rw,g=r,o=r dblas.c
 
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 <<\!!!
      45     284    2558 READ.ME
      53     154    1169 makefile
     314    1449   10369 dcgdrv.c
     136     856    5250 dcg.c
      52     376    2427 dcg.h
      50     185    1252 dbatimes.c
      63     267    1780 dbmsolve.c
     124     461    3281 dcgutil.c
     427    1640   10047 dblas.c
    1264    5672   38133 total
!!!
wc  READ.ME makefile dcgdrv.c dcg.c dcg.h dbatimes.c dbmsolve.c dcgutil.c dblas.c | 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