#include "BSprivate.h"

/*@ BSfor_solve1 - Forward triangular matrix solution on a 
                  single vector

    Input Parameters:
.   A - The sparse matrix
.   x - The rhs
.   comm - The communication structure for A
.   procinfo - the usual processor information

    Output Parameters:
.   x - on exit contains the solution vector

    Returns:
    void

	Notes:
    We assume that A has no i-nodes or cliques

 @*/
void BSfor_solve1(BSpar_mat *A, FLOAT *x, BScomm *comm, BSprocinfo *procinfo)
{
	BMphase *to_phase, *from_phase;
	BMmsg *msg;
	int	i, j, k;
	int	cl_ind, in_ind, symmetric;
	int	count, size, ind;
	int *row;
	FLOAT *nz;
	BScl_2_inode *clique2inode;
	BSnumbering *color2clique;
	BSinode *inodes;
	int	*data_ptr, msg_len;
	FLOAT *msg_buf;
	FLOAT	t;
	int *gnum, *iperm;

	/* Is the symmetric data structure used? */
	symmetric = A->icc_storage;

	color2clique = A->color2clique;
	clique2inode = A->clique2inode;
	inodes = A->inodes->list;
	gnum = A->global_row_num->numbers;
	iperm = A->inv_perm->perm;

	/* post for all messages */
	BMinit_comp_msg(comm->from_msg,procinfo); CHKERR(0);

	/* now do this phase by phase */
	for (i=0;i<color2clique->length-1;i++) {
		if(symmetric) {
			/* find my portion of the solution using the cliques on the diagonal */
			for (cl_ind=color2clique->numbers[i];
				cl_ind<color2clique->numbers[i+1];cl_ind++) {
				if (procinfo->my_id == clique2inode->proc[cl_ind]) {
					/* first, multiply the clique */
					/* the clique is stored, inverted, in the upper triangle */
					ind = clique2inode->d_mats[cl_ind].local_ind;
					x[ind] *= *(clique2inode->d_mats[cl_ind].matrix);
				}
			}
		}

		/* now send my messages */
		to_phase = BMget_phase(comm->to_msg,i); CHKERR(0);
		msg = NULL;
		while ((msg = BMnext_msg(to_phase,msg)) != NULL) {
			CHKERR(0);
			msg_buf = (FLOAT *) BMget_msg_ptr(msg); CHKERR(0);
			data_ptr = BMget_user(msg,&msg_len); CHKERR(0);
			for (j=0;j<msg_len;j++) {
				msg_buf[j] = x[data_ptr[j]];
			}
			BMsendf_msg(msg,procinfo); CHKERR(0);
		}
		CHKERR(0);

		/* do some local work */
		for (cl_ind=color2clique->numbers[i];
			cl_ind<color2clique->numbers[i+1];cl_ind++) {
			if (procinfo->my_id == clique2inode->proc[cl_ind]) {
				/* multiply the inodes */
				in_ind=clique2inode->inode_index[cl_ind];
				size = inodes[in_ind].length;
				if (size > 0) {
					ind = clique2inode->d_mats[cl_ind].local_ind;
					row = inodes[in_ind].row_num;
					nz = inodes[in_ind].nz;
					t = x[ind];
					if(symmetric) {
						for (k=0;k<size;k++) x[row[k]] -= t*nz[k];
					} else {
						for (k=0;k<size;k++) {
							if (gnum[iperm[row[k]]] > inodes[in_ind].gcol_num)
								x[row[k]] -= nz[k]*t;
						}
					}
				}
			}
		}

		/* receive my messages and do non-local work */
		from_phase = BMget_phase(comm->from_msg,i); CHKERR(0);
		while ((msg = BMrecv_msg(from_phase)) != NULL) {
			CHKERR(0);
			msg_buf = (FLOAT *) BMget_msg_ptr(msg); CHKERR(0);
			data_ptr = BMget_user(msg,&msg_len); CHKERR(0);
			count = 0;
			for (cl_ind=data_ptr[0];cl_ind<=data_ptr[1];cl_ind++) {
				in_ind=clique2inode->inode_index[cl_ind];
				size = inodes[in_ind].length;
				if (size > 0) {
					row = inodes[in_ind].row_num;
					nz = inodes[in_ind].nz;
					t = msg_buf[count];
					if(symmetric) {
						for (k=0;k<size;k++) x[row[k]] -= t*nz[k];
					} else {
						for (k=0;k<size;k++) {
							if (gnum[iperm[row[k]]] > inodes[in_ind].gcol_num)
								x[row[k]] -= nz[k]*t;
						}
					}
				}
				count++;
			}
			BMfree_msg(msg); CHKERR(0);
		}
		CHKERR(0);
	}
	/* wait for all of the sent messages to finish */
	BMfinish_comp_msg(comm->to_msg,procinfo); CHKERR(0);
	MLOG_flop((2*A->local_nnz));
}