home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / matrix.h < prev    next >
C/C++ Source or Header  |  1994-01-13  |  20KB  |  655 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.         Type definitions for general purpose maths package
  29. */
  30.  
  31. #ifndef    MATRIXH
  32.  
  33. /* RCS id: $Id: matrix.h,v 1.15 1994/01/13 05:30:42 des Exp $ */
  34.  
  35. #define    MATRIXH    
  36.  
  37. #include    "machine.h"
  38. #include        "err.h"
  39. #include     "meminfo.h"
  40.  
  41. /* unsigned integer type */
  42. #ifndef U_INT_DEF
  43. typedef    unsigned int    u_int;
  44. #define U_INT_DEF
  45. #endif
  46.  
  47. /* vector definition */
  48. typedef    struct    {
  49.         u_int    dim, max_dim;
  50.         Real    *ve;
  51.         } VEC;
  52.  
  53. /* matrix definition */
  54. typedef    struct    {
  55.         u_int    m, n;
  56.         u_int    max_m, max_n, max_size;
  57.         Real    **me,*base;    /* base is base of alloc'd mem */
  58.         } MAT;
  59.  
  60. /* band matrix definition */
  61. typedef struct {
  62.                MAT   *mat;       /* matrix */
  63.                int   lb,ub;    /* lower and upper bandwidth */
  64.                } BAND;
  65.  
  66.  
  67. /* permutation definition */
  68. typedef    struct    {
  69.         u_int    size, max_size, *pe;
  70.         } PERM;
  71.  
  72. /* integer vector definition */
  73. typedef struct    {
  74.         u_int    dim, max_dim;
  75.         int    *ive;
  76.             } IVEC;
  77.  
  78.  
  79. #ifndef MALLOCDECL
  80. #ifndef ANSI_C
  81. extern    char    *malloc(), *calloc(), *realloc();
  82. #else
  83. extern    void    *malloc(size_t),
  84.         *calloc(size_t,size_t),
  85.         *realloc(void *,size_t);
  86. #endif
  87. #endif
  88.  
  89. #ifndef ANSI_C
  90. /* allocate one object of given type */
  91. #define    NEW(type)    ((type *)calloc(1,sizeof(type)))
  92.  
  93. /* allocate num objects of given type */
  94. #define    NEW_A(num,type)    ((type *)calloc((unsigned)(num),sizeof(type)))
  95.  
  96.  /* re-allocate arry to have num objects of the given type */
  97. #define    RENEW(var,num,type) \
  98.     ((var)=(type *)((var) ? \
  99.             realloc((char *)(var),(unsigned)(num)*sizeof(type)) : \
  100.             calloc((unsigned)(num),sizeof(type))))
  101.  
  102. #define    MEMCOPY(from,to,n_items,type) \
  103.     MEM_COPY((char *)(from),(char *)(to),(unsigned)(n_items)*sizeof(type))
  104.  
  105. #else
  106. /* allocate one object of given type */
  107. #define    NEW(type)    ((type *)calloc((size_t)1,(size_t)sizeof(type)))
  108.  
  109. /* allocate num objects of given type */
  110. #define    NEW_A(num,type)    ((type *)calloc((size_t)(num),(size_t)sizeof(type)))
  111.  
  112.  /* re-allocate arry to have num objects of the given type */
  113. #define    RENEW(var,num,type) \
  114.     ((var)=(type *)((var) ? \
  115.             realloc((char *)(var),(size_t)((num)*sizeof(type))) : \
  116.             calloc((size_t)(num),(size_t)sizeof(type))))
  117.  
  118. #define    MEMCOPY(from,to,n_items,type) \
  119.  MEM_COPY((char *)(from),(char *)(to),(unsigned)(n_items)*sizeof(type))
  120.  
  121. #endif
  122.  
  123. /* type independent min and max operations */
  124. #ifndef max
  125. #define    max(a,b)    ((a) > (b) ? (a) : (b))
  126. #endif
  127. #ifndef min
  128. #define    min(a,b)    ((a) > (b) ? (b) : (a))
  129. #endif
  130.  
  131.  
  132. #undef TRUE
  133. #define    TRUE    1
  134. #undef FALSE
  135. #define    FALSE    0
  136.  
  137.  
  138. /* for input routines */
  139. #define MAXLINE 81
  140.  
  141.  
  142. /* Dynamic memory allocation */
  143.  
  144. /* Should use M_FREE/V_FREE/PX_FREE in programs instead of m/v/px_free()
  145.    as this is considerably safer -- also provides a simple type check ! */
  146.  
  147. #ifndef ANSI_C
  148.  
  149. extern    VEC *v_get(), *v_resize();
  150. extern    MAT *m_get(), *m_resize();
  151. extern    PERM *px_get(), *px_resize();
  152. extern    IVEC *iv_get(), *iv_resize();
  153. extern    int m_free(),v_free();
  154. extern  int px_free();
  155. extern  int iv_free();
  156. extern  BAND *bd_get(), *bd_resize();
  157. extern  int bd_free();
  158.  
  159. #else
  160.  
  161. /* get/resize vector to given dimension */
  162. extern    VEC *v_get(int), *v_resize(VEC *,int);
  163. /* get/resize matrix to be m x n */
  164. extern    MAT *m_get(int,int), *m_resize(MAT *,int,int);
  165. /* get/resize permutation to have the given size */
  166. extern    PERM *px_get(int), *px_resize(PERM *,int);
  167. /* get/resize an integer vector to given dimension */
  168. extern    IVEC *iv_get(int), *iv_resize(IVEC *,int);
  169. /* get/resize a band matrix to given dimension */
  170. extern  BAND *bd_get(int,int,int), *bd_resize(BAND *,int,int,int);
  171.  
  172. /* free (de-allocate) (band) matrices, vectors, permutations and 
  173.    integer vectors */
  174. extern  int iv_free(IVEC *);
  175. extern    m_free(MAT *),v_free(VEC *),px_free(PERM *);
  176. extern   int bd_free(BAND *);
  177.  
  178. #endif
  179.  
  180.  
  181. /* MACROS */
  182.  
  183. /* macros that also check types and sets pointers to NULL */
  184. #define    M_FREE(mat)    ( m_free(mat),    (mat)=(MAT *)NULL )
  185. #define V_FREE(vec)    ( v_free(vec),    (vec)=(VEC *)NULL )
  186. #define    PX_FREE(px)    ( px_free(px),    (px)=(PERM *)NULL )
  187. #define    IV_FREE(iv)    ( iv_free(iv),    (iv)=(IVEC *)NULL )
  188.  
  189. #define MAXDIM      2001
  190.  
  191.  
  192. /* Entry level access to data structures */
  193. #ifdef DEBUG
  194.  
  195. /* returns x[i] */
  196. #define    v_entry(x,i)    (((i) < 0 || (i) >= (x)->dim) ? \
  197.              error(E_BOUNDS,"v_entry"), 0.0 : (x)->ve[i] )
  198.  
  199. /* x[i] <- val */
  200. #define    v_set_val(x,i,val) ((x)->ve[i] = ((i) < 0 || (i) >= (x)->dim) ? \
  201.                 error(E_BOUNDS,"v_set_val"), 0.0 : (val))
  202.  
  203. /* x[i] <- x[i] + val */
  204. #define    v_add_val(x,i,val) ((x)->ve[i] += ((i) < 0 || (i) >= (x)->dim) ? \
  205.                 error(E_BOUNDS,"v_set_val"), 0.0 : (val))
  206.  
  207. /* x[i] <- x[i] - val */
  208. #define    v_sub_val(x,i,val) ((x)->ve[i] -= ((i) < 0 || (i) >= (x)->dim) ? \
  209.                 error(E_BOUNDS,"v_set_val"), 0.0 : (val))
  210.  
  211. /* returns A[i][j] */
  212. #define    m_entry(A,i,j)    (((i) < 0 || (i) >= (A)->m || \
  213.               (j) < 0 || (j) >= (A)->n) ? \
  214.              error(E_BOUNDS,"m_entry"), 0.0 : (A)->me[i][j] )
  215.  
  216. /* A[i][j] <- val */
  217. #define    m_set_val(A,i,j,val) ((A)->me[i][j] = ((i) < 0 || (i) >= (A)->m || \
  218.                            (j) < 0 || (j) >= (A)->n) ? \
  219.                   error(E_BOUNDS,"m_set_val"), 0.0 : (val) )
  220.  
  221. /* A[i][j] <- A[i][j] + val */
  222. #define    m_add_val(A,i,j,val) ((A)->me[i][j] += ((i) < 0 || (i) >= (A)->m || \
  223.                         (j) < 0 || (j) >= (A)->n) ? \
  224.                   error(E_BOUNDS,"m_set_val"), 0.0 : (val) )
  225.  
  226. /* A[i][j] <- A[i][j] - val */
  227. #define    m_sub_val(A,i,j,val) ((A)->me[i][j] -= ((i) < 0 || (i) >= (A)->m || \
  228.                         (j) < 0 || (j) >= (A)->n) ? \
  229.                   error(E_BOUNDS,"m_set_val"), 0.0 : (val) )
  230. #else
  231.  
  232. /* returns x[i] */
  233. #define    v_entry(x,i)        ((x)->ve[i])
  234.  
  235. /* x[i] <- val */
  236. #define    v_set_val(x,i,val)    ((x)->ve[i]  = (val))
  237.  
  238. /* x[i] <- x[i] + val */
  239. #define    v_add_val(x,i,val)    ((x)->ve[i] += (val))
  240.  
  241.  /* x[i] <- x[i] - val */
  242. #define    v_sub_val(x,i,val)    ((x)->ve[i] -= (val))
  243.  
  244. /* returns A[i][j] */
  245. #define    m_entry(A,i,j)        ((A)->me[i][j])
  246.  
  247. /* A[i][j] <- val */
  248. #define    m_set_val(A,i,j,val)    ((A)->me[i][j]  = (val) )
  249.  
  250. /* A[i][j] <- A[i][j] + val */
  251. #define    m_add_val(A,i,j,val)    ((A)->me[i][j] += (val) )
  252.  
  253. /* A[i][j] <- A[i][j] - val */
  254. #define    m_sub_val(A,i,j,val)    ((A)->me[i][j] -= (val) )
  255.  
  256. #endif
  257.  
  258.  
  259. /* I/O routines */
  260. #ifndef ANSI_C
  261.  
  262. extern    void v_foutput(),m_foutput(),px_foutput();
  263. extern  void iv_foutput();
  264. extern    VEC *v_finput();
  265. extern    MAT *m_finput();
  266. extern    PERM *px_finput();
  267. extern    IVEC *iv_finput();
  268. extern    int fy_or_n(), fin_int(), yn_dflt(), skipjunk();
  269. extern    double fin_double();
  270.  
  271. #else
  272.  
  273. /* print x on file fp */
  274. void v_foutput(FILE *fp,VEC *x),
  275.        /* print A on file fp */
  276.     m_foutput(FILE *fp,MAT *A),
  277.        /* print px on file fp */
  278.     px_foutput(FILE *fp,PERM *px);
  279. /* print ix on file fp */
  280. void iv_foutput(FILE *fp,IVEC *ix);
  281.  
  282. /* Note: if out is NULL, then returned object is newly allocated;
  283.         Also: if out is not NULL, then that size is assumed */
  284.  
  285. /* read in vector from fp */
  286. VEC *v_finput(FILE *fp,VEC *out);
  287. /* read in matrix from fp */
  288. MAT *m_finput(FILE *fp,MAT *out);
  289. /* read in permutation from fp */
  290. PERM *px_finput(FILE *fp,PERM *out);
  291. /* read in int vector from fp */
  292. IVEC *iv_finput(FILE *fp,IVEC *out);
  293.  
  294. /* fy_or_n -- yes-or-no to question in string s
  295.         -- question written to stderr, input from fp 
  296.         -- if fp is NOT a tty then return y_n_dflt */
  297. int fy_or_n(FILE *fp,char *s);
  298.  
  299. /* yn_dflt -- sets the value of y_n_dflt to val */
  300. int yn_dflt(int val);
  301.  
  302. /* fin_int -- return integer read from file/stream fp
  303.         -- prompt s on stderr if fp is a tty
  304.         -- check that x lies between low and high: re-prompt if
  305.                 fp is a tty, error exit otherwise
  306.         -- ignore check if low > high           */
  307. int fin_int(FILE *fp,char *s,int low,int high);
  308.  
  309. /* fin_double -- return double read from file/stream fp
  310.         -- prompt s on stderr if fp is a tty
  311.         -- check that x lies between low and high: re-prompt if
  312.                 fp is a tty, error exit otherwise
  313.         -- ignore check if low > high           */
  314. double fin_double(FILE *fp,char *s,double low,double high);
  315.  
  316. /* it skips white spaces and strings of the form #....\n
  317.    Here .... is a comment string */
  318. int skipjunk(FILE *fp);
  319.  
  320. #endif
  321.  
  322.  
  323. /* MACROS */
  324.  
  325. /* macros to use stdout and stdin instead of explicit fp */
  326. #define    v_output(vec)    v_foutput(stdout,vec)
  327. #define    v_input(vec)    v_finput(stdin,vec)
  328. #define    m_output(mat)    m_foutput(stdout,mat)
  329. #define    m_input(mat)    m_finput(stdin,mat)
  330. #define    px_output(px)    px_foutput(stdout,px)
  331. #define    px_input(px)    px_finput(stdin,px)
  332. #define    iv_output(iv)    iv_foutput(stdout,iv)
  333. #define    iv_input(iv)    iv_finput(stdin,iv)
  334.  
  335. /* general purpose input routine; skips comments # ... \n */
  336. #define    finput(fp,prompt,fmt,var) \
  337.     ( ( isatty(fileno(fp)) ? fprintf(stderr,prompt) : skipjunk(fp) ), \
  338.                             fscanf(fp,fmt,var) )
  339. #define    input(prompt,fmt,var)    finput(stdin,prompt,fmt,var)
  340. #define    fprompter(fp,prompt) \
  341.     ( isatty(fileno(fp)) ? fprintf(stderr,prompt) : skipjunk(fp) )
  342. #define    prompter(prompt)    fprompter(stdin,prompt)
  343. #define    y_or_n(s)    fy_or_n(stdin,s)
  344. #define    in_int(s,lo,hi)    fin_int(stdin,s,lo,hi)
  345. #define    in_double(s,lo,hi)    fin_double(stdin,s,lo,hi)
  346.  
  347. /* Copying routines */
  348. #ifndef ANSI_C
  349. extern    MAT    *_m_copy(), *m_move(), *vm_move();
  350. extern    VEC    *_v_copy(), *v_move(), *mv_move();
  351. extern    PERM    *px_copy();
  352. extern    IVEC    *iv_copy(), *iv_move();
  353. extern  BAND    *bd_copy();
  354.  
  355. #else
  356.  
  357. /* copy in to out starting at out[i0][j0] */
  358. extern    MAT    *_m_copy(MAT *in,MAT *out,u_int i0,u_int j0),
  359.         * m_move(MAT *in, int, int, int, int, MAT *out, int, int),
  360.         *vm_move(VEC *in, int, MAT *out, int, int, int, int);
  361. /* copy in to out starting at out[i0] */
  362. extern    VEC    *_v_copy(VEC *in,VEC *out,u_int i0),
  363.         * v_move(VEC *in, int, int, VEC *out, int),
  364.         *mv_move(MAT *in, int, int, int, int, VEC *out, int);
  365. extern    PERM    *px_copy(PERM *in,PERM *out);
  366. extern    IVEC    *iv_copy(IVEC *in,IVEC *out),
  367.         *iv_move(IVEC *in, int, int, IVEC *out, int);
  368. extern  BAND    *bd_copy(BAND *in,BAND *out);
  369.  
  370. #endif
  371.  
  372.  
  373. /* MACROS */
  374. #define    m_copy(in,out)    _m_copy(in,out,0,0)
  375. #define    v_copy(in,out)    _v_copy(in,out,0)
  376.  
  377.  
  378. /* Initialisation routines -- to be zero, ones, random or identity */
  379. #ifndef ANSI_C
  380. extern    VEC     *v_zero(), *v_rand(), *v_ones();
  381. extern    MAT     *m_zero(), *m_ident(), *m_rand(), *m_ones();
  382. extern    PERM    *px_ident();
  383. extern  IVEC    *iv_zero();
  384. #else
  385. extern    VEC     *v_zero(VEC *), *v_rand(VEC *), *v_ones(VEC *);
  386. extern    MAT     *m_zero(MAT *), *m_ident(MAT *), *m_rand(MAT *),
  387.                         *m_ones(MAT *);
  388. extern    PERM    *px_ident(PERM *);
  389. extern  IVEC    *iv_zero(IVEC *);
  390. #endif
  391.  
  392. /* Basic vector operations */
  393. #ifndef ANSI_C
  394. extern    VEC *sv_mlt(), *mv_mlt(), *vm_mlt(), *v_add(), *v_sub(),
  395.         *px_vec(), *pxinv_vec(), *v_mltadd(), *v_map(), *_v_map(),
  396.         *v_lincomb(), *v_linlist();
  397. extern    double    v_min(), v_max(), v_sum();
  398. extern    VEC    *v_star(), *v_slash(), *v_sort();
  399. extern    double _in_prod(), __ip__();
  400. extern    void    __mltadd__(), __add__(), __sub__(), 
  401.                 __smlt__(), __zero__();
  402. #else
  403.  
  404. extern    VEC    *sv_mlt(double,VEC *,VEC *),    /* out <- s.x */
  405.         *mv_mlt(MAT *,VEC *,VEC *),    /* out <- A.x */
  406.         *vm_mlt(MAT *,VEC *,VEC *),    /* out^T <- x^T.A */
  407.         *v_add(VEC *,VEC *,VEC *),     /* out <- x + y */
  408.                 *v_sub(VEC *,VEC *,VEC *),    /* out <- x - y */
  409.         *px_vec(PERM *,VEC *,VEC *),    /* out <- P.x */
  410.         *pxinv_vec(PERM *,VEC *,VEC *),      /* out <- P^{-1}.x */
  411.         *v_mltadd(VEC *,VEC *,double,VEC *),   /* out <- x + s.y */
  412.         *v_map(double (*f)(double),VEC *,VEC *),  
  413.                                                  /* out[i] <- f(x[i]) */
  414.         *_v_map(double (*f)(void *,double),void *,VEC *,VEC *),
  415.         *v_lincomb(int,VEC **,Real *,VEC *),   
  416.                                                  /* out <- sum_i s[i].x[i] */
  417.                 *v_linlist(VEC *out,VEC *v1,double a1,...);
  418.                                               /* out <- s1.x1 + s2.x2 + ... */
  419.  
  420. /* returns min_j x[j] (== x[i]) */
  421. extern    double    v_min(VEC *, int *), 
  422.      /* returns max_j x[j] (== x[i]) */        
  423.         v_max(VEC *, int *), 
  424.         /* returns sum_i x[i] */
  425.         v_sum(VEC *);
  426.  
  427. /* Hadamard product: out[i] <- x[i].y[i] */
  428. extern    VEC    *v_star(VEC *, VEC *, VEC *),
  429.                  /* out[i] <- x[i] / y[i] */
  430.         *v_slash(VEC *, VEC *, VEC *),
  431.                /* sorts x, and sets order so that sorted x[i] = x[order[i]] */ 
  432.         *v_sort(VEC *, PERM *);
  433.  
  434. /* returns inner product starting at component i0 */
  435. extern    double    _in_prod(VEC *x,VEC *y,u_int i0), 
  436.                 /* returns sum_{i=0}^{len-1} x[i].y[i] */
  437.                 __ip__(Real *,Real *,int);
  438.  
  439. /* see v_mltadd(), v_add(), v_sub() and v_zero() */
  440. extern    void    __mltadd__(Real *,Real *,double,int),
  441.         __add__(Real *,Real *,Real *,int),
  442.         __sub__(Real *,Real *,Real *,int),
  443.                 __smlt__(Real *,double,Real *,int),
  444.         __zero__(Real *,int);
  445.  
  446. #endif
  447.  
  448.  
  449. /* MACRO */
  450. /* usual way of computing the inner product */
  451. #define    in_prod(a,b)    _in_prod(a,b,0)
  452.  
  453. /* Norms */
  454. /* scaled vector norms -- scale == NULL implies unscaled */
  455. #ifndef ANSI_C
  456.  
  457. extern    double    _v_norm1(), _v_norm2(), _v_norm_inf(),
  458.         m_norm1(), m_norm_inf(), m_norm_frob();
  459.  
  460. #else
  461.                /* returns sum_i |x[i]/scale[i]| */
  462. extern    double    _v_norm1(VEC *x,VEC *scale),   
  463.                /* returns (scaled) Euclidean norm */
  464.                 _v_norm2(VEC *x,VEC *scale),
  465.                /* returns max_i |x[i]/scale[i]| */
  466.         _v_norm_inf(VEC *x,VEC *scale);
  467.  
  468. /* unscaled matrix norms */
  469. extern double m_norm1(MAT *A), m_norm_inf(MAT *A), m_norm_frob(MAT *A);
  470.  
  471. #endif
  472.  
  473.  
  474. /* MACROS */
  475. /* unscaled vector norms */
  476. #define    v_norm1(x)    _v_norm1(x,VNULL)
  477. #define    v_norm2(x)    _v_norm2(x,VNULL)
  478. #define    v_norm_inf(x)    _v_norm_inf(x,VNULL)
  479.  
  480. /* Basic matrix operations */
  481. #ifndef ANSI_C
  482.  
  483. extern    MAT *sm_mlt(), *m_mlt(), *mmtr_mlt(), *mtrm_mlt(), *m_add(), *m_sub(),
  484.         *sub_mat(), *m_transp(), *ms_mltadd();
  485.  
  486. extern   BAND *bd_transp();
  487. extern    MAT *px_rows(), *px_cols(), *swap_rows(), *swap_cols(),
  488.              *_set_row(), *_set_col();
  489. extern    VEC *get_row(), *get_col(), *sub_vec(),
  490.         *mv_mltadd(), *vm_mltadd();
  491.  
  492. #else
  493.  
  494. extern    MAT    *sm_mlt(double s,MAT *A,MAT *out),     /* out <- s.A */
  495.         *m_mlt(MAT *A,MAT *B,MAT *out),    /* out <- A.B */
  496.         *mmtr_mlt(MAT *A,MAT *B,MAT *out),    /* out <- A.B^T */
  497.         *mtrm_mlt(MAT *A,MAT *B,MAT *out),    /* out <- A^T.B */
  498.         *m_add(MAT *A,MAT *B,MAT *out),    /* out <- A + B */
  499.         *m_sub(MAT *A,MAT *B,MAT *out),    /* out <- A - B */
  500.         *sub_mat(MAT *A,u_int,u_int,u_int,u_int,MAT *out),
  501.         *m_transp(MAT *A,MAT *out),        /* out <- A^T */
  502.                 /* out <- A + s.B */ 
  503.         *ms_mltadd(MAT *A,MAT *B,double s,MAT *out);   
  504.  
  505.  
  506. extern  BAND    *bd_transp(BAND *in, BAND *out);   /* out <- A^T */
  507. extern    MAT    *px_rows(PERM *px,MAT *A,MAT *out),    /* out <- P.A */
  508.         *px_cols(PERM *px,MAT *A,MAT *out),    /* out <- A.P^T */
  509.         *swap_rows(MAT *,int,int,int,int),
  510.         *swap_cols(MAT *,int,int,int,int),
  511.                  /* A[i][j] <- out[j], j >= j0 */
  512.         *_set_col(MAT *A,u_int i,VEC *out,u_int j0),
  513.                  /* A[i][j] <- out[i], i >= i0 */
  514.         *_set_row(MAT *A,u_int j,VEC *out,u_int i0);
  515.  
  516. extern    VEC    *get_row(MAT *,u_int,VEC *),
  517.         *get_col(MAT *,u_int,VEC *),
  518.         *sub_vec(VEC *,int,int,VEC *),
  519.                    /* out <- x + s.A.y */
  520.         *mv_mltadd(VEC *x,VEC *y,MAT *A,double s,VEC *out),
  521.                   /* out^T <- x^T + s.y^T.A */
  522.         *vm_mltadd(VEC *x,VEC *y,MAT *A,double s,VEC *out);
  523. #endif
  524.  
  525.  
  526. /* MACROS */
  527. /* row i of A <- vec */
  528. #define    set_row(mat,row,vec)    _set_row(mat,row,vec,0) 
  529. /* col j of A <- vec */
  530. #define    set_col(mat,col,vec)    _set_col(mat,col,vec,0)
  531.  
  532.  
  533. /* Basic permutation operations */
  534. #ifndef ANSI_C
  535.  
  536. extern    PERM *px_mlt(), *px_inv(), *px_transp();
  537. extern    int  px_sign();
  538.  
  539. #else
  540.  
  541. extern    PERM    *px_mlt(PERM *px1,PERM *px2,PERM *out),    /* out <- px1.px2 */
  542.         *px_inv(PERM *px,PERM *out),    /* out <- px^{-1} */
  543.                  /* swap px[i] and px[j] */
  544.         *px_transp(PERM *px,u_int i,u_int j);
  545.  
  546.      /* returns sign(px) = +1 if px product of even # transpositions
  547.                            -1 if ps product of odd  # transpositions */
  548. extern    int    px_sign(PERM *);
  549.  
  550. #endif
  551.  
  552.  
  553. /* Basic integer vector operations */
  554. #ifndef ANSI_C
  555.  
  556. extern    IVEC    *iv_add(), *iv_sub(), *iv_sort();
  557.  
  558. #else
  559.  
  560. extern    IVEC    *iv_add(IVEC *ix,IVEC *iy,IVEC *out),  /* out <- ix + iy */
  561.         *iv_sub(IVEC *ix,IVEC *iy,IVEC *out),  /* out <- ix - iy */
  562.         /* sorts ix & sets order so that sorted ix[i] = old ix[order[i]] */
  563.         *iv_sort(IVEC *ix, PERM *order);
  564.  
  565. #endif
  566.  
  567.  
  568. /* miscellaneous functions */
  569.  
  570. #ifndef ANSI_C
  571.  
  572. extern    double    square(), cube(), mrand();
  573. extern    void    smrand(), mrandlist();
  574. extern  void    m_dump(), px_dump(), v_dump(), iv_dump();
  575. extern MAT *band2mat();
  576. extern BAND *mat2band();
  577.  
  578. #else
  579.  
  580. double    square(double x),     /* returns x^2 */
  581.   cube(double x),         /* returns x^3 */
  582.   mrand(void);                  /* returns random # in [0,1) */
  583.  
  584. void    smrand(int seed),            /* seeds mrand() */
  585.   mrandlist(Real *x, int len);       /* generates len random numbers */
  586.  
  587. void    m_dump(FILE *fp,MAT *a), px_dump(FILE *,PERM *px),
  588.         v_dump(FILE *fp,VEC *x), iv_dump(FILE *fp, IVEC *ix);
  589.  
  590. MAT *band2mat(BAND *bA, MAT *A);
  591. BAND *mat2band(MAT *A, int lb,int ub, BAND *bA);
  592.  
  593. #endif
  594.  
  595.  
  596. /* miscellaneous constants */
  597. #define    VNULL    ((VEC *)NULL)
  598. #define    MNULL    ((MAT *)NULL)
  599. #define    PNULL    ((PERM *)NULL)
  600. #define    IVNULL    ((IVEC *)NULL)
  601. #define BDNULL  ((BAND *)NULL)
  602.  
  603.  
  604.  
  605. /* varying number of arguments */
  606.  
  607. #ifdef ANSI_C
  608. #include <stdarg.h>
  609.  
  610. /* prototypes */
  611.  
  612. int v_get_vars(int dim,...);
  613. int iv_get_vars(int dim,...);
  614. int m_get_vars(int m,int n,...);
  615. int px_get_vars(int dim,...);
  616.  
  617. int v_resize_vars(int new_dim,...);
  618. int iv_resize_vars(int new_dim,...);
  619. int m_resize_vars(int m,int n,...);
  620. int px_resize_vars(int new_dim,...);
  621.  
  622. int v_free_vars(VEC **,...);
  623. int iv_free_vars(IVEC **,...);
  624. int px_free_vars(PERM **,...);
  625. int m_free_vars(MAT **,...);
  626.  
  627. #elif VARARGS
  628. /* old varargs is used */
  629.  
  630. #include  <varargs.h>
  631.  
  632. /* prototypes */
  633.  
  634. int v_get_vars();
  635. int iv_get_vars();
  636. int m_get_vars();
  637. int px_get_vars();
  638.  
  639. int v_resize_vars();
  640. int iv_resize_vars();
  641. int m_resize_vars();
  642. int px_resize_vars();
  643.  
  644. int v_free_vars();
  645. int iv_free_vars();
  646. int px_free_vars();
  647. int m_free_vars();
  648.  
  649. #endif
  650.  
  651.  
  652. #endif
  653.  
  654.  
  655.