/*
memory.c
alloc, set, copy, free
*/

#include "oth.h"
#include "complex.h"

void memalloc(void)
{
	// mesh factor

	KRXn = (double *)malloc((Nx + 1) * sizeof(double));
	KRYn = (double *)malloc((Ny + 1) * sizeof(double));
	KRZn = (double *)malloc((Nz + 1) * sizeof(double));
	KRXc = (double *)malloc((Nx + 0) * sizeof(double));
	KRYc = (double *)malloc((Ny + 0) * sizeof(double));
	KRZc = (double *)malloc((Nz + 0) * sizeof(double));

	KDXn = (double *)malloc((Nx + 1) * sizeof(double));
	KDYn = (double *)malloc((Ny + 1) * sizeof(double));
	KDZn = (double *)malloc((Nz + 1) * sizeof(double));
	KDXc = (double *)malloc((Nx + 0) * sizeof(double));
	KDYc = (double *)malloc((Ny + 0) * sizeof(double));
	KDZc = (double *)malloc((Nz + 0) * sizeof(double));

	KRXp = (double *)malloc((Nx + 1) * sizeof(double));
	KRYp = (double *)malloc((Ny + 1) * sizeof(double));
	KRZp = (double *)malloc((Nz + 1) * sizeof(double));
	KRXm = (double *)malloc((Nx + 1) * sizeof(double));
	KRYm = (double *)malloc((Ny + 1) * sizeof(double));
	KRZm = (double *)malloc((Nz + 1) * sizeof(double));

	// iteration

	Residual = (double **)malloc(NFrequency * sizeof(double *));
	Eaverage = (double **)malloc(NFrequency * sizeof(double *));
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		size_t size = (Solver.maxiter + 1) * sizeof(double);
		Residual[ifreq] = (double *)malloc(size);
		Eaverage[ifreq] = (double *)malloc(size);
		memset(Residual[ifreq], 0, size);
		memset(Eaverage[ifreq], 0, size);
	}

	NIter = (int *)malloc(NFrequency * sizeof(int));
	memset(NIter, 0, NFrequency * sizeof(int));

	// Ex, Ey, Ez

	Ex = (f_complex_t ****)malloc(NFrequency * sizeof(f_complex_t ***));
	Ey = (f_complex_t ****)malloc(NFrequency * sizeof(f_complex_t ***));
	Ez = (f_complex_t ****)malloc(NFrequency * sizeof(f_complex_t ***));

	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		// Ex
		Ex[ifreq] = (f_complex_t ***)malloc((Nx + 0) * sizeof(f_complex_t **));
		for (int i = 0; i < Nx + 0; i++) {
			Ex[ifreq][i] = (f_complex_t **)malloc((Ny + 1) * sizeof(f_complex_t *));
			for (int j = 0; j < Ny + 1; j++) {
				Ex[ifreq][i][j] = (f_complex_t *)malloc((Nz + 1) * sizeof(f_complex_t));
				memset(Ex[ifreq][i][j], 0, (Nz + 1) * sizeof(f_complex_t));
			}
		}
		// Ey
		Ey[ifreq] = (f_complex_t ***)malloc((Nx + 1) * sizeof(f_complex_t **));
		for (int i = 0; i < Nx + 1; i++) {
			Ey[ifreq][i] = (f_complex_t **)malloc((Ny + 0) * sizeof(f_complex_t *));
			for (int j = 0; j < Ny + 0; j++) {
				Ey[ifreq][i][j] = (f_complex_t *)malloc((Nz + 1) * sizeof(f_complex_t));
				memset(Ey[ifreq][i][j], 0, (Nz + 1) * sizeof(f_complex_t));
			}
		}
		// Ez
		Ez[ifreq] = (f_complex_t ***)malloc((Nx + 1) * sizeof(f_complex_t **));
		for (int i = 0; i < Nx + 1; i++) {
			Ez[ifreq][i] = (f_complex_t **)malloc((Ny + 1) * sizeof(f_complex_t *));
			for (int j = 0; j < Ny + 1; j++) {
				Ez[ifreq][i][j] = (f_complex_t *)malloc((Nz + 0) * sizeof(f_complex_t));
				memset(Ez[ifreq][i][j], 0, (Nz + 0) * sizeof(f_complex_t));
			}
		}
	}
}

// copy E to Ex+Ey+Ez
void copyEtoE(int ifreq)
{
	// Ex
	for (int i = 0; i < Nx + 0; i++) {
	for (int j = 0; j < Ny + 1; j++) {
	for (int k = 0; k < Nz + 1; k++) {
		Ex[ifreq][i][j][k] = f_cast(Ecg[NEX(i, j, k)]);
	}
	}
	}

	// Ey
	for (int i = 0; i < Nx + 1; i++) {
	for (int j = 0; j < Ny + 0; j++) {
	for (int k = 0; k < Nz + 1; k++) {
		Ey[ifreq][i][j][k] = f_cast(Ecg[NEY(i, j, k)]);
	}
	}
	}

	// Ez
	for (int i = 0; i < Nx + 1; i++) {
	for (int j = 0; j < Ny + 1; j++) {
	for (int k = 0; k < Nz + 0; k++) {
		Ez[ifreq][i][j][k] = f_cast(Ecg[NEZ(i, j, k)]);
	}
	}
	}
}

void memalloc_cg(void)
{
	size_t size = NE * sizeof(d_complex_t);

	Ecg = (d_complex_t *)malloc(size);
	Bcg = (d_complex_t *)malloc(size);
	Rcg = (d_complex_t *)malloc(size);
	Pcg = (d_complex_t *)malloc(size);
	Qcg = (d_complex_t *)malloc(size);
	Tcg = (d_complex_t *)malloc(size);
	Ucg = (d_complex_t *)malloc(size);
	//R0cg = (d_complex_t *)malloc(size);

	if (MATRIX) {
		Acg = (d_complex_t *)malloc(NNZ * size);
	}
}

void memfree_cg(void)
{
	free(Ecg);
	free(Bcg);
	free(Rcg);
	free(Pcg);
	free(Qcg);
	free(Tcg);
	free(Ucg);

	if (MATRIX) {
		free(Acg);
	}
}
