| /*BHEADER********************************************************************** |
| * (c) 1997 The Regents of the University of California |
| * |
| * See the file COPYRIGHT_and_DISCLAIMER for a complete copyright |
| * notice, contact person, and disclaimer. |
| * |
| * $Revision$ |
| *********************************************************************EHEADER*/ |
| /****************************************************************************** |
| * |
| * Member functions for hypre_StructMatrix class. |
| * |
| *****************************************************************************/ |
| |
| #include "headers.h" |
| |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixExtractPointerByIndex |
| * Returns pointer to data for stencil entry coresponding to |
| * `index' in `matrix'. If the index does not exist in the matrix's |
| * stencil, the NULL pointer is returned. |
| *--------------------------------------------------------------------------*/ |
| |
| double * |
| hypre_StructMatrixExtractPointerByIndex( hypre_StructMatrix *matrix, |
| int b, |
| hypre_Index index ) |
| { |
| hypre_StructStencil *stencil; |
| int rank; |
| |
| stencil = hypre_StructMatrixStencil(matrix); |
| rank = hypre_StructStencilElementRank( stencil, index ); |
| |
| if ( rank >= 0 ) |
| return hypre_StructMatrixBoxData(matrix, b, rank); |
| else |
| return NULL; /* error - invalid index */ |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixCreate |
| *--------------------------------------------------------------------------*/ |
| |
| hypre_StructMatrix * |
| hypre_StructMatrixCreate( MPI_Comm comm, |
| hypre_StructGrid *grid, |
| hypre_StructStencil *user_stencil ) |
| { |
| hypre_StructMatrix *matrix; |
| |
| int i; |
| |
| matrix = hypre_CTAlloc(hypre_StructMatrix, 1); |
| |
| hypre_StructMatrixComm(matrix) = comm; |
| hypre_StructGridRef(grid, &hypre_StructMatrixGrid(matrix)); |
| hypre_StructMatrixUserStencil(matrix) = |
| hypre_StructStencilRef(user_stencil); |
| hypre_StructMatrixDataAlloced(matrix) = 1; |
| hypre_StructMatrixRefCount(matrix) = 1; |
| |
| /* set defaults */ |
| hypre_StructMatrixSymmetric(matrix) = 0; |
| for (i = 0; i < 6; i++) |
| hypre_StructMatrixNumGhost(matrix)[i] = 0; |
| |
| return matrix; |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixRef |
| *--------------------------------------------------------------------------*/ |
| |
| hypre_StructMatrix * |
| hypre_StructMatrixRef( hypre_StructMatrix *matrix ) |
| { |
| hypre_StructMatrixRefCount(matrix) ++; |
| |
| return matrix; |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixDestroy |
| *--------------------------------------------------------------------------*/ |
| |
| int |
| hypre_StructMatrixDestroy( hypre_StructMatrix *matrix ) |
| { |
| int i; |
| int ierr = 0; |
| |
| if (matrix) |
| { |
| hypre_StructMatrixRefCount(matrix) --; |
| if (hypre_StructMatrixRefCount(matrix) == 0) |
| { |
| if (hypre_StructMatrixDataAlloced(matrix)) |
| { |
| hypre_SharedTFree(hypre_StructMatrixData(matrix)); |
| } |
| hypre_CommPkgDestroy(hypre_StructMatrixCommPkg(matrix)); |
| |
| hypre_ForBoxI(i, hypre_StructMatrixDataSpace(matrix)) |
| hypre_TFree(hypre_StructMatrixDataIndices(matrix)[i]); |
| hypre_TFree(hypre_StructMatrixDataIndices(matrix)); |
| |
| hypre_BoxArrayDestroy(hypre_StructMatrixDataSpace(matrix)); |
| |
| hypre_TFree(hypre_StructMatrixSymmElements(matrix)); |
| hypre_StructStencilDestroy(hypre_StructMatrixUserStencil(matrix)); |
| hypre_StructStencilDestroy(hypre_StructMatrixStencil(matrix)); |
| hypre_StructGridDestroy(hypre_StructMatrixGrid(matrix)); |
| |
| hypre_TFree(matrix); |
| } |
| } |
| |
| return ierr; |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixInitializeShell |
| *--------------------------------------------------------------------------*/ |
| |
| int |
| hypre_StructMatrixInitializeShell( hypre_StructMatrix *matrix ) |
| { |
| int ierr = 0; |
| |
| hypre_StructGrid *grid; |
| |
| hypre_StructStencil *user_stencil; |
| hypre_StructStencil *stencil; |
| hypre_Index *stencil_shape; |
| int stencil_size; |
| int num_values; |
| int *symm_elements; |
| |
| int *num_ghost; |
| int extra_ghost[] = {0, 0, 0, 0, 0, 0}; |
| |
| hypre_BoxArray *data_space; |
| hypre_BoxArray *boxes; |
| hypre_Box *box; |
| hypre_Box *data_box; |
| |
| int **data_indices; |
| int data_size; |
| int data_box_volume; |
| |
| int i, j, d; |
| |
| grid = hypre_StructMatrixGrid(matrix); |
| |
| /*----------------------------------------------------------------------- |
| * Set up stencil and num_values: |
| * The stencil is a "symmetrized" version of the user's stencil |
| * as computed by hypre_StructStencilSymmetrize. |
| * |
| * The `symm_elements' array is used to determine what data is |
| * explicitely stored (symm_elements[i] < 0) and what data does is |
| * not explicitely stored (symm_elements[i] >= 0), but is instead |
| * stored as the transpose coefficient at a neighboring grid point. |
| *-----------------------------------------------------------------------*/ |
| |
| if (hypre_StructMatrixStencil(matrix) == NULL) |
| { |
| user_stencil = hypre_StructMatrixUserStencil(matrix); |
| |
| hypre_StructStencilSymmetrize(user_stencil, &stencil, &symm_elements); |
| |
| stencil_shape = hypre_StructStencilShape(stencil); |
| stencil_size = hypre_StructStencilSize(stencil); |
| |
| if (!hypre_StructMatrixSymmetric(matrix)) |
| { |
| /* store all element data */ |
| for (i = 0; i < stencil_size; i++) |
| symm_elements[i] = -1; |
| num_values = stencil_size; |
| } |
| else |
| { |
| num_values = (stencil_size + 1) / 2; |
| } |
| |
| hypre_StructMatrixStencil(matrix) = stencil; |
| hypre_StructMatrixSymmElements(matrix) = symm_elements; |
| hypre_StructMatrixNumValues(matrix) = num_values; |
| } |
| |
| /*----------------------------------------------------------------------- |
| * Set ghost-layer size for symmetric storage |
| * - All stencil coeffs are to be available at each point in the |
| * grid, as well as in the user-specified ghost layer. |
| *-----------------------------------------------------------------------*/ |
| |
| num_ghost = hypre_StructMatrixNumGhost(matrix); |
| stencil = hypre_StructMatrixStencil(matrix); |
| stencil_shape = hypre_StructStencilShape(stencil); |
| stencil_size = hypre_StructStencilSize(stencil); |
| symm_elements = hypre_StructMatrixSymmElements(matrix); |
| |
| for (i = 0; i < stencil_size; i++) |
| { |
| if (symm_elements[i] >= 0) |
| { |
| for (d = 0; d < 3; d++) |
| { |
| extra_ghost[2*d] = |
| hypre_max(extra_ghost[2*d], -hypre_IndexD(stencil_shape[i], d)); |
| extra_ghost[2*d + 1] = |
| hypre_max(extra_ghost[2*d + 1], hypre_IndexD(stencil_shape[i], d)); |
| } |
| } |
| } |
| |
| for (d = 0; d < 3; d++) |
| { |
| num_ghost[2*d] += extra_ghost[2*d]; |
| num_ghost[2*d + 1] += extra_ghost[2*d + 1]; |
| } |
| |
| /*----------------------------------------------------------------------- |
| * Set up data_space |
| *-----------------------------------------------------------------------*/ |
| |
| if (hypre_StructMatrixDataSpace(matrix) == NULL) |
| { |
| boxes = hypre_StructGridBoxes(grid); |
| data_space = hypre_BoxArrayCreate(hypre_BoxArraySize(boxes)); |
| |
| hypre_ForBoxI(i, boxes) |
| { |
| box = hypre_BoxArrayBox(boxes, i); |
| data_box = hypre_BoxArrayBox(data_space, i); |
| |
| hypre_CopyBox(box, data_box); |
| for (d = 0; d < 3; d++) |
| { |
| hypre_BoxIMinD(data_box, d) -= num_ghost[2*d]; |
| hypre_BoxIMaxD(data_box, d) += num_ghost[2*d + 1]; |
| } |
| } |
| |
| hypre_StructMatrixDataSpace(matrix) = data_space; |
| } |
| |
| /*----------------------------------------------------------------------- |
| * Set up data_indices array and data-size |
| *-----------------------------------------------------------------------*/ |
| |
| if (hypre_StructMatrixDataIndices(matrix) == NULL) |
| { |
| data_space = hypre_StructMatrixDataSpace(matrix); |
| data_indices = hypre_CTAlloc(int *, hypre_BoxArraySize(data_space)); |
| |
| data_size = 0; |
| hypre_ForBoxI(i, data_space) |
| { |
| data_box = hypre_BoxArrayBox(data_space, i); |
| data_box_volume = hypre_BoxVolume(data_box); |
| |
| data_indices[i] = hypre_CTAlloc(int, stencil_size); |
| |
| /* set pointers for "stored" coefficients */ |
| for (j = 0; j < stencil_size; j++) |
| { |
| if (symm_elements[j] < 0) |
| { |
| data_indices[i][j] = data_size; |
| data_size += data_box_volume; |
| } |
| } |
| |
| /* set pointers for "symmetric" coefficients */ |
| for (j = 0; j < stencil_size; j++) |
| { |
| if (symm_elements[j] >= 0) |
| { |
| data_indices[i][j] = data_indices[i][symm_elements[j]] + |
| hypre_BoxOffsetDistance(data_box, stencil_shape[j]); |
| } |
| } |
| } |
| |
| hypre_StructMatrixDataIndices(matrix) = data_indices; |
| hypre_StructMatrixDataSize(matrix) = data_size; |
| } |
| |
| /*----------------------------------------------------------------------- |
| * Set total number of nonzero coefficients |
| *-----------------------------------------------------------------------*/ |
| |
| hypre_StructMatrixGlobalSize(matrix) = |
| hypre_StructGridGlobalSize(grid) * stencil_size; |
| |
| /*----------------------------------------------------------------------- |
| * Return |
| *-----------------------------------------------------------------------*/ |
| |
| return ierr; |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixInitializeData |
| *--------------------------------------------------------------------------*/ |
| |
| int |
| hypre_StructMatrixInitializeData( hypre_StructMatrix *matrix, |
| double *data ) |
| { |
| int ierr = 0; |
| |
| hypre_BoxArray *data_boxes; |
| hypre_Box *data_box; |
| hypre_Index loop_size; |
| hypre_Index index; |
| hypre_IndexRef start; |
| hypre_Index stride; |
| double *datap; |
| int datai; |
| int i; |
| int loopi, loopj, loopk; |
| |
| hypre_StructMatrixData(matrix) = data; |
| hypre_StructMatrixDataAlloced(matrix) = 0; |
| |
| /*------------------------------------------------- |
| * If the matrix has a diagonal, set these entries |
| * to 1 everywhere. This reduces the complexity of |
| * many computations by eliminating divide-by-zero |
| * in the ghost region. |
| *-------------------------------------------------*/ |
| |
| hypre_SetIndex(index, 0, 0, 0); |
| hypre_SetIndex(stride, 1, 1, 1); |
| |
| data_boxes = hypre_StructMatrixDataSpace(matrix); |
| hypre_ForBoxI(i, data_boxes) |
| { |
| datap = hypre_StructMatrixExtractPointerByIndex(matrix, i, index); |
| |
| if (datap) |
| { |
| data_box = hypre_BoxArrayBox(data_boxes, i); |
| start = hypre_BoxIMin(data_box); |
| |
| hypre_BoxGetSize(data_box, loop_size); |
| |
| hypre_BoxLoop1Begin(loop_size, |
| data_box, start, stride, datai); |
| #define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,datai |
| #include "hypre_box_smp_forloop.h" |
| hypre_BoxLoop1For(loopi, loopj, loopk, datai) |
| { |
| datap[datai] = 1.0; |
| } |
| hypre_BoxLoop1End(datai); |
| } |
| } |
| |
| return ierr; |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixInitialize |
| *--------------------------------------------------------------------------*/ |
| int |
| hypre_StructMatrixInitialize( hypre_StructMatrix *matrix ) |
| { |
| int ierr = 0; |
| |
| double *data; |
| |
| ierr = hypre_StructMatrixInitializeShell(matrix); |
| |
| data = hypre_SharedCTAlloc(double, hypre_StructMatrixDataSize(matrix)); |
| hypre_StructMatrixInitializeData(matrix, data); |
| hypre_StructMatrixDataAlloced(matrix) = 1; |
| |
| return ierr; |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixSetValues |
| *--------------------------------------------------------------------------*/ |
| |
| int |
| hypre_StructMatrixSetValues( hypre_StructMatrix *matrix, |
| hypre_Index grid_index, |
| int num_stencil_indices, |
| int *stencil_indices, |
| double *values, |
| int add_to ) |
| { |
| int ierr = 0; |
| |
| hypre_BoxArray *boxes; |
| hypre_Box *box; |
| |
| double *matp; |
| |
| int i, s; |
| |
| boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(matrix)); |
| |
| hypre_ForBoxI(i, boxes) |
| { |
| box = hypre_BoxArrayBox(boxes, i); |
| |
| if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) && |
| (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) && |
| (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) && |
| (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) && |
| (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) && |
| (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) ) |
| { |
| if (add_to) |
| { |
| for (s = 0; s < num_stencil_indices; s++) |
| { |
| matp = hypre_StructMatrixBoxDataValue(matrix, i, |
| stencil_indices[s], |
| grid_index); |
| *matp += values[s]; |
| } |
| } |
| else |
| { |
| for (s = 0; s < num_stencil_indices; s++) |
| { |
| matp = hypre_StructMatrixBoxDataValue(matrix, i, |
| stencil_indices[s], |
| grid_index); |
| *matp = values[s]; |
| } |
| } |
| } |
| } |
| |
| return(ierr); |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixSetBoxValues |
| *--------------------------------------------------------------------------*/ |
| |
| int |
| hypre_StructMatrixSetBoxValues( hypre_StructMatrix *matrix, |
| hypre_Box *value_box, |
| int num_stencil_indices, |
| int *stencil_indices, |
| double *values, |
| int add_to ) |
| { |
| int ierr = 0; |
| |
| hypre_BoxArray *grid_boxes; |
| hypre_Box *grid_box; |
| hypre_BoxArray *box_array; |
| hypre_Box *box; |
| |
| hypre_BoxArray *data_space; |
| hypre_Box *data_box; |
| hypre_IndexRef data_start; |
| hypre_Index data_stride; |
| int datai; |
| double *datap; |
| |
| hypre_Box *dval_box; |
| hypre_Index dval_start; |
| hypre_Index dval_stride; |
| int dvali; |
| |
| hypre_Index loop_size; |
| |
| int i, s; |
| int loopi, loopj, loopk; |
| |
| /*----------------------------------------------------------------------- |
| * Set up `box_array' by intersecting `box' with the grid boxes |
| *-----------------------------------------------------------------------*/ |
| |
| grid_boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(matrix)); |
| box_array = hypre_BoxArrayCreate(hypre_BoxArraySize(grid_boxes)); |
| box = hypre_BoxCreate(); |
| hypre_ForBoxI(i, grid_boxes) |
| { |
| grid_box = hypre_BoxArrayBox(grid_boxes, i); |
| hypre_IntersectBoxes(value_box, grid_box, box); |
| hypre_CopyBox(box, hypre_BoxArrayBox(box_array, i)); |
| } |
| hypre_BoxDestroy(box); |
| |
| /*----------------------------------------------------------------------- |
| * Set the matrix coefficients |
| *-----------------------------------------------------------------------*/ |
| |
| if (box_array) |
| { |
| data_space = hypre_StructMatrixDataSpace(matrix); |
| hypre_SetIndex(data_stride, 1, 1, 1); |
| |
| dval_box = hypre_BoxDuplicate(value_box); |
| hypre_BoxIMinD(dval_box, 0) *= num_stencil_indices; |
| hypre_BoxIMaxD(dval_box, 0) *= num_stencil_indices; |
| hypre_BoxIMaxD(dval_box, 0) += num_stencil_indices - 1; |
| hypre_SetIndex(dval_stride, num_stencil_indices, 1, 1); |
| |
| hypre_ForBoxI(i, box_array) |
| { |
| box = hypre_BoxArrayBox(box_array, i); |
| data_box = hypre_BoxArrayBox(data_space, i); |
| |
| /* if there was an intersection */ |
| if (box) |
| { |
| data_start = hypre_BoxIMin(box); |
| hypre_CopyIndex(data_start, dval_start); |
| hypre_IndexD(dval_start, 0) *= num_stencil_indices; |
| |
| for (s = 0; s < num_stencil_indices; s++) |
| { |
| datap = hypre_StructMatrixBoxData(matrix, i, |
| stencil_indices[s]); |
| |
| hypre_BoxGetSize(box, loop_size); |
| |
| if (add_to) |
| { |
| hypre_BoxLoop2Begin(loop_size, |
| data_box,data_start,data_stride,datai, |
| dval_box,dval_start,dval_stride,dvali); |
| #define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,datai,dvali |
| #include "hypre_box_smp_forloop.h" |
| hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali) |
| { |
| datap[datai] += values[dvali]; |
| } |
| hypre_BoxLoop2End(datai, dvali); |
| } |
| else |
| { |
| hypre_BoxLoop2Begin(loop_size, |
| data_box,data_start,data_stride,datai, |
| dval_box,dval_start,dval_stride,dvali); |
| #define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,datai,dvali |
| #include "hypre_box_smp_forloop.h" |
| hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali) |
| { |
| datap[datai] = values[dvali]; |
| } |
| hypre_BoxLoop2End(datai, dvali); |
| } |
| |
| hypre_IndexD(dval_start, 0) ++; |
| } |
| } |
| } |
| |
| hypre_BoxDestroy(dval_box); |
| } |
| |
| hypre_BoxArrayDestroy(box_array); |
| |
| return(ierr); |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixAssemble |
| *--------------------------------------------------------------------------*/ |
| |
| int |
| hypre_StructMatrixAssemble( hypre_StructMatrix *matrix ) |
| { |
| int ierr = 0; |
| |
| int *num_ghost = hypre_StructMatrixNumGhost(matrix); |
| |
| hypre_BoxArrayArray *send_boxes; |
| hypre_BoxArrayArray *recv_boxes; |
| int **send_processes; |
| int **recv_processes; |
| |
| hypre_Index unit_stride; |
| |
| hypre_CommPkg *comm_pkg; |
| |
| hypre_CommHandle *comm_handle; |
| |
| /*----------------------------------------------------------------------- |
| * If the CommPkg has not been set up, set it up |
| *-----------------------------------------------------------------------*/ |
| |
| comm_pkg = hypre_StructMatrixCommPkg(matrix); |
| |
| if (!comm_pkg) |
| { |
| hypre_SetIndex(unit_stride, 1, 1, 1); |
| |
| hypre_CreateCommInfoFromNumGhost(hypre_StructMatrixGrid(matrix), |
| num_ghost, |
| &send_boxes, &recv_boxes, |
| &send_processes, &recv_processes); |
| |
| comm_pkg = hypre_CommPkgCreate(send_boxes, recv_boxes, |
| unit_stride, unit_stride, |
| hypre_StructMatrixDataSpace(matrix), |
| hypre_StructMatrixDataSpace(matrix), |
| send_processes, recv_processes, |
| hypre_StructMatrixNumValues(matrix), |
| hypre_StructMatrixComm(matrix), |
| hypre_StructGridPeriodic( |
| hypre_StructMatrixGrid(matrix))); |
| |
| hypre_StructMatrixCommPkg(matrix) = comm_pkg; |
| } |
| |
| /*----------------------------------------------------------------------- |
| * Update the ghost data |
| *-----------------------------------------------------------------------*/ |
| |
| hypre_InitializeCommunication(comm_pkg, |
| hypre_StructMatrixData(matrix), |
| hypre_StructMatrixData(matrix), |
| &comm_handle); |
| hypre_FinalizeCommunication(comm_handle); |
| |
| return(ierr); |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixSetNumGhost |
| *--------------------------------------------------------------------------*/ |
| |
| int |
| hypre_StructMatrixSetNumGhost( hypre_StructMatrix *matrix, |
| int *num_ghost ) |
| { |
| int ierr = 0; |
| int i; |
| |
| for (i = 0; i < 6; i++) |
| hypre_StructMatrixNumGhost(matrix)[i] = num_ghost[i]; |
| |
| return ierr; |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixPrint |
| *--------------------------------------------------------------------------*/ |
| |
| int |
| hypre_StructMatrixPrint( char *filename, |
| hypre_StructMatrix *matrix, |
| int all ) |
| { |
| int ierr = 0; |
| |
| FILE *file; |
| char new_filename[255]; |
| |
| hypre_StructGrid *grid; |
| hypre_BoxArray *boxes; |
| |
| hypre_StructStencil *stencil; |
| hypre_Index *stencil_shape; |
| |
| int num_values; |
| |
| hypre_BoxArray *data_space; |
| |
| int *symm_elements; |
| |
| int i, j; |
| |
| int myid; |
| |
| /*---------------------------------------- |
| * Open file |
| *----------------------------------------*/ |
| |
| #ifdef HYPRE_USE_PTHREADS |
| #if MPI_Comm_rank == hypre_thread_MPI_Comm_rank |
| #undef MPI_Comm_rank |
| #endif |
| #endif |
| |
| MPI_Comm_rank(hypre_StructMatrixComm(matrix), &myid); |
| |
| sprintf(new_filename, "%s.%05d", filename, myid); |
| |
| if ((file = fopen(new_filename, "w")) == NULL) |
| { |
| printf("Error: can't open output file %s\n", new_filename); |
| exit(1); |
| } |
| |
| /*---------------------------------------- |
| * Print header info |
| *----------------------------------------*/ |
| |
| fprintf(file, "StructMatrix\n"); |
| |
| fprintf(file, "\nSymmetric: %d\n", hypre_StructMatrixSymmetric(matrix)); |
| |
| /* print grid info */ |
| fprintf(file, "\nGrid:\n"); |
| grid = hypre_StructMatrixGrid(matrix); |
| hypre_StructGridPrint(file, grid); |
| |
| /* print stencil info */ |
| fprintf(file, "\nStencil:\n"); |
| stencil = hypre_StructMatrixStencil(matrix); |
| stencil_shape = hypre_StructStencilShape(stencil); |
| |
| num_values = hypre_StructMatrixNumValues(matrix); |
| symm_elements = hypre_StructMatrixSymmElements(matrix); |
| fprintf(file, "%d\n", num_values); |
| j = 0; |
| for (i = 0; i < hypre_StructStencilSize(stencil); i++) |
| { |
| if (symm_elements[i] < 0) |
| { |
| fprintf(file, "%d: %d %d %d\n", j++, |
| hypre_IndexX(stencil_shape[i]), |
| hypre_IndexY(stencil_shape[i]), |
| hypre_IndexZ(stencil_shape[i])); |
| } |
| } |
| |
| /*---------------------------------------- |
| * Print data |
| *----------------------------------------*/ |
| |
| data_space = hypre_StructMatrixDataSpace(matrix); |
| |
| if (all) |
| boxes = data_space; |
| else |
| boxes = hypre_StructGridBoxes(grid); |
| |
| fprintf(file, "\nData:\n"); |
| hypre_PrintBoxArrayData(file, boxes, data_space, num_values, |
| hypre_StructMatrixData(matrix)); |
| |
| /*---------------------------------------- |
| * Close file |
| *----------------------------------------*/ |
| |
| fflush(file); |
| fclose(file); |
| |
| return ierr; |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixMigrate |
| *--------------------------------------------------------------------------*/ |
| |
| int |
| hypre_StructMatrixMigrate( hypre_StructMatrix *from_matrix, |
| hypre_StructMatrix *to_matrix ) |
| { |
| hypre_BoxArrayArray *send_boxes; |
| hypre_BoxArrayArray *recv_boxes; |
| int **send_processes; |
| int **recv_processes; |
| |
| hypre_Index unit_stride; |
| |
| hypre_CommPkg *comm_pkg; |
| hypre_CommHandle *comm_handle; |
| |
| int ierr = 0; |
| |
| /*------------------------------------------------------ |
| * Set up hypre_CommPkg |
| *------------------------------------------------------*/ |
| |
| hypre_SetIndex(unit_stride, 1, 1, 1); |
| |
| hypre_CreateCommInfoFromGrids(hypre_StructMatrixGrid(from_matrix), |
| hypre_StructMatrixGrid(to_matrix), |
| &send_boxes, &recv_boxes, |
| &send_processes, &recv_processes); |
| |
| comm_pkg = hypre_CommPkgCreate(send_boxes, recv_boxes, |
| unit_stride, unit_stride, |
| hypre_StructMatrixDataSpace(from_matrix), |
| hypre_StructMatrixDataSpace(to_matrix), |
| send_processes, recv_processes, |
| hypre_StructMatrixNumValues(from_matrix), |
| hypre_StructMatrixComm(from_matrix), |
| hypre_StructGridPeriodic( |
| hypre_StructMatrixGrid(from_matrix))); |
| /* is this correct for periodic? */ |
| |
| /*----------------------------------------------------------------------- |
| * Migrate the matrix data |
| *-----------------------------------------------------------------------*/ |
| |
| hypre_InitializeCommunication(comm_pkg, |
| hypre_StructMatrixData(from_matrix), |
| hypre_StructMatrixData(to_matrix), |
| &comm_handle); |
| hypre_FinalizeCommunication(comm_handle); |
| |
| return ierr; |
| } |
| |
| /*-------------------------------------------------------------------------- |
| * hypre_StructMatrixRead |
| *--------------------------------------------------------------------------*/ |
| |
| hypre_StructMatrix * |
| hypre_StructMatrixRead( MPI_Comm comm, |
| char *filename, |
| int *num_ghost ) |
| { |
| FILE *file; |
| char new_filename[255]; |
| |
| hypre_StructMatrix *matrix; |
| |
| hypre_StructGrid *grid; |
| hypre_BoxArray *boxes; |
| int dim; |
| |
| hypre_StructStencil *stencil; |
| hypre_Index *stencil_shape; |
| int stencil_size; |
| |
| int num_values; |
| |
| hypre_BoxArray *data_space; |
| |
| int symmetric; |
| |
| int i, idummy; |
| |
| int myid; |
| |
| /*---------------------------------------- |
| * Open file |
| *----------------------------------------*/ |
| |
| #ifdef HYPRE_USE_PTHREADS |
| #if MPI_Comm_rank == hypre_thread_MPI_Comm_rank |
| #undef MPI_Comm_rank |
| #endif |
| #endif |
| |
| MPI_Comm_rank(comm, &myid ); |
| |
| sprintf(new_filename, "%s.%05d", filename, myid); |
| |
| if ((file = fopen(new_filename, "r")) == NULL) |
| { |
| printf("Error: can't open output file %s\n", new_filename); |
| exit(1); |
| } |
| |
| /*---------------------------------------- |
| * Read header info |
| *----------------------------------------*/ |
| |
| fscanf(file, "StructMatrix\n"); |
| |
| fscanf(file, "\nSymmetric: %d\n", &symmetric); |
| |
| /* read grid info */ |
| fscanf(file, "\nGrid:\n"); |
| hypre_StructGridRead(comm,file,&grid); |
| |
| /* read stencil info */ |
| fscanf(file, "\nStencil:\n"); |
| dim = hypre_StructGridDim(grid); |
| fscanf(file, "%d\n", &stencil_size); |
| stencil_shape = hypre_CTAlloc(hypre_Index, stencil_size); |
| for (i = 0; i < stencil_size; i++) |
| { |
| fscanf(file, "%d: %d %d %d\n", &idummy, |
| &hypre_IndexX(stencil_shape[i]), |
| &hypre_IndexY(stencil_shape[i]), |
| &hypre_IndexZ(stencil_shape[i])); |
| } |
| stencil = hypre_StructStencilCreate(dim, stencil_size, stencil_shape); |
| |
| /*---------------------------------------- |
| * Initialize the matrix |
| *----------------------------------------*/ |
| |
| matrix = hypre_StructMatrixCreate(comm, grid, stencil); |
| hypre_StructMatrixSymmetric(matrix) = symmetric; |
| hypre_StructMatrixSetNumGhost(matrix, num_ghost); |
| hypre_StructMatrixInitialize(matrix); |
| |
| /*---------------------------------------- |
| * Read data |
| *----------------------------------------*/ |
| |
| boxes = hypre_StructGridBoxes(grid); |
| data_space = hypre_StructMatrixDataSpace(matrix); |
| num_values = hypre_StructMatrixNumValues(matrix); |
| |
| fscanf(file, "\nData:\n"); |
| hypre_ReadBoxArrayData(file, boxes, data_space, num_values, |
| hypre_StructMatrixData(matrix)); |
| |
| /*---------------------------------------- |
| * Assemble the matrix |
| *----------------------------------------*/ |
| |
| hypre_StructMatrixAssemble(matrix); |
| |
| /*---------------------------------------- |
| * Close file |
| *----------------------------------------*/ |
| |
| fclose(file); |
| |
| return matrix; |
| } |
| |
