az_examples.c

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

/*====================================================================
 * ------------------------
 * | CVS File Information |
 * ------------------------
 *
 * $RCSfile: az_examples.c,v $
 *
 * $Author: tuminaro $
 *
 * $Date: 1998/12/21 19:17:26 $
 *
 * $Revision: 1.10 $
 *
 * $Name:  $
 *====================================================================*/
#ifndef lint
static char rcsid[] = "$Id: az_examples.c,v 1.10 1998/12/21 19:17:26 tuminaro Exp $";
#endif


/*******************************************************************************
 * Copyright 1995, Sandia Corporation.  The United States Government retains a *
 * nonexclusive license in this software as prescribed in AL 88-1 and AL 91-7. *
 * Export of this program may require a license from the United States         *
 * Government.                                                                 *
 ******************************************************************************/


/******************************************************************************
 *
 * These routines correspond to the matrix setup phase for a number of sample
 * PDE problems using AZTEC. All of the examples are discretizations of
 * the Poisson equation (2D structured using 5pt operator, 2D structured using
 * 9pt operator, 2D unstructured using finite elements, 3D structured using
 * 7pt operator). All the examples are for MSR matrices excluding one VBR
 * example.  The examples are discussed in further detail in the paper:
 *
 *    Tuminaro, R.S., Shadid, J.N., and Hutchinson, H.A., Parallel Sparse
 *    Matrix Vector Multiply Software for Matrix with Data Locality",
 *    Sandia Technical Report, 1995.
 *
 *****************************************************************************/

extern int N_grid_pts, num_PDE_eqns;

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "az_aztec.h"

/* external function declarations */

extern void create_msr_matrix(int update[], double **val, int **bindx, 
                              int N_update);
extern void add_row_3D(int row, int i, double val[], int bindx[], int *n);
extern void add_row_5pt(int row, int location, double val[], int bindx[], 
			int *n);
extern void add_row_9pt(int row, int location, int *next_nonzero, double val[],
                 int bindx[], int *n);
extern void create_vbr_matrix(int update_blks[], double **val, int **indx,
                       int num_update_blks, int **rnptr, int **bnptr,
                       int **bindx);
extern void check_memory_space(int bnptr_next, int total_blks,
                               int indx_next, int total_nz);
extern void add_vbr_row(int row, int location, double val[], int indx[],
                        int rnptr[], int bnptr[], int bindx[]);
extern int get_block_size(int );
extern void new_block(int *next , int row_blk_size, int col_blk_size,
                      int blk_col, double number, int bindx[], int indx[],
                      double val[]);

extern void read_triangles(int proc, int N_update);
extern void init_msr(double **val, int **bindx, int N_update);
extern void compress_matrix(double val[], int bindx[], int N_update);
extern void add_to_element(int row, int column, double element, double val[],
                           int bindx[], int diag);
extern void setup_Ke(double Ke[][3], double x1, double , double ,
                     double y2, double x3, double y3);
extern void read_coordinates(double **x, double **y, int update_indx[],int,int);
extern void add_to_element(int row, int column, double element, double val[],
                           int bindx[], int diag);
extern void create_fe_matrix(int update[], int proc, int **bindx, double **val,
                      int N_update);
extern void fill_fe_matrix(double val[], int bindx[], int update[], 
	int update_indx[], int external[], int extern_indx[], int data_org[],
	int proc_config[]);
extern void add_row_matrix_free(int row, int location, int *next_nonzero, 
        double val[], int bindx[]);
extern void matrix_vector_mult(int N_update, int update[], int update_indx[],
                      int N_external, int external[], int extern_indx[],
                      double x[], double y[], int data_org[],int proc_config[]);





/*****************************************************************************
 *                                                                           *
 *                             MSR                                           *
 *                                                                           *
 *****************************************************************************/

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

