cwplib.doc

seismic software,very useful

Notes:
	Adapted from LINPACK FORTRAN.


Prototypes:

int isamax (int n, float *sx, int incx);
float sasum (int n, float *sx, int incx);
void saxpy (int n, float sa, float *sx, int incx, float *sy, int incy);
void scopy (int n, float *sx, int incx, float *sy, int incy);
float sdot (int n, float *sx, int incx, float *sy, int incy);
float snrm2 (int n, float *sx, int incx);
void sscal (int n, float sa, float *sx, int incx);
void sswap (int n, float *sx, int incx, float *sy, int incy);

------------------------------------------------------------------------
14. Double precision BLAS (basic linear algebra subroutines)

Contents:
	a. idamax - return index of element with maximum absolute value
	b. dasum - return sum of absolute values
	c. daxpy - compute y[i] = a*x[i]+y[i]
	d. dcopy - copy x[i] to y[i] (i.e., set y[i] = x[i])
	e. ddot - return sum of x[i]*y[i] (i.e., dot product of x and y)
	f. dnrm2 - return square root of sum of squares of x[i]
	g. dscal - compute x[i] = a*x[i]
	h. dswap - swap x[i] and y[i]

Notes:
	Adapted from LINPACK FORTRAN.


Prototypes:

int idamax (int n, double *sx, int incx);
double dasum (int n, double *sx, int incx);
void daxpy (int n, double sa, double *sx, int incx, double *sy, int incy);
void dcopy (int n, double *sx, int incx, double *sy, int incy);
double ddot (int n, double *sx, int incx, double *sy, int incy);
double dnrm2 (int n, double *sx, int incx);
void dscal (int n, double sa, double *sx, int incx);
void dswap (int n, double *sx, int incx, double *sy, int incy);

------------------------------------------------------------------------
15. Single precision LINPACK general matrix functions

Contents:
	a. sgefa - LU factorization
	b. sgeco - LU factorization and condition number
	c. sgesl - Solve linear system Ax = b or A'x = b after LU factorization.


Prototypes:

void sgeco (float **a, int n, int *ipvt, float *rcond, float *z);
void sgefa (float **a, int n, int *ipvt, int *info);
void sgesl (float **a, int n, int *ipvt, float *b, int job);

------------------------------------------------------------------------
16. Double precision LINPACK general matrix functions

Contents:
	a. dgefa - LU factorization
	b. dgeco - LU factorization and condition number
	c. dgesl - Solve linear system Ax = b or A'x = b after LU factorization.


Prototypes:

void dgeco (double **a, int n, int *ipvt, double *rcond, double *z);
void dgefa (double **a, int n, int *ipvt, int *info);
void dgesl (double **a, int n, int *ipvt, double *b, int job);

------------------------------------------------------------------------
17. Single precision LINPACK QR decomposition functions

Contents:
	a. sqrdc - QR matrix decomposition

	b. sqrsl - Use QR decomposition to solve for coordinate transformations,
		projections, and least squares solutions.

	c. sqrst - Solve under-determined or over-determined least squares
		problems, with a user-specified tolerance.


Prototypes:

void sqrdc (float **x, int n, int p, float *qraux, int *jpvt,
	float *work, int job);
void sqrsl (float **x, int n, int k, float *qraux,
	float *y, float *qy, float *qty,
	float *b, float *rsd, float *xb, int job, int *info);
void sqrst (float **x, int n, int p, float *y, float tol,
	float *b, float *rsd, int *k,
	int *jpvt, float *qraux, float *work);

------------------------------------------------------------------------
18. Other Linear System Solving routines

	a. stoepf - symmetric Toeplitz linear system of equations Rf=g
		for f (single precision)

Notes:
	This routine does NOT solve the case when the main diagonal is
	zero, it just silently returns.


Prototype and parameters:

void stoepd (int n, double r[], double g[], double f[], double a[]);
n	i dimension of system
r	i array of top row of Toeplitz matrix
g	i array of right-hand-side column vector
f	o array of solution (left-hand-side) column vector
a	o array of solution to Ra=v (spiking decon filter)

-----------------------------------------
	b. stoepd - symmetric Toeplitz linear system of equations Rf=g
		for f (double precision)

