From oleg@nrcbsa.bio.nrc.ca Thu May 23 08:51:21 1991
Return-Path: <oleg@nrcbsa.bio.nrc.ca>
Received: from nrcnet0.nrc.ca by CS.UTK.EDU with SMTP (5.61++/2.5.1s-UTK)
	id AA25974; Thu, 23 May 91 08:51:05 -0400
Message-Id: <9105231250.AA28661@ nrcbsa.bio.nrc.ca>
Date: Thu, 23 May 91 08:50:38 EDT
From: Dr. Oleg Keselyov <oleg@nrcbsa.bio.nrc.ca>
To: dongarra@cs.utk.edu
Subject: serv.shar
Status: RO

#--------------------------------CUT HERE-------------------------------------
#! /bin/sh
#
# This is a shell archive. Save this into a file, edit it
# and delete all lines above this comment.  Then give this
# file to sh by executing the command 'sh file'. The files
# will be extracted into the current directory owned by
# you with default permissions.

# This arcive contains h-files and a set of service functions that provide 
# a friendly environment, in particular, for Numerical 
# programming in C. These files are needed to compile the
# programs in packages
# brent.shar vector.shar fft.shar
# Verification programs and data are included as well.
 
#			Contents of the arcive
#
# H-file needed for compilation. They are special files not to be confused
# with the general purpose files, say, /usr/include/math.h
# They are 'included' as "math.h" rather than <math.h>.
#	math.h		Contains machine-specific parameters (EPSILON)
#			as well as declaration of all the my functions.
#	assert.h	Declaration of assure and assert macros
#	stdio.h		Constants and functions used by the test drivers
#
# Set of the service math and general purpose functions
#	servio.c	Only _error() function is used in this package
#	servmath.c	Especially powi functions
#	epsilon.c	Obtain the machine epsilon
#
#	vservmath.c	Validate service mathematical functions
#	vservmath.dat
#	epsilon.dat


echo 'x - math.h'
sed 's/^X//' << '________This_Is_The_END________' > math.h
X/* Revision history
X *	EPSILON, SQRT_EPSILON added		28-apr-1990	
X *	zeroin function       added		 4-may-1990	
X *      fmin, fmax, 
X *      iabs, imin, imax,
X *      powi  functions	      added		 8-may-1990
X *	nrandom		      added		10-may-1990 
X * 	fminbr		      added		11-may-1990
X *	alis		      added		18-jun-1990
X *	PI		      added		18-jun-1990	
X *	minvhol		      added		10-sep-1990
X *	hjmin		      added		15-oct-1990
X */
X
X#pragma once
X
X#include </usr/include/math.h>		/* Include the standard math.h file	*/
X
X#if CC$gfloat
X# define HUGE_VAL 8.988465674311578540726371186585e+307
X#else
X# define EPSILON	2.22045e-16
X# define SQRT_EPSILON	1.49012e-08
X# define PI		3.14159265358979300
X#endif
X			/* Service mathematical functions		*/
Xint iabs(int);    		/* Absolute value of a integer   	*/
Xint imax(int,int);		/* Maximum of two integers		*/
Xint imin(int,int);		/* Minimum of two integers		*/
X
Xdouble fmax(double,double);  	/* Maximum of two doubles		*/
Xdouble fmin(double,double);  	/* Minimum of two doubles		*/
X
Xdouble powi(double r,int exp);	/* r up exp				*/
X
Xdouble nrandom(void);		/* Random number distributed as N(0,1)	*/
X
X				/* Obtain a zero of function f		*/
X				/* over the range [ax,bx] with the 	*/
X				/* accuracy tol. tol may be specified as 0.0*/
Xdouble zeroin(double ax, double bx, double (*f)(), double tol);
X
X				/* Find a minimum of function f		*/
X				/* over the range [a,b] with the	*/
X				/* accuracy tol > 0.			*/
X				/* Returns an approx. to the min location*/
Xdouble fminbr(double a, double b, double (*f)(), double tol);
X
X#ifdef V$elem_type
X				/* Aitken-Lagrange interpolation to the	*/
X				/* point x over the table with even mesh*/
X				/* x[i] = x0 + s*(i-Lwb(y)),		*/
X				/* i = Lwb(y)..Upb(y)			*/
Xextern double alis(double q, double x0, double s, VECTOR y);
Xextern double alin(double q, VECTOR x, VECTOR y);
X
X				/* Inverse a positive definite symmetric*/
X				/* matrix, return condition number	*/
Xextern double minvhol(float *A,const int N);
X
X				/* Find the local minimum of a given	*/
X				/* function by the Hook-Jeevse method	*/
Xdouble hjmin(				/* Return the function value at min  */
X	VECTOR b,			/* On input - initial guess to min loc*/
X					/* On output - loc for the min found  */
X 	const VECTOR h0,		/* Initial values for the steps along */
X 					/* each directions		*/
X 	double (*f)(const VECTOR x)	/* Procedure to compute a function   */	
X 					/* value at the specified point	*/
X	    );
X#endif
________This_Is_The_END________
if test `wc -l < math.h` -ne 71; then
echo 'shar: math.h was damaged during transit (should have had 71 lines)'
fi


