az_examples.c

来自「并行解法器,功能强大」· C语言 代码 · 共 1,282 行 · 第 1/3 页

  int avg_blk_size     = num_PDE_eqns*num_PDE_eqns;
  int i, total_nz, total_blks;

  /* --------------------- execution begins ----------------------------------*/

  total_blks = num_update_blks*avg_blks_per_row;
  total_nz   = total_blks*avg_blk_size + 1;

  *indx  = (int *)    AZ_allocate( (total_blks+1)*sizeof(int));
  *bindx = (int *)    AZ_allocate( (total_blks+1)*sizeof(int));
  *rnptr = (int *)    AZ_allocate( (num_update_blks+1)*sizeof(int));
  *bnptr = (int *)    AZ_allocate( (num_update_blks+1)*sizeof(int));
  *val   = (double *) AZ_allocate( (total_nz)*sizeof(double));

  if (*val == NULL) {
    (void) fprintf(stderr, "Error: Not enough space to create matrix\n"
                   "      Try reducing avg_blks_per_row,avg_blk_size\n");
    exit(1);
  }

  (*bnptr)[0] = 0;
  (*rnptr)[0] = 0;
  (*indx )[0] = 0;
  for (i = 0; i < num_update_blks; i++) {
    add_vbr_row(update_blks[i], i, *val, *indx, *rnptr, *bnptr, *bindx);

    /* do a series of memory checks to make sure we have not exceeded space */

    check_memory_space((*bnptr)[i+1], total_blks, (*indx)[(*bnptr)[i+1]],
                       total_nz);
  }

} /* create_vbr_matrix */

/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/

void add_vbr_row(int row, int location, double val[], int indx[], int rnptr[],
                 int bnptr[], int bindx[])

/*
 * Create one block row of a VBR matrix corresponding to block row 'row'
 * of a series of m uncoupled 2D Poisson equations.
 *
 * Author: Ray Tuminaro, Div 1422
 * Date:   3/1/95
 *
 * Parameter list:
 *   row           == Global index of new block row to be added.
 *   location      == Place (local index) in rnptr,bnptr where the
 *                    new row is to be added.
 *   val           == On output, contains matrix nonzeros for new row.
 *   indx          == On output, contains new pointers to blocks of nonzeros
 *                    in val[] corresponding to new row. In particular, the
 *                    pointers to the nonzero blocks are in
 *                      indx[bnptr[location], ... , bnptr[location+1]-1]
 *   rnptr         == On output, rnptr[location+1] - rnptr[location] gives the
 *                    row dimension of the new block.
 *   bnptr         == On output, bnptr[location] points to
 *                       1) the global block column list of nonzero blocks in
 *                          bindx[] of the new block row.
 *                       2) the list of pointers to nonzero matrix values
 *                          in indx[] of the new block row.
 *   bindx         == On output, bindx[bnptr[location],...,bnptr[location+1]-1]
 *                    lists the global block column indices of the nonzero
 *                    blocks in the new row.
 */

{

  int    i, j, k, stride;
  static int n = -1;
  int        row_blk_size, col_blk_size;

  /* -------------------------- execution begins -----------------------------*/

  if (n == -1) {
    n = N_grid_pts;
    n = (int ) pow( ((double) n )+0.01,0.33334);
    if (n <= 1) {
      (void) fprintf(stderr, "grid size too small\n");
      exit(1);
    }

    if ( n*n*n != N_grid_pts) {
      (void) fprintf(stderr, "Error: the total number of points (%d) needs\n",
                     N_grid_pts);
      (void) fprintf(stderr, "       to equal k^3 where k is an integer\n");
      exit(1);
    }
  }

  /* figure out where we are in the global grid */

  i = row%n;
  k = (row-i)/n;
  j = k%n;
  k = (k-j)/n;

 /* put in the diagonal block */

  row_blk_size = get_block_size(row);
  bnptr[location+1] = bnptr[location];
  rnptr[location+1] = rnptr[location] + row_blk_size;
  new_block(&bnptr[location+1], row_blk_size, row_blk_size, row, 6.0,
            bindx, indx, val);

  /* check in each of the directions and add a new block column */

  for (stride = 1; stride <= n*n; stride *= n) {
    if (i != 0) {
      col_blk_size = get_block_size(row-stride);
      new_block(&bnptr[location+1], row_blk_size, col_blk_size, row - stride,
                -1.0, bindx, indx, val);
    }

    if (i != n-1) {
      col_blk_size = get_block_size(row+stride);
      new_block(&bnptr[location+1], row_blk_size, col_blk_size, row + stride,
                -1.0, bindx, indx, val);
    }
    if (stride == 1) i = j;
    else i = k;
  }

} /* add_vbr_row */

