/*
blas1c.c

BLAS functions, level-1, double complex
*/

#include <stdint.h>
#include <math.h>

typedef struct {double r, i;} d_complex_t;

// [y] = [x]
void Zcopy(int64_t n, const d_complex_t *x, d_complex_t *y)
{
	int64_t i;

#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (i = 0; i < n; i++) {
		y[i].r = x[i].r;
		y[i].i = x[i].i;
	}
}
/*
// [x] * [y]
d_complex_t Zdotu(int64_t n, const d_complex_t *x, const d_complex_t *y)
{
	d_complex_t ret;
	double sum_r = 0, sum_i = 0;
	int64_t i;
#ifdef _OPENMP
#pragma omp parallel for reduction (+:sum_r, sum_i)
#endif
	for (i = 0; i < n; i++) {
		sum_r += (x[i].r * y[i].r) - (x[i].i * y[i].i);
		sum_i += (x[i].r * y[i].i) + (x[i].i * y[i].r);
	}
	ret.r = sum_r;
	ret.i = sum_i;

	return ret;
}
*/
// conj[x] * [y]
d_complex_t Zdotc(int64_t n, const d_complex_t *x, const d_complex_t *y)
{
	d_complex_t ret;
	double sum_r = 0, sum_i = 0;
	int64_t i;
#ifdef _OPENMP
#pragma omp parallel for reduction (+:sum_r, sum_i)
#endif
	for (i = 0; i < n; i++) {
		sum_r += (x[i].r * y[i].r) + (x[i].i * y[i].i);
		sum_i += (x[i].r * y[i].i) - (x[i].i * y[i].r);
	}
	ret.r = sum_r;
	ret.i = sum_i;

	return ret;
}

// y[i] = a * x[i] + y[i]
void Zaxpy(int64_t n, d_complex_t alpha, const d_complex_t *x, d_complex_t *y)
{
	int64_t i;
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (i = 0; i < n; i++) {
		y[i].r += (alpha.r * x[i].r) - (alpha.i * x[i].i);
		y[i].i += (alpha.r * x[i].i) + (alpha.i * x[i].r);
	}
}

// |x[i]|
double Dznrm2(int64_t n, const d_complex_t *x)
{
	double sum = 0;
	int64_t i;
#ifdef _OPENMP
#pragma omp parallel for reduction (+:sum)
#endif
	for (i = 0; i < n; i++) {
		sum += (x[i].r * x[i].r) + (x[i].i * x[i].i);
	}

	return sqrt(sum);
}

// x[i] = a * x[i] + b * y[i] + c * z[i] (custom)
void Zaxbycz(int64_t n, d_complex_t a, d_complex_t b, d_complex_t c, d_complex_t *x, const d_complex_t *y, const d_complex_t *z)
{
	int64_t i;
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (i = 0; i < n; i++) {
		double x_r = ((a.r * x[i].r) - (a.i * x[i].i)) + ((b.r * y[i].r) - (b.i * y[i].i)) + ((c.r * z[i].r) - (c.i * z[i].i));
		double x_i = ((a.r * x[i].i) + (a.i * x[i].r)) + ((b.r * y[i].i) + (b.i * y[i].r)) + ((c.r * z[i].i) + (c.i * z[i].r));
		x[i].r = x_r;
		x[i].i = x_i;
	}
}
