tutorial

来自「并行解法器,功能强大」· 代码 · 共 1,185 行 · 第 1/4 页

Replace

 68b      printf("Enter grid size (n, for n x n grid)\n");

by 

 68b      printf("Enter grid size (n, n horizontal lines in grid)\n");

Replace

 73    AZ_read_update(&N_update,&update,proc_config,nrow,1,AZ_linear);

by

 73    total = 0;   /* total is set to the total number of grid points */
 73a   for (i = 0 ; i < n ; i++ ) { total += (right(i) - left(i)+1);}
 73b   AZ_read_update(&N_update,&update,proc_config,total,1,AZ_linear);

Note: total needs to be declared as well.
Replace

 86    for (i = 0 ; i < N_update ; i++ ) {
 87       create_matrix_row(update[i],i,val,bindx);
 88    }

by 

 86    yval     = 0;   yline_head = 0;
 87    next     = right(0) - left(0) + 1;
 88    for (i = 0 ; i < N_update ; i++ ) {
 88a      while ( next <= update[i] ) {
 88b         yline_head = next;   yval++;
 88c         next = yline_head + right(yval) - left(yval) + 1;
 88d      }
 88e      create_matrix_row(update[i],i,val,bindx,yline_head,yval,n);
 88f   }

and add the subroutine create_matrix_row():

void create_matrix_row(int row,int location,double val[],int bindx[],
             int yline_head, int yval,int n)
{
/* Add one row to an MSR matrix corresponding to a discrete approximation
 * on an irregular 2D domain. Specifically, the grid contains n horizontal
 * lines. The x values on line j start at left(j) and end at right(j).
 *
 * Parameters:
 *    row          == global row number of the new row to be added.
 *    location     == local row where 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.
 *    yline_head   == global row number of the first grid point in the
 *                    horizontal line containing 'row'.
 *    yval         == horizontal line y = yval contains 'row'.
 *    n            == total number of horizontal lines in grid.
 */
   int b_neigh = 0, t_neigh = 0, l_neigh = 0, r_neigh = 0;
   int lval0 = 0, rval1 = 0, lval2 = 0, rval2 = 0, lval3 = 0, rval3 = 0;
   extern int left(), right();
   int normalized, k;
 
   normalized = yline_head - left(yval);
   l_neigh    = row-1;   lval0 = normalized +  left(yval) + 1;
   r_neigh    = row+1;   rval1 = normalized + right(yval) - 1;
 
   if (yval != 0) {
      b_neigh = row+left(yval)-right(yval-1)-1;
      lval2   = normalized +  left(yval-1);
      rval2   = normalized + right(yval-1);
   }
   else rval2 = lval2 - 1;
   if (yval != n-1) {
      t_neigh = row-left(yval+1)+right(yval)+1;
      lval3   = normalized +  left(yval+1);
      rval3   = normalized + right(yval+1);
   }
   else rval3 = lval3 - 1;
 
   /* check neighbors in each direction and add nonzero if neighbor exits */

   k = bindx[location];
   bindx[k] = l_neigh; if  (row >= lval0)                    val[k++] = -1.;
   bindx[k] = r_neigh; if  (row <= rval1)                    val[k++] = -1.;
   bindx[k] = b_neigh; if ((row >= lval2) && (row <= rval2)) val[k++] = -1.;
   bindx[k] = t_neigh; if ((row >= lval3) && (row <= rval3)) val[k++] = -1.;
   bindx[location+1] = k;  val[location]     = 4.; /* matrix diagonal */
}
  
Problem 10:
==========
Change line 68b in Problem 4 program to 

 68b      printf("Enter grid size (n, for n x n x n grid)\n");

Change line 72 in Problem 4 program to 

 72    nrow=n*n*n;

Change 
       #define MAX_NZ_ROW 5
to 
       #define MAX_NZ_ROW 7

so that the proper amount of memory is allocated in lines 80 & 81.

