tutorial

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

#====================================================================
# ------------------------
# | CVS File Information |
# ------------------------
# 
# $RCSfile: az_tutorial,v $
#
# $Author: tuminaro $
#
# $Date: 2000/06/02 16:49:21 $
#
# $Revision: 1.6 $
#
# $Name:  $
#====================================================================

                          AZTEC Tutorial

A sample AZTEC application is provided in the file 'az_tutorial.c'. The 11 
exercises below are performed by continually modifying the sample program. 
By doing the exercises, the new user will be introduced to many of AZTEC's
features.  Unless otherwise stated, problem x will require modification to 
the program used for problem x-1. Complete answers to the exercises are
given at the end of this document.

To compile and run the original sample program, the OBJ line in the 
application Makefile must be changed to : OBJ = az_tutorial.o. 'az_tutorial.c' 
corresponds to setting up and solving a 2D Poisson approximation on a 6 x 6 
grid where the right hand side corresponds to a delta function in the lower 
left corner of the grid. Below we illustrate the program body:

line no.
-------

    62    /* get number of processors and the name of this processor */
    63  #ifdef AZ_MPI
    64    MPI_Init(&argc,&argv);
    65    AZ_set_proc_config(proc_config, MPI_COMM_WORLD);
    66  #else
    67    AZ_set_proc_config(proc_config, AZ_NOT_MPI);
    68 #endif
    69
    70    /* Define partitioning:  matrix rows (ascending order) owned by this node */
    71    
    72    nrow = n*n;
    73    AZ_read_update(&N_update, &update, proc_config, nrow, 1, AZ_linear);
    74
    75    /*
    76     * Create the matrix: each processor creates only rows appearing in update[]
    77     * (using global col. numbers).
    78     */
    79
    80    bindx = (int    *) calloc(N_update*MAX_NZ_ROW+1,sizeof(int));
    81    val   = (double *) calloc(N_update*MAX_NZ_ROW+1,sizeof(double));
    82    if (val == NULL) perror("Error: Not enough space to create matrix");
    83
    84    bindx[0] = N_update+1;
    85
    86    for (i = 0; i < N_update; i++) {
    87      create_matrix_row(update[i], i, val, bindx);
    88    }
    89  
    90    /* convert matrix to a local distributed matrix */
    91  
    92    AZ_transform(proc_config, &external, bindx, val, update, &update_index,
    93                 &extern_index, &data_org, N_update, NULL, NULL, NULL, NULL,
    94                 AZ_MSR_MATRIX);
    95  
    96    /* initialize AZTEC options */
    97  
    98    AZ_defaults(options, params);
    99  
   100    /* Set rhs (delta function at lower left corner) and initialize guess */
   101  
   102    b = (double *) calloc(N_update, sizeof(double));
   103    x = (double *) calloc(N_update + data_org[AZ_N_external], sizeof(double));
   104    if ((x == NULL) && (i != 0)) perror("Not enough space in rhs");
   105  
   106    for (i = 0; i < N_update; i++) {
   107      x[update_index[i]] = 0.0;
   108      b[update_index[i]] = 0.0;
   109      if (update[i] == 0) b[update_index[i]] = 1.0;
   110    }
   111  
   112  
   113    /* solve the system of equations using b  as the right hand side */
   114  
   115    AZ_solve(x, b, options, params, NULL, bindx, NULL, NULL, NULL, val, data_org,
   116             status, proc_config);
   117  
   118    /* Free allocated memory */
   119  
   120  
   121    free((void *) update);   free((void *) update_index);
   122    free((void *) external); free((void *) extern_index);
   123  
   124    free((void *) x);    free((void *) b);       free((void *) bindx);
   125    free((void *) val);  free((void *) data_org);
   126  
   127  
   128  #ifdef AZ_MPI
   129    MPI_Finalize();
   130  #endif
   131  
   132  } /* main */
   133  
   134  
   135  
   136  
   137  
   138  
   139  
   140  /***************************************************************************/
   141  /***************************************************************************/
   142  /***************************************************************************/
   143  
   144  void create_matrix_row(int row, int location, double val[], int bindx[])
   145  
   146  /* Add one row to an MSR matrix corresponding to a discrete approximation to the
   147   * 2D Poisson operator on an n x n square.
   148   *
   149   * Parameters:
   150   *    row          == global row number of the new row to be added.
   151   *    location     == local row where diagonal of the new row will be stored.
   152   *    val,bindx    == (see user's guide). On output, val[] and bindx[]
   153   *                    are appended such that the new row has been added.
   154   */
   155  
   156  {
   157    int k;
   158  
   159    /* check neighbors in each direction and add nonzero if neighbor exits */
   160  
   161    k = bindx[location];
   162    bindx[k]  = row + 1; if ((row  )%n != n-1) val[k++] = -1.;
   163    bindx[k]  = row - 1; if ((row  )%n !=   0) val[k++] = -1.;
   164    bindx[k]  = row + n; if ((row/n)%n != n-1) val[k++] = -1.;
   165    bindx[k]  = row - n; if ((row/n)%n !=   0) val[k++] = -1.;
   166  
   167    bindx[location+1] = k;  val[location]     = 4.; /* matrix diagonal */
   168  }



