#include "BSprivate.h"

/*@ BSpar_solve - General solver of a system of equations preconditioned 
                  by one of several preconditioners and using one of several
                  possible methods.  The rhs can be a block of vectors.

    Input Parameters:
.   A - a sparse matrix
.   fact_A - the incomplete factored version of A, if any
.   comm_A - the communication structure for A
.   rhs - the contiguous block of vectors forming the rhs
.   residual - the final computed residual

    Output Parameters:
.   x - the contiguous block of vectors containing the solution (can
	contain an initial guess if BSctx_set_guess() is set.

Now specified in procinfo context:
.   pre_option - the preconditioner to use:
$           PRE_ICC: incomplete Cholesky factorization
$           PRE_ILU: incomplete LU factorization
$           PRE_SSOR: Successive over relaxation
$           PRE_BJACOBI: Block Jacobi
.   err_tol - the tolerance to which to solve the problem
              stop if the estimated norm of the residual divided by
              the norm of the rhs is less than err_tol
.   max_iter - the maximum number of iterations to take
.   guess - if TRUE, then initialize out_x to 0, otherwise the program
            assumes that out_x contains an initial guess
.   procinfo - the usual processor stuff

    Returns:
    The number of iterations or a negative number indicating the number
    of iterations prior to finding that the matrix (or preconditioner) 
    is not positive definite.

    Notes:
    The preconditioners must be computed prior to calling BSpar_solve.
    For more information on the preconditioners, see the manual or
	BSctx_set_pre().

 @*/
int BSpar_solve(BSpar_mat *A, BSpar_mat *fact_A, BScomm *comm_A,
	FLOAT *rhs, FLOAT *x, FLOAT *residual, BSprocinfo *procinfo)
{
	int num_iter;

	/* Do some error checking */
	if((A->icc_storage)&&(!procinfo->single)) {
		if(comm_A->num_rhs!=procinfo->num_rhs) {
			MY_SETERRCN(PAR_SOLVE_ERROR,"num_rhs in context does not agree num_rhs in comm\n");
		}
	}
	if((procinfo->scaling)&&(A->scale_diag==NULL)) {
		MY_SETERRCN(SCALING_ERROR,"Scaling set in context however matrix not scaled\n");
	}

	if ((procinfo->print) && (PSISROOT(procinfo))) {
		BSctx_print(procinfo);
	}

	if(procinfo->method==GMRES) {
		num_iter = BSpar_gmres(procinfo->num_rhs,A,fact_A,comm_A,rhs,x,
			procinfo->preconditioner,procinfo->restart,
			procinfo->residual_tol,procinfo->max_iterations,
			residual,procinfo->guess,procinfo);
	} else {
		num_iter = BSpar_sym_solve(procinfo->num_rhs,A,fact_A,
			comm_A,rhs,x,procinfo->preconditioner,
			procinfo->residual_tol,procinfo->max_iterations,
			residual,procinfo->guess,procinfo);
	}
	return(num_iter);
}