echo 'x - assert.h'
sed 's/^X//' << '________This_Is_The_END________' > assert.h
X/*	ASSERT - V3.0	*/
X
X#include </usr/include/assert.h>	/* Include the standard assert	*/
Xextern void _error( const char * message,... ); /* Print a message and abort*/
X# define assure(expr,message) 					\
X	if	(expr) ;				\
X	else _error(message);
________This_Is_The_END________
if test `wc -l < assert.h` -ne 7; then
echo 'shar: assert.h was damaged during transit (should have had 7 lines)'
fi


echo 'x - stdio.h'
sed 's/^X//' << '________This_Is_The_END________' > stdio.h
X/*
X ************************************************************************
X *
X *		     Standard input/output
X *			Service functions
X */
X
X#include </usr/include/stdio.h>
X
X				/* Strings of symbols			*/
X				/* They may be used as a delimiting lines*/
Xextern const char _Minuses [];
Xextern const char _Asteriscs [];
Xextern const char _Equals [];
X
X
X				/*	Safe file open function		*/
X				/* Null pointer is never returned, so   */
X				/* the user may be safe			*/
XFILE * fopenc(
X	const char * name, 		/* File name                     */
X	const char * mode  		/* Open mode                     */
X	     );
X
X
X				/* Print an error message at stderr and	*/
X				/* abort				*/
Xvoid _error(
X	const char * message,		/* Message to be printed	*/
X	...                             /* Additional args to printf	*/
X	   );
X
X				/* Print a message at stderr 		*/
Xvoid message(
X	const char * text,		/* Message to be printed	*/
X	...                             /* Additional args to printf	*/
X	   );
X
X				/* PASCAL-like reading functions	*/
Xvoid nl(				/* Skip all characters until	*/
X	const FILE *fp			/* the end-of-line		*/
X       );
X
X               			/* Format read of the input stream	*/
Xvoid readf(
X	const FILE * fp, 		/* File pointer			*/
X	const char * fstring,		/* Format string as in scanf	*/
X	...);                           /* Other args - ptr to vars to	*/
X                                        /* be read			*/
X
X
X               			/* Same as above excepting i/o pointer	*/
X				/* being moved to the beginning of the	*/
X				/* next record				*/
Xvoid readln(
X	const FILE * fp, 		/* File pointer			*/
X	const char * fstring,		/* Format string as in scanf	*/
X	...);                           /* Other args - ptr to vars to	*/
X                                        /* be read			*/
X
________This_Is_The_END________
if test `wc -l < stdio.h` -ne 60; then
echo 'shar: stdio.h was damaged during transit (should have had 60 lines)'
fi