Notes:
	This routine does NOT solve the case when the main diagonal is
	zero, it just silently returns.


Prototype and parameters:

void stoepf (int n, double r[], double g[], double f[], double a[]);
n	i dimension of system
r	i array of top row of Toeplitz matrix
g	i array of right-hand-side column vector
f	o array of solution (left-hand-side) column vector
a	o array of solution to Ra=v (spiking decon filter)

-----------------------------------------
	c. vanded - solve Vandermonde system of equations Vx=b (double version)

Notes:
	The arrays b and x may be equivalenced.

	Adapted from Algorithm 5.6-2 in Golub, G. H., and Van Loan, C.  F.,
	1983, Matrix Computations, John-Hopkins University Press.


Prototype and parameters:

void vanded (int n, double v[], double b[], double x[]);
n	i dimension of system
v	i array of 2nd row of Vandermonde matrix (1st row is all ones)
b	i array of right-hand-side column vector
x	o array of solution column vector

-----------------------------------------
	d. vandef - solve Vandermonde system of equations Vx=b (single version)

Notes:
	The arrays b and x may be equivalenced.

	Adapted from Algorithm 5.6-2 in Golub, G. H., and Van Loan, C.  F.,
	1983, Matrix Computations, John-Hopkins University Press.


Prototype and parameters:

void vanded (int n, double v[], double b[], double x[]);
n	i dimension of system
v	i array of 2nd row of Vandermonde matrix (1st row is all ones)
b	i array of right-hand-side column vector
x	o array of solution column vector

-----------------------------------------
	e. tridif - solve a tridiagonal linear system of equations Tu=r
		(single version)

Notes:
	For example, a tridiagonal system of dimension 4 is specified as:

	    |b[0]    c[0]     0       0  | |u[0]|     |r[0]|
	    |a[1]    b[1]    c[1]     0  | |u[1]|  =  |r[1]|
	    | 0      a[2]    b[2]    c[2]| |u[2]|     |r[2]|
	    | 0       0      a[3]    b[3]| |u[3]|     |r[3]|

	The tridiagonal matrix is assumed to be non-singular.


Prototype and parameters:

void tridif (int n, float a[], float b[], float c[], float r[], float u[]);
n	i dimension of system
a	i coefficients of lower sub-diagonal of T (a[0] ignored)
b	i coefficients of diagonal of T
c	i coefficients of upper super-diagonal of T (c[n-1] ignored)
r	i array of right-hand-side column vector
u	o array of solution (left-hand-side) column vector

------------------------------------------------------------------------
19. Graphics Utility routines

	a. rfwtva - rasterize a float array as wiggle-trace-variable-area


Notes:
	The raster coordinate of the (top,left) bit in the image is
	(0,0).  In other words, x increases downward and y increases to
	the right.  Raster scan lines run from left to right, and from
	top to bottom.  Therefore, xfirst, xlast, yzmin, and yzmax
	should not be less than 0.  Likewise, yzmin and yzmax should
	not be greater than nbpr*8-1, and care should be taken to
	ensure that xfirst and xlast do not cause bits to be set
	outside (off the bottom) of the image.

	Variable area fill is performed on the right-hand (increasing
	y) side of the wiggle.  If yzmin is greater than yzmax, then z
	values between zmin will be plotted to the right of zmax, and z
	values between zbase and zmin are filled.  Swapping yzmin and
	yzmax is an easy way to reverse the polarity of a wiggle.


Prototype and parameters:

void rfwtva (int n, float z[], float zmin, float zmax, float zbase,
		int yzmin, int yzmax, int xfirst, int xlast,
		int wiggle, int nbpr, unsigned char *bits);
n	i number of samples in array to rasterize
z	i array to rasterize
zmin	i z values below zmin will be clipped
zmax	i z values above zmax will be clipped
zbase	i z values between zbase and zmax will be filled (see notes)
yzmin	i horizontal raster coordinate corresponding to zmin
yzmax	i horizontal raster coordinate corresponding to zmax
xfirst	i vertical raster coordinate of z[0] (see notes)
xlast	i vertical raster coordinate of z[n-1] (see notes)
wiggle	i =0 for no wiggle (VA only); =1 for wiggle (with VA)
nbpr	i number of bytes per row of bits
bits	b pointer to first (top,left) byte in image

