#include "BSprivate.h"

/*@ BStri_solve - Multiply the matrix A(-1) by a block of vectors

    Input Parameters:
.   A - a sparse matrix
.   fact_A - a possibly incompletely factored version of A
.   comm_A - the communication structure for A (fact_A)
.   vec - the block of vectors to multiply by
.   pre_option - the preconditioner to select
                 PRE_ICC: incomplete Cholesky factorization
                 PRE_ILU: incomplete LU factorization
                 PRE_SSOR: Successive over relaxation
                 PRE_BJACOBI: Block Jacobi
.   BS - the number of vectors in vec
.   procinfo - the usual processor stuff

    Output Parameters:
.   vec - the resulting block of vectors

    Returns:
    void

    Notes:
    Different code is used to multiply a single vector than is used
    to multiply a block of vectors (this improves efficiency).

 @*/
void BStri_solve(BSpar_mat *A, BSpar_mat *fact_A, BScomm *comm_A, FLOAT *vec,
		int pre_option, int BS, BSprocinfo *procinfo)
{
	int	i;

	if (pre_option == PRE_DIAG) {
		return;
	} else if ((pre_option==PRE_ICC)||(pre_option==PRE_ILU)) {
		if((pre_option==PRE_ICC)&&(!fact_A->icc_storage)) {
			MY_SETERRC(TRI_SOLVE_ERROR,"ICC preconditioner requested with ILU storage\n");
		} else if ((pre_option==PRE_ILU)&&(fact_A->icc_storage)) {
			MY_SETERRC(TRI_SOLVE_ERROR,"ILU preconditioner requested with ICC storage\n");
		}
		if (BS == 1) {
			MLOG_ELM(procinfo->procset);
			if (procinfo->single) {
				BSfor_solve1(fact_A,vec,comm_A,procinfo); CHKERR(0);
			} else {
				BSfor_solve(fact_A,vec,comm_A,procinfo); CHKERR(0);
			}
			MLOG_ACC(MS_FORWARD);
			MLOG_ELM(procinfo->procset);
			if (procinfo->single) {
				BSback_solve1(fact_A,vec,comm_A,procinfo); CHKERR(0);
			} else {
				BSback_solve(fact_A,vec,comm_A,procinfo); CHKERR(0);
			}
			MLOG_ACC(MS_BACKWARD);
		} else {
			MLOG_ELM(procinfo->procset);
			BSb_for_solve(fact_A,vec,comm_A,BS,procinfo); CHKERR(0);
			MLOG_ACC(BMS_FORWARD);
			MLOG_ELM(procinfo->procset);
			BSb_back_solve(fact_A,vec,comm_A,BS,procinfo); CHKERR(0);
			MLOG_ACC(BMS_BACKWARD);
		}
	} else if (pre_option == PRE_SSOR) {
		if (BS == 1) {
			MLOG_ELM(procinfo->procset);
			if (procinfo->single) {
				BSfor_solve1(A,vec,comm_A,procinfo); CHKERR(0);
			} else {
				BSfor_solve(A,vec,comm_A,procinfo); CHKERR(0);
			}
			MLOG_ACC(MS_FORWARD);
			MLOG_ELM(procinfo->procset);
			if (procinfo->single) {
				BSback_solve1(A,vec,comm_A,procinfo); CHKERR(0);
			} else {
				BSback_solve(A,vec,comm_A,procinfo); CHKERR(0);
			}
			MLOG_ACC(MS_BACKWARD);
		} else {
			MLOG_ELM(procinfo->procset);
			BSb_for_solve(A,vec,comm_A,BS,procinfo); CHKERR(0);
			MLOG_ACC(BMS_FORWARD);
			MLOG_ELM(procinfo->procset);
			BSb_back_solve(A,vec,comm_A,BS,procinfo); CHKERR(0);
			MLOG_ACC(BMS_BACKWARD);
		}
	} else if (pre_option == PRE_BJACOBI) {
		if (BS == 1) {
			MLOG_ELM(procinfo->procset);
			BSbjacobi(A,vec,comm_A,procinfo); CHKERR(0);
			MLOG_ACC(BJ_FORWARD);
		} else {
			for (i=0;i<BS;i++) {
				MLOG_ELM(procinfo->procset);
				BSbjacobi(A,&(vec[A->num_rows*i]),comm_A,procinfo); CHKERR(0);
				MLOG_ACC(BJ_FORWARD);
			}
		}
	}
}