echo 'x - servio.c'
sed 's/^X//' << '________This_Is_The_END________' > servio.c
X/*
X ************************************************************************
X *			Service C functions 
X *	       concerning the input/output facilities
X */
X
X#include "stdio.h"
X#include "assert.h"
X#include <stdarg.h>
X#include <malloc.h>
X#include <string.h>
X
X/*
X *-----------------------------------------------------------------------
X *		Some global constant pertaining to input/output
X */
X
Xconst char _Minuses [] = "\
X-------------------------------------------------------------------------------";
X
Xconst char _Asteriscs [] = "\
X*******************************************************************************";
X
Xconst char _Equals [] = "\
X===============================================================================";
X
X
X
X/*
X *------------------------------------------------------------------------
X *			Safe file open function
X * Synopsis
X *	FILE *fp = fopenc(const char * file_name, const char * mode)
X *
X * Opens the file specified like 'fopen' and returns ptr to the FILE
X * structure corresponding to the file.
X * On open failure an error message is typed and program is aborted.
X * Null pointer is never returned, so the user may be safe using fp in
X * his own.
X */
X
XFILE * fopenc(name,mode)
Xconst char * name;			/* File name                     */
Xconst char * mode;			/* Open mode                     */
X{
X  register FILE * fp;
X  if( (fp = fopen(name,mode)) == NULL )
X    {
X      fprintf(stderr,"\n\nCan't open a file '%s'\n",name);
X      perror("");
X      exit(2);
X    }
X  return fp;
X}
X
X
X/*
X *------------------------------------------------------------------------
X *                     Open the level II file
X *              with checking of successfullness
X *         and explicit specified file name extension
X */
X
XFILE * fopencex(name,ext,mode)
Xconst char * name;			/* File name                     */
Xconst char * ext;			/* File name extension           */
Xconst char * mode;			/* Open mode                     */
X{
X  char fname[65];                      /* The full file name            */
X  register int i;
X
X			/* Copy the name of file
X  strncpy(fname,name,sizeof(fname));
X  for(i=0; i<sizeof(fname); i++)
X  {
X    if( name[i] == '\0' )
X      break;
X  }
X/*  strmfe(fname,name,ext); */
X  return fopenc(fname,mode);
X}
X
X
X
X/*
X *-----------------------------------------------------------------------
X *			PASCAL-like read routines
X */
X
X			/* Flush all data until the end-of-line		*/
Xvoid nl(fp)
Xconst FILE * fp;
X{
X  register char c;
X
X  while( c=getc(fp), c != '\n' && c != EOF )
X    ;
X}
X
X
X                        /* Format read of the input stream		*/
Xvoid readf(
X	const FILE *fp,			/* File being read		*/
X	const char * fstring,		/* Format string, as scanf	*/
X	...
X)
X{
X  register int no;
X  register int no_read;
X  va_list args;
X  char *p1, *p2, *p3, *p4, *p5, *p6, *p7, *p8, *p9;
X  register char *p = fstring;
X  char scratch [40];			/* Scratch buffer		*/
X
X  va_start(args,fstring); 		/* Init 'args' to the beginning of */
X  p1 = (char *)va_arg(args,char *);
X  p2 = (char *)va_arg(args,char *);
X  p3 = (char *)va_arg(args,char *);
X  p4 = (char *)va_arg(args,char *);
X  p5 = (char *)va_arg(args,char *);
X  p6 = (char *)va_arg(args,char *);
X  p7 = (char *)va_arg(args,char *);
X  p8 = (char *)va_arg(args,char *);
X  p9 = (char *)va_arg(args,char *);
X
X  for(no=0; (p=strchr(p,'%')) != 0; no++)
X    p++;				/* Estimate no. of data item to be*/
X					/* read ( = no. of % signs )	*/
X
X  assure( no < 9, 
X	"No more than 9 data items are allowed to read with one read call");
X
X  if( (no_read = fscanf(fp,fstring,p1,p2,p3,p4,p5,p6,p7,p8,p9)) != no )
X   _error("Only %d items have been read of %d required\ntail string '%s'",
X          no_read, no, fgets(scratch,sizeof(scratch)-1,fp) );
X}
X
X
X				/* Same as above but i/o ptr		*/
X				/* is moved to the begiining of next rec*/
Xvoid readln(
X	const FILE * fp,    		/* File being read		*/
X	const char * fstring,		/* Format string, as scanf	*/
X	...
X)
X{
X  va_list args;
X  char *p1, *p2, *p3, *p4, *p5, *p6, *p7, *p8, *p9;
X
X  va_start(args,fstring); 		/* Init 'args' to the beginning of */
X  p1 = (char *)va_arg(args,char *);
X  p2 = (char *)va_arg(args,char *);
X  p3 = (char *)va_arg(args,char *);
X  p4 = (char *)va_arg(args,char *);
X  p5 = (char *)va_arg(args,char *);
X  p6 = (char *)va_arg(args,char *);
X  p7 = (char *)va_arg(args,char *);
X  p8 = (char *)va_arg(args,char *);
X  p9 = (char *)va_arg(args,char *);
X
X  readf(fp,fstring,p1,p2,p3,p4,p5,p6,p7,p8,p9);
X  nl(fp);
X}
X
X
X/*
X *------------------------------------------------------------------------
X *                      Write the "string of" to a file
X * File is assumed to be opened in binary mode and set at appropriate
X * position. Integer descriptor (string length) preceds to string
X * characters.
X *------------------------------------------------------------------------
X */
X
Xvoid write_string(str,fp)
Xchar * str;
XFILE * fp;
X{
X  int len;                             /* String length                 */
X
X  len = ( str == 0L ? 0 : strlen(str) );
X/*  writei(len,fp);   */
X  if( len != 0 )
X    assert( fwrite(str,len,1,fp) == 1 ); /* Write a string symbols     */
X}
X
X
X
X/*
X *------------------------------------------------------------------------
X *                     Read the "string of" from a file
X *                         and return the ptr to it.
X * File is assumed to be opened in binary mode and set at appropriate
X * position. Integer descriptor (string length) presedes to string
X * characters.
X *
X */
X
Xchar * read_string(fp)
XFILE * fp;
X{
X  int len;                             /* String length                 */
X  char *p;
X
X/*  readi(len,fp);  */                     /* Read the string length        */
X  if( len == 0 )
X    return 0L;                         /* Null string has been read     */
X  assert( len > 0 && len < 25*80*2 );
X  p = malloc(len+1);                   /* Plus null char at end of str  */
X  assert( fread(p,len,1,fp) == 1 );    /* Read string symbols           */
X  p[len] = '\0';                       /* Append the null char          */
X  return p;
X}
X
X/*
X *------------------------------------------------------------------------
X *	        Print an error message at stderr and abort
X * Synopsis
X *	void _error(const char * message,... );
X *	Message may contain format control sequences %x. Additional 
X *	arguments for them a passed next to the message.
X */
X
X
Xvoid _error(
X	const char * message,
X	...
X)
X{
X  va_list args;
X  va_start(args,message);		/* Init 'args' to the beginning of */
X					/* the variable length list of args*/
X  fprintf(stderr,"\n_error:\n"); 	
X  vfprintf(stderr,message,args);
X  fputs("\n",stderr);
X  abort();
X}
X
X
X/*
X *------------------------------------------------------------------------
X *	       		 Print a message at stderr
X * Synopsis
X *	void message(const char * text,... );
X *	Message may contain format control sequences %x. Additional 
X *	arguments for them a passed next to the message.
X */
X
Xvoid message(
X	const char * text,
X	...
X)
X{
X  va_list args;
X  va_start(args,text);		/* Init 'args' to the beginning of */
X					/* the variable length list of args*/
X  vfprintf(stderr,text,args);
X}
________This_Is_The_END________
if test `wc -l < servio.c` -ne 258; then
echo 'shar: servio.c was damaged during transit (should have had 258 lines)'
fi


