init.c

/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**			     Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/


/*
	This is a file of routines for zero-ing, and initialising
	vectors, matrices and permutations.
	This is to be included in the matrix.a library
*/

static	char	rcsid[] = "$Id: init.c,v 1.6 1994/01/13 05:36:58 des Exp $";

#include	<stdio.h>
#include	"matrix.h"

/* v_zero -- zero the vector x */
#ifndef ANSI_C
VEC	*v_zero(x)
VEC	*x;
#else
VEC	*v_zero(VEC *x)
#endif
{
	if ( x == VNULL )
		error(E_NULL,"v_zero");

	__zero__(x->ve,x->dim);
	/* for ( i = 0; i < x->dim; i++ )
		x->ve[i] = 0.0; */

	return x;
}


/* iv_zero -- zero the vector ix */
#ifndef ANSI_C
IVEC	*iv_zero(ix)
IVEC	*ix;
#else
IVEC	*iv_zero(IVEC *ix)
#endif
{
   int i;
   
   if ( ix == IVNULL )
     error(E_NULL,"iv_zero");
   
   for ( i = 0; i < ix->dim; i++ )
     ix->ive[i] = 0; 
   
   return ix;
}


/* m_zero -- zero the matrix A */
#ifndef ANSI_C
MAT	*m_zero(A)
MAT	*A;
#else
MAT	*m_zero(MAT *A)
#endif
{
	int	i, A_m, A_n;
	Real	**A_me;

	if ( A == MNULL )
		error(E_NULL,"m_zero");

	A_m = A->m;	A_n = A->n;	A_me = A->me;
	for ( i = 0; i < A_m; i++ )
		__zero__(A_me[i],A_n);
		/* for ( j = 0; j < A_n; j++ )
			A_me[i][j] = 0.0; */

	return A;
}

/* mat_id -- set A to being closest to identity matrix as possible
	-- i.e. A[i][j] == 1 if i == j and 0 otherwise */
#ifndef ANSI_C
MAT	*m_ident(A)
MAT	*A;
#else
MAT	*m_ident(MAT *A)
#endif
{
	int	i, size;

	if ( A == MNULL )
		error(E_NULL,"m_ident");

	m_zero(A);
	size = min(A->m,A->n);
	for ( i = 0; i < size; i++ )
		A->me[i][i] = 1.0;

	return A;
}
	
/* px_ident -- set px to identity permutation */
#ifndef ANSI_C
PERM	*px_ident(px)
PERM	*px;
#else
PERM	*px_ident(PERM *px)
#endif
{
	int	i, px_size;
	unsigned int	*px_pe;

	if ( px == PNULL )
		error(E_NULL,"px_ident");

	px_size = px->size;	px_pe = px->pe;
	for ( i = 0; i < px_size; i++ )
		px_pe[i] = i;

	return px;
}

/* Pseudo random number generator data structures */
/* Knuth's lagged Fibonacci-based generator: See "Seminumerical Algorithms:
   The Art of Computer Programming" sections 3.2-3.3 */

#ifdef ANSI_C
#ifndef LONG_MAX
#include	<limits.h>
#endif
#endif

#ifdef LONG_MAX
#define MODULUS	LONG_MAX
#else
#define MODULUS	1000000000L	/* assuming long's at least 32 bits long */
#endif
#define MZ	0L

static long mrand_list[56];
static int  started = FALSE;
static int  inext = 0, inextp = 31;


/* mrand -- pseudo-random number generator */
#ifdef ANSI_C
double mrand(void)
#else
double mrand()
#endif
{
    long	lval;
    static Real  factor = 1.0/((Real)MODULUS);
    
    if ( ! started )
	smrand(3127);
    
    inext = (inext >= 54) ? 0 : inext+1;
    inextp = (inextp >= 54) ? 0 : inextp+1;

    lval = mrand_list[inext]-mrand_list[inextp];
    if ( lval < 0L )
	lval += MODULUS;
    mrand_list[inext] = lval;
    
    return (double)lval*factor;
}

