home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / init.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  6KB  |  299 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*
  28.     This is a file of routines for zero-ing, and initialising
  29.     vectors, matrices and permutations.
  30.     This is to be included in the matrix.a library
  31. */
  32.  
  33. static    char    rcsid[] = "$Id: init.c,v 1.6 1994/01/13 05:36:58 des Exp $";
  34.  
  35. #include    <stdio.h>
  36. #include    "matrix.h"
  37.  
  38. /* v_zero -- zero the vector x */
  39. VEC    *v_zero(x)
  40. VEC    *x;
  41. {
  42.     if ( x == VNULL )
  43.         error(E_NULL,"v_zero");
  44.  
  45.     __zero__(x->ve,x->dim);
  46.     /* for ( i = 0; i < x->dim; i++ )
  47.         x->ve[i] = 0.0; */
  48.  
  49.     return x;
  50. }
  51.  
  52.  
  53. /* iv_zero -- zero the vector ix */
  54. IVEC    *iv_zero(ix)
  55. IVEC    *ix;
  56. {
  57.    int i;
  58.    
  59.    if ( ix == IVNULL )
  60.      error(E_NULL,"iv_zero");
  61.    
  62.    for ( i = 0; i < ix->dim; i++ )
  63.      ix->ive[i] = 0; 
  64.    
  65.    return ix;
  66. }
  67.  
  68.  
  69. /* m_zero -- zero the matrix A */
  70. MAT    *m_zero(A)
  71. MAT    *A;
  72. {
  73.     int    i, A_m, A_n;
  74.     Real    **A_me;
  75.  
  76.     if ( A == MNULL )
  77.         error(E_NULL,"m_zero");
  78.  
  79.     A_m = A->m;    A_n = A->n;    A_me = A->me;
  80.     for ( i = 0; i < A_m; i++ )
  81.         __zero__(A_me[i],A_n);
  82.         /* for ( j = 0; j < A_n; j++ )
  83.             A_me[i][j] = 0.0; */
  84.  
  85.     return A;
  86. }
  87.  
  88. /* mat_id -- set A to being closest to identity matrix as possible
  89.     -- i.e. A[i][j] == 1 if i == j and 0 otherwise */
  90. MAT    *m_ident(A)
  91. MAT    *A;
  92. {
  93.     int    i, size;
  94.  
  95.     if ( A == MNULL )
  96.         error(E_NULL,"m_ident");
  97.  
  98.     m_zero(A);
  99.     size = min(A->m,A->n);
  100.     for ( i = 0; i < size; i++ )
  101.         A->me[i][i] = 1.0;
  102.  
  103.     return A;
  104. }
  105.     
  106. /* px_ident -- set px to identity permutation */
  107. PERM    *px_ident(px)
  108. PERM    *px;
  109. {
  110.     int    i, px_size;
  111.     u_int    *px_pe;
  112.  
  113.     if ( px == PNULL )
  114.         error(E_NULL,"px_ident");
  115.  
  116.     px_size = px->size;    px_pe = px->pe;
  117.     for ( i = 0; i < px_size; i++ )
  118.         px_pe[i] = i;
  119.  
  120.     return px;
  121. }
  122.  
  123. /* Pseudo random number generator data structures */
  124. /* Knuth's lagged Fibonacci-based generator: See "Seminumerical Algorithms:
  125.    The Art of Computer Programming" sections 3.2-3.3 */
  126.  
  127. #ifdef ANSI_C
  128. #ifndef LONG_MAX
  129. #include    <limits.h>
  130. #endif
  131. #endif
  132.  
  133. #ifdef LONG_MAX
  134. #define MODULUS    LONG_MAX
  135. #else
  136. #define MODULUS    1000000000L    /* assuming long's at least 32 bits long */
  137. #endif
  138. #define MZ    0L
  139.  
  140. static long mrand_list[56];
  141. static int  started = FALSE;
  142. static int  inext = 0, inextp = 31;
  143.  
  144.  
  145. /* mrand -- pseudo-random number generator */
  146. #ifdef ANSI_C
  147. double mrand(void)
  148. #else
  149. double mrand()
  150. #endif
  151. {
  152.     long    lval;
  153.     static Real  factor = 1.0/((Real)MODULUS);
  154.     
  155.     if ( ! started )
  156.     smrand(3127);
  157.     
  158.     inext = (inext >= 54) ? 0 : inext+1;
  159.     inextp = (inextp >= 54) ? 0 : inextp+1;
  160.  
  161.     lval = mrand_list[inext]-mrand_list[inextp];
  162.     if ( lval < 0L )
  163.     lval += MODULUS;
  164.     mrand_list[inext] = lval;
  165.     
  166.     return (double)lval*factor;
  167. }
  168.  
  169. /* mrandlist -- fills the array a[] with len random numbers */
  170. void    mrandlist(a, len)
  171. Real    a[];
  172. int    len;
  173. {
  174.     int        i;
  175.     long    lval;
  176.     static Real  factor = 1.0/((Real)MODULUS);
  177.     
  178.     if ( ! started )
  179.     smrand(3127);
  180.     
  181.     for ( i = 0; i < len; i++ )
  182.     {
  183.     inext = (inext >= 54) ? 0 : inext+1;
  184.     inextp = (inextp >= 54) ? 0 : inextp+1;
  185.     
  186.     lval = mrand_list[inext]-mrand_list[inextp];
  187.     if ( lval < 0L )
  188.         lval += MODULUS;
  189.     mrand_list[inext] = lval;
  190.     
  191.     a[i] = (Real)lval*factor;
  192.     }
  193. }
  194.  
  195.  
  196. /* smrand -- set seed for mrand() */
  197. void smrand(seed)
  198. int    seed;
  199. {
  200.     int        i;
  201.  
  202.     mrand_list[0] = (123413*seed) % MODULUS;
  203.     for ( i = 1; i < 55; i++ )
  204.     mrand_list[i] = (123413*mrand_list[i-1]) % MODULUS;
  205.  
  206.     started = TRUE;
  207.  
  208.     /* run mrand() through the list sufficient times to
  209.        thoroughly randomise the array */
  210.     for ( i = 0; i < 55*55; i++ )
  211.     mrand();
  212. }
  213. #undef MODULUS
  214. #undef MZ
  215. #undef FAC
  216.  
  217. /* v_rand -- initialises x to be a random vector, components
  218.     independently & uniformly ditributed between 0 and 1 */
  219. VEC    *v_rand(x)
  220. VEC    *x;
  221. {
  222.     /* int    i; */
  223.  
  224.     if ( ! x )
  225.         error(E_NULL,"v_rand");
  226.  
  227.     /* for ( i = 0; i < x->dim; i++ ) */
  228.         /* x->ve[i] = rand()/((Real)MAX_RAND); */
  229.         /* x->ve[i] = mrand(); */
  230.     mrandlist(x->ve,x->dim);
  231.  
  232.     return x;
  233. }
  234.  
  235. /* m_rand -- initialises A to be a random vector, components
  236.     independently & uniformly distributed between 0 and 1 */
  237. MAT    *m_rand(A)
  238. MAT    *A;
  239. {
  240.     int    i /* , j */;
  241.  
  242.     if ( ! A )
  243.         error(E_NULL,"m_rand");
  244.  
  245.     for ( i = 0; i < A->m; i++ )
  246.         /* for ( j = 0; j < A->n; j++ ) */
  247.             /* A->me[i][j] = rand()/((Real)MAX_RAND); */
  248.             /* A->me[i][j] = mrand(); */
  249.         mrandlist(A->me[i],A->n);
  250.  
  251.     return A;
  252. }
  253.  
  254. /* v_ones -- fills x with one's */
  255. VEC    *v_ones(x)
  256. VEC    *x;
  257. {
  258.     int    i;
  259.  
  260.     if ( ! x )
  261.         error(E_NULL,"v_ones");
  262.  
  263.     for ( i = 0; i < x->dim; i++ )
  264.         x->ve[i] = 1.0;
  265.  
  266.     return x;
  267. }
  268.  
  269. /* m_ones -- fills matrix with one's */
  270. MAT    *m_ones(A)
  271. MAT    *A;
  272. {
  273.     int    i, j;
  274.  
  275.     if ( ! A )
  276.         error(E_NULL,"m_ones");
  277.  
  278.     for ( i = 0; i < A->m; i++ )
  279.         for ( j = 0; j < A->n; j++ )
  280.             A->me[i][j] = 1.0;
  281.  
  282.     return A;
  283. }
  284.  
  285. /* v_count -- initialises x so that x->ve[i] == i */
  286. VEC    *v_count(x)
  287. VEC    *x;
  288. {
  289.     int    i;
  290.  
  291.     if ( ! x )
  292.         error(E_NULL,"v_count");
  293.  
  294.     for ( i = 0; i < x->dim; i++ )
  295.         x->ve[i] = (Real)i;
  296.  
  297.     return x;
  298. }
  299.