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