/* mrandlist -- fills the array a[] with len random numbers */
#ifndef ANSI_C
void	mrandlist(a, len)
Real	a[];
int	len;
#else
void	mrandlist(Real a[], int len)
#endif
{
    int		i;
    long	lval;
    static Real  factor = 1.0/((Real)MODULUS);
    
    if ( ! started )
	smrand(3127);
    
    for ( i = 0; i < len; i++ )
    {
	inext = (inext >= 54) ? 0 : inext+1;
	inextp = (inextp >= 54) ? 0 : inextp+1;
	
	lval = mrand_list[inext]-mrand_list[inextp];
	if ( lval < 0L )
	    lval += MODULUS;
	mrand_list[inext] = lval;
	
	a[i] = (Real)lval*factor;
    }
}


/* smrand -- set seed for mrand() */
#ifndef ANSI_C
void smrand(seed)
int	seed;
#else
void smrand(int seed)
#endif
{
    int		i;

    mrand_list[0] = (123413*seed) % MODULUS;
    for ( i = 1; i < 55; i++ )
	mrand_list[i] = (123413*mrand_list[i-1]) % MODULUS;

    started = TRUE;

    /* run mrand() through the list sufficient times to
       thoroughly randomise the array */
    for ( i = 0; i < 55*55; i++ )
	mrand();
}
#undef MODULUS
#undef MZ
#undef FAC

/* v_rand -- initialises x to be a random vector, components
	independently & uniformly ditributed between 0 and 1 */
#ifndef ANSI_C
VEC	*v_rand(x)
VEC	*x;
#else
VEC	*v_rand(VEC *x)
#endif
{
	/* int	i; */

	if ( ! x )
		error(E_NULL,"v_rand");

	/* for ( i = 0; i < x->dim; i++ ) */
	    /* x->ve[i] = rand()/((Real)MAX_RAND); */
	    /* x->ve[i] = mrand(); */
	mrandlist(x->ve,x->dim);

	return x;
}

/* m_rand -- initialises A to be a random vector, components
	independently & uniformly distributed between 0 and 1 */
#ifndef ANSI_C
MAT	*m_rand(A)
MAT	*A;
#else
MAT	*m_rand(MAT *A)
#endif
{
	int	i /* , j */;

	if ( ! A )
		error(E_NULL,"m_rand");

	for ( i = 0; i < A->m; i++ )
		/* for ( j = 0; j < A->n; j++ ) */
		    /* A->me[i][j] = rand()/((Real)MAX_RAND); */
		    /* A->me[i][j] = mrand(); */
	    mrandlist(A->me[i],A->n);

	return A;
}

/* v_ones -- fills x with one's */
#ifndef ANSI_C
VEC	*v_ones(x)
VEC	*x;
#else
VEC	*v_ones(VEC *x)
#endif
{
	int	i;

	if ( ! x )
		error(E_NULL,"v_ones");

	for ( i = 0; i < x->dim; i++ )
		x->ve[i] = 1.0;

	return x;
}

/* m_ones -- fills matrix with one's */
#ifndef ANSI_C
MAT	*m_ones(A)
MAT	*A;
#else
MAT	*m_ones(MAT *A)
#endif
{
	int	i, j;

	if ( ! A )
		error(E_NULL,"m_ones");

	for ( i = 0; i < A->m; i++ )
		for ( j = 0; j < A->n; j++ )
		    A->me[i][j] = 1.0;

	return A;
}

/* v_count -- initialises x so that x->ve[i] == i */
#ifndef ANSI_C
VEC	*v_count(x)
VEC	*x;
#else
VEC	*v_count(VEC *x)
#endif
{
	int	i;

	if ( ! x )
	    error(E_NULL,"v_count");

	for ( i = 0; i < x->dim; i++ )
	    x->ve[i] = (Real)i;

	return x;
}