#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 *i_ptr;
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 *) MALLOC(sizeof(int));
i_ptr = (int *) A->map->vlocal2global; /* pointer to mapping data */
*(i_ptr) = grid->offset;
A->map->flocal2global = BSloc2glob; /* the mapping function */
A->map->free_l2g = BSfreel2g; /* 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 *) MALLOC(sizeof(int));
i_ptr = (int *) A->map->vglobal2local; /* pointer to mapping data */
*(i_ptr) = grid->offset;
A->map->fglobal2local = BSglob2loc; /* the mapping function */
A->map->free_g2l = BSfreeg2l; /* 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 *)
BSmake_off_map(grid->offset,procinfo,A->global_num_rows); CHKERRN(0);
A->map->fglobal2proc = BSglob2proc; /* the mapping function */
A->map->free_g2p = BSfree_off_map; /* 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 */
row->col[row->length] = points[i+1][j][k].num;
row->nz[row->length] = -(1.0/8.0);
(row->length)++;
/* west */
row->col[row->length] = points[i-1][j][k].num;
row->nz[row->length] = -(1.0/8.0);
(row->length)++;
/* north */
row->col[row->length] = points[i][j+1][k].num;
row->nz[row->length] = -(1.0/8.0);
(row->length)++;
/* south */
row->col[row->length] = points[i][j-1][k].num;
row->nz[row->length] = -(1.0/8.0);
(row->length)++;
/* up */
row->col[row->length] = points[i][j][k+1].num;
row->nz[row->length] = -(1.0/8.0);
(row->length)++;
/* down */
row->col[row->length] = points[i][j][k-1].num;
row->nz[row->length] = -(1.0/8.0);
(row->length)++;
/* east, north, up */
if ((ptype == 27) || (points[i+1][j+1][k+1].type == 27)) {
row->col[row->length] = points[i+1][j+1][k+1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* east, south, up */
if ((ptype == 27) || (points[i+1][j-1][k+1].type == 27)) {
row->col[row->length] = points[i+1][j-1][k+1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* east, north, even */
if ((ptype == 27) || (points[i+1][j+1][k].type == 27)) {
row->col[row->length] = points[i+1][j+1][k].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* east, south, even */
if ((ptype == 27) || (points[i+1][j-1][k].type == 27)) {
row->col[row->length] = points[i+1][j-1][k].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* east, north, down */
if ((ptype == 27) || (points[i+1][j+1][k-1].type == 27)) {
row->col[row->length] = points[i+1][j+1][k-1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* east, south, down */
if ((ptype == 27) || (points[i+1][j-1][k-1].type == 27)) {
row->col[row->length] = points[i+1][j-1][k-1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* east, even, up */
if ((ptype == 27) || (points[i+1][j][k+1].type == 27)) {
row->col[row->length] = points[i+1][j][k+1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* east, even, down */
if ((ptype == 27) || (points[i+1][j][k-1].type == 27)) {
row->col[row->length] = points[i+1][j][k-1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* even, north, up */
if ((ptype == 27) || (points[i][j+1][k+1].type == 27)) {
row->col[row->length] = points[i][j+1][k+1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* even, south, up */
if ((ptype == 27) || (points[i][j-1][k+1].type == 27)) {
row->col[row->length] = points[i][j-1][k+1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* even, north, down */
if ((ptype == 27) || (points[i][j+1][k-1].type == 27)) {
row->col[row->length] = points[i][j+1][k-1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* even, south, down */
if ((ptype == 27) || (points[i][j-1][k-1].type == 27)) {
row->col[row->length] = points[i][j-1][k-1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* west, north, up */
if ((ptype == 27) || (points[i-1][j+1][k+1].type == 27)) {
row->col[row->length] = points[i-1][j+1][k+1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* west, south, up */
if ((ptype == 27) || (points[i-1][j-1][k+1].type == 27)) {
row->col[row->length] = points[i-1][j-1][k+1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* west, north, even */
if ((ptype == 27) || (points[i-1][j+1][k].type == 27)) {
row->col[row->length] = points[i-1][j+1][k].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* west, south, even */
if ((ptype == 27) || (points[i-1][j-1][k].type == 27)) {
row->col[row->length] = points[i-1][j-1][k].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* west, north, down */
if ((ptype == 27) || (points[i-1][j+1][k-1].type == 27)) {
row->col[row->length] = points[i-1][j+1][k-1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* west, south, down */
if ((ptype == 27) || (points[i-1][j-1][k-1].type == 27)) {
row->col[row->length] = points[i-1][j-1][k-1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* west, even, up */
if ((ptype == 27) || (points[i-1][j][k+1].type == 27)) {
row->col[row->length] = points[i-1][j][k+1].num;
row->nz[row->length] = -(1.0/80.0);
(row->length)++;
}
/* west, even, down */
if ((ptype == 27) || (points[i-1][j][k-1].type == 27)) {
row->col[row->length] = points[i-1][j][k-1].num;
row->nz[row->length] = -(1.0/80.0);
(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);
}
syntax highlighted by Code2HTML, v. 0.9.1