Change line 147 to 

 147  * to the 3D Poisson operator on an n x n x n square.

After line 165, add the following lines:

165a   bindx[k]  = row + n*n; if ((row/(n*n))%n != n-1) val[k++] = -1.;
165b   bindx[k]  = row - n*n; if ((row/(n*n))%n !=   0) val[k++] = -1.;

and finally change the matrix diagonal in line 167 to 

117    bindx[location+1] = k;  val[location]     = 6.; /* matrix diagonal */

Problem 11:
==========
Our new program is as follows:

/* This software was developed at Sandia National Labs under US Energy Dept.
 * contract DE-AC-4-76DP00789 and is copyrighted by Sandia Corporation.  */

#include <stdio.h>
#include <stdlib.h>
#ifdef AZ_MPI
#include <mpi.h>
#endif
#include "az_aztec.h"
#define perror(str) { fprintf(stderr,"%s\n",str);   exit(-1); }

int m, n = 6;      /* EQUATIONS WILL BE SOLVED ON an n x n x n GRID. */

void main(int argc, char *argv[])

/* Set up a coupled PDE test problem and solve it with AZTEC.                */

{
   double *b,*x;                    /* rhs and approximate solution          */
   int    i, j, nrow;
   extern void create_matrix_row(int row,int location, int rpntr[],int bpntr[],
			int bindx[], int indx[],double val[]);

            /* See Aztec User's Guide for the variables that follow:         */
                               
   int    proc_config[AZ_PROC_SIZE];/* Processor information.                */
   int    options[AZ_OPTIONS_SIZE]; /* Array used to select solver options.  */
   double params[AZ_PARAMS_SIZE];   /* User selected solver paramters.       */
   int    *data_org;                /* Array to specify data layout          */
   double status[AZ_STATUS_SIZE];   /* Information returned from AZ_solve(). */
   int    *update,                  /* vector elements updated on this node. */
          *external;                /* vector elements needed by this node.  */
   int    *update_index;            /* ordering of update[] and external[]   */
   int    *extern_index;            /* locally on this processor.            */
   int    *bindx;                   /* Sparse matrix to be solved is stored  */
   double *val;                     /* in these MSR arrays.                  */
   int    *indx,*cpntr,*bpntr,*rpntr;
   int    N_update;                 /* # of unknowns updated on this node    */

#ifdef AZ_MPI
   MPI_Init(&argc,&argv);
#endif

   /* get number of processors and the name of this processor */

   AZ_processor_info(proc_config);

   if (proc_config[AZ_node] == 0) {
      printf("Enter grid size (n, for n x n x n grid)\n");
      scanf("%d",&n);
      printf("Enter number of PDEs\n");
      scanf("%d",&m);
   }
   AZ_broadcast((char *) &n, sizeof(int), proc_config, AZ_PACK);
   AZ_broadcast((char *) &m, sizeof(int), proc_config, AZ_PACK);
   AZ_broadcast(NULL       ,           0, proc_config, AZ_SEND);

   /* Define partitioning:matrix rows (ascending order) owned by this node */

   nrow = n*n*n;
   AZ_read_update(&N_update,&update,proc_config,nrow,1,AZ_linear);

   /* create the matrix: each processor creates only rows  */
   /* appearing in update[] (using global col. numbers).   */

   rpntr   = (int    *) calloc(N_update+1,sizeof(int));
   bpntr   = (int    *) calloc(N_update+1,sizeof(int));
   bindx   = (int    *) calloc(N_update*7+1,sizeof(int));
   indx    = (int    *) calloc(N_update*7+1,sizeof(int));
   val     = (double *) calloc(m*m*N_update*7+1,sizeof(double));
   if (val == NULL) perror("Error: Not enough space to create matrix")

   rpntr[0] = 0; bpntr[0] = 0; indx[0]  = 0;
   for (i = 0 ; i < N_update ; i++ ) {
      create_matrix_row(update[i],i,rpntr,bpntr,bindx,indx,val);
   }

   /* convert matrix to a local distributed matrix */

   AZ_transform(proc_config,&external,bindx,val,update,&update_index,
        &extern_index,&data_org,N_update,indx,bpntr,rpntr,&cpntr,AZ_VBR_MATRIX);

   /* initialize AZTEC options */

   AZ_defaults(options, params);
   options[AZ_solver] = AZ_tfqmr;
   params[AZ_tol]     = 1.0e-5;

  /* Set rhs (delta function at lower left corner) and initialize guess */

  b = (double *) calloc(m*N_update,sizeof(double));
  x = (double *) calloc(m*N_update + data_org[AZ_N_external],sizeof(double));
  if ((x == NULL) && (i != 0)) perror("Not enough space in rhs")

  for (i=0 ; i < N_update; i++ ) { 
    for (j = 0 ; j < m ; j++ ) {
       x[m*update_index[i]+j] = 0.0;
       b[m*update_index[i]+j] = 0.0;
    }
    if (update[i] == 0 ) b[m*update_index[i]] = 1.0;
  }

  /* solve the system of equations using b  as the right hand side */

   AZ_solve(x,b, options, params, indx,bindx,rpntr,cpntr,bpntr,val,
	    data_org, status, proc_config);

#ifdef AZ_MPI
   MPI_Finalize();
#endif
}

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