/******************************************************************************/
/******************************************************************************/
/******************************************************************************/

void check_memory_space(int indx_next, int total_blks, int val_next,
                        int total_nz)

/*
 * Check that we have not used more space then we allocated in creating
 * the vbr matrix.
 *
 * Author: Ray Tuminaro, Div 1422
 * Date:   3/1/95
 *
 * Parameters:
 *
 *    indx_next        == pointer to next available location in indx[].
 *    total_blks       == total number of blocks for which we have allocated
 *                        space.
 *    val_next         == point to next available location in val[].
 *    total_nz         == total number of nonzeros allocated for matrix.
 */

{

  if (indx_next > total_blks) {
    (void) fprintf(stderr, "Error: total_blks not large enough to accomodate "
                   "nonzeros\n       Try increasing the variable "
                   "'avg_blks_per_row'\n");
    exit(1);
  }

  if (val_next > total_nz) {
    (void) fprintf(stderr,
                   "Error: total_nz not large enough to accomodate nonzeros\n"
                   "      Try increasing the variable 'avg_blk_size'\n");
    exit(1);
  }

}

/******************************************************************************/
/******************************************************************************/
/******************************************************************************/

int get_block_size(int location)

/*
 * Return the row block size of the global block row 'location'.
 */

{
  if (location < 0)  {
     (void) fprintf(stderr,"Error: Improper grid location (%d) ",location);
     (void) fprintf(stderr,"inside get_block_size()\n");
     exit(-1);
  }

  return num_PDE_eqns;

}

/******************************************************************************/
/******************************************************************************/
/******************************************************************************/

void new_block(int *next, int row_blk_size, int col_blk_size, int blk_col,
               double number, int bindx[], int indx[], double val[])

/*
 * Add a new block to the VBR matrix, initialize all the offdiagonals to zeros,
 * and initialize all the diagonal elements to 'number'.
 *
 * Author: Ray Tuminaro
 * Date:   3/15/95
 * Parameters:
 *   next          == On input, indx[*next] points to the new block.
 *                    On output, *next points to the next free location
 *                    in indx[] and bindx[].
 *   row_blk_size  == number of rows in the new block to be added.
 *   col_blk_size  == number of columns in the new block to be added.
 *   blk_col       == Block column index of new block to be added.
 *   number        == On output, each element of the 'diagonal' of the new
 *                    block contains 'number'
 *   bindx,indx,val== See file 'parameters'. On output, bindx[],indx[],val[]
 *                    have been appended such that the new block is added
 *                    to the VBR matrix.
 */

{

  int k, kk, start;

  bindx[*next]  = blk_col;
  indx[*next+1] = indx[*next] + row_blk_size*col_blk_size;
  start = indx[*next];
  for (k = 0; k < row_blk_size; k++) {
    for (kk = 0; kk < col_blk_size; kk++) {
      if (k == kk) val[start + kk*row_blk_size+k] = number;
      else val[start + kk*row_blk_size+k] = 0.0;
    }
  }

  (*next)++;

}


/*****************************************************************************
 *                                                                           *
 *                             FINITE ELEMENT MSR                            *
 *                                                                           *
 *****************************************************************************/

/******************************************************************************/
/******************************************************************************/

#define NOT_FOUND -1
FILE *fp;
int  *T[3];
int  N_triangle;

void create_fe_matrix(int update[], int proc, int **bindx, double **val,
                      int N_update)

/*
 * Create a matrix in MSR format corresponding to a finite element
 * discretization of the 2D Poisson operator. The finite element grid
 * is contained in the files 'fe_grid_#' (where each file contains only
 * the grid information for processor #).
 *
 * NOTE: this routine does not actually compute the nonzero entries
 * of the matrix. It only computes the nonzero pattern (i.e. bindx[]).
 * The subroutine fill_fe_matrix() computes the nonzero entries and
 * must be called after the communication pattern has been set.
 *
 *
 *      Author:         Lydie Prevost Div 1422 SNL
 *      Date:           11/25/1994
 *
 * Parameter List:
 *
 *  update  == list of unknowns updated by this processor
 *  bindx   == On output, contains nonzero locations in matrix where
 *
 *                 bindx[i],  0 <=i< N_update: ptr to offdiagonal nonzeros
 *                                             for row 'update[i]'.
 *                 bindx[i],  N_update > i     contains column numbers
 *                                             corresponding to nonzero val[i].
 *  N_update== Number of unknowns updated by this processor.
 *  proc    == Node number of this processor.
 */

