init.c

Go to the documentation of this file.

#include <../../nrnconf.h>

/**************************************************************************
**
** 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[] = "init.c,v 1.1 1997/12/04 17:55:25 hines Exp";

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

/* v_zero -- zero the vector x */
VEC *v_zero(x)
VEC *x;
{
    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 */
IVEC    *iv_zero(ix)
IVEC    *ix;
{
   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 */
MAT *m_zero(A)
MAT *A;
{
    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 */
MAT *m_ident(A)
MAT *A;
{
    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 */
PERM    *px_ident(px)
PERM    *px;
{
    int i, px_size;
    u_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 */
void    mrandlist(a, len)
Real    a[];
int len;
{
    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() */
void smrand(seed)
int seed;
{
    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 */
VEC *v_rand(x)
VEC *x;
{
    /* 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 */
MAT *m_rand(A)
MAT *A;
{
    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 */
VEC *v_ones(x)
VEC *x;
{
    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 */
MAT *m_ones(A)
MAT *A;
{
    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 */
VEC *v_count(x)
VEC *x;
{
    int i;

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

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

    return x;
}