If you run this program, the iteration output should approximately be the 
following:

                *******************************************************
                ***** Preconditioned GMRES solution
                ***** No preconditioning
                ***** No scaling
                *******************************************************
 
                iter:    0              residual = 1.000000e+00
                iter:    1              residual = 3.333333e-01
                iter:    2              residual = 1.549193e-01
                iter:    3              residual = 8.498208e-02
                iter:    4              residual = 5.178089e-02
                iter:    5              residual = 3.392084e-02
                iter:    6              residual = 2.343361e-02
                iter:    7              residual = 1.673219e-02
                iter:    8              residual = 1.183812e-02
                iter:    9              residual = 7.732077e-03
                iter:   10              residual = 4.398328e-03
                iter:   11              residual = 1.241776e-03
                iter:   12              residual = 5.528073e-04
                iter:   13              residual = 1.918082e-04
                iter:   14              residual = 6.826307e-05
                iter:   15              residual = 1.470372e-05
                iter:   16              residual = 3.001987e-06
                iter:   17              residual = 2.111772e-07
 
---

IMPORTANT: SAVE A COPY OF THIS PROGRAM IN THE FILE 'original.c'.

Problem 1: 2D Poisson PARTITIONING (6 x 6)
=========
To partition or distribute problems over processors, each processor
must set 'N_update' to the number of matrix rows it will own and the 
array 'update' to the specific matrix row numbers (in ascending order)
owned by this processor. In our sample program, a linear partitioning is 
used via the function AZ_read_update() (line 73).

Change the AZ_read_update() call so that it instead reads the paritioning
from the file '.update'. Initialize this file so that the row r (0 <= r < 36)
is assigned to processor number r%P (mod function in "C") where P is
the number of processors used with this '.update' file. 

Note: AZ_read_update() is not particularly fast in that it reads a 
single file and thus distributes rows serially. For large data sets
it is generally better to use some kind of parallel input/output.
See the AZ_read_update() description in the Aztec User's Guide for 
more information.

The output of this program should be the same as that of the original
sample problem as we have only changed the assignment of rows to processors.  

Problem 2: DENSE MSR matrix (4 x 4)
=========
IMPORTANT: COPY THE FILE 'original.c' INTO 'az_tutorial.c'. YOU MIGHT WANT 
TO SAVE 'az_tutorial.c' BEFORE DOING THIS.

Change this program to solve the matrix application

   /              \        / \
   |  7  3  2  1  |       | 1 |  
   |  5  8  3  1  |  x  = | 0 |
   |  5  3  9  2  |       | 0 |
   |  1  2  1  6  |       | 0 |
   \              /       \   /

Using x = 0 as an initial guess.

Hint: Change in line 72 to indicate the number of matrix rows (nrow=4). 
Other than this, only the subroutine create_matrix_row() needs to 
be changed. The "Data Formats" and "High Level Data Interface" sections 
of the Aztec User's Guide explain the MSR data format used to represent 
the matrix. Essentially, to append a row to a matrix already containing 
j rows (rows are locally numbered from 0 to j-1), we need to do the 
following:

   a) set index, k, to point to the free space for the offdiagonal elements:
      k = bindx[j]

   b) Store each off diagonal element as follows:

         val[k] = 1st off-diagonal, bindx[k] = 1st off-diagonal column number
         k = k + 1
         val[k] = 2nd off-diagonal, bindx[k] = 2nd off-diagonal column number
         etc

   c) put the diagonal element of the new row in  val[j]

   d) store free space pointer in bindx[j+1] (i.e. bindx[j+1] = k).


This program should produce approximately the following iteration output:

 
                *******************************************************
                ***** Preconditioned GMRES solution
                ***** No preconditioning
                ***** No scaling
                *******************************************************
 
                iter:    0              residual = 1.000000e+00
                iter:    1              residual = 7.141428e-01
                iter:    2              residual = 1.687498e-01
                iter:    3              residual = 2.119830e-02
                iter:    4              residual = 1.928932e-16

---


Problem 3: Aztec OPTIONS (solver) - 2D Poisson (6 x 6)
=========
IMPORTANT: COPY THE FILE 'original.c' INTO 'az_tutorial.c'. YOU MIGHT WANT 
TO SAVE 'az_tutorial.c' BEFORE DOING THIS.

The original sample program uses the GMRES algorithm as this is the default 
AZTEC method (indicated by setting options[AZ_solver] in AZ_defaults). Change
the sample program so that we use the TFQMR method and so that the algorithm
terminates when the initial residual is reduced by 10^5 (instead of the 
default 10^6).  See the "Aztec Options" section of the Aztec User's Guide for 
more information.
Note: this example is just for demonstration purposes as we should really use
the conjugate gradient algorithm for this symmetric problem.

The  resulting iteration output of the program should now approximately be:


                *******************************************************
                ***** Preconditioned TFQMR solution
                ***** No preconditioning
                ***** No scaling
                *******************************************************
 
                iter:    0              residual = 1.000000e+00
                iter:    1              residual = 3.078276e-01
                iter:    2              residual = 1.352989e-01
                iter:    3              residual = 7.202456e-02
                iter:    4              residual = 4.273759e-02
                iter:    5              residual = 2.571415e-02
                iter:    6              residual = 1.391762e-02
                iter:    7              residual = 6.317909e-03
                iter:    8              residual = 2.292283e-03
                iter:    9              residual = 6.648821e-04
                iter:   10              residual = 1.506405e-04
                iter:   11              residual = 1.590716e-05
                iter:   12              residual = 9.021697e-07
 
---
 
Problem 4: BROADCASTING n ( grid-size parameter ) - 2D Poisson (n x n)
=========
Currently, the grid size parameter, n (for an n x n grid), is fixed.