@ -0,0 +1,45 @@ | |||||
#include "cplx.h" | |||||
#include <gsl/gsl_blas.h> | |||||
#include <gsl/gsl_complex_math.h> | |||||
gsl_vector_complex *cvec_dot_gslcplx(gsl_vector_complex *v, gsl_complex x) { | |||||
gsl_vector_complex *out = gsl_vector_complex_alloc(v->size); | |||||
for (size_t i=0; i < v->size; i++) { | |||||
gsl_vector_complex_set(out, i, gsl_complex_mul(x, gsl_vector_complex_get(v, i))); | |||||
} | |||||
return out; | |||||
} | |||||
gsl_vector_complex *vec_dot_gslcplx(gsl_vector *v, gsl_complex x) { | |||||
gsl_vector_complex *out = gsl_vector_complex_alloc(v->size); | |||||
for (size_t i=0; i < v->size; i++) { | |||||
gsl_vector_complex_set(out, i, gsl_complex_mul_real(x, gsl_vector_get(v, i))); | |||||
} | |||||
return out; | |||||
} | |||||
gsl_vector_complex *cvec_dot_c(gsl_vector_complex *v, complex double x) { | |||||
return cvec_dot_gslcplx(v, gsl_cplx_from_c99(x)); | |||||
} | |||||
gsl_vector_complex *vec_dot_c(gsl_vector *v, complex double x) { | |||||
return vec_dot_gslcplx(v, gsl_cplx_from_c99(x)); | |||||
} | |||||
complex double ddot(complex double x, complex double y) { | |||||
return x * y; | |||||
} | |||||
void gsl_vector_complex_print(gsl_vector_complex *v) { | |||||
for (size_t i = 0; i < v->size; i++) { | |||||
gsl_complex x = gsl_vector_complex_get(v, i); | |||||
printf("%4g + %4gi%c", GSL_REAL(x), GSL_IMAG(x), i < v->size-1 ? '\t' : '\n'); | |||||
} | |||||
} |
@ -0,0 +1,31 @@ | |||||
#include <complex.h> | |||||
#include <gsl/gsl_vector.h> | |||||
gsl_vector_complex *cvec_dot_gslcplx(gsl_vector_complex *v, gsl_complex x); | |||||
gsl_vector_complex *vec_dot_gslcplx(gsl_vector *v, gsl_complex x); | |||||
gsl_vector_complex *cvec_dot_c(gsl_vector_complex *v, complex double x); | |||||
gsl_vector_complex *vec_dot_c(gsl_vector *v, complex double x); | |||||
void gsl_vector_complex_print(gsl_vector_complex *v); | |||||
#define gsl_cplx_from_c99(x) (gsl_complex){.dat = {creal(x), cimag(x)}} | |||||
complex double ddot(complex double x, complex double y); | |||||
#define dot(x, y) _Generic((x), \ | |||||
gsl_vector* : dot_given_vec(y), \ | |||||
gsl_vector_complex* : dot_given_cplx_vec(y), \ | |||||
default : ddot) \ | |||||
((x), (y)) | |||||
#define dot_given_vec(y) _Generic ((y), \ | |||||
gsl_complex : vec_dot_gslcplx, \ | |||||
default : vec_dot_c) | |||||
#define dot_given_cplx_vec(y) _Generic((y), \ | |||||
gsl_complex : cvec_dot_gslcplx, \ | |||||
default : cvec_dot_c) |
@ -0,0 +1,31 @@ | |||||
#include <stdio.h> | |||||
#include "cplx.h" | |||||
int main() { | |||||
int complex a = 1 + 2I; | |||||
complex double b = 2 + I; | |||||
gsl_complex c = gsl_cplx_from_c99(a); | |||||
gsl_vector *v = gsl_vector_alloc(8); | |||||
for (size_t i = 0; i < v->size; i++) { | |||||
gsl_vector_set(v, i, i/8.); | |||||
} | |||||
complex double adotb = dot(a, b); | |||||
printf("(1+2i) dot (2+i): %g+%gi\n", creal(adotb), cimag(adotb)); | |||||
printf("v dot 2:\n"); | |||||
double d = 2; | |||||
gsl_vector_complex_print(dot(v, d)); | |||||
printf("v dot (1+2i)\n"); | |||||
gsl_vector_complex *vc = dot(v, a); | |||||
gsl_vector_complex_print(vc); | |||||
printf("v dot (1+2i) again\n"); | |||||
gsl_vector_complex_print(dot(v, c)); | |||||
} |