echo 'x - servmath.c'
sed 's/^X//' << '________This_Is_The_END________' > servmath.c
X/*
X ************************************************************************
X *			Service C functions 
X *	      concerning the mathematical computations
X */
X
X#include "math.h"
X
Xdouble fmin(x1,x2) 	/* Minimum of two doubles			*/
Xconst double x1,x2;
X{
X  if( x1 > x2 )
X    return x2;
X  else
X    return x1;
X}
X
X
Xdouble fmax(x1,x2) 	/* Maximum of two doubles			*/
Xconst double x1,x2;
X{
X  if( x1 > x2 )
X    return x1;
X  else
X    return x2;
X}
X
Xint iabs(x) 		/* Absolute value of a integer			*/
Xint x;
X{
X  if( x < 0 )
X    x = -x;
X  return x;
X}
X
Xint imin(x1,x2) 	/* Minimum of two integers			*/
Xconst int x1,x2;
X{
X  if( x1 > x2 )
X    return x2;
X  else
X    return x1;
X}
X
X
Xint imax(x1,x2) 	/* Maximum of two integers			*/
Xconst int x1,x2;
X{
X  if( x1 > x2 )
X    return x1;
X  else
X    return x2;
X}
X
X
Xdouble powi(x,e)	/* x up integer power				*/
Xconst double x;
Xconst int e;
X{
X  register int i;
X  register double r;
X  
X  switch(e)
X  {
X	case 0:
X		return 1.0;
X
X	case 2:
X		return x*x;
X
X	case 3:
X		return x*x*x;
X
X	case -1:
X		return 1/x;
X
X	case -2:
X		return 1/x/x;
X  }
X
X  {
X    register double xn = x;
X    register int en = e;
X    if( e < 0 )
X    {
X      xn = 1/x;
X      en = -e;
X    }
X
X    for(r=xn,i=1; 2*i<=en; i*=2)
X      r *= r;
X    for(; i<en; i++)
X      r *= xn;
X
X    return r;
X  }
X
X}
X
X
X#include <stdlib.h>
X
Xdouble nrandom()	/* Random number NORMAL distributed as N(0,1)	*/
X{
X  register int i;
X  register double accum = 0;
X
X  for(i=0; i<12; i++)
X     accum += rand();
X
X  return accum/(double)RAND_MAX - 6.;
X}
________This_Is_The_END________
if test `wc -l < servmath.c` -ne 112; then
echo 'shar: servmath.c was damaged during transit (should have had 112 lines)'
fi