{

  int i, j,k;
  int row;

  /**************************** execution begins ******************************/

  /* read data */

  read_triangles(proc, N_update);

  /* initialize msr matrix */

  init_msr(val, bindx, N_update);

  for (k = 0; k < N_triangle; k++) {
    for (i = 0; i < 3; i++) {
      row = AZ_find_index(T[i][k], update, N_update);
      if (row != NOT_FOUND ) {
        for (j = 0; j < 3; j++)
          add_to_element(row, T[j][k], 0.0, *val, *bindx, i==j);
      }
    }
  }

  compress_matrix(*val, *bindx, N_update);

} /* create_fe_matrix */

/******************************************************************************/
/******************************************************************************/
/******************************************************************************/

void fill_fe_matrix(double val[], int bindx[], int update[], int update_indx[],
                    int external[], int extern_indx[], int data_org[],
		    int proc_config[])

/*
 * Fill in the nonzero entries of an MSR matrix corresponding to a finite
 * element discretization of the 2D Poisson operator. It is assumed that
 * the routine create_fe_matrix() has already been called.
 *
 *      Author:         Lydie Prevost Div 1422 SNL
 *      Date:           11/25/1994
 *
 * Parameter List:
 *
 *  update   == list of unknowns updated by this processor
 *  val      == On output, contains matrix nonzeros where
 *
 *               val[i],0 <= i < N_update: diagonal for row update[i].
 *               val[bindx[i]-->val[bindx[i+1]-1]: nonzeros for row 'update[i]'.
 *
 *  N_update == Number of unknowns updated by this processor.
 *  proc     == Node number of this processor.
 *
 *      Author:         Lydie Prevost Div 1422 SNL
 *      Date:           11/25/1994
 *
 */

{

  /* local variables */

  int    i_triangle, i, j;
  double Ke[3][3];
  double *x,*y;
  int    row;
  int    N_update,N_external;

  /**************************** execution begins ******************************/

  N_update   = data_org[AZ_N_internal] + data_org[AZ_N_border];
  N_external = data_org[AZ_N_external];

  read_coordinates(&x, &y, update_indx, N_update, N_external);

  AZ_exchange_bdry(x, data_org, proc_config);
  AZ_exchange_bdry(y, data_org, proc_config);

  /* relabel the triangle with the appropriate local index number */

  for (i_triangle = 0; i_triangle < N_triangle; i_triangle++) {
    for (i = 0; i < 3; i++) {
      row = AZ_find_index(T[i][i_triangle], update, N_update);

      if (row == NOT_FOUND) {
        row = AZ_find_index(T[i][i_triangle], external, N_external);
        T[i][i_triangle] = extern_indx[row];
      }
      else T[i][i_triangle] = update_indx[row];
    }
  }

  /* fill in the discretization */

  for (i_triangle = 0; i_triangle < N_triangle; i_triangle++) {

    /* initialize submatrix Ke containing contributions for this triangle */

    setup_Ke(Ke, x[T[0][i_triangle]], y[T[0][i_triangle]], x[T[1][i_triangle]],
             y[T[1][i_triangle]], x[T[2][i_triangle]], y[T[2][i_triangle]]);

    /* store submatrix Ke elements in a[] */

    for (i = 0; i < 3; i++) {
      if (T[i][i_triangle] < N_update) {
        for (j = 0; j < 3; j++) {
          add_to_element(T[i][i_triangle], T[j][i_triangle], Ke[i][j], val,
                         bindx, i==j);
        }
      }
    }
  }

} /* fill_fe_matrix */

/******************************************************************************/
/******************************************************************************/
/******************************************************************************/

void add_to_element(int row, int column, double element, double val[],
                    int bindx[], int diag)

/*
 * Add the value 'element' to the matrix entry A(row,column) and store the
 * result back in the matrix.
 *
 * Note: A is encoded in the MSR arrays val[] and bindx[].
 */

{

  int i, start, end;

  /* (row,column) corresponds to a diagonal element of the matrix */

  if (diag) val[row] += element;
  else {
    start = bindx[row];
    end   = bindx[row+1] - 1;

    /* search for the column */

    for (i = start; i <= end; i++) {
      if ( (bindx[i] == -1) || (bindx[i] == column)) break;