/*
bicgstab.c
*/

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

void bicgstab(int ifreq, int64_t n, d_complex_t *a, d_complex_t *x, d_complex_t *b, d_complex_t *r, d_complex_t *p, d_complex_t *q, d_complex_t *t, d_complex_t *u, d_complex_t *r0, int *niter, double *residual, double *eaverage, FILE *fp, int iout)
{
	const char msg1[] = "  --- converged ---";
	const char msg2[] = "  --- diverged ---";
	const char msg3[] = "  --- not converged ---";
	const char msg4[] = "  --- devided by zero ---";
	const double inf = 1e14;
	char str[BUFSIZ];

	// RHS norm
	double bnrm = Dznrm2(n, b);

	// == initial value ==

	// q = Ax
	copyV(x);
	prodmv(0, ifreq, a, x, q, NULL);

	// r = b - Ax = b - q
	Zcopy(n, b, r);
	Zaxpy(n, d_complex(-1, 0), q, r);

	// r0* <- b
	r0 = b;

	// r0*
	if (!Iseed) {
		// r0* = r
		Zcopy(n, r, r0);
	}
	else {
		// r0* = p (work array)
		Zcopy(n, p, r0);
	}

	// p = r
	Zcopy(n, r, p);

	// dot1 = (r0*, r)
	d_complex_t dot1 = Zdotc(n, r0, r);

	// == iteration ==

	int iter;
	for (iter = 0; iter <= Solver.maxiter; iter++) {
		// update

		// q = Ap
		copyV(p);
		prodmv(0, ifreq, a, p, q, NULL);

		// a = (r0*, r) / (r0*, Ap) = dot1 / (r0*, q)
		d_complex_t dot2 = Zdotc(n, r0, q);
		if (d_abs(dot2) == 0) {
			if (iout) {
				sprintf(str, "%s (1) (%d %e %e)", msg4, iter, d_abs(dot1), d_abs(dot2));
				monitor1(fp, str);
			}
			break;
		}
		d_complex_t acg = d_div(dot1, dot2);

		// t = r - aAp = r - aq
		Zaxpy(n, d_neg(acg), q, r);
		Zcopy(n, r, t);

		// u = At
		copyV(t);
		prodmv(0, ifreq, a, t, u, NULL);

		// z = (At, t) / (At, At) = (u, t) / (u, u)
		d_complex_t dot3 = Zdotc(n, u, t);
		d_complex_t dot4 = Zdotc(n, u, u);
		if (d_abs(dot4) == 0) {
			if (iout) {
				sprintf(str, "%s (2) (%d %e %e)", msg4, iter, d_abs(dot3), d_abs(dot4));
				monitor1(fp, str);
			}
			break;
		}
		d_complex_t zcg = d_div(dot3, dot4);

		// x = x + ap + zt
/*
		Zaxpy(n, acg, p, x);
		Zaxpy(n, zcg, t, x);
*/
		Zaxbycz(n, d_complex(1, 0), acg, zcg, x, p, t);

		// r = t - zAt = t - zu
		Zaxpy(n, d_neg(zcg), u, t);
		Zcopy(n, t, r);

		// residual

		double xnrm = Dznrm2(n, x);
		eaverage[iter] = xnrm;

		double rnrm = Dznrm2(n, r);
		residual[iter] = rnrm / bnrm;

		// monitor

		if ((iter % Solver.nout) == 0 || (residual[iter] < Solver.converg)) {
			if (iout) {
				monitor6(fp, iter, residual[iter], eaverage[iter]);
			}
		}

		// converged / diverged / not converged

		if ((iter > 0) && (residual[iter] < Solver.converg)) {
			if (iout) {
				monitor1(fp, msg1);
			}
			break;
		}
		if (residual[iter] > inf) {
			if (iout) {
				monitor1(fp, msg2);
			}
			break;
		}
		if (iter == Solver.maxiter) {
			if (iout) {
				monitor1(fp, msg3);
			}
			break;
		}

		// update

		// b = (a / z) (r0*, r[k+1]) / (r0*, r[k]) = (a / z) (r0*, r) / dot1
		d_complex_t dot5 = Zdotc(n, r0, r);
		if (d_abs(d_mul(zcg, dot1)) == 0) {
			if (iout) {
				sprintf(str, "%s (3) (%d %.3e %.3e %.3e %.3e)", msg4, iter, d_abs(acg), d_abs(zcg), d_abs(dot1), d_abs(dot5));
				monitor1(fp, str);
			}
			break;
		}
		d_complex_t bcg = d_div(d_mul(acg, dot5), d_mul(zcg, dot1));
		//d_complex_t bcg = d_mul(d_div(acg, dot1), d_div(dot5, zcg));
		dot1 = dot5;

		// p = r + b(p - zAp) = r + b(p - zq)
/*
		Zaxpy(n, d_neg(zcg), q, p);
		Zcopy(n, r, q);
		Zaxpy(n, bcg, p, q);
		Zcopy(n, q, p);
*/
		Zaxbycz(n, bcg, d_complex(1, 0), d_neg(d_mul(bcg, zcg)), p, r, q);
	}

	// iterations
	niter[ifreq] = iter + 1;
}
