#include "pargrid.h"

/*+ get_mat - Generate a sparse matrix from the given grid

     Input Parameters:
     grid - the given grid
     procinfo - the processor information (in BlockSolve format)

     Returns:
     the sparse matrix

 +*/

BSspmat *get_mat(par_grid *grid, BSprocinfo *procinfo)
{
	BSspmat *A;
	int	i, j, k, l;
	int	ptype, pnum;
	int	count;
	BSsprow *row;
	point	***points = grid->points;

	grid->rp = NULL;
	grid->cval = NULL;
	grid->aval = NULL;

	/* ****************************************************** */
	/* set up the structure and mapping for the sparse matrix */
	/* ****************************************************** */
	/* allocate the pointer to A */
	A = (BSspmat *) MALLOC(sizeof(BSspmat));

	/* set the number of local rows */
	A->num_rows = grid->local_total;

	/* set the number of global rows */
	A->global_num_rows = grid->global_total;

	/* allocate the array of rows, and the space in each row */
	/* allow for a max length in each row, but set the current length to 0 */
	A->rows = (BSsprow **) MALLOC(sizeof(BSsprow *)*A->num_rows);
	for (i=0;i<A->num_rows;i++) {
		A->rows[i] = (BSsprow *) MALLOC(sizeof(BSsprow));
		A->rows[i]->length = 0;
		A->rows[i]->col = (int *) MALLOC(sizeof(int)*MAX_LEN);
		A->rows[i]->nz = (FLOAT *) MALLOC(sizeof(FLOAT)*MAX_LEN);
	}

	/* allocate a pointer to a mapping structure */
	A->map = (BSmapping *) MALLOC(sizeof(BSmapping));

	/* set up the local to global mapping */
	/* all we need for this is the beginning number (offset) of */
	/* the local rows in the global numbering (see BSloc2glob) */
	A->map->vlocal2global = (void *) grid;
	A->map->flocal2global = ex3_loc2glob; /* the mapping function */
	A->map->free_l2g = ex3_freel2g; /* the routine to free the mapping */
		
	/* set up the global to local mapping */
	/* all we need for this is the beginning number (offset) of */
	/* the local rows in the global numbering (see BSglob2loc) */
	A->map->vglobal2local = (void *) grid;
	A->map->fglobal2local = ex3_glob2loc; /* the mapping function */
	A->map->free_g2l = ex3_freeg2l;	/* the routine to free the mapping */
		
	/* set up the global to processor number mapping */
	/* we call the routine BSmake_off_map to create the mapping data */
	/* the local rows in the global numbering (see BSglob2proc) */
	A->map->vglobal2proc = (void *) grid; 
	A->map->fglobal2proc = ex3_glob2proc;	/* the mapping function */
	A->map->free_g2p = ex3_freeg2l;	/* the routine to free the mapping */

	/* ****************************************************** */
	/* now put values in the matrix                           */
	/* ****************************************************** */
	count = 0;
	for (i=1;i<grid->l_num_x+1;i++) {
		for (j=1;j<grid->l_num_y+1;j++) {
			for (k=1;k<grid->l_num_z+1;k++) {
				/* get a pointer to the correct row */
				row = A->rows[count];
				pnum = points[i][j][k].num;
				ptype = points[i][j][k].type;

				/* diagonal, set the column number, the nonzero values and */
				/* increment the length of the row */
				row->col[row->length] = pnum;
				row->nz[row->length] = 1.0;
				(row->length)++;

				/* east */
				if (points[i+1][j][k].num >= 0) {
					if (ptype == 5) {
						row->nz[row->length] = -(1.0/4.0);
					} else {
						row->nz[row->length] = -(1.0/5.0);
					}
					row->col[row->length] = points[i+1][j][k].num;
					(row->length)++;
				}

				/* west */
				if (points[i-1][j][k].num >= 0) {
					if (ptype == 5) {
						row->nz[row->length] = -(1.0/4.0);
					} else {
						row->nz[row->length] = -(1.0/5.0);
					}
					row->col[row->length] = points[i-1][j][k].num;
					(row->length)++;
				}

				/* north */
				if (points[i][j+1][k].num >= 0) {
					if (ptype == 5) {
						row->nz[row->length] = -(1.0/4.0);
					} else {
						row->nz[row->length] = -(1.0/5.0);
					}
					row->col[row->length] = points[i][j+1][k].num;
					(row->length)++;
				}

				/* south */
				if (points[i][j-1][k].num >= 0) {
					if (ptype == 5) {
						row->nz[row->length] = -(1.0/4.0);
					} else {
						row->nz[row->length] = -(1.0/5.0);
					}
					row->col[row->length] = points[i][j-1][k].num;
					(row->length)++;
				}

				/* east, north, even */
				if (points[i+1][j+1][k].num >= 0) {
					if (ptype == 9) {
						row->nz[row->length] = -(1.0/20.0);
						row->col[row->length] = points[i+1][j+1][k].num;
						(row->length)++;
					}
				}

				/* east, south, even */
				if (points[i+1][j-1][k].num >= 0) {
					if (ptype == 9) {
						row->nz[row->length] = -(1.0/20.0);
						row->col[row->length] = points[i+1][j-1][k].num;
						(row->length)++;
					}
				}

				/* west, north, even */
				if (points[i-1][j+1][k].num >= 0) {
					if (ptype == 9) {
						row->nz[row->length] = -(1.0/20.0);
						row->col[row->length] = points[i-1][j+1][k].num;
						(row->length)++;
					}
				}

				/* west, south, even */
				if (points[i-1][j-1][k].num >= 0) {
					if (ptype == 9) {
						row->nz[row->length] = -(1.0/20.0);
						row->col[row->length] = points[i-1][j-1][k].num;
						(row->length)++;
					}
				}

				/* sort the column values in the row */
				BSheap_sort1d(row->length,row->col,row->nz); CHKERRN(0);
				/* now find the index of the diagonal in the row */
				row->diag_ind = -1;
				for (l=0;l<row->length;l++) {
					if (row->col[l] == pnum) {
						row->diag_ind = l;
						break;
					}
				}
				if (row->diag_ind < 0) {
					printf("Diagonal problem\n");
				}
				count++;
			}
		}
	}
	return(A);
}