#include #define RVEC(A) int A##n, double*A##p #define RMAT(A) int A##r, int A##c, double* A##p #define KRVEC(A) int A##n, const double*A##p #define KRMAT(A) int A##r, int A##c, const double* A##p #define CVEC(A) int A##n, gsl_complex*A##p #define CMAT(A) int A##r, int A##c, gsl_complex* A##p #define KCVEC(A) int A##n, const gsl_complex*A##p #define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p #define FVEC(A) int A##n, float*A##p #define FMAT(A) int A##r, int A##c, float* A##p #define KFVEC(A) int A##n, const float*A##p #define KFMAT(A) int A##r, int A##c, const float* A##p #define QVEC(A) int A##n, gsl_complex_float*A##p #define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p #define KQVEC(A) int A##n, const gsl_complex_float*A##p #define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p #include #include #include #include #include #include #define MACRO(B) do {B} while (0) #define ERROR(CODE) MACRO(return CODE;) #define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) #define OK return 0; #define MIN(A,B) ((A)<(B)?(A):(B)) #define MAX(A,B) ((A)>(B)?(A):(B)) #ifdef DBG #define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M); #else #define DEBUGMSG(M) #endif #define CHECK(RES,CODE) MACRO(if(RES) return CODE;) #ifdef DBG #define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n"); #else #define DEBUGMAT(MSG,X) #endif #ifdef DBG #define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n"); #else #define DEBUGVEC(MSG,X) #endif #define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n) #define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c) #define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n) #define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c) #define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n) #define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c) #define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n) #define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c) #define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n) #define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c) #define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n) #define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c) #define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n) #define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c) #define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n) #define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c) #define V(a) (&a.vector) #define M(a) (&a.matrix) #define GCVEC(A) int A##n, gsl_complex*A##p #define KGCVEC(A) int A##n, const gsl_complex*A##p #define GQVEC(A) int A##n, gsl_complex_float*A##p #define KGQVEC(A) int A##n, const gsl_complex_float*A##p #define BAD_SIZE 2000 #define BAD_CODE 2001 #define MEM 2002 #define BAD_FILE 2003 int sumF(KFVEC(x),FVEC(r)) { DEBUGMSG("sumF"); REQUIRES(rn==1,BAD_SIZE); int i; float res = 0; for (i = 0; i < xn; i++) res += xp[i]; rp[0] = res; OK } int sumR(KRVEC(x),RVEC(r)) { DEBUGMSG("sumR"); REQUIRES(rn==1,BAD_SIZE); int i; double res = 0; for (i = 0; i < xn; i++) res += xp[i]; rp[0] = res; OK } int sumQ(KQVEC(x),QVEC(r)) { DEBUGMSG("sumQ"); REQUIRES(rn==1,BAD_SIZE); int i; gsl_complex_float res; res.dat[0] = 0; res.dat[1] = 0; for (i = 0; i < xn; i++) { res.dat[0] += xp[i].dat[0]; res.dat[1] += xp[i].dat[1]; } rp[0] = res; OK } int sumC(KCVEC(x),CVEC(r)) { DEBUGMSG("sumC"); REQUIRES(rn==1,BAD_SIZE); int i; gsl_complex res; res.dat[0] = 0; res.dat[1] = 0; for (i = 0; i < xn; i++) { res.dat[0] += xp[i].dat[0]; res.dat[1] += xp[i].dat[1]; } rp[0] = res; OK } int prodF(KFVEC(x),FVEC(r)) { DEBUGMSG("prodF"); REQUIRES(rn==1,BAD_SIZE); int i; float res = 1; for (i = 0; i < xn; i++) res *= xp[i]; rp[0] = res; OK } int prodR(KRVEC(x),RVEC(r)) { DEBUGMSG("prodR"); REQUIRES(rn==1,BAD_SIZE); int i; double res = 1; for (i = 0; i < xn; i++) res *= xp[i]; rp[0] = res; OK } int prodQ(KQVEC(x),QVEC(r)) { DEBUGMSG("prodQ"); REQUIRES(rn==1,BAD_SIZE); int i; gsl_complex_float res; float temp; res.dat[0] = 1; res.dat[1] = 0; for (i = 0; i < xn; i++) { temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; res.dat[0] = temp; } rp[0] = res; OK } int prodC(KCVEC(x),CVEC(r)) { DEBUGMSG("prodC"); REQUIRES(rn==1,BAD_SIZE); int i; gsl_complex res; double temp; res.dat[0] = 1; res.dat[1] = 0; for (i = 0; i < xn; i++) { temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; res.dat[0] = temp; } rp[0] = res; OK } int toScalarR(int code, KRVEC(x), RVEC(r)) { REQUIRES(rn==1,BAD_SIZE); DEBUGMSG("toScalarR"); KDVVIEW(x); double res; switch(code) { case 0: { res = gsl_blas_dnrm2(V(x)); break; } case 1: { res = gsl_blas_dasum(V(x)); break; } case 2: { res = gsl_vector_max_index(V(x)); break; } case 3: { res = gsl_vector_max(V(x)); break; } case 4: { res = gsl_vector_min_index(V(x)); break; } case 5: { res = gsl_vector_min(V(x)); break; } default: ERROR(BAD_CODE); } rp[0] = res; OK } int toScalarF(int code, KFVEC(x), FVEC(r)) { REQUIRES(rn==1,BAD_SIZE); DEBUGMSG("toScalarF"); KFVVIEW(x); float res; switch(code) { case 0: { res = gsl_blas_snrm2(V(x)); break; } case 1: { res = gsl_blas_sasum(V(x)); break; } case 2: { res = gsl_vector_float_max_index(V(x)); break; } case 3: { res = gsl_vector_float_max(V(x)); break; } case 4: { res = gsl_vector_float_min_index(V(x)); break; } case 5: { res = gsl_vector_float_min(V(x)); break; } default: ERROR(BAD_CODE); } rp[0] = res; OK } int toScalarC(int code, KCVEC(x), RVEC(r)) { REQUIRES(rn==1,BAD_SIZE); DEBUGMSG("toScalarC"); KCVVIEW(x); double res; switch(code) { case 0: { res = gsl_blas_dznrm2(V(x)); break; } case 1: { res = gsl_blas_dzasum(V(x)); break; } default: ERROR(BAD_CODE); } rp[0] = res; OK } int toScalarQ(int code, KQVEC(x), FVEC(r)) { REQUIRES(rn==1,BAD_SIZE); DEBUGMSG("toScalarQ"); KQVVIEW(x); float res; switch(code) { case 0: { res = gsl_blas_scnrm2(V(x)); break; } case 1: { res = gsl_blas_scasum(V(x)); break; } default: ERROR(BAD_CODE); } rp[0] = res; OK } inline double sign(double x) { if(x>0) { return +1.0; } else if (x<0) { return -1.0; } else { return 0.0; } } inline float float_sign(float x) { if(x>0) { return +1.0; } else if (x<0) { return -1.0; } else { return 0.0; } } inline gsl_complex complex_abs(gsl_complex z) { gsl_complex r; r.dat[0] = gsl_complex_abs(z); r.dat[1] = 0; return r; } inline gsl_complex complex_signum(gsl_complex z) { gsl_complex r; double mag; if (z.dat[0] == 0 && z.dat[1] == 0) { r.dat[0] = 0; r.dat[1] = 0; } else { mag = gsl_complex_abs(z); r.dat[0] = z.dat[0]/mag; r.dat[1] = z.dat[1]/mag; } return r; } #define OP(C,F) case C: { for(k=0;k