-----------------------------------------
	b. scaxis - compute a readable scale for use in plotting axes

Notes:
	scaxis attempts to honor the user-specified nxnum.  However,
	nxnum will be modified if necessary for readability.  Also,
	fxnum and nxnum will be adjusted to compensate for roundoff
	error; in particular, fxnum will not be less than xmin-eps, and
	fxnum+(nxnum-1)*dxnum will not be greater than xmax+eps, where
	eps = 0.0001*(xmax-xmin).  xmin is the minimum of x1 and x2.
	xmax is the maximum of x1 and x2.


Prototype and parameters:

void scaxis (float x1, float x2, int *nxnum, float *dxnum, float *fxnum)
x1	i first x value
x2	i second x value
nxnum	i desired number of numbered values
nxnum	o number of numbered values
dxnum	o increment between numbered values (dxnum>0.0)
fxnum	o first numbered value

-----------------------------------------
	c. yclip - clip a uniformly sampled function y(x)

Notes:
	y(x) is defined by linear interpolation of the uniformly
	sampled values:  y(fx), y(fx+dx), ..., y(fx+(nx-1)*dx).

	Returns the number of samples in the clipped function.

	The output arrays xc and yc should contain space 2*nx values,
	which is the maximum possible number (nc) of xc and yc
	returned.

PARAMETERS:
int yclip (int nx, float dx, float fx, float y[], float ymin, float ymax,
		float xc[], float yc[]);
nc	f returned number of samples in output arrays xc and yc
nx	i number of x (and y) values
dx	i x sampling interval
fx	i first x
y	i array of uniformly sampled y(x) values
ymin	i minimum y value; must not be greater than ymax
ymax	i maximum y value; must not be less than ymin
xc	o array of x values for clipped y(x)
yc	o array of y values for clipped y(x)

------------------------------------------------------------------------
20. Special Function routines

	a. airya - return approximation to the Airy function Ai(x)

Notes:
	The approximation is derived from tables and formulas in
	Abramowitz and Stegun, p. 475-477.


Prototype and parameter:

float airya (float x);
x	i value at which to evaluate Ai(x)

-----------------------------------------
	b. airyb - return approximation to the Airy function Bi(x)

Notes:
	The approximation is derived from tables and formulas in
	Abramowitz and Stegun, p. 475-477.


Prototype and parameter:

float airyb (float x);
x	i value at which to evaluate Bi(x)

------------------------------------------------------------------------
21. Timing routines
	a. cpusec - return cpu time (UNIX user time) in seconds

Prototype:
float cpusec (void);

-----------------------------------------
	b. wallsec - return elapsed time (wall clock time) in seconds

Prototype:
float wallsec (void);

-----------------------------------------
	c. cputime - return cpu time (UNIX user time) in seconds using
		ANSI C built-ins

Prototype:
float cputime (void);

-----------------------------------------
	b. walltime - return elapsed time (wall clock time) in seconds using
		ANSI C built-ins

Notes:
	return value will be an integral number of seconds since t1 and t2,
	as returned by the time() intrinsic, are the number of seconds
	since the epoch


Prototype:
float walltime (void);

------------------------------------------------------------------------
22. Random Number Generating routines

	a. franuni - return a uniformly distributed random float
	b. sranuni - seed random number generator for uniform distribution

Notes:
	functions to generate a pseudo-random float uniformly
	distributed on [0,1)

	Adapted from subroutine uni in "Numerical Methods and
	Software", D. Kahaner, C. Moler, S. Nash, Prentice Hall, 1988.
	This book references, Marsaglia G., "Comments on the perfect
	uniform random number generator", Unpublished notes, Wash S. U.
	According to the reference, this random number generator
	"passes all known tests and has a period that is ...
	approximately 10^19".


Prototypes:

void sranuni (int seed);
float franuni (void);

-----------------------------------------
	c. frannor - return a normally distributed random float
	d. srannor - seed random number generator for normal distribution

Notes:
	functions to generate a pseudo-random float normally
	distributed with N(0,1); i.e., with zero mea