/************************************************************************** ** ** 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. ** ***************************************************************************/ /* Type definitions for general purpose maths package */ #ifndef MATRIXH /* RCS id: $Id: matrix.h,v 1.18 1994/04/16 00:33:37 des Exp $ */ #define MATRIXH #include "machine.h" #include "err.h" #include "meminfo.h" /* unsigned integer type */ #ifndef U_INT_DEF typedef unsigned int u_int; #define U_INT_DEF #endif /* vector definition */ typedef struct { u_int dim, max_dim; Real *ve; } VEC; /* matrix definition */ typedef struct { u_int m, n; u_int max_m, max_n, max_size; Real **me,*base; /* base is base of alloc'd mem */ } MAT; /* band matrix definition */ typedef struct { MAT *mat; /* matrix */ int lb,ub; /* lower and upper bandwidth */ } BAND; /* permutation definition */ typedef struct { u_int size, max_size, *pe; } PERM; /* integer vector definition */ typedef struct { u_int dim, max_dim; int *ive; } IVEC; #ifndef MALLOCDECL #ifndef ANSI_C extern char *malloc(), *calloc(), *realloc(); #else extern void *malloc(size_t), *calloc(size_t,size_t), *realloc(void *,size_t); #endif #endif #ifndef ANSI_C extern void m_version(); #else void m_version( void ); #endif #ifndef ANSI_C /* allocate one object of given type */ #define NEW(type) ((type *)calloc(1,sizeof(type))) /* allocate num objects of given type */ #define NEW_A(num,type) ((type *)calloc((unsigned)(num),sizeof(type))) /* re-allocate arry to have num objects of the given type */ #define RENEW(var,num,type) \ ((var)=(type *)((var) ? \ realloc((char *)(var),(unsigned)(num)*sizeof(type)) : \ calloc((unsigned)(num),sizeof(type)))) #define MEMCOPY(from,to,n_items,type) \ MEM_COPY((char *)(from),(char *)(to),(unsigned)(n_items)*sizeof(type)) #else /* allocate one object of given type */ #define NEW(type) ((type *)calloc((size_t)1,(size_t)sizeof(type))) /* allocate num objects of given type */ #define NEW_A(num,type) ((type *)calloc((size_t)(num),(size_t)sizeof(type))) /* re-allocate arry to have num objects of the given type */ #define RENEW(var,num,type) \ ((var)=(type *)((var) ? \ realloc((char *)(var),(size_t)((num)*sizeof(type))) : \ calloc((size_t)(num),(size_t)sizeof(type)))) #define MEMCOPY(from,to,n_items,type) \ MEM_COPY((char *)(from),(char *)(to),(unsigned)(n_items)*sizeof(type)) #endif /* type independent min and max operations */ #ifndef max #define max(a,b) ((a) > (b) ? (a) : (b)) #endif #ifndef min #define min(a,b) ((a) > (b) ? (b) : (a)) #endif #undef TRUE #define TRUE 1 #undef FALSE #define FALSE 0 /* for input routines */ #define MAXLINE 81 /* Dynamic memory allocation */ /* Should use M_FREE/V_FREE/PX_FREE in programs instead of m/v/px_free() as this is considerably safer -- also provides a simple type check ! */ #ifndef ANSI_C extern VEC *v_get(), *v_resize(); extern MAT *m_get(), *m_resize(); extern PERM *px_get(), *px_resize(); extern IVEC *iv_get(), *iv_resize(); extern int m_free(),v_free(); extern int px_free(); extern int iv_free(); extern BAND *bd_get(), *bd_resize(); extern int bd_free(); #else /* get/resize vector to given dimension */ extern VEC *v_get(int), *v_resize(VEC *,int); /* get/resize matrix to be m x n */ extern MAT *m_get(int,int), *m_resize(MAT *,int,int); /* get/resize permutation to have the given size */ extern PERM *px_get(int), *px_resize(PERM *,int); /* get/resize an integer vector to given dimension */ extern IVEC *iv_get(int), *iv_resize(IVEC *,int); /* get/resize a band matrix to given dimension */ extern BAND *bd_get(int,int,int), *bd_resize(BAND *,int,int,int); /* free (de-allocate) (band) matrices, vectors, permutations and integer vectors */ extern int iv_free(IVEC *); extern m_free(MAT *),v_free(VEC *),px_free(PERM *); extern int bd_free(BAND *); #endif /* MACROS */ /* macros that also check types and sets pointers to NULL */ #define M_FREE(mat) ( m_free(mat), (mat)=(MAT *)NULL ) #define V_FREE(vec) ( v_free(vec), (vec)=(VEC *)NULL ) #define PX_FREE(px) ( px_free(px), (px)=(PERM *)NULL ) #define IV_FREE(iv) ( iv_free(iv), (iv)=(IVEC *)NULL ) #define MAXDIM 2001 /* Entry level access to data structures */ #ifdef DEBUG /* returns x[i] */ #define v_entry(x,i) (((i) < 0 || (i) >= (x)->dim) ? \ error(E_BOUNDS,"v_entry"), 0.0 : (x)->ve[i] ) /* x[i] <- val */ #define v_set_val(x,i,val) ((x)->ve[i] = ((i) < 0 || (i) >= (x)->dim) ? \ error(E_BOUNDS,"v_set_val"), 0.0 : (val)) /* x[i] <- x[i] + val */ #define v_add_val(x,i,val) ((x)->ve[i] += ((i) < 0 || (i) >= (x)->dim) ? \ error(E_BOUNDS,"v_add_val"), 0.0 : (val)) /* x[i] <- x[i] - val */ #define v_sub_val(x,i,val) ((x)->ve[i] -= ((i) < 0 || (i) >= (x)->dim) ? \ error(E_BOUNDS,"v_sub_val"), 0.0 : (val)) /* returns A[i][j] */ #define m_entry(A,i,j) (((i) < 0 || (i) >= (A)->m || \ (j) < 0 || (j) >= (A)->n) ? \ error(E_BOUNDS,"m_entry"), 0.0 : (A)->me[i][j] ) /* A[i][j] <- val */ #define m_set_val(A,i,j,val) ((A)->me[i][j] = ((i) < 0 || (i) >= (A)->m || \ (j) < 0 || (j) >= (A)->n) ? \ error(E_BOUNDS,"m_set_val"), 0.0 : (val) ) /* A[i][j] <- A[i][j] + val */ #define m_add_val(A,i,j,val) ((A)->me[i][j] += ((i) < 0 || (i) >= (A)->m || \ (j) < 0 || (j) >= (A)->n) ? \ error(E_BOUNDS,"m_add_val"), 0.0 : (val) ) /* A[i][j] <- A[i][j] - val */ #define m_sub_val(A,i,j,val) ((A)->me[i][j] -= ((i) < 0 || (i) >= (A)->m || \ (j) < 0 || (j) >= (A)->n) ? \ error(E_BOUNDS,"m_sub_val"), 0.0 : (val) ) #else /* returns x[i] */ #define v_entry(x,i) ((x)->ve[i]) /* x[i] <- val */ #define v_set_val(x,i,val) ((x)->ve[i] = (val)) /* x[i] <- x[i] + val */ #define v_add_val(x,i,val) ((x)->ve[i] += (val)) /* x[i] <- x[i] - val */ #define v_sub_val(x,i,val) ((x)->ve[i] -= (val)) /* returns A[i][j] */ #define m_entry(A,i,j) ((A)->me[i][j]) /* A[i][j] <- val */ #define m_set_val(A,i,j,val) ((A)->me[i][j] = (val) ) /* A[i][j] <- A[i][j] + val */ #define m_add_val(A,i,j,val) ((A)->me[i][j] += (val) ) /* A[i][j] <- A[i][j] - val */ #define m_sub_val(A,i,j,val) ((A)->me[i][j] -= (val) ) #endif /* I/O routines */ #ifndef ANSI_C extern void v_foutput(),m_foutput(),px_foutput(); extern void iv_foutput(); extern VEC *v_finput(); extern MAT *m_finput(); extern PERM *px_finput(); extern IVEC *iv_finput(); extern int fy_or_n(), fin_int(), yn_dflt(), skipjunk(); extern double fin_double(); #else /* print x on file fp */ void v_foutput(FILE *fp,VEC *x), /* print A on file fp */ m_foutput(FILE *fp,MAT *A), /* print px on file fp */ px_foutput(FILE *fp,PERM *px); /* print ix on file fp */ void iv_foutput(FILE *fp,IVEC *ix); /* Note: if out is NULL, then returned object is newly allocated; Also: if out is not NULL, then that size is assumed */ /* read in vector from fp */ VEC *v_finput(FILE *fp,VEC *out); /* read in matrix from fp */ MAT *m_finput(FILE *fp,MAT *out); /* read in permutation from fp */ PERM *px_finput(FILE *fp,PERM *out); /* read in int vector from fp */ IVEC *iv_finput(FILE *fp,IVEC *out); /* fy_or_n -- yes-or-no to question in string s -- question written to stderr, input from fp -- if fp is NOT a tty then return y_n_dflt */ int fy_or_n(FILE *fp,char *s); /* yn_dflt -- sets the value of y_n_dflt to val */ int yn_dflt(int val); /* fin_int -- return integer read from file/stream fp -- prompt s on stderr if fp is a tty -- check that x lies between low and high: re-prompt if fp is a tty, error exit otherwise -- ignore check if low > high */ int fin_int(FILE *fp,char *s,int low,int high); /* fin_double -- return double read from file/stream fp -- prompt s on stderr if fp is a tty -- check that x lies between low and high: re-prompt if fp is a tty, error exit otherwise -- ignore check if low > high */ double fin_double(FILE *fp,char *s,double low,double high); /* it skips white spaces and strings of the form #....\n Here .... is a comment string */ int skipjunk(FILE *fp); #endif /* MACROS */ /* macros to use stdout and stdin instead of explicit fp */ #define v_output(vec) v_foutput(stdout,vec) #define v_input(vec) v_finput(stdin,vec) #define m_output(mat) m_foutput(stdout,mat) #define m_input(mat) m_finput(stdin,mat) #define px_output(px) px_foutput(stdout,px) #define px_input(px) px_finput(stdin,px) #define iv_output(iv) iv_foutput(stdout,iv) #define iv_input(iv) iv_finput(stdin,iv) /* general purpose input routine; skips comments # ... \n */ #define finput(fp,prompt,fmt,var) \ ( ( isatty(fileno(fp)) ? fprintf(stderr,prompt) : skipjunk(fp) ), \ fscanf(fp,fmt,var) ) #define input(prompt,fmt,var) finput(stdin,prompt,fmt,var) #define fprompter(fp,prompt) \ ( isatty(fileno(fp)) ? fprintf(stderr,prompt) : skipjunk(fp) ) #define prompter(prompt) fprompter(stdin,prompt) #define y_or_n(s) fy_or_n(stdin,s) #define in_int(s,lo,hi) fin_int(stdin,s,lo,hi) #define in_double(s,lo,hi) fin_double(stdin,s,lo,hi) /* Copying routines */ #ifndef ANSI_C extern MAT *_m_copy(), *m_move(), *vm_move(); extern VEC *_v_copy(), *v_move(), *mv_move(); extern PERM *px_copy(); extern IVEC *iv_copy(), *iv_move(); extern BAND *bd_copy(); #else /* copy in to out starting at out[i0][j0] */ extern MAT *_m_copy(MAT *in,MAT *out,u_int i0,u_int j0), * m_move(MAT *in, int, int, int, int, MAT *out, int, int), *vm_move(VEC *in, int, MAT *out, int, int, int, int); /* copy in to out starting at out[i0] */ extern VEC *_v_copy(VEC *in,VEC *out,u_int i0), * v_move(VEC *in, int, int, VEC *out, int), *mv_move(MAT *in, int, int, int, int, VEC *out, int); extern PERM *px_copy(PERM *in,PERM *out); extern IVEC *iv_copy(IVEC *in,IVEC *out), *iv_move(IVEC *in, int, int, IVEC *out, int); extern BAND *bd_copy(BAND *in,BAND *out); #endif /* MACROS */ #define m_copy(in,out) _m_copy(in,out,0,0) #define v_copy(in,out) _v_copy(in,out,0) /* Initialisation routines -- to be zero, ones, random or identity */ #ifndef ANSI_C extern VEC *v_zero(), *v_rand(), *v_ones(); extern MAT *m_zero(), *m_ident(), *m_rand(), *m_ones(); extern PERM *px_ident(); extern IVEC *iv_zero(); #else extern VEC *v_zero(VEC *), *v_rand(VEC *), *v_ones(VEC *); extern MAT *m_zero(MAT *), *m_ident(MAT *), *m_rand(MAT *), *m_ones(MAT *); extern PERM *px_ident(PERM *); extern IVEC *iv_zero(IVEC *); #endif /* Basic vector operations */ #ifndef ANSI_C extern VEC *sv_mlt(), *mv_mlt(), *vm_mlt(), *v_add(), *v_sub(), *px_vec(), *pxinv_vec(), *v_mltadd(), *v_map(), *_v_map(), *v_lincomb(), *v_linlist(); extern double v_min(), v_max(), v_sum(); extern VEC *v_star(), *v_slash(), *v_sort(); extern double _in_prod(), __ip__(); extern void __mltadd__(), __add__(), __sub__(), __smlt__(), __zero__(); #else extern VEC *sv_mlt(double,VEC *,VEC *), /* out <- s.x */ *mv_mlt(MAT *,VEC *,VEC *), /* out <- A.x */ *vm_mlt(MAT *,VEC *,VEC *), /* out^T <- x^T.A */ *v_add(VEC *,VEC *,VEC *), /* out <- x + y */ *v_sub(VEC *,VEC *,VEC *), /* out <- x - y */ *px_vec(PERM *,VEC *,VEC *), /* out <- P.x */ *pxinv_vec(PERM *,VEC *,VEC *), /* out <- P^{-1}.x */ *v_mltadd(VEC *,VEC *,double,VEC *), /* out <- x + s.y */ #ifdef PROTOTYPES_IN_STRUCT *v_map(double (*f)(double),VEC *,VEC *), /* out[i] <- f(x[i]) */ *_v_map(double (*f)(void *,double),void *,VEC *,VEC *), #else *v_map(double (*f)(),VEC *,VEC *), /* out[i] <- f(x[i]) */ *_v_map(double (*f)(),void *,VEC *,VEC *), #endif *v_lincomb(int,VEC **,Real *,VEC *), /* out <- sum_i s[i].x[i] */ *v_linlist(VEC *out,VEC *v1,double a1,...); /* out <- s1.x1 + s2.x2 + ... */ /* returns min_j x[j] (== x[i]) */ extern double v_min(VEC *, int *), /* returns max_j x[j] (== x[i]) */ v_max(VEC *, int *), /* returns sum_i x[i] */ v_sum(VEC *); /* Hadamard product: out[i] <- x[i].y[i] */ extern VEC *v_star(VEC *, VEC *, VEC *), /* out[i] <- x[i] / y[i] */ *v_slash(VEC *, VEC *, VEC *), /* sorts x, and sets order so that sorted x[i] = x[order[i]] */ *v_sort(VEC *, PERM *); /* returns inner product starting at component i0 */ extern double _in_prod(VEC *x,VEC *y,u_int i0), /* returns sum_{i=0}^{len-1} x[i].y[i] */ __ip__(Real *,Real *,int); /* see v_mltadd(), v_add(), v_sub() and v_zero() */ extern void __mltadd__(Real *,Real *,double,int), __add__(Real *,Real *,Real *,int), __sub__(Real *,Real *,Real *,int), __smlt__(Real *,double,Real *,int), __zero__(Real *,int); #endif /* MACRO */ /* usual way of computing the inner product */ #define in_prod(a,b) _in_prod(a,b,0) /* Norms */ /* scaled vector norms -- scale == NULL implies unscaled */ #ifndef ANSI_C extern double _v_norm1(), _v_norm2(), _v_norm_inf(), m_norm1(), m_norm_inf(), m_norm_frob(); #else /* returns sum_i |x[i]/scale[i]| */ extern double _v_norm1(VEC *x,VEC *scale), /* returns (scaled) Euclidean norm */ _v_norm2(VEC *x,VEC *scale), /* returns max_i |x[i]/scale[i]| */ _v_norm_inf(VEC *x,VEC *scale); /* unscaled matrix norms */ extern double m_norm1(MAT *A), m_norm_inf(MAT *A), m_norm_frob(MAT *A); #endif /* MACROS */ /* unscaled vector norms */ #define v_norm1(x) _v_norm1(x,VNULL) #define v_norm2(x) _v_norm2(x,VNULL) #define v_norm_inf(x) _v_norm_inf(x,VNULL) /* Basic matrix operations */ #ifndef ANSI_C extern MAT *sm_mlt(), *m_mlt(), *mmtr_mlt(), *mtrm_mlt(), *m_add(), *m_sub(), *sub_mat(), *m_transp(), *ms_mltadd(); extern BAND *bd_transp(); extern MAT *px_rows(), *px_cols(), *swap_rows(), *swap_cols(), *_set_row(), *_set_col(); extern VEC *get_row(), *get_col(), *sub_vec(), *mv_mltadd(), *vm_mltadd(); #else extern MAT *sm_mlt(double s,MAT *A,MAT *out), /* out <- s.A */ *m_mlt(MAT *A,MAT *B,MAT *out), /* out <- A.B */ *mmtr_mlt(MAT *A,MAT *B,MAT *out), /* out <- A.B^T */ *mtrm_mlt(MAT *A,MAT *B,MAT *out), /* out <- A^T.B */ *m_add(MAT *A,MAT *B,MAT *out), /* out <- A + B */ *m_sub(MAT *A,MAT *B,MAT *out), /* out <- A - B */ *sub_mat(MAT *A,u_int,u_int,u_int,u_int,MAT *out), *m_transp(MAT *A,MAT *out), /* out <- A^T */ /* out <- A + s.B */ *ms_mltadd(MAT *A,MAT *B,double s,MAT *out); extern BAND *bd_transp(BAND *in, BAND *out); /* out <- A^T */ extern MAT *px_rows(PERM *px,MAT *A,MAT *out), /* out <- P.A */ *px_cols(PERM *px,MAT *A,MAT *out), /* out <- A.P^T */ *swap_rows(MAT *,int,int,int,int), *swap_cols(MAT *,int,int,int,int), /* A[i][j] <- out[j], j >= j0 */ *_set_col(MAT *A,u_int i,VEC *out,u_int j0), /* A[i][j] <- out[i], i >= i0 */ *_set_row(MAT *A,u_int j,VEC *out,u_int i0); extern VEC *get_row(MAT *,u_int,VEC *), *get_col(MAT *,u_int,VEC *), *sub_vec(VEC *,int,int,VEC *), /* out <- x + s.A.y */ *mv_mltadd(VEC *x,VEC *y,MAT *A,double s,VEC *out), /* out^T <- x^T + s.y^T.A */ *vm_mltadd(VEC *x,VEC *y,MAT *A,double s,VEC *out); #endif /* MACROS */ /* row i of A <- vec */ #define set_row(mat,row,vec) _set_row(mat,row,vec,0) /* col j of A <- vec */ #define set_col(mat,col,vec) _set_col(mat,col,vec,0) /* Basic permutation operations */ #ifndef ANSI_C extern PERM *px_mlt(), *px_inv(), *px_transp(); extern int px_sign(); #else extern PERM *px_mlt(PERM *px1,PERM *px2,PERM *out), /* out <- px1.px2 */ *px_inv(PERM *px,PERM *out), /* out <- px^{-1} */ /* swap px[i] and px[j] */ *px_transp(PERM *px,u_int i,u_int j); /* returns sign(px) = +1 if px product of even # transpositions -1 if ps product of odd # transpositions */ extern int px_sign(PERM *); #endif /* Basic integer vector operations */ #ifndef ANSI_C extern IVEC *iv_add(), *iv_sub(), *iv_sort(); #else extern IVEC *iv_add(IVEC *ix,IVEC *iy,IVEC *out), /* out <- ix + iy */ *iv_sub(IVEC *ix,IVEC *iy,IVEC *out), /* out <- ix - iy */ /* sorts ix & sets order so that sorted ix[i] = old ix[order[i]] */ *iv_sort(IVEC *ix, PERM *order); #endif /* miscellaneous functions */ #ifndef ANSI_C extern double square(), cube(), mrand(); extern void smrand(), mrandlist(); extern void m_dump(), px_dump(), v_dump(), iv_dump(); extern MAT *band2mat(); extern BAND *mat2band(); #else double square(double x), /* returns x^2 */ cube(double x), /* returns x^3 */ mrand(void); /* returns random # in [0,1) */ void smrand(int seed), /* seeds mrand() */ mrandlist(Real *x, int len); /* generates len random numbers */ void m_dump(FILE *fp,MAT *a), px_dump(FILE *,PERM *px), v_dump(FILE *fp,VEC *x), iv_dump(FILE *fp, IVEC *ix); MAT *band2mat(BAND *bA, MAT *A); BAND *mat2band(MAT *A, int lb,int ub, BAND *bA); #endif /* miscellaneous constants */ #define VNULL ((VEC *)NULL) #define MNULL ((MAT *)NULL) #define PNULL ((PERM *)NULL) #define IVNULL ((IVEC *)NULL) #define BDNULL ((BAND *)NULL) /* varying number of arguments */ #ifdef ANSI_C #include /* prototypes */ int v_get_vars(int dim,...); int iv_get_vars(int dim,...); int m_get_vars(int m,int n,...); int px_get_vars(int dim,...); int v_resize_vars(int new_dim,...); int iv_resize_vars(int new_dim,...); int m_resize_vars(int m,int n,...); int px_resize_vars(int new_dim,...); int v_free_vars(VEC **,...); int iv_free_vars(IVEC **,...); int px_free_vars(PERM **,...); int m_free_vars(MAT **,...); #elif VARARGS /* old varargs is used */ #include /* prototypes */ int v_get_vars(); int iv_get_vars(); int m_get_vars(); int px_get_vars(); int v_resize_vars(); int iv_resize_vars(); int m_resize_vars(); int px_resize_vars(); int v_free_vars(); int iv_free_vars(); int px_free_vars(); int m_free_vars(); #endif #endif