/*
solve.c
*/

#include "oth.h"
#include "oth_prototype.h"

void solve(FILE *fp, int iout)
{
	const size_t size = NE * sizeof(d_complex_t);

	// alloc
	memalloc();
	memalloc_cg();

	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		// monitor
		if (iout) {
			monitor5(fp, ifreq);
		}

		// setup
		setup2(ifreq);

		// setup RHS vector
		memset(Bcg, 0, size);
		rhs(Bcg, ifreq);

		// (setup matrix and) scaling RHS vector
		if (MATRIX) memset(Acg, 0, NNZ * size);
		prodmv(1, ifreq, Acg, NULL, NULL, Bcg);

		// clear
		memset(Rcg, 0, size);
		memset(Pcg, 0, size);
		memset(Qcg, 0, size);
		memset(Tcg, 0, size);
		memset(Ucg, 0, size);

		// initial value
		if (!Iseed) {
			memcpy(Ecg, Bcg, size);
		}
		else {
			memset(Ecg, 0, size);

			// r0* : random vector
			//srand((unsigned int) time(NULL));
			srand(Iseed);
			for (int64_t n = 0; n < NE; n++) {
				Pcg[n].r = rand() / (double)RAND_MAX;
				Pcg[n].i = rand() / (double)RAND_MAX;
			}
		}

		// solve
		bicgstab(ifreq, NE, Acg, Ecg, Bcg, Rcg, Pcg, Qcg, Tcg, Ucg, R0cg, NIter, Residual[ifreq], Eaverage[ifreq], fp, iout);

		// copy E to Ex+Ey+Ez
		copyE(ifreq);
	}

	// free
	memfree_cg();
}
