diff --git a/.gitignore b/.gitignore index 9a59314..d441a80 100644 --- a/.gitignore +++ b/.gitignore @@ -81,6 +81,8 @@ target/ # svn .svn +# visual studio conf +.vscode # coverage .coverage diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..783ced5 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "files.associations": { + "meminfo.h": "c" + } +} \ No newline at end of file diff --git a/src/cpp/meschach/mesch12a/src/bdfactor.c b/src/cpp/meschach/mesch12a/src/bdfactor.c index 52279f3..5cd0974 100644 --- a/src/cpp/meschach/mesch12a/src/bdfactor.c +++ b/src/cpp/meschach/mesch12a/src/bdfactor.c @@ -49,8 +49,7 @@ static char rcsid[] = "$Id: "; in one column are now also in one column */ -BAND *bd_get(lb,ub,n) -int lb, ub, n; +BAND *bd_get(int lb, int ub, int n) { BAND *A; @@ -70,8 +69,7 @@ int lb, ub, n; return A; } -int bd_free(A) -BAND *A; +int bd_free(BAND *A) { if ( A == (BAND *)NULL || A->lb < 0 || A->ub < 0 ) /* don't trust it */ @@ -91,9 +89,7 @@ BAND *A; /* resize band matrix */ -BAND *bd_resize(A,new_lb,new_ub,new_n) -BAND *A; -int new_lb,new_ub,new_n; +BAND *bd_resize(BAND *A, int new_lb, int new_ub, int new_n) { int lb,ub,i,j,l,shift,umin; Real **Av; @@ -152,8 +148,7 @@ int new_lb,new_ub,new_n; -BAND *bd_copy(A,B) -BAND *A,*B; +BAND *bd_copy(BAND *A, BAND *B) { int lb,ub,i,j,n; @@ -183,9 +178,7 @@ BAND *A,*B; /* copy band matrix to a square matrix */ -MAT *band2mat(bA,A) -BAND *bA; -MAT *A; +MAT *band2mat(BAND *bA, MAT *A) { int i,j,l,n,n1; int lb, ub; @@ -214,10 +207,7 @@ MAT *A; /* copy a square matrix to a band matrix with lb subdiagonals and ub superdiagonals */ -BAND *mat2band(A,lb,ub,bA) -BAND *bA; -MAT *A; -int lb, ub; +BAND *mat2band(MAT *A, int lb, int ub, BAND *bA) { int i, j, l, n1; Real **bmat; @@ -249,8 +239,7 @@ int lb, ub; can be done in situ */ -BAND *bd_transp(in,out) -BAND *in, *out; +BAND *bd_transp(BAND *in, BAND *out) { int i, j, jj, l, k, lb, ub, lub, n, n1; int in_situ; @@ -362,9 +351,7 @@ BAND *in, *out; Matrix U is permuted, whereas L is not permuted !!! Therefore we save some memory. */ -BAND *bdLUfactor(bA,pivot) -BAND *bA; -PERM *pivot; +BAND *bdLUfactor(BAND *bA, PERM *pivot) { int i, j, k, l, n, n1, lb, ub, lub, k_end, k_lub; int i_max, shift; @@ -445,10 +432,7 @@ PERM *pivot; /* bdLUsolve -- given an LU factorisation in bA, solve bA*x=b */ /* pivot is changed upon return */ -VEC *bdLUsolve(bA,pivot,b,x) -BAND *bA; -PERM *pivot; -VEC *b,*x; +VEC *bdLUsolve(BAND *bA, PERM *pivot, VEC *b, VEC *x) { int i,j,l,n,n1,pi,lb,ub,jmin, maxj; Real c; @@ -503,8 +487,7 @@ VEC *b,*x; so it is possible to set A->ub = 0 */ -BAND *bdLDLfactor(A) -BAND *A; +BAND *bdLDLfactor(BAND *A) { int i,j,k,n,n1,lb,ki,jk,ji,lbkm,lbkp; Real **Av; @@ -550,9 +533,7 @@ BAND *A; /* solve A*x = b, where A is factorized by Choleski LDL^T factorization */ -VEC *bdLDLsolve(A,b,x) -BAND *A; -VEC *b, *x; +VEC *bdLDLsolve(BAND *A, VEC *b, VEC *x) { int i,j,l,n,n1,lb,ilb; Real **Av, *Avlb; @@ -609,9 +590,7 @@ VEC *b, *x; * may not work in situ (x != out) */ -VEC *bd_mv_mlt(A, x, out) -BAND *A; -VEC *x, *out; +VEC *bd_mv_mlt(BAND *A, VEC *x, VEC *out) { int i, j, j_end, k; int start_idx, end_idx; diff --git a/src/cpp/meschach/mesch12a/src/bkpfacto.c b/src/cpp/meschach/mesch12a/src/bkpfacto.c index a06062a..df6b1db 100644 --- a/src/cpp/meschach/mesch12a/src/bkpfacto.c +++ b/src/cpp/meschach/mesch12a/src/bkpfacto.c @@ -42,14 +42,13 @@ static char rcsid[] = "$Id: bkpfacto.c,v 1.7 1994/01/13 05:45:50 des Exp $"; #define alpha 0.6403882032022076 /* = (1+sqrt(17))/8 */ /* sqr -- returns square of x -- utility function */ -double sqr(x) -double x; +double sqr(double x) { return x*x; } /* interchange -- a row/column swap routine */ -static void interchange(A,i,j) -MAT *A; /* assumed != NULL & also SQUARE */ -int i, j; /* assumed in range */ +static void interchange(MAT *A, int i, int j) +//MAT *A; /* assumed != NULL & also SQUARE */ +//int i, j; /* assumed in range */ { Real **A_me, tmp; int k, n; @@ -99,9 +98,7 @@ int i, j; /* assumed in range */ P is a permutation matrix, M lower triangular and D is block diagonal with blocks of size 1 or 2 -- P is stored in pivot; blocks[i]==i iff D[i][i] is a block */ -MAT *BKPfactor(A,pivot,blocks) -MAT *A; -PERM *pivot, *blocks; +MAT *BKPfactor(MAT *A,PERM *pivot, PERM *blocks) { int i, j, k, n, onebyone, r; Real **A_me, aii, aip1, aip1i, lambda, sigma, tmp; @@ -224,10 +221,7 @@ PERM *pivot, *blocks; /* BKPsolve -- solves A.x = b where A has been factored a la BKPfactor() -- returns x, which is created if NULL */ -VEC *BKPsolve(A,pivot,block,b,x) -MAT *A; -PERM *pivot, *block; -VEC *b, *x; +VEC *BKPsolve(MAT *A, PERM *pivot, PERM *block, VEC *b, VEC *x) { static VEC *tmp=VNULL; /* dummy storage needed */ int i, j, n, onebyone; diff --git a/src/cpp/meschach/mesch12a/src/fmacheps.c b/src/cpp/meschach/mesch12a/src/fmacheps.c index df1bb08..9ed6c2a 100644 --- a/src/cpp/meschach/mesch12a/src/fmacheps.c +++ b/src/cpp/meschach/mesch12a/src/fmacheps.c @@ -26,8 +26,7 @@ #include -double fclean(x) -double x; +double fclean(double x) { static float y; y = x; diff --git a/src/cpp/meschach/mesch12a/src/init.c b/src/cpp/meschach/mesch12a/src/init.c index 2d96b16..b2e543c 100644 --- a/src/cpp/meschach/mesch12a/src/init.c +++ b/src/cpp/meschach/mesch12a/src/init.c @@ -36,8 +36,7 @@ static char rcsid[] = "$Id: init.c,v 1.6 1994/01/13 05:36:58 des Exp $"; #include "matrix.h" /* v_zero -- zero the vector x */ -VEC *v_zero(x) -VEC *x; +VEC *v_zero(VEC *x) { if ( x == VNULL ) error(E_NULL,"v_zero"); @@ -51,8 +50,7 @@ VEC *x; /* iv_zero -- zero the vector ix */ -IVEC *iv_zero(ix) -IVEC *ix; +IVEC *iv_zero(IVEC *ix) { int i; @@ -67,8 +65,7 @@ IVEC *ix; /* m_zero -- zero the matrix A */ -MAT *m_zero(A) -MAT *A; +MAT *m_zero(MAT *A) { int i, A_m, A_n; Real **A_me; @@ -87,8 +84,7 @@ MAT *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; +MAT *m_ident(MAT *A) { int i, size; @@ -104,8 +100,7 @@ MAT *A; } /* px_ident -- set px to identity permutation */ -PERM *px_ident(px) -PERM *px; +PERM *px_ident(PERM *px) { int i, px_size; u_int *px_pe; @@ -167,9 +162,7 @@ double mrand() } /* mrandlist -- fills the array a[] with len random numbers */ -void mrandlist(a, len) -Real a[]; -int len; +void mrandlist(Real a[], int len) { int i; long lval; @@ -194,8 +187,7 @@ int len; /* smrand -- set seed for mrand() */ -void smrand(seed) -int seed; +void smrand(int seed) { int i; @@ -216,8 +208,7 @@ int seed; /* v_rand -- initialises x to be a random vector, components independently & uniformly ditributed between 0 and 1 */ -VEC *v_rand(x) -VEC *x; +VEC *v_rand(VEC *x) { /* int i; */ @@ -234,8 +225,7 @@ VEC *x; /* m_rand -- initialises A to be a random vector, components independently & uniformly distributed between 0 and 1 */ -MAT *m_rand(A) -MAT *A; +MAT *m_rand(MAT *A) { int i /* , j */; @@ -252,8 +242,7 @@ MAT *A; } /* v_ones -- fills x with one's */ -VEC *v_ones(x) -VEC *x; +VEC *v_ones(VEC *x) { int i; @@ -267,8 +256,7 @@ VEC *x; } /* m_ones -- fills matrix with one's */ -MAT *m_ones(A) -MAT *A; +MAT *m_ones(MAT *A) { int i, j; @@ -283,8 +271,7 @@ MAT *A; } /* v_count -- initialises x so that x->ve[i] == i */ -VEC *v_count(x) -VEC *x; +VEC *v_count(VEC *x) { int i; diff --git a/src/cpp/meschach/mesch12a/src/itersym.c b/src/cpp/meschach/mesch12a/src/itersym.c index e40463f..6ed4e06 100644 --- a/src/cpp/meschach/mesch12a/src/itersym.c +++ b/src/cpp/meschach/mesch12a/src/itersym.c @@ -64,11 +64,7 @@ VEC *trieig(); x = iter_spcg(A,LLT,b,eps,VNULL,limit,steps); In the second case the solution vector is created. */ -VEC *iter_spcg(A,LLT,b,eps,x,limit,steps) -SPMAT *A, *LLT; -VEC *b, *x; -double eps; -int *steps, limit; +VEC *iter_spcg(SPMAT *A, SPMAT *LLT, VEC *b, double eps, VEC *x, int limit, int *steps) { ITER *ip; @@ -93,8 +89,7 @@ int *steps, limit; /* Conjugate gradients method; */ -VEC *iter_cg(ip) -ITER *ip; +VEC *iter_cg(ITER *ip) { static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL; Real alpha, beta, inner, old_inner, nres; @@ -181,11 +176,7 @@ ITER *ip; -- creates T matrix of size == m, but no larger than before beta_k == 0 -- uses passed routine to do matrix-vector multiplies */ -void iter_lanczos(ip,a,b,beta2,Q) -ITER *ip; -VEC *a, *b; -Real *beta2; -MAT *Q; +void iter_lanczos(ITER *ip, VEC *a, VEC *b, Real *beta2, MAT *Q) { int j; static VEC *v = VNULL, *w = VNULL, *tmp = VNULL; @@ -253,12 +244,7 @@ MAT *Q; } /* iter_splanczos -- version that uses sparse matrix data structure */ -void iter_splanczos(A,m,x0,a,b,beta2,Q) -SPMAT *A; -int m; -VEC *x0, *a, *b; -Real *beta2; -MAT *Q; +void iter_splanczos(SPMAT *A, int m, VEC *x0, VEC *a, VEC *b, Real *beta2, MAT *Q) { ITER *ip; @@ -278,10 +264,7 @@ extern double frexp(), ldexp(); /* product -- returns the product of a long list of numbers -- answer stored in mant (mantissa) and expt (exponent) */ -static double product(a,offset,expt) -VEC *a; -double offset; -int *expt; +static double product(VEC *a, double offset, int *expt) { Real mant, tmp_fctr; int i, tmp_expt; @@ -325,10 +308,8 @@ int *expt; /* product2 -- returns the product of a long list of numbers -- answer stored in mant (mantissa) and expt (exponent) */ -static double product2(a,k,expt) -VEC *a; -int k; /* entry of a to leave out */ -int *expt; +static double product2(VEC *a, int k, int *expt) +//int k; /* entry of a to leave out */ { Real mant, mu, tmp_fctr; int i, tmp_expt; @@ -362,8 +343,7 @@ int *expt; } /* dbl_cmp -- comparison function to pass to qsort() */ -static int dbl_cmp(x,y) -Real *x, *y; +static int dbl_cmp(Real *x, Real *y) { Real tmp; @@ -374,11 +354,12 @@ Real *x, *y; /* iter_lanczos2 -- lanczos + error estimate for every e-val -- uses Cullum & Willoughby approach, Sparse Matrix Proc. 1978 -- returns multiple e-vals where multiple e-vals may not exist - -- returns evals vector */ -VEC *iter_lanczos2(ip,evals,err_est) -ITER *ip; /* ITER structure */ -VEC *evals; /* eigenvalue vector */ -VEC *err_est; /* error estimates of eigenvalues */ + -- returns evals vector + ITER *ip; ITER structure + VEC *evals; eigenvalue vector + VEC *err_est; error estimates of eigenvalues + */ +VEC *iter_lanczos2(ITER *ip, VEC *evals, VEC *err_est) { VEC *a; static VEC *b=VNULL, *a2=VNULL, *b2=VNULL; @@ -467,14 +448,11 @@ VEC *err_est; /* error estimates of eigenvalues */ } /* iter_splanczos2 -- version of iter_lanczos2() that uses sparse matrix data - structure */ - -VEC *iter_splanczos2(A,m,x0,evals,err_est) -SPMAT *A; -int m; -VEC *x0; /* initial vector */ -VEC *evals; /* eigenvalue vector */ -VEC *err_est; /* error estimates of eigenvalues */ + structure + VEC *x0; initial vector + VEC *evals; eigenvalue vector + VEC *err_est; error estimates of eigenvalues */ +VEC *iter_splanczos2(SPMAT *A, int m, VEC *x0, VEC *evals, VEC *err_est) { ITER *ip; VEC *a; @@ -498,8 +476,7 @@ VEC *err_est; /* error estimates of eigenvalues */ Another variant - mainly for testing */ -VEC *iter_cg1(ip) -ITER *ip; +VEC *iter_cg1(ITER *ip) { static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL; Real alpha; diff --git a/src/cpp/meschach/mesch12a/src/lufactor.c b/src/cpp/meschach/mesch12a/src/lufactor.c index b732487..6554430 100644 --- a/src/cpp/meschach/mesch12a/src/lufactor.c +++ b/src/cpp/meschach/mesch12a/src/lufactor.c @@ -42,9 +42,7 @@ static char rcsid[] = "$Id: lufactor.c,v 1.9 1995/04/20 19:21:54 des Exp $"; /* LUfactor -- gaussian elimination with scaled partial pivoting -- Note: returns LU matrix which is A */ -MAT *LUfactor(A,pivot) -MAT *A; -PERM *pivot; +MAT *LUfactor(MAT *A, PERM *pivot) { u_int i, j, k, k_max, m, n; int i_max; @@ -136,10 +134,7 @@ PERM *pivot; /* LUsolve -- given an LU factorisation in A, solve Ax=b */ -VEC *LUsolve(A,pivot,b,x) -MAT *A; -PERM *pivot; -VEC *b,*x; +VEC *LUsolve(MAT *A, PERM *pivot, VEC *b, VEC *x) { if ( A==(MAT *)NULL || b==(VEC *)NULL || pivot==(PERM *)NULL ) error(E_NULL,"LUsolve"); @@ -155,10 +150,7 @@ VEC *b,*x; } /* LUTsolve -- given an LU factorisation in A, solve A^T.x=b */ -VEC *LUTsolve(LU,pivot,b,x) -MAT *LU; -PERM *pivot; -VEC *b,*x; +VEC *LUTsolve(MAT *LU, PERM *pivot, VEC *b, VEC *x) { if ( ! LU || ! b || ! pivot ) error(E_NULL,"LUTsolve"); @@ -175,8 +167,7 @@ VEC *b,*x; /* m_inverse -- returns inverse of A, provided A is not too rank deficient -- uses LU factorisation */ -MAT *m_inverse(A,out) -MAT *A, *out; +MAT *m_inverse(MAT *A, MAT *out) { int i; static VEC *tmp = VNULL, *tmp2 = VNULL; @@ -213,9 +204,7 @@ MAT *A, *out; /* LUcondest -- returns an estimate of the condition number of LU given the LU factorisation in compact form */ -double LUcondest(LU,pivot) -MAT *LU; -PERM *pivot; +double LUcondest(MAT *LU, PERM *pivot) { static VEC *y = VNULL, *z = VNULL; Real cond_est, L_norm, U_norm, sum, tiny; diff --git a/src/cpp/meschach/mesch12a/src/memory.c b/src/cpp/meschach/mesch12a/src/memory.c index 3fc1594..f613906 100644 --- a/src/cpp/meschach/mesch12a/src/memory.c +++ b/src/cpp/meschach/mesch12a/src/memory.c @@ -32,8 +32,7 @@ static char rcsid[] = "$Id: memory.c,v 1.13 1994/04/05 02:10:37 des Exp $"; /* m_get -- gets an mxn matrix (in MAT form) by dynamic memory allocation */ -MAT *m_get(m,n) -int m,n; +MAT *m_get(int m, int n) { MAT *matrix; int i; @@ -90,8 +89,7 @@ int m,n; /* px_get -- gets a PERM of given 'size' by dynamic memory allocation -- Note: initialized to the identity permutation */ -PERM *px_get(size) -int size; +PERM *px_get(int size) { PERM *permute; int i; @@ -121,8 +119,7 @@ int size; /* v_get -- gets a VEC of dimension 'dim' -- Note: initialized to zero */ -VEC *v_get(size) -int size; +VEC *v_get(int size) { VEC *vector; @@ -150,8 +147,7 @@ int size; } /* m_free -- returns MAT & asoociated memory back to memory heap */ -int m_free(mat) -MAT *mat; +int m_free(MAT *mat) { #ifdef SEGMENTED int i; @@ -197,8 +193,7 @@ MAT *mat; /* px_free -- returns PERM & asoociated memory back to memory heap */ -int px_free(px) -PERM *px; +int px_free(PERM *px) { if ( px==(PERM *)NULL || (int)(px->size) < 0 ) /* don't trust it */ @@ -227,8 +222,7 @@ PERM *px; /* v_free -- returns VEC & asoociated memory back to memory heap */ -int v_free(vec) -VEC *vec; +int v_free(VEC *vec) { if ( vec==(VEC *)NULL || (int)(vec->dim) < 0 ) /* don't trust it */ @@ -258,9 +252,7 @@ VEC *vec; /* m_resize -- returns the matrix A of size new_m x new_n; A is zeroed -- if A == NULL on entry then the effect is equivalent to m_get() */ -MAT *m_resize(A,new_m,new_n) -MAT *A; -int new_m, new_n; +MAT *m_resize(MAT *A, int new_m, int new_n) { int i; int new_max_m, new_max_n, new_size, old_m, old_n; @@ -395,9 +387,7 @@ int new_m, new_n; /* px_resize -- returns the permutation px with size new_size -- px is set to the identity permutation */ -PERM *px_resize(px,new_size) -PERM *px; -int new_size; +PERM *px_resize(PERM *px, int new_size) { int i; @@ -437,9 +427,7 @@ int new_size; /* v_resize -- returns the vector x with dim new_dim -- x is set to the zero vector */ -VEC *v_resize(x,new_dim) -VEC *x; -int new_dim; +VEC *v_resize(VEC *x, int new_dim) { if (new_dim < 0) diff --git a/src/cpp/meschach/mesch12a/src/memstat.c b/src/cpp/meschach/mesch12a/src/memstat.c index a5a975e..775d3f8 100644 --- a/src/cpp/meschach/mesch12a/src/memstat.c +++ b/src/cpp/meschach/mesch12a/src/memstat.c @@ -87,8 +87,7 @@ static unsigned int mem_hash(void **ptr) /* look for a place in mem_stat_var */ -static int mem_lookup(var) -void **var; +static int mem_lookup(void **var) { int k, j; @@ -148,9 +147,7 @@ void **var; returned value >= 0 --> registered with this mark; */ -int mem_stat_reg_list(var,type,list) -void **var; -int type,list; +int mem_stat_reg_list(void **var, int type, int list) { int n; @@ -188,8 +185,7 @@ int type,list; -1 if mark is negative. */ -int mem_stat_mark(mark) -int mark; +int mem_stat_mark(int mark) { if (mark < 0) { mem_stat_mark_curr = 0; @@ -217,8 +213,7 @@ int mark; 0 if mark == 0; */ -int mem_stat_free_list(mark,list) -int mark,list; +int mem_stat_free_list(int mark, int list) { u_int i,j; int (*free_fn)(); @@ -273,9 +268,7 @@ int mark,list; /* only for diagnostic purposes */ -void mem_stat_dump(fp,list) -FILE *fp; -int list; +void mem_stat_dump(FILE *fp, int list) { u_int i,j,k=1; diff --git a/src/cpp/meschach/mesch12a/src/norm.c b/src/cpp/meschach/mesch12a/src/norm.c index aaf2db7..045541f 100644 --- a/src/cpp/meschach/mesch12a/src/norm.c +++ b/src/cpp/meschach/mesch12a/src/norm.c @@ -35,8 +35,7 @@ static char rcsid[] = "$Id: norm.c,v 1.6 1994/01/13 05:34:35 des Exp $"; /* _v_norm1 -- computes (scaled) 1-norms of vectors */ -double _v_norm1(x,scale) -VEC *x, *scale; +double _v_norm1(VEC *x, VEC *scale) { int i, dim; Real s, sum; @@ -61,18 +60,15 @@ VEC *x, *scale; } /* square -- returns x^2 */ -double square(x) -double x; +double square(double x) { return x*x; } /* cube -- returns x^3 */ -double cube(x) -double x; +double cube(double x) { return x*x*x; } /* _v_norm2 -- computes (scaled) 2-norm (Euclidean norm) of vectors */ -double _v_norm2(x,scale) -VEC *x, *scale; +double _v_norm2(VEC *x, VEC *scale) { int i, dim; Real s, sum; @@ -100,8 +96,7 @@ VEC *x, *scale; #define max(a,b) ((a) > (b) ? (a) : (b)) /* _v_norm_inf -- computes (scaled) infinity-norm (supremum norm) of vectors */ -double _v_norm_inf(x,scale) -VEC *x, *scale; +double _v_norm_inf(VEC *x, VEC *scale) { int i, dim; Real s, maxval, tmp; @@ -129,8 +124,7 @@ VEC *x, *scale; } /* m_norm1 -- compute matrix 1-norm -- unscaled */ -double m_norm1(A) -MAT *A; +double m_norm1(MAT *A) { int i, j, m, n; Real maxval, sum; @@ -153,8 +147,7 @@ MAT *A; } /* m_norm_inf -- compute matrix infinity-norm -- unscaled */ -double m_norm_inf(A) -MAT *A; +double m_norm_inf(MAT *A) { int i, j, m, n; Real maxval, sum; @@ -177,8 +170,7 @@ MAT *A; } /* m_norm_frob -- compute matrix frobenius-norm -- unscaled */ -double m_norm_frob(A) -MAT *A; +double m_norm_frob(MAT *A) { int i, j, m, n; Real sum; diff --git a/src/cpp/meschach/mesch12a/src/qrfactor.c b/src/cpp/meschach/mesch12a/src/qrfactor.c index 716acf9..2c1dbff 100644 --- a/src/cpp/meschach/mesch12a/src/qrfactor.c +++ b/src/cpp/meschach/mesch12a/src/qrfactor.c @@ -61,9 +61,8 @@ extern VEC *Usolve(); /* See matrix2.h */ /* QRfactor -- forms the QR factorisation of A -- factorisation stored in compact form as described above ( not quite standard format ) */ /* MAT *QRfactor(A,diag,beta) */ -MAT *QRfactor(A,diag) -MAT *A; -VEC *diag /* ,*beta */; +MAT *QRfactor(MAT *A, VEC *diag) +//VEC *diag /* ,*beta */; { u_int k,limit; Real beta; @@ -98,10 +97,8 @@ VEC *diag /* ,*beta */; -- factorisation stored in compact form as described above ( not quite standard format ) */ /* MAT *QRCPfactor(A,diag,beta,px) */ -MAT *QRCPfactor(A,diag,px) -MAT *A; -VEC *diag /* , *beta */; -PERM *px; +MAT *QRCPfactor(MAT *A, VEC *diag, PERM *px) +//VEC *diag /* , *beta */; { u_int i, i_max, j, k, limit; static VEC *gamma=VNULL, *tmp1=VNULL, *tmp2=VNULL; @@ -181,9 +178,7 @@ PERM *px; /* Qsolve -- solves Qx = b, Q is an orthogonal matrix stored in compact form a la QRfactor() -- may be in-situ */ /* VEC *_Qsolve(QR,diag,beta,b,x,tmp) */ -VEC *_Qsolve(QR,diag,b,x,tmp) -MAT *QR; -VEC *diag /* ,*beta */ , *b, *x, *tmp; +VEC *_Qsolve(MAT *QR, VEC *diag, VEC *b, VEC *x, VEC *tmp) { u_int dynamic; int k, limit; @@ -222,9 +217,7 @@ VEC *diag /* ,*beta */ , *b, *x, *tmp; /* makeQ -- constructs orthogonal matrix from Householder vectors stored in compact QR form */ /* MAT *makeQ(QR,diag,beta,Qout) */ -MAT *makeQ(QR,diag,Qout) -MAT *QR,*Qout; -VEC *diag /* , *beta */; +MAT *makeQ(MAT *QR, VEC *diag, MAT *Qout) { static VEC *tmp1=VNULL,*tmp2=VNULL; u_int i, limit; @@ -272,8 +265,7 @@ VEC *diag /* , *beta */; /* makeR -- constructs upper triangular matrix from QR (compact form) -- may be in-situ (all it does is zero the lower 1/2) */ -MAT *makeR(QR,Rout) -MAT *QR,*Rout; +MAT *makeR(MAT *QR, MAT *Rout) { u_int i,j; @@ -291,9 +283,7 @@ MAT *QR,*Rout; /* QRsolve -- solves the system Q.R.x=b where Q & R are stored in compact form -- returns x, which is created if necessary */ /* VEC *QRsolve(QR,diag,beta,b,x) */ -VEC *QRsolve(QR,diag,b,x) -MAT *QR; -VEC *diag /* , *beta */ , *b, *x; +VEC *QRsolve(MAT *QR, VEC *diag, VEC *b, VEC *x) { int limit; static VEC *tmp = VNULL; @@ -317,11 +307,7 @@ VEC *diag /* , *beta */ , *b, *x; /* QRCPsolve -- solves A.x = b where A is factored by QRCPfactor() -- assumes that A is in the compact factored form */ /* VEC *QRCPsolve(QR,diag,beta,pivot,b,x) */ -VEC *QRCPsolve(QR,diag,pivot,b,x) -MAT *QR; -VEC *diag /* , *beta */; -PERM *pivot; -VEC *b, *x; +VEC *QRCPsolve(MAT *QR, VEC *diag, PERM *pivot, VEC *b, VEC *x) { static VEC *tmp=VNULL; @@ -339,9 +325,7 @@ VEC *b, *x; /* Umlt -- compute out = upper_triang(U).x -- may be in situ */ -static VEC *Umlt(U,x,out) -MAT *U; -VEC *x, *out; +static VEC *Umlt(MAT *U, VEC *x, VEC *out) { int i, limit; @@ -359,9 +343,7 @@ VEC *x, *out; } /* UTmlt -- returns out = upper_triang(U)^T.x */ -static VEC *UTmlt(U,x,out) -MAT *U; -VEC *x, *out; +static VEC *UTmlt(MAT *U, VEC *x, VEC *out) { Real sum; int i, j, limit; @@ -386,9 +368,7 @@ VEC *x, *out; compact form -- returns sc -- original due to Mike Osborne modified Wed 09th Dec 1992 */ -VEC *QRTsolve(A,diag,c,sc) -MAT *A; -VEC *diag, *c, *sc; +VEC *QRTsolve(MAT *A, VEC *diag, VEC *c, VEC *sc) { int i, j, k, n, p; Real beta, r_ii, s, tmp_val; @@ -445,8 +425,7 @@ VEC *diag, *c, *sc; -- if the matrix is exactly singular, HUGE is returned -- note that QRcondest() is likely to be more reliable for matrices factored using QRCPfactor() */ -double QRcondest(QR) -MAT *QR; +double QRcondest(MAT *QR) { static VEC *y=VNULL; Real norm1, norm2, sum, tmp1, tmp2; diff --git a/src/cpp/meschach/mesch12a/src/sprow.c b/src/cpp/meschach/mesch12a/src/sprow.c index 0a0a0d8..b5aa0d3 100644 --- a/src/cpp/meschach/mesch12a/src/sprow.c +++ b/src/cpp/meschach/mesch12a/src/sprow.c @@ -41,9 +41,7 @@ static char rcsid[] = "$Id: sprow.c,v 1.1 1994/01/13 05:35:36 des Exp $"; /* sprow_dump - prints relevant information about the sparse row r */ -void sprow_dump(fp,r) -FILE *fp; -SPROW *r; +void sprow_dump(FILE *fp, SPROW *r) { int j_idx; row_elt *elts; @@ -71,9 +69,7 @@ SPROW *r; /* sprow_idx -- get index into row for a given column in a given row -- return -1 on error -- return -(idx+2) where idx is index to insertion point */ -int sprow_idx(r,col) -SPROW *r; -int col; +int sprow_idx(SPROW *r, int col) { int lo, hi, mid; int tmp; @@ -118,8 +114,7 @@ int col; /* sprow_get -- gets, initialises and returns a SPROW structure -- max. length is maxlen */ -SPROW *sprow_get(maxlen) -int maxlen; +SPROW *sprow_get(int maxlen) { SPROW *r; @@ -151,9 +146,7 @@ int maxlen; -- type must be TYPE_SPMAT if r is a row of a SPMAT structure, otherwise it must be TYPE_SPROW -- returns r */ -SPROW *sprow_xpd(r,n,type) -SPROW *r; -int n,type; +SPROW *sprow_xpd(SPROW *r, int n, int type) { int newlen; @@ -211,9 +204,7 @@ int n,type; /* sprow_resize -- resize a SPROW variable by means of realloc() -- n is a new size -- returns r */ -SPROW *sprow_resize(r,n,type) -SPROW *r; -int n,type; +SPROW *sprow_resize(SPROW *r, int n, int type) { if (n < 0) error(E_NEG,"sprow_resize"); @@ -255,8 +246,7 @@ int n,type; /* release a row of a matrix */ -int sprow_free(r) -SPROW *r; +int sprow_free(SPROW *r) { if ( ! r ) return -1; @@ -284,9 +274,7 @@ SPROW *r; whether r_out is a row of a SPMAT structure or a SPROW variable -- returns r_out */ -SPROW *sprow_merge(r1,r2,r_out,type) -SPROW *r1, *r2, *r_out; -int type; +SPROW *sprow_merge(SPROW *r1, SPROW *r2, SPROW *r_out, int type) { int idx1, idx2, idx_out, len1, len2, len_out; row_elt *elt1, *elt2, *elt_out; @@ -339,9 +327,7 @@ int type; whether r_out is a row of a SPMAT structure or a SPROW variable -- returns r_out */ -SPROW *sprow_copy(r1,r2,r_out,type) -SPROW *r1, *r2, *r_out; -int type; +SPROW *sprow_copy(SPROW *r1, SPROW *r2, SPROW *r_out, int type) { int idx1, idx2, idx_out, len1, len2, len_out; row_elt *elt1, *elt2, *elt_out; @@ -395,10 +381,7 @@ int type; whether r_out is a row of a SPMAT structure or a SPROW variable -- returns r_out */ -SPROW *sprow_mltadd(r1,r2,alpha,j0,r_out,type) -SPROW *r1, *r2, *r_out; -double alpha; -int j0, type; +SPROW *sprow_mltadd(SPROW *r1, SPROW *r2, double alpha, int j0, SPROW * r_out, int type) { int idx1, idx2, idx_out, len1, len2, len_out; row_elt *elt1, *elt2, *elt_out; @@ -465,9 +448,7 @@ int j0, type; whether r_out is a row of a SPMAT structure or a SPROW variable -- returns r_out */ -SPROW *sprow_add(r1,r2,j0,r_out,type) -SPROW *r1, *r2, *r_out; -int j0, type; +SPROW *sprow_add(SPROW *r1, SPROW *r2, int j0, SPROW *r_out, int type) { int idx1, idx2, idx_out, len1, len2, len_out; row_elt *elt1, *elt2, *elt_out; @@ -534,9 +515,7 @@ int j0, type; whether r_out is a row of a SPMAT structure or a SPROW variable -- returns r_out */ -SPROW *sprow_sub(r1,r2,j0,r_out,type) -SPROW *r1, *r2, *r_out; -int j0, type; +SPROW *sprow_sub(SPROW *r1, SPROW *r2, int j0, SPROW *r_out, int type) { int idx1, idx2, idx_out, len1, len2, len_out; row_elt *elt1, *elt2, *elt_out; @@ -601,10 +580,7 @@ int j0, type; -- can be in situ -- only for columns j0, j0+1, ... -- returns r_out */ -SPROW *sprow_smlt(r1,alpha,j0,r_out,type) -SPROW *r1, *r_out; -double alpha; -int j0, type; +SPROW *sprow_smlt(SPROW *r1, double alpha, int j0, SPROW *r_out,int type) { int idx1, idx_out, len1; row_elt *elt1, *elt_out; @@ -640,9 +616,7 @@ int j0, type; /* sprow_foutput -- print a representation of r on stream fp */ -void sprow_foutput(fp,r) -FILE *fp; -SPROW *r; +void sprow_foutput(FILE *fp, SPROW *r) { int i, len; row_elt *e; @@ -662,10 +636,7 @@ SPROW *r; /* sprow_set_val -- sets the j-th column entry of the sparse row r -- Note: destroys the usual column & row access paths */ -double sprow_set_val(r,j,val) -SPROW *r; -int j; -double val; +double sprow_set_val(SPROW *r, int j, double val) { int idx, idx2, new_len; diff --git a/src/cpp/meschach/mesch12a/src/spswap.c b/src/cpp/meschach/mesch12a/src/spswap.c index 1b11655..99fd5d6 100644 --- a/src/cpp/meschach/mesch12a/src/spswap.c +++ b/src/cpp/meschach/mesch12a/src/spswap.c @@ -41,10 +41,7 @@ static char rcsid[] = "$Id: spswap.c,v 1.3 1994/01/13 05:44:43 des Exp $"; /* scan_to -- updates scan (int) vectors to point to the last row in each column with row # <= max_row, if any */ -void scan_to(A, scan_row, scan_idx, col_list, max_row) -SPMAT *A; -IVEC *scan_row, *scan_idx, *col_list; -int max_row; +void scan_to(SPMAT *A, IVEC *scan_row, IVEC *scan_idx, IVEC *col_list, int max_row) { int col, idx, j_idx, row_num; SPROW *r; @@ -103,9 +100,7 @@ int max_row; } /* patch_col -- patches column access paths for fill-in */ -void patch_col(A, col, old_row, old_idx, row_num, idx) -SPMAT *A; -int col, old_row, old_idx, row_num, idx; +void patch_col(SPMAT *A, int col, int old_row, int old_idx, int row_num, int idx) { SPROW *r; row_elt *e; @@ -130,9 +125,7 @@ int col, old_row, old_idx, row_num, idx; -- row_num is returned; idx is also set by this routine -- assumes that the column access paths (possibly without the nxt_idx fields) are set up */ -row_elt *chase_col(A, col, row_num, idx, max_row) -SPMAT *A; -int col, *row_num, *idx, max_row; +row_elt *chase_col(SPMAT *A, int col, int *row_num, int *idx, int max_row) { int old_idx, old_row, tmp_idx, tmp_row; SPROW *r; @@ -206,9 +199,7 @@ int col, *row_num, *idx, max_row; /* chase_past -- as for chase_col except that we want the first row whose row # >= min_row; -1 indicates no such row */ -row_elt *chase_past(A, col, row_num, idx, min_row) -SPMAT *A; -int col, *row_num, *idx, min_row; +row_elt *chase_past(SPMAT *A, int col, int *row_num, int *idx, int min_row) { SPROW *r; row_elt *e; @@ -255,9 +246,7 @@ int col, *row_num, *idx, min_row; /* bump_col -- move along to next nonzero entry in column col after row_num -- update row_num and idx */ -row_elt *bump_col(A, col, row_num, idx) -SPMAT *A; -int col, *row_num, *idx; +row_elt *bump_col(SPMAT *A, int col, int *row_num, int *idx) { SPROW *r; row_elt *e; diff --git a/src/cpp/meschach/mesch12a/src/submat.c b/src/cpp/meschach/mesch12a/src/submat.c index 56c58bb..ffbcdb4 100644 --- a/src/cpp/meschach/mesch12a/src/submat.c +++ b/src/cpp/meschach/mesch12a/src/submat.c @@ -33,10 +33,7 @@ static char rcsid[] = "$Id: submat.c,v 1.2 1994/01/13 05:28:12 des Exp $"; /* get_col -- gets a specified column of a matrix and retruns it as a vector */ -VEC *get_col(mat,col,vec) -u_int col; -MAT *mat; -VEC *vec; +VEC *get_col(MAT *mat, u_int col, VEC *vec) { u_int i; @@ -54,10 +51,7 @@ VEC *vec; } /* get_row -- gets a specified row of a matrix and retruns it as a vector */ -VEC *get_row(mat,row,vec) -u_int row; -MAT *mat; -VEC *vec; +VEC *get_row(MAT *mat, u_int row, VEC *vec) { u_int i; @@ -75,10 +69,7 @@ VEC *vec; } /* _set_col -- sets column of matrix to values given in vec (in situ) */ -MAT *_set_col(mat,col,vec,i0) -MAT *mat; -VEC *vec; -u_int col,i0; +MAT *_set_col(MAT *mat, u_int col, VEC *vec, u_int i0) { u_int i,lim; @@ -94,10 +85,7 @@ u_int col,i0; } /* _set_row -- sets row of matrix to values given in vec (in situ) */ -MAT *_set_row(mat,row,vec,j0) -MAT *mat; -VEC *vec; -u_int row,j0; +MAT *_set_row(MAT *mat, u_int row, VEC *vec, u_int j0) { u_int j,lim; @@ -116,9 +104,7 @@ u_int row,j0; from (row1,col1) to (row2,col2) -- Note: storage is shared so that altering the "new" matrix will alter the "old" matrix */ -MAT *sub_mat(old,row1,col1,row2,col2,new) -MAT *old,*new; -u_int row1,col1,row2,col2; +MAT *sub_mat(MAT *old, u_int row1, u_int col1, u_int row2, u_int col2, MAT *new) { u_int i; @@ -153,9 +139,7 @@ u_int row1,col1,row2,col2; /* sub_vec -- returns sub-vector which is formed by the elements i1 to i2 -- as for sub_mat, storage is shared */ -VEC *sub_vec(old,i1,i2,new) -VEC *old, *new; -int i1, i2; +VEC *sub_vec(VEC *old, int i1, int i2, VEC *new) { if ( old == (VEC *)NULL ) error(E_NULL,"sub_vec"); diff --git a/src/cpp/meschach/mesch12a/src/vecop.c b/src/cpp/meschach/mesch12a/src/vecop.c index 69eaa4e..1e1f6d3 100644 --- a/src/cpp/meschach/mesch12a/src/vecop.c +++ b/src/cpp/meschach/mesch12a/src/vecop.c @@ -33,9 +33,7 @@ static char rcsid[] = "$Id: vecop.c,v 1.4 1994/03/08 05:50:39 des Exp $"; /* _in_prod -- inner product of two vectors from i0 downwards */ -double _in_prod(a,b,i0) -VEC *a,*b; -u_int i0; +double _in_prod(VEC *a, VEC *b, u_int i0) { u_int limit; /* Real *a_v, *b_v; */ @@ -59,9 +57,7 @@ u_int i0; } /* sv_mlt -- scalar-vector multiply -- may be in-situ */ -VEC *sv_mlt(scalar,vector,out) -double scalar; -VEC *vector,*out; +VEC *sv_mlt(double scalar,VEC *vector, VEC *out) { /* u_int dim, i; */ /* Real *out_ve, *vec_ve; */ @@ -87,8 +83,7 @@ VEC *vector,*out; } /* v_add -- vector addition -- may be in-situ */ -VEC *v_add(vec1,vec2,out) -VEC *vec1,*vec2,*out; +VEC *v_add(VEC *vec1, VEC *vec2, VEC *out) { u_int dim; /* Real *out_ve, *vec1_ve, *vec2_ve; */ @@ -113,9 +108,7 @@ VEC *vec1,*vec2,*out; /* v_mltadd -- scalar/vector multiplication and addition -- out = v1 + scale.v2 */ -VEC *v_mltadd(v1,v2,scale,out) -VEC *v1,*v2,*out; -double scale; +VEC *v_mltadd(VEC *v1, VEC *v2, double scale, VEC *out) { /* register u_int dim, i; */ /* Real *out_ve, *v1_ve, *v2_ve; */ @@ -152,8 +145,7 @@ double scale; } /* v_sub -- vector subtraction -- may be in-situ */ -VEC *v_sub(vec1,vec2,out) -VEC *vec1,*vec2,*out; +VEC *v_sub(VEC *vec1, VEC *vec2, VEC *out) { /* u_int i, dim; */ /* Real *out_ve, *vec1_ve, *vec2_ve; */ @@ -179,13 +171,7 @@ VEC *vec1,*vec2,*out; /* v_map -- maps function f over components of x: out[i] = f(x[i]) -- _v_map sets out[i] = f(params,x[i]) */ -VEC *v_map(f,x,out) -#ifdef PROTOTYPES_IN_STRUCT -double (*f)(double); -#else -double (*f)(); -#endif -VEC *x, *out; +VEC *v_map(double (*f)(double), VEC *x, VEC *out) { Real *x_ve, *out_ve; int i, dim; @@ -202,14 +188,7 @@ VEC *x, *out; return out; } -VEC *_v_map(f,params,x,out) -#ifdef PROTOTYPES_IN_STRUCT -double (*f)(void *,double); -#else -double (*f)(); -#endif -VEC *x, *out; -void *params; +VEC *_v_map(double (*f)(void *, double), void *params, VEC *x, VEC *out) { Real *x_ve, *out_ve; int i, dim; @@ -227,10 +206,10 @@ void *params; } /* v_lincomb -- returns sum_i a[i].v[i], a[i] real, v[i] vectors */ -VEC *v_lincomb(n,v,a,out) -int n; /* number of a's and v's */ -Real a[]; -VEC *v[], *out; +VEC *v_lincomb(int n, VEC *v[], Real a[], VEC *out) +//int n; /* number of a's and v's */ +//Real a[]; +//VEC *v[], *out; { int i; @@ -350,8 +329,7 @@ VEC *v_linlist(va_alist) va_dcl /* v_star -- computes componentwise (Hadamard) product of x1 and x2 -- result out is returned */ -VEC *v_star(x1, x2, out) -VEC *x1, *x2, *out; +VEC *v_star(VEC *x1, VEC *x2, VEC *out) { int i; @@ -371,8 +349,7 @@ VEC *x1, *x2, *out; -- out[i] = x2[i] / x1[i] -- if x1[i] == 0 for some i, then raise E_SING error -- result out is returned */ -VEC *v_slash(x1, x2, out) -VEC *x1, *x2, *out; +VEC *v_slash(VEC *x1, VEC *x2, VEC *out) { int i; Real tmp; @@ -396,9 +373,7 @@ VEC *x1, *x2, *out; /* v_min -- computes minimum component of x, which is returned -- also sets min_idx to the index of this minimum */ -double v_min(x, min_idx) -VEC *x; -int *min_idx; +double v_min(VEC *x, int *min_idx) { int i, i_min; Real min_val, tmp; @@ -426,9 +401,7 @@ int *min_idx; /* v_max -- computes maximum component of x, which is returned -- also sets max_idx to the index of this maximum */ -double v_max(x, max_idx) -VEC *x; -int *max_idx; +double v_max(VEC *x, int *max_idx) { int i, i_max; Real max_val, tmp; @@ -462,9 +435,7 @@ int *max_idx; the permutation is order = [2, 0, 1]. -- if order is NULL on entry then it is ignored -- the sorted vector x is returned */ -VEC *v_sort(x, order) -VEC *x; -PERM *order; +VEC *v_sort(VEC *x, PERM *order) { Real *x_ve, tmp, v; /* int *order_pe; */ @@ -541,8 +512,7 @@ PERM *order; } /* v_sum -- returns sum of entries of a vector */ -double v_sum(x) -VEC *x; +double v_sum(VEC *x) { int i; Real sum; @@ -558,8 +528,7 @@ VEC *x; } /* v_conv -- computes convolution product of two vectors */ -VEC *v_conv(x1, x2, out) -VEC *x1, *x2, *out; +VEC *v_conv(VEC *x1, VEC *x2, VEC *out) { int i; @@ -580,8 +549,7 @@ VEC *x1, *x2, *out; /* v_pconv -- computes a periodic convolution product -- the period is the dimension of x2 */ -VEC *v_pconv(x1, x2, out) -VEC *x1, *x2, *out; +VEC *v_pconv(VEC *x1, VEC *x2, VEC *out) { int i;