void create_matrix_row(int row,int location, int rpntr[],int bpntr[],int bindx[],
	int indx[],double val[])
{
/* Add one row to an MSR matrix corresponding to a discrete approximation
 * to the 3D Poisson operator on an n x n x n square.
 *
 * Parameters:
 *    row          == global row number of the new row to be added.
 *    location     == local row where 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.
 */
   int k;
   extern void VBR_block(int indx[],int *k,double val[],double factor);

   /* check neighbors in each direction and add nonzero if neighbor exits */

   rpntr[location+1] = rpntr[location] + m;

   k = bpntr[location];
   bindx[k]  = row;                                 VBR_block(indx,&k,val, 6.);
   bindx[k]  = row + 1;   if ((row  )%n != n-1)     VBR_block(indx,&k,val,-1.);
   bindx[k]  = row - 1;   if ((row  )%n !=   0)     VBR_block(indx,&k,val,-1.);
   bindx[k]  = row + n;   if ((row/n)%n != n-1)     VBR_block(indx,&k,val,-1.);
   bindx[k]  = row - n;   if ((row/n)%n !=   0)     VBR_block(indx,&k,val,-1.);
   bindx[k]  = row + n*n; if ((row/(n*n))%n != n-1) VBR_block(indx,&k,val,-1.);
   bindx[k]  = row - n*n; if ((row/(n*n))%n !=   0) VBR_block(indx,&k,val,-1.);

   bpntr[location+1] = k; 
}

void VBR_block(int indx[],int *k,double val[],double factor)
{
   int i,j,indx_ptr;

   indx_ptr = indx[*k];
   indx[*k+1] = indx_ptr + m*m;  
   for (i = 0 ; i < m ; i++ ) {
     for (j = 0 ; j < m ; j++ ) {
        if (i == j) val[indx_ptr++] = factor;
	else val[indx_ptr++] = factor * .1;
     }
   }
   (*k)++;
}

All the necessary changes include:
  1) new declarations (m,j,*indx,*cpntr,*bpntr,*rpntr)
  2) changes in the create_matrix_row(), AZ_transform(), and AZ_solve() parameters
  3) adding 3 lines to input m and broadcast it.
  4) storage allocation for val,rpntr,bpntr,indx,val,b and x
  5) initializing rpntr[0], bpntr[0], indx[0] to 0
  6) changing the right hand side and initial guess for VBR format
  7) changing create_matrix_row() to call VBR_block instead of setting val[].
  8) adding the subroutine VBR_block()