diff --git a/Makefile.dspl b/Makefile.dspl index 27f8129..534ff84 100644 --- a/Makefile.dspl +++ b/Makefile.dspl @@ -14,7 +14,7 @@ all: $(RELEASE_DIR)/$(LIB_NAME) $(RELEASE_DIR)/$(LIB_NAME): $(DSPL_OBJ_FILES) $(BLAS_LIB_NAME) $(LAPACK_LIB_NAME) - $(CC) -shared -o $(RELEASE_DIR)/$(LIB_NAME) $(DSPL_OBJ_FILES) -lm -L$(BLAS_LIB_DIR) -lblas -L$(LAPACK_LIB_DIR) -llapack + $(CC) -shared -o $(RELEASE_DIR)/$(LIB_NAME) $(DSPL_OBJ_FILES) -lm -L$(LAPACK_LIB_DIR) -llapack -L$(BLAS_LIB_DIR) -lblas -lgfortran -lquadmath $(DSPL_OBJ_DIR)/%.o:$(DSPL_SRC_DIR)/%.c @@ -32,4 +32,4 @@ clean: rm -f $(BLAS_OBJ_DIR)/*.o rm -f $(BLAS_OBJ_DIR)/*.a rm -f $(RELEASE_DIR)/$(LIB_NAME) - \ No newline at end of file + diff --git a/dspl.project.win.geany b/dspl.project.win.geany index 3329281..fed0662 100644 --- a/dspl.project.win.geany +++ b/dspl.project.win.geany @@ -28,21 +28,23 @@ long_line_behaviour=1 long_line_column=72 [files] -current_page=7 +current_page=14 FILE_NAME_0=360;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Ccheby_poly1_test.c;0;2 FILE_NAME_1=507;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Ciir_test.c;0;2 FILE_NAME_2=80;None;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Cbin%5Cgnuplot%5Ciir_test.plt;0;2 FILE_NAME_3=1097;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Cru%5Ccontent%5Cfourier_series%5Cc%5Cfourier_series_dirichlet_ex.c;0;2 FILE_NAME_4=2672;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Cru%5Ccontent%5Cfourier_series%5Cc%5Cfourier_series_pimp_spectrum.c;0;2 FILE_NAME_5=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Cru%5Ccontent%5Cfourier_series%5Cc%5Cfourier_series_rec.c;0;2 -FILE_NAME_6=12388;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Carray.c;0;2 -FILE_NAME_7=424;C++;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cblas.h;0;2 -FILE_NAME_8=1737;F77;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cblas_src%5Cddot.f;0;2 -FILE_NAME_9=51351;C++;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cinclude%5Cdspl.h;0;2 -FILE_NAME_10=13061;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cinclude%5Cdspl.c;0;2 -FILE_NAME_11=323;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Carray_test.c;0;2 -FILE_NAME_12=1454;F77;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cblas_src%5Cdasum.f;0;2 -FILE_NAME_13=419;F77;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cblas_src%5Cdaxpy.f;0;2 +FILE_NAME_6=12457;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Carray.c;0;2 +FILE_NAME_7=28861;C++;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cblas.h;0;2 +FILE_NAME_8=42098;C++;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cinclude%5Cdspl.h;0;2 +FILE_NAME_9=12508;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cinclude%5Cdspl.c;0;2 +FILE_NAME_10=323;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Carray_test.c;0;2 +FILE_NAME_11=1012;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cmatrix.c;0;2 +FILE_NAME_12=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cmatrix_print.c;0;2 +FILE_NAME_13=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cmatrix_transpose.c;0;2 +FILE_NAME_14=131;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cmatrix_mul.c;0;2 +FILE_NAME_15=523;Make;0;EUTF-8;1;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5CMakefile.dspl;0;2 [build-menu] NF_00_LB=_Собрать diff --git a/dspl/dspl_src/blas.h b/dspl/dspl_src/blas.h index 4506ac9..b0b67d6 100644 --- a/dspl/dspl_src/blas.h +++ b/dspl/dspl_src/blas.h @@ -3,7 +3,6 @@ #ifndef BLAS_H #define BLAS_H - #define FORTRAN_FUNC(FUNC) FUNC##_ int FORTRAN_FUNC(xerbla)(const char*, int*info, int); @@ -361,7 +360,7 @@ int FORTRAN_FUNC(zher2m)(const char*, const char*, const char*, const int*, cons int FORTRAN_FUNC(xher2m)(const char*, const char*, const char*, const int*, const int*, const double*, const double*, const int*, const double*, const int*, const double*, double*, const int*); - +void FORTRAN_FUNC(zgees)(const char*, const char*, void*, int*, complex_t*, int*, int*, complex_t*, complex_t*, int*, complex_t*, int*, double*, int*, int*); #endif diff --git a/dspl/dspl_src/dspl_internal.h b/dspl/dspl_src/dspl_internal.h index 1388504..bfee1e4 100644 --- a/dspl/dspl_src/dspl_internal.h +++ b/dspl/dspl_src/dspl_internal.h @@ -30,10 +30,6 @@ // sqrt(2^31) #define FFT_COMPOSITE_MAX 46340 -void transpose(double* a, int n, int m, double* b); -void transpose_cmplx(complex_t* a, int n, int m, complex_t* b); -void transpose_hermite(complex_t* a, int n, int m, complex_t* b); - int fft_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr); diff --git a/dspl/dspl_src/fft.c b/dspl/dspl_src/fft.c index ad0c02b..0f50b0e 100644 --- a/dspl/dspl_src/fft.c +++ b/dspl/dspl_src/fft.c @@ -150,7 +150,7 @@ label_size: if(n2>1) { memcpy(t1, t0, n*sizeof(complex_t)); - transpose_cmplx(t1, n2, n1, t0); + matrix_transpose_cmplx(t1, n2, n1, t0); } if(n1 == 16) @@ -186,13 +186,13 @@ label_size: IM(t0[k]) = CMIM(t1[k], pw[k]); } - transpose_cmplx(t0, n1, n2, t1); + matrix_transpose_cmplx(t0, n1, n2, t1); for(k = 0; k < n1; k++) { fft_krn(t1+k*n2, t0+k*n2, p, n2, addr+n); } - transpose_cmplx(t0, n2, n1, t1); + matrix_transpose_cmplx(t0, n2, n1, t1); } } return RES_OK; diff --git a/dspl/dspl_src/matrix.c b/dspl/dspl_src/matrix.c index 11249f7..9807565 100644 --- a/dspl/dspl_src/matrix.c +++ b/dspl/dspl_src/matrix.c @@ -27,94 +27,93 @@ -/******************************************************************************* -matrix_create -*******************************************************************************/ -int DSPL_API matrix_create(matrix_t* a, int n, int m, int type) -{ - if(!a) - return ERROR_PTR; - if(n < 1 || m < 1) - return ERROR_MATRIX_SIZE; - if(a->dat) +int DSPL_API matrix_eig_cmplx(complex_t* a, int n, complex_t* v, int* info) +{ + int err; + int sdim = 0; + int ldvs = 1; + int lwork = 2*n; + if(!a || !v) + return ERROR_PTR; + + if(n<1) + return ERROR_MATRIX_SIZE; + + complex_t *work=(complex_t*)malloc(lwork*sizeof(complex_t)); + double *rwork = (double*)malloc(n*sizeof(double)); + + zgees_("N", "N", NULL, &n, a, &n, &sdim, v, NULL, &ldvs, work, &lwork, + rwork, NULL, &err); + + if(err!=0) { - a->dat = (type & DAT_MASK) ? - (void*) realloc(a->dat, n*m*sizeof(complex_t)): - (void*) realloc(a->dat, n*m*sizeof(double)); + if(info) + *info = err; + err = ERROR_LAPACK; } else - { - a->dat = (type & DAT_MASK) ? - (void*) malloc(n*m*sizeof(complex_t)): - (void*) malloc(n*m*sizeof(double)); - } - - a->n = n; - a->m = m; - a->type = type; - - return RES_OK; + err = RES_OK; + + free(work); + free(rwork); + return err; } -/******************************************************************************* -matrix_create eye -*******************************************************************************/ -int DSPL_API matrix_create_eye(matrix_t* a, int n, int type) -{ - double *pr; - complex_t *pc; - int err, m, k; - - err = matrix_create(a, n, n, type); - if(err != RES_OK) - return RES_OK; - - k = 0; - if((a->type & DAT_MASK) == DAT_DOUBLE) - { - pr = (double*) a->dat; - memset(pr, 0, n*n*sizeof(double)); - for(m = 0; m < n; m++) - { - pr[k] = 1.0; - k += n+1; - } - } - - if((a->type & DAT_MASK) == DAT_COMPLEX) - { - pc = (complex_t*) a->dat; - memset(pc, 0, n*n*sizeof(complex_t)); - for(m = 0; m < n; m++) - { - RE(pc[k]) = 1.0; - k += n+1; - } - } - return RES_OK; -} - - /******************************************************************************* -matrix_free +Real matrix eye *******************************************************************************/ -void DSPL_API matrix_free(matrix_t* a) +int DSPL_API matrix_eye(double* a, int n, int m) { + int p, k; if(!a) - return; - if(a->dat) - free(a->dat); - a->n = a->m = a->type = 0; + return ERROR_PTR; + if (n < 1 || m < 1) + return ERROR_MATRIX_SIZE; + + k = 0; + memset(a, 0, n*m*sizeof(double)); + for(p = 0; p < m; p++) + { + a[k] = 1.0; + k += n+1; + } + + return RES_OK; } +/******************************************************************************* +Complex matrix eye +*******************************************************************************/ +int DSPL_API matrix_eye_cmplx(complex_t* a, int n, int m) +{ + int p, k; + if(!a) + return ERROR_PTR; + if (n < 1 || m < 1) + return ERROR_MATRIX_SIZE; + + k = 0; + memset(a, 0, n*m*sizeof(complex_t)); + for(p = 0; p < m; p++) + { + RE(a[k]) = 1.0; + k += n+1; + } + + return RES_OK; +} + + + + /******************************************************************************* matrix LU decomposition -*******************************************************************************/ +****************************************************************************** int DSPL_API matrix_lu(matrix_t* a, matrix_t* L, matrix_t* U, matrix_t* P) { int err, k, n, m, N, ind; @@ -174,159 +173,99 @@ int DSPL_API matrix_lu(matrix_t* a, matrix_t* L, matrix_t* U, matrix_t* P) return RES_OK; } +*/ + + +/******************************************************************************* +real matrix multiplication +*******************************************************************************/ +int DSPL_API matrix_mul(double* a, int na, int ma, + double* b, int nb, int mb, + double* c) +{ + + double alpha = 1; + double beta = 0.0; + + if(!a || !b || !c) + return ERROR_PTR; + if(na < 1 || ma < 1 || nb < 1 || mb < 1 || ma != nb) + return ERROR_MATRIX_SIZE; + + /* BLAS DGEMM */ + dgemm_("N", "N", &na, &mb, &ma, &alpha, a, &na, b, &nb, &beta, c, &na); + + return RES_OK; +} /******************************************************************************* -matrix transposition +real matrix print *******************************************************************************/ -int DSPL_API matrix_print(matrix_t* a, const char* name, const char* format) +int DSPL_API matrix_print(double* a, int n, int m, + const char* name, const char* format) { - int n,m; + int p,q; + if(!a) return ERROR_PTR; - if(!a->dat) - return ERROR_PTR; - - if((a->type & DAT_MASK) == DAT_DOUBLE) + if(n < 1 || m < 1) + return ERROR_SIZE; + + printf("\n%s = [ %% size [%d x %d] type: real", name, n, m); + + for(p = 0; p < n; p++) { - printf("\nMatrix %s size [%d x %d] type: real\n", - name, a->n, a->m); - double* p = (double*)(a->dat); - for(n = 0; n < a->n; n++) + printf("\n"); + for(q = 0; q < m; q++) { - for(m = 0; m < a->m; m++) - { - printf(format, p[m*a->n + n]); - } - printf("\n"); + printf(format, a[q*n + p]); + if(q == m-1) + printf(";"); + else + printf(", "); } } + printf("];\n"); + + return RES_OK; +} - if((a->type & DAT_MASK) == DAT_COMPLEX) +/******************************************************************************* +complex matrix print +*******************************************************************************/ +int DSPL_API matrix_print_cmplx(complex_t* a, int n, int m, + const char* name, const char* format) +{ + int p,q; + + if(!a) + return ERROR_PTR; + if(n < 1 || m < 1) + return ERROR_MATRIX_SIZE; + + if(!a) + return ERROR_PTR; + if(n < 1 || m < 1) + return ERROR_SIZE; + + printf("\n%s = [ %% size [%d x %d] type: complex", name, n, m); + + for(p = 0; p < n; p++) { - printf("\nMatrix %s size [%d x %d] type: complex\n", - name, a->n, a->m); - complex_t* p = (complex_t*)(a->dat); - for(n = 0; n < a->n; n++) + printf("\n"); + for(q = 0; q < m; q++) { - for(m = 0; m < a->m; m++) - { - printf(format, RE(p[m*a->n + n]), IM(p[m*a->n + n])); - } - printf("\n"); + printf(format, RE(a[q*n + p]), IM(a[q*n + p])); + if(q == m-1) + printf(";"); + else + printf(", "); } } - return RES_OK; -} - - -/******************************************************************************* -matrix swap 2 elements -*******************************************************************************/ -int DSPL_API matrix_swap(matrix_t* a, int r0, int c0, int r1, int c1) -{ - double tr; - complex_t tc; - double *pr; - complex_t *pc; - - if(!a) - return ERROR_PTR; - if(r0 >= a->n || r1 >= a->n || c0 >= a->m || c1 >= a->m) - return ERROR_MATRIX_INDEX; - - if((a->type & DAT_MASK) == DAT_DOUBLE) - { - pr = (double*)(a->dat); - tr = pr[r0 + c0 * a->n]; - pr[r0 + c0 * a->n] = pr[r1 + c1 * a->n]; - pr[r1 + c1 * a->n] = tr; - } - - if((a->type & DAT_MASK) == DAT_COMPLEX) - { - pc = (complex_t*)(a->dat); - RE(tc) = RE(pc[r0 + c0 * a->n]); - IM(tc) = IM(pc[r0 + c0 * a->n]); - - RE(pc[r0 + c0 * a->n]) = RE(pc[r1 + c1 * a->n]); - IM(pc[r0 + c0 * a->n]) = IM(pc[r1 + c1 * a->n]); - - RE(pc[r1 + c1 * a->n]) = RE(tc); - IM(pc[r1 + c1 * a->n]) = IM(tc); - } - return RES_OK; -} - - - -/******************************************************************************* -matrix swap 2 rows -*******************************************************************************/ -int DSPL_API matrix_swap_rows(matrix_t* a, int r0, int r1) -{ - int c, err; - if(!a) - return ERROR_PTR; - if(r0 >= a->n || r1 >= a->n) - return ERROR_MATRIX_INDEX; - - for(c = 0; c < a->m; c++) - { - err = matrix_swap(a, r0, c, r1, c); - if(err != RES_OK) - break; - } - return err; -} - - - - -/******************************************************************************* -matrix transposition -*******************************************************************************/ -int DSPL_API matrix_transpose(matrix_t* a, matrix_t* b) -{ - int err; - if(!a || !b) - return ERROR_PTR; - - err = matrix_create(b, a->m, a->n, a->type); - if(err != RES_OK) - return err; - - if((a->type & DAT_MASK) == DAT_DOUBLE) - transpose((double*)(a->dat), a->n, a->m, (double*)(b->dat)); - if((a->type & DAT_MASK) == DAT_COMPLEX) - transpose_cmplx((complex_t*)(a->dat), a->n, a->m, (complex_t*)(b->dat)); - - return RES_OK; -} - - - - - -/******************************************************************************* -matrix Hermite transposition -*******************************************************************************/ -int DSPL_API matrix_transpose_hermite(matrix_t* a, matrix_t* b) -{ - int err; - if(!a || !b) - return ERROR_PTR; - - err = matrix_create(b, a->m, a->n, a->type); - if(err != RES_OK) - return err; - - if((a->type & DAT_MASK) == DAT_DOUBLE) - transpose((double*)(a->dat), a->n, a->m, (double*)(b->dat)); - if((a->type & DAT_MASK) == DAT_COMPLEX) - transpose_hermite((complex_t*)(a->dat), a->n, a->m, (complex_t*)(b->dat)); + printf("];\n"); return RES_OK; } @@ -336,14 +275,18 @@ int DSPL_API matrix_transpose_hermite(matrix_t* a, matrix_t* b) - /******************************************************************************* -Real matrx transpose +Real matrix transpose *******************************************************************************/ -void transpose(double* a, int n, int m, double* b) +int DSPL_API matrix_transpose(double* a, int n, int m, double* b) { int p, q, i, j, aind, bind; - + if(!a || !b) + return ERROR_PTR; + if(n < 1 || m < 1) + return ERROR_MATRIX_SIZE; + + for(p = 0; p < n - DSPL_MATRIX_BLOCK; p+=DSPL_MATRIX_BLOCK) { for(q = 0; q < m - DSPL_MATRIX_BLOCK; q+=DSPL_MATRIX_BLOCK) @@ -366,7 +309,8 @@ void transpose(double* a, int n, int m, double* b) for(i = 0; i < p; i++) for(j = q; j < m; j++) b[i*m + j] = a[j*n+i]; - + + return RES_OK; } @@ -374,12 +318,17 @@ void transpose(double* a, int n, int m, double* b) /******************************************************************************* -Complex matrx transpose +Complex matrix transpose *******************************************************************************/ -void transpose_cmplx(complex_t* a, int n, int m, complex_t* b) +int DSPL_API matrix_transpose_cmplx(complex_t* a, int n, int m, complex_t* b) { int p, q, i, j, aind, bind; - + + if(!a || !b) + return ERROR_PTR; + if(n < 1 || m < 1) + return ERROR_MATRIX_SIZE; + for(p = 0; p < n - DSPL_MATRIX_BLOCK; p+=DSPL_MATRIX_BLOCK) { for(q = 0; q < m - DSPL_MATRIX_BLOCK; q+=DSPL_MATRIX_BLOCK) @@ -413,6 +362,7 @@ void transpose_cmplx(complex_t* a, int n, int m, complex_t* b) IM(b[i*m + j]) = IM(a[j*n+i]); } } + return RES_OK; } @@ -420,12 +370,17 @@ void transpose_cmplx(complex_t* a, int n, int m, complex_t* b) /******************************************************************************* -Hermite matrx transpose +Hermite matrix transpose *******************************************************************************/ -void transpose_hermite(complex_t* a, int n, int m, complex_t* b) +int DSPL_API matrix_transpose_hermite(complex_t* a, int n, int m, complex_t* b) { int p, q, i, j, aind, bind; - + + if(!a || !b) + return ERROR_PTR; + if(n < 1 || m < 1) + return ERROR_MATRIX_SIZE; + for(p = 0; p < n - DSPL_MATRIX_BLOCK; p+=DSPL_MATRIX_BLOCK) { for(q = 0; q < m - DSPL_MATRIX_BLOCK; q+=DSPL_MATRIX_BLOCK) @@ -459,5 +414,7 @@ void transpose_hermite(complex_t* a, int n, int m, complex_t* b) IM(b[i*m + j]) = -IM(a[j*n+i]); } } + + return RES_OK; } diff --git a/examples/src/matrix_eig.c b/examples/src/matrix_eig.c new file mode 100644 index 0000000..917282a --- /dev/null +++ b/examples/src/matrix_eig.c @@ -0,0 +1,32 @@ +#include +#include +#include +#include "dspl.h" + +#define N 3 + +int main() +{ + void* handle; /* DSPL handle */ + handle = dspl_load(); /* Load DSPL function */ + + complex_t a[N*N] = {{1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}, + {2.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}, + {3.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}; + + complex_t v[N]; + int err, info; + matrix_print_cmplx(a, N, N, "A", "%8.2f%+8.2fi"); + + err = matrix_eig_cmplx(a, N, v, &info); + if(err!=RES_OK) + { + printf("ERROR CODE: 0x%.8x, info = %d\n", err, info); + } + + matrix_print_cmplx(v, N, 1, "v", "%10.6f%+10.6fi"); + + dspl_free(handle); /* free dspl handle */ + return 0; +} + diff --git a/examples/src/matrix_mul.c b/examples/src/matrix_mul.c new file mode 100644 index 0000000..52d9958 --- /dev/null +++ b/examples/src/matrix_mul.c @@ -0,0 +1,31 @@ +#include +#include +#include +#include "dspl.h" + +#define N 4 +#define M 3 +#define K 10 + +int main() +{ + void* handle; /* DSPL handle */ + handle = dspl_load(); /* Load DSPL function */ + + double a[N*K], b[K*M], c[N*M]; + + + linspace(0, N*K, N*K, DSPL_PERIODIC, a); + matrix_print(a, N, K, "A", "%8.2f"); + + linspace(0, K*M, K*M, DSPL_PERIODIC, b); + matrix_print(b, K, M, "B", "%8.2f"); + + matrix_mul(a, N, K, b, K, M, c); + matrix_print(c, N, M, "C", "%8.2f"); + + + dspl_free(handle); /* free dspl handle */ + return 0; +} + diff --git a/examples/src/matrix_print.c b/examples/src/matrix_print.c index 5df607f..6df81d2 100644 --- a/examples/src/matrix_print.c +++ b/examples/src/matrix_print.c @@ -3,31 +3,30 @@ #include #include "dspl.h" +#define N 4 +#define M 3 + int main() { - void* handle; // DSPL handle - handle = dspl_load(); // Load DSPL function - - matrix_t r, c; - - memset(&r, 0, sizeof(matrix_t)); - memset(&c, 0, sizeof(matrix_t)); - - matrix_create(&r, 4, 3, DAT_DOUBLE); - linspace(0, 12, 12, DSPL_PERIODIC, (double*)r.dat); - matrix_print(&r, "R", "%8.2f"); - - matrix_create(&c, 2, 3, DAT_COMPLEX); - linspace(0, 12, 12, DSPL_PERIODIC, (double*)c.dat); - matrix_print(&c, "C", "%8.2f%+8.2fj "); - - - matrix_free(&r); - matrix_free(&c); - dspl_free(handle); // free dspl handle - - // выполнить скрипт GNUPLOT для построения графиков - // по рассчитанным данным - return system("gnuplot gnuplot/sinc_test.plt");; + void* handle; /* DSPL handle */ + handle = dspl_load(); /* Load DSPL function */ + + double r[N*M]; + complex_t c[N*M]; + + linspace(0, N*M, N*M, DSPL_PERIODIC, r); + matrix_print(r, N, M, "R", "%8.2f"); + + linspace(0, N*M*2, N*M*2, DSPL_PERIODIC, (double*)c); + matrix_print_cmplx(c, N, M, "C", "%8.2f%+8.2fi"); + + matrix_eye(r, N, M); + matrix_print(r, N, M, "I", "%8.2f"); + + matrix_eye_cmplx(c, N, M); + matrix_print_cmplx(c, N, M, "I", "%8.2f%+8.2fi"); + + dspl_free(handle); /* free dspl handle */ + return 0; } diff --git a/examples/src/matrix_transpose.c b/examples/src/matrix_transpose.c index d3535e0..0e9219f 100644 --- a/examples/src/matrix_transpose.c +++ b/examples/src/matrix_transpose.c @@ -3,48 +3,37 @@ #include #include "dspl.h" +#define N 4 +#define M 3 + int main() { - void* handle; // DSPL handle - handle = dspl_load(); // Load DSPL function + + void* handle; /* DSPL handle */ + handle = dspl_load(); /* Load DSPL function */ + + double r[N*M], p[N*M]; + complex_t c[N*M], d[N*M]; + + linspace(0, N*M, N*M, DSPL_PERIODIC, r); + matrix_print(r, N, M, "R", "%8.2f"); - matrix_t a, b, c; - - memset(&a, 0, sizeof(matrix_t)); - memset(&b, 0, sizeof(matrix_t)); - memset(&c, 0, sizeof(matrix_t)); - - matrix_create(&a, 4, 3, DAT_DOUBLE); - linspace(0, 12, 12, DSPL_PERIODIC, (double*)a.dat); - matrix_print(&a, "A", "%8.2f"); - - matrix_transpose(&a, &b); - matrix_print(&b, "B=A^T", "%8.2f"); - - matrix_transpose_hermite(&a, &b); - matrix_print(&b, "B=A^H", "%8.2f"); - - - - matrix_create(&c, 2, 3, DAT_COMPLEX); - linspace(0, 12, 12, DSPL_PERIODIC, (double*)c.dat); - matrix_print(&c, "C", "%8.2f%+8.2fj "); - - matrix_transpose(&c, &b); - matrix_print(&b, "B=C^T", "%8.2f%+8.2fj "); - - - matrix_transpose_hermite(&c, &b); - matrix_print(&b, "B=C^H", "%8.2f%+8.2fj "); - - - matrix_free(&a); - matrix_free(&b); - matrix_free(&c); - dspl_free(handle); // free dspl handle - - // выполнить скрипт GNUPLOT для построения графиков - // по рассчитанным данным - return system("gnuplot gnuplot/sinc_test.plt");; + matrix_transpose(r, N, M, p); + matrix_print(p, M, N, "P", "%8.2f"); + + linspace(0, N*M*2, N*M*2, DSPL_PERIODIC, (double*)c); + matrix_print_cmplx(c, N, M, "C", "%8.2f%+8.2fi"); + + matrix_transpose_cmplx(c, N, M, d); + matrix_print_cmplx(d, M, N, "D", "%8.2f%+8.2fi"); + + matrix_transpose_hermite(c, N, M, d); + matrix_print_cmplx(d, M, N, "D", "%8.2f%+8.2fi"); + + + dspl_free(handle); /* free dspl handle */ + + return 0; + } diff --git a/include/dspl.c b/include/dspl.c index c2f8c34..fccfe67 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -114,14 +114,18 @@ p_low2bp low2bp ; p_low2bs low2bs ; p_low2high low2high ; p_low2low low2low ; -p_matrix_create matrix_create ; -p_matrix_create_eye matrix_create_eye ; -p_matrix_free matrix_free ; + +p_matrix_eig_cmplx matrix_eig_cmplx ; +p_matrix_eye matrix_eye ; +p_matrix_eye_cmplx matrix_eye_cmplx ; +p_matrix_mul matrix_mul ; p_matrix_print matrix_print ; -p_matrix_swap matrix_swap ; -p_matrix_swap_rows matrix_swap_rows ; +p_matrix_print_cmplx matrix_print_cmplx ; p_matrix_transpose matrix_transpose ; +p_matrix_transpose_cmplx matrix_transpose_cmplx ; p_matrix_transpose_hermite matrix_transpose_hermite ; + + p_minmax minmax ; p_poly_z2a_cmplx poly_z2a_cmplx ; p_polyval polyval ; @@ -282,13 +286,14 @@ void* dspl_load() LOAD_FUNC(low2bs); LOAD_FUNC(low2high); LOAD_FUNC(low2low); - LOAD_FUNC(matrix_create); - LOAD_FUNC(matrix_create_eye); - LOAD_FUNC(matrix_free); + LOAD_FUNC(matrix_eig_cmplx); + LOAD_FUNC(matrix_eye); + LOAD_FUNC(matrix_eye_cmplx); + LOAD_FUNC(matrix_mul); LOAD_FUNC(matrix_print); - LOAD_FUNC(matrix_swap); - LOAD_FUNC(matrix_swap_rows); + LOAD_FUNC(matrix_print_cmplx); LOAD_FUNC(matrix_transpose); + LOAD_FUNC(matrix_transpose_cmplx); LOAD_FUNC(matrix_transpose_hermite); LOAD_FUNC(minmax); LOAD_FUNC(poly_z2a_cmplx); diff --git a/include/dspl.h b/include/dspl.h index 2bfa4fe..1ee2ea9 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -126,6 +126,7 @@ typedef struct /* J 0x10xxxxxx*/ /* K 0x11xxxxxx*/ /* L 0x12xxxxxx*/ +#define ERROR_LAPACK 0x12011601 /* M 0x13xxxxxx*/ #define ERROR_MATRIX_INDEX 0x13010914 #define ERROR_MATRIX_SINGULAR 0x13011914 @@ -715,41 +716,53 @@ DECLARE_FUNC(int, low2low, double* b COMMA double* beta COMMA double* alpha); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_create, matrix_t* a +DECLARE_FUNC(int, matrix_eig_cmplx, complex_t* a COMMA int n - COMMA int m - COMMA int type); + COMMA complex_t* v + COMMA int* info); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_create_eye, matrix_t* a +DECLARE_FUNC(int, matrix_eye, double* a COMMA int n - COMMA int type); + COMMA int m); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(void, matrix_free, matrix_t* a); +DECLARE_FUNC(int, matrix_eye_cmplx, complex_t* a + COMMA int n + COMMA int m); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_lu, matrix_t* a - COMMA matrix_t* L - COMMA matrix_t* U - COMMA matrix_t* P); +DECLARE_FUNC(int, matrix_mul, double* a + COMMA int na + COMMA int ma + COMMA double* b + COMMA int nb + COMMA int mb + COMMA double* c); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_print, matrix_t* a +DECLARE_FUNC(int, matrix_print, double* a + COMMA int n + COMMA int m COMMA const char* name COMMA const char* format); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_swap, matrix_t* a - COMMA int r0 - COMMA int c0 - COMMA int r1 - COMMA int c1); +DECLARE_FUNC(int, matrix_print_cmplx, complex_t* a + COMMA int n + COMMA int m + COMMA const char* name + COMMA const char* format); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_swap_rows, matrix_t* a - COMMA int r0 - COMMA int r1); +DECLARE_FUNC(int, matrix_transpose, double* a + COMMA int n + COMMA int m + COMMA double* b); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_transpose, matrix_t* a - COMMA matrix_t* b); +DECLARE_FUNC(int, matrix_transpose_cmplx, complex_t* a + COMMA int n + COMMA int m + COMMA complex_t* b); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_transpose_hermite, matrix_t* a - COMMA matrix_t* b); +DECLARE_FUNC(int, matrix_transpose_hermite, complex_t* a + COMMA int n + COMMA int m + COMMA complex_t* b); /*----------------------------------------------------------------------------*/ DECLARE_FUNC(int, minmax, double* x COMMA int n diff --git a/release/include/dspl.c b/release/include/dspl.c index c2f8c34..fccfe67 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -114,14 +114,18 @@ p_low2bp low2bp ; p_low2bs low2bs ; p_low2high low2high ; p_low2low low2low ; -p_matrix_create matrix_create ; -p_matrix_create_eye matrix_create_eye ; -p_matrix_free matrix_free ; + +p_matrix_eig_cmplx matrix_eig_cmplx ; +p_matrix_eye matrix_eye ; +p_matrix_eye_cmplx matrix_eye_cmplx ; +p_matrix_mul matrix_mul ; p_matrix_print matrix_print ; -p_matrix_swap matrix_swap ; -p_matrix_swap_rows matrix_swap_rows ; +p_matrix_print_cmplx matrix_print_cmplx ; p_matrix_transpose matrix_transpose ; +p_matrix_transpose_cmplx matrix_transpose_cmplx ; p_matrix_transpose_hermite matrix_transpose_hermite ; + + p_minmax minmax ; p_poly_z2a_cmplx poly_z2a_cmplx ; p_polyval polyval ; @@ -282,13 +286,14 @@ void* dspl_load() LOAD_FUNC(low2bs); LOAD_FUNC(low2high); LOAD_FUNC(low2low); - LOAD_FUNC(matrix_create); - LOAD_FUNC(matrix_create_eye); - LOAD_FUNC(matrix_free); + LOAD_FUNC(matrix_eig_cmplx); + LOAD_FUNC(matrix_eye); + LOAD_FUNC(matrix_eye_cmplx); + LOAD_FUNC(matrix_mul); LOAD_FUNC(matrix_print); - LOAD_FUNC(matrix_swap); - LOAD_FUNC(matrix_swap_rows); + LOAD_FUNC(matrix_print_cmplx); LOAD_FUNC(matrix_transpose); + LOAD_FUNC(matrix_transpose_cmplx); LOAD_FUNC(matrix_transpose_hermite); LOAD_FUNC(minmax); LOAD_FUNC(poly_z2a_cmplx); diff --git a/release/include/dspl.h b/release/include/dspl.h index 2bfa4fe..1ee2ea9 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -126,6 +126,7 @@ typedef struct /* J 0x10xxxxxx*/ /* K 0x11xxxxxx*/ /* L 0x12xxxxxx*/ +#define ERROR_LAPACK 0x12011601 /* M 0x13xxxxxx*/ #define ERROR_MATRIX_INDEX 0x13010914 #define ERROR_MATRIX_SINGULAR 0x13011914 @@ -715,41 +716,53 @@ DECLARE_FUNC(int, low2low, double* b COMMA double* beta COMMA double* alpha); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_create, matrix_t* a +DECLARE_FUNC(int, matrix_eig_cmplx, complex_t* a COMMA int n - COMMA int m - COMMA int type); + COMMA complex_t* v + COMMA int* info); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_create_eye, matrix_t* a +DECLARE_FUNC(int, matrix_eye, double* a COMMA int n - COMMA int type); + COMMA int m); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(void, matrix_free, matrix_t* a); +DECLARE_FUNC(int, matrix_eye_cmplx, complex_t* a + COMMA int n + COMMA int m); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_lu, matrix_t* a - COMMA matrix_t* L - COMMA matrix_t* U - COMMA matrix_t* P); +DECLARE_FUNC(int, matrix_mul, double* a + COMMA int na + COMMA int ma + COMMA double* b + COMMA int nb + COMMA int mb + COMMA double* c); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_print, matrix_t* a +DECLARE_FUNC(int, matrix_print, double* a + COMMA int n + COMMA int m COMMA const char* name COMMA const char* format); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_swap, matrix_t* a - COMMA int r0 - COMMA int c0 - COMMA int r1 - COMMA int c1); +DECLARE_FUNC(int, matrix_print_cmplx, complex_t* a + COMMA int n + COMMA int m + COMMA const char* name + COMMA const char* format); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_swap_rows, matrix_t* a - COMMA int r0 - COMMA int r1); +DECLARE_FUNC(int, matrix_transpose, double* a + COMMA int n + COMMA int m + COMMA double* b); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_transpose, matrix_t* a - COMMA matrix_t* b); +DECLARE_FUNC(int, matrix_transpose_cmplx, complex_t* a + COMMA int n + COMMA int m + COMMA complex_t* b); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, matrix_transpose_hermite, matrix_t* a - COMMA matrix_t* b); +DECLARE_FUNC(int, matrix_transpose_hermite, complex_t* a + COMMA int n + COMMA int m + COMMA complex_t* b); /*----------------------------------------------------------------------------*/ DECLARE_FUNC(int, minmax, double* x COMMA int n