{
  int i,total_nz, n = -1;
  int avg_nonzeros_per_row = 9;


  total_nz = N_update*avg_nonzeros_per_row + 1;
  *bindx   = (int *) AZ_allocate(total_nz*sizeof(int));
  *val     = (double *) AZ_allocate(total_nz*sizeof(double));
  if ((*val == NULL) && (total_nz != 0) ) {
    (void) fprintf(stderr, "Error: Not enough space to create matrix\n");
    (void) fprintf(stderr,
                   "      Try reducing the variable 'avg_nonzeros_per_row'\n");
    exit(1);
  }
  for (i = 0 ; i < total_nz ; i++ ) (*bindx)[i] = 0;

  (*bindx)[0] = N_update+1;
  for (i = 0; i < N_update; i++) {
    add_row_3D(update[i], i, *val, *bindx,&n);
    if ( (*bindx)[i+1] > total_nz) {
      (void) fprintf(stderr, "Error:total_nz not large enough to accomodate");
      (void) fprintf(stderr, " nonzeros\n       Try increasing the variable");
      (void) fprintf(stderr, " 'avg_nonzeros_per_row'\n       Finished \n");
      (void) fprintf(stderr, "first %d rows out of %d\n", i, N_update);
      exit(1);
    }
  }

} /* create_msr_matrix */

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

void add_row_5pt(int row, int location, double val[], int bindx[], int *n)

/*
 * Add one row to an MSR matrix corresponding to a 5pt discrete approximation
 * to the 2D Poisson operator on an n x n square.
 *
 * Author: Ray Tuminaro, Div 1422
 * Date:   3/15/95
 *
 * Parameters:
 *    row          == the global row number of the new row to be added.
 *    location     == the local row number where the diagonal of the new row
 *                    will be stored.
 *    val,bindx    == (see user's guide). On output, val[] and bindx[]
 *                    are appended such that the new row has been added.
 */

{

/*
  static int n = -1;
*/
  int m;
  int        k;
  int        NP;

  /* determine grid dimensions */

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

      /* special case */

      val[0] = 1.0;
      bindx[1] = bindx[0];
      return;
    }
  }
  *n = m;

  k  = bindx[location];
  NP = num_PDE_eqns;

  /*
   * Check neighboring points in each direction and add nonzero entry if
   * neighbor exists.
   */

  bindx[k] = row + NP;   if ((row/NP)%m !=     m-1) val[k++] = -1.00;
  bindx[k] = row - NP;   if ((row/NP)%m !=       0) val[k++] = -1.00;
  bindx[k] = row + m*NP; if ((row/(NP*m))%m != m-1) val[k++] = -1.00;
  bindx[k] = row - m*NP; if ((row/(NP*m))%m !=   0) val[k++] = -1.00;

  bindx[location+1] = k;
  val[location]     = 4.0;

} /* add_row_5pt */


/***************************************************************
 *                                                             *
 *           Other examples                                    *
 *                                                             *
 **************************************************************/


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

void add_row_3D(int row, int location, double val[], int bindx[], int *n)

/*
 * Add one row to an MSR matrix corresponding to a 7pt discrete approximation
 * to the 3D Poisson operator on an n x n x n square.
 *
 * Author: Ray Tuminaro, Div 1422
 * Date:   3/15/95
 *
 * Parameters:
 *    row          == the global row number of the new row to be added.
 *    location     == the local row number where the diagonal of the new row
 *                    will be stored.
 *    val,bindx    == (see Aztec User's guide). On output, val[] and bindx[]
 *                    are appended such that the new row has been added.
 */

{

  int        m;
  int        k, NP;

  /* determine grid dimensions */

  m = *n; 
  if (m == -1) {
    m = N_grid_pts;
    m = (int ) pow( ((double) m )+0.01, 0.33334);
    if (m == 1) {

      /* special case */

      val[0] = 1.0;
      bindx[1] = bindx[0];
      return;
    }

    if (m*m*m != 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);
    }
  }
  *n = m;
  k = bindx[location];
  NP = num_PDE_eqns;

  /*
   * Check neighboring points in each direction and add nonzero entry if
   * neighbor exists.
   */

  bindx[k] = row + NP;     if ((row/NP)%m      != m-1) val[k++] = -1.0;
  bindx[k] = row - NP;     if ((row/NP)%m      != 0  ) val[k++] = -1.0;
  bindx[k] = row + m*NP;   if ((row/(NP*m))%m  != m-1) val[k++] = -1.0;
  bindx[k] = row - m*NP;   if ((row/(NP*m))%m  != 0  ) val[k++] = -1.0;
  bindx[k] = row + m*m*NP; if (bindx[k]    < m*m*m*NP) val[k++] = -1.0;
  bindx[k] = row - m*m*NP; if (bindx[k]          >= 0) val[k++] = -1.0;

  bindx[location+1] = k;
  val[location]     = 6.0;

} /* add_row_3D */

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

void add_row_9pt(int row, int location, int *next_nonzero, double val[],
                 int bindx[], int *n)