echo 'x - vservmath.c'
sed 's/^X//' << '________This_Is_The_END________' > vservmath.c
X/*		Verify SERVMATH functions		*/
X
X#include "math.h"
X
Xvoid timin(n1,n2)			/* Test imin function		*/
Xint n1,n2;
X{ printf("\nimin(%d,%d)=%d",n1,n2,imin(n1,n2)); }
X
Xvoid timax(n1,n2)		/* Test imax function		*/
Xint n1,n2;
X{ printf("\nimax(%d,%d)=%d",n1,n2,imax(n1,n2)); }
X
Xvoid tfmin(n1,n2)		/* Test fmin function		*/
Xdouble n1,n2;
X{ printf("\nfmin(%g,%g)=%g",n1,n2,fmin(n1,n2)); }
X
Xvoid tfmax(n1,n2)		/* Test fmax function		*/
Xdouble n1,n2;
X{ printf("\nfmax(%g,%g)=%g",n1,n2,fmax(n1,n2)); }
X
X
Xvoid tiabs(n1)   		/* Test iabs function		*/
Xint n1;      
X{ printf("\niabs(%d)=%d",n1,iabs(n1)); }
X
X
Xvoid tnrandom()			/* Test nrandom() function		*/
X{               		/* Making use of the fact		*/
X  register int i;               /* Chi2(n-1) = SUM[ nrandom() ], i=1,n	*/
X  register int trial_counter;
X
X  printf("\n\nTest the normal random number generator\n\n");
X  printf("No. of freedoms\t\tChi^2\n");
X
X  for(trial_counter=1; trial_counter < 25; trial_counter++ )
X  {
X    double chi2 = 0;	
X    for(i=0; i<trial_counter; i++)
X    {
X       register double t = nrandom();
X       chi2 += t*t;
X    }
X    printf("%d\t\t\t%g\n",trial_counter-1,chi2);
X  }
X  printf("\nEnd of test nrandom()\n\n");
X}
X
X
Xmain()
X{
X  register int i;
X
X  printf("\n\n\t\tTest SERVMATH functions\n\n");
X
X  tiabs(0);  tiabs(2);  tiabs(-2);
X
X  timin(1,0); timin(0,32767); timin(-32767,0); timin(4,55); timin(-6,10);
X  timax(1,0); timax(0,32767); timax(-32767,0); timax(4,55); timax(-6,10);
X
X  tfmin(1.0,0); tfmin(0,HUGE_VAL);
X  tfmin(-HUGE_VAL,0); tfmin(4.,55e20); tfmin(-6e-5,-4e-3);
X  tfmax(1.0,0); tfmax(0,HUGE_VAL);  /* Note! tfmax(4,55e20) is bad call! */
X  tfmax(-HUGE_VAL,0); tfmax(4.,55e20); tfmax(-6e-5,-4e-3);
X
X  tnrandom();
X
X  for(i= -7; i<=10; i++)
X     printf("\npowi(2,%d)=%g",i,powi(2,i));
X}
________This_Is_The_END________
if test `wc -l < vservmath.c` -ne 69; then
echo 'shar: vservmath.c was damaged during transit (should have had 69 lines)'
fi


echo 'x - vservmath.dat'
sed 's/^X//' << '________This_Is_The_END________' > vservmath.dat
X
X
X		Test SERVMATH functions
X
X
Xiabs(0)=0
Xiabs(2)=2
Xiabs(-2)=2
Ximin(1,0)=0
Ximin(0,32767)=0
Ximin(-32767,0)=-32767
Ximin(4,55)=4
Ximin(-6,10)=-6
Ximax(1,0)=1
Ximax(0,32767)=32767
Ximax(-32767,0)=0
Ximax(4,55)=55
Ximax(-6,10)=10
Xfmin(1,0)=0
Xfmin(0,0)=0
Xfmin(-Infinity,0)=-Infinity
Xfmin(4,5.5e+21)=4
Xfmin(-6e-05,-0.004)=-0.004
Xfmax(1,-0.004)=1
Xfmax(1,-0.004)=1
Xfmax(-Infinity,-0.004)=-0.004
Xfmax(4,5.5e+21)=5.5e+21
Xfmax(-6e-05,-0.004)=-6e-05
X
XTest the normal random number generator
X
XNo. of freedoms		Chi^2
X0			1.75892
X1			6.18685
X2			0.0880468
X3			6.95459
X4			2.3511
X5			9.1004
X6			15.0925
X7			12.8436
X8			8.80412
X9			10.6615
X10			9.68159
X11			12.4396
X12			25.9775
X13			13.9597
X14			13.9225
X15			12.5468
X16			33.985
X17			8.75921
X18			22.0895
X19			20.0303
X20			15.2557
X21			19.0718
X22			13.7847
X23			26.6989
X
XEnd of test nrandom()
X
X
Xpowi(2,-7)=0.0078125
Xpowi(2,-6)=0.015625
Xpowi(2,-5)=0.03125
Xpowi(2,-4)=0.0625
Xpowi(2,-3)=0.125
Xpowi(2,-2)=0.25
Xpowi(2,-1)=0.5
Xpowi(2,0)=1
Xpowi(2,1)=2
Xpowi(2,2)=4
Xpowi(2,3)=8
Xpowi(2,4)=16
Xpowi(2,5)=32
Xpowi(2,6)=64
Xpowi(2,7)=128
Xpowi(2,8)=256
Xpowi(2,9)=512
Xpowi(2,10)=1024
________This_Is_The_END________
if test `wc -l < vservmath.dat` -ne 78; then
echo 'shar: vservmath.dat was damaged during transit (should have had 78 lines)'
fi

echo 'x - epsilon.c'
sed 's/^X//' << '________This_Is_The_END________' > epsilon.c
X/*	  
X *   Obtain the machine EPSILON
X *  i.e. the smallest positive number which, been added to 1., yields 
X *  the result other than 1.
X */
X#include <math.h>
X
Xmain()
X{
X    double eps;
X    eps = 1.0;
X    while( eps + 1.0 != 1.0 )
X	eps /= 2.0;
X    eps *= 2.0;
X    
X    printf("FPU epsilon is\t %g\nsqrt(eps) is\t %g\n",
X	    eps, sqrt(eps) );
X}
________This_Is_The_END________
if test `wc -l < epsilon.c` -ne 18; then
echo 'shar: epsilon.c was damaged during transit (should have had 18 lines)'
fi


echo 'x - epsilon.dat'
sed 's/^X//' << '________This_Is_The_END________' > epsilon.dat
XFPU epsilon is	 2.22045e-16
Xsqrt(eps) is	 1.49012e-08
________This_Is_The_END________
if test `wc -l < epsilon.dat` -ne 2; then
echo 'shar: epsilon.dat was damaged during transit (should have had 2 lines)'
fi