/*
 * Add one row to an MSR matrix corresponding to a 9pt discrete approximation
 * to the 2D Poisson operator on an n x n square.
 *
 * Author: Ray Tuminaro, Div 1422
 * Date:   3/15/95
 *
 * Parameters:
 *    row          == the global row number of the new row to be added.
 *    location     == the local row number where the diagonal of the new row
 *                    will be stored.
 *    next_nonzero == points to the next free storage location in val[] and
 *                    bindx[].
 *                    On input, next_nonzero points to location where off
 *                    diagonals in new row will be stored.
 *                    On output, next_nonzero points to the next free location
 *                    (i.e. after the new row that was added).
 *    val,bindx    == (see file 'parameters'). On output, val[] and bindx[]
 *                    are appended such that the new row has been added.
 */

{
  int        i, j, point, old_ptr, stride, m;

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

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

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

  point = row/num_PDE_eqns;
  i     = point%m;
  j     = (point-i)/m;

  old_ptr = *next_nonzero;
  val[location] = 0.0;

  for (stride = 1; stride <= m; stride = stride*m) {

    /*
     * If we are on the bottom boundary, create entries for 2nd order
     * discretization. Otherwise create high order discretization.
     */

    if ((i != 0) && (i != m-1)) {
      if (i != 1) {
        val[*next_nonzero]   = 1.0/12.0;
        bindx[*next_nonzero] = row - stride*2*num_PDE_eqns;
        *next_nonzero        = *next_nonzero + 1;
      }
      val[*next_nonzero]   = -16.0/12.0;
      bindx[*next_nonzero] = row - stride*num_PDE_eqns;
      *next_nonzero        = *next_nonzero + 1;
      if (i != m-2) {
        val[*next_nonzero]   =   1.0/12.0;
        bindx[*next_nonzero] = row + stride*2*num_PDE_eqns;
        *next_nonzero        = *next_nonzero + 1;
      }
      val[*next_nonzero]   = -16.0/12.0;
      bindx[*next_nonzero] = row + stride*num_PDE_eqns;
      *next_nonzero        = *next_nonzero + 1;
    }
    else if (i != m-1) {
      val[*next_nonzero]   = -1.0;
      bindx[*next_nonzero] = row + stride*num_PDE_eqns;
      *next_nonzero        = *next_nonzero + 1;
    }
    else if (i != 0) {
      val[*next_nonzero]   = -1.0;
      bindx[*next_nonzero] = row - stride*num_PDE_eqns;
      *next_nonzero        = *next_nonzero + 1;
    }
    if ((i != 0) && (i != m-1)) val[location] += 30.0/12.0;
    else val[location] += 2.0;
    i = j;
  }
  bindx[location+1] = bindx[location] + (*next_nonzero - old_ptr);

} /* add_row_9pt */


/*****************************************************************************
 *                                                                           *
 *                             VBR                                           *
 *                                                                           *
 *****************************************************************************/


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

void create_vbr_matrix(int update_blks[], double **val, int **indx,
                       int num_update_blks, int **rnptr, int **bnptr,
                       int **bindx)
/*
 * Create a matrix in VBR format (containing 'num_update_blks' block rows).
 *
 * Author: Ray Tuminaro, Div 1422
 * Date:   3/1/95
 *
 * Parameter list:
 *
 *  update_blks     == list of blocks updated by this processor
 *  val             == On output, contains matrix nonzeros
 *  indx            == On output, contains points to blocks of nonzeros
 *                     in val[]. IN particular, the pointers to the nonzero
 *                     blocks in block row i are
 *                        indx[bnptr[i], ... , bnptr[i+1]-1 ]
 *  num_update_blks == Number of blocks updated by this processor.
 *  rnptr           == On output, rnptr[0] = 0 and rnptr[i+1]-rnptr[i] gives
 *                     the row dimension of global block row 'update_blks[i]'.
 *  bnptr           == On output, bnptr[i] points to
 *                        1) the list of nonzero blocks in bindx[] of the
 *                           global row 'update_blks[i]'.
 *                        2) the list of pointers to nonzero matrix values
 *                           in indx[] of the global row 'update_blks[i]'.
 *  bindx           == On output, bindx[bnptr[i], ..., bnptr[i+1]-1] lists
 *                     the global block column indices of the nonzero
 *                     blocks in global row 'update_blks[i]'.
 */

{

  int avg_blks_per_row = 7;