tutorial

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

                *******************************************************
                ***** Preconditioned TFQMR solution
                ***** No preconditioning
                ***** No scaling
                *******************************************************
 
                iter:    0              residual = 1.000000e+00
                iter:    1              residual = 2.193551e-01
                iter:    2              residual = 7.539287e-02
                iter:    3              residual = 3.318523e-02
                iter:    4              residual = 1.690645e-02
                iter:    5              residual = 9.507550e-03
                iter:    6              residual = 5.724587e-03
                iter:    7              residual = 3.588558e-03
                iter:    8              residual = 2.265129e-03
                iter:    9              residual = 1.380014e-03
                iter:   10              residual = 7.732288e-04
                iter:   11              residual = 3.832266e-04
                iter:   12              residual = 1.652267e-04
                iter:   13              residual = 6.181623e-05
                iter:   14              residual = 2.009490e-05
                iter:   15              residual = 1.346214e-06


Problem 11: 3D VBR example (n x n x n)
==========
Change the 3D Poisson operator to approximate the system of 
m coupled PDE equations on an n x n x n grid:

Pois( u^(0) ) + .1    sum    Pois( u^(i) )  = f^0
                  0 <= i  < m 
                    i != 0

Pois( u^(1) ) + .1    sum    Pois( u^(i) )  = f^1
                  0 <= i  < m 
                    i != 1


                       .
                       .
                       .
Pois( u^(m-1) ) + .1    sum    Pois( u^(i) )  = f^(m-1)
                  0 <= i  < m 
                    i != m-1

where Pois(v) = v_xx + v_yy + v_zz.

Use the VBR format and let the value of m be a global variable
input by the user. Run this program with the same initial guess
(x = 0) and the same combined right hand side (B = e_1, that is 
f^0 = e_1, f^1 = 0, f^2 = 0, ..., f^(m-1) = 0).

Hint: The m x m diagonal blocks of this matrix should consist of 
the value 6. on the diagonal entries and the value .6 on all the 
off-diagonals entries.  The m x m off-diagonal blocks of this 
matrix should consist of the value -1. on the diagonal entries and 
the value -.1 on all the off diagonals entries.

If this program is run with 7 and 3 as input (i.e. a system of
3 coupled PDEs on a 7 x 7 x 7 grid), the iteration output should
approximately be:

                *******************************************************
                ***** Preconditioned TFQMR solution
                ***** No preconditioning
                ***** No scaling
                *******************************************************
 
                iter:    0              residual = 1.000000e+00
                iter:    1              residual = 2.885746e-01
                iter:    2              residual = 1.241962e-01
                iter:    3              residual = 6.121277e-02
                iter:    4              residual = 3.385888e-02
                iter:    5              residual = 2.037280e-02
                iter:    6              residual = 1.293223e-02
                iter:    7              residual = 8.521449e-03
                iter:    8              residual = 5.685374e-03
                iter:    9              residual = 3.711614e-03
                iter:   10              residual = 2.284157e-03
                iter:   11              residual = 1.284914e-03
                iter:   12              residual = 6.570774e-04
                iter:   13              residual = 3.102746e-04
                iter:   14              residual = 1.369291e-04
                iter:   15              residual = 5.682475e-05
                iter:   16              residual = 6.471154e-06


+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
                     Answers to exercises:
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   
Problem 1:
=========
Change line 73 to 
 73    AZ_read_update(&N_update,&update,proc_config,nrow,1,AZ_file);

Create the file .update. To run on 4 processors, this file should
contain the following information:

   9
 3  7 11 15 19 23 27 31 35
   9
 2  6 10 14 18 22 26 30 34
   9
 1  5  9 13 17 21 25 29 33
   9
 0  4  8 12 16 20 24 28 32


To run on 1 processor, this file should contain:

   36
 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 
22 23 24 25 26 27 28 29 30 31 32 33 34 35

Problem 2:
=========
Change line 72 to 
 72    nrow=4;

The uncommented new create_matrix_row() routine should resemble: 

void create_matrix_row(int row,int location,double val[], int bindx[])
{
   int k;
 
   k = bindx[location];
   if (row == 0) {
      bindx[k] = 1;   val[k++] = 3;
      bindx[k] = 2;   val[k++] = 2;
      bindx[k] = 3;   val[k++] = 1;
                      val[location] = 7.;
   }
   else if (row == 1) {
      bindx[k] = 0;   val[k++] = 5;
      bindx[k] = 2;   val[k++] = 3;
      bindx[k] = 3;   val[k++] = 1;
                      val[location] = 8.;
   }
   else if (row == 2) {
      bindx[k] = 0;   val[k++] = 5;
      bindx[k] = 1;   val[k++] = 3;
      bindx[k] = 3;   val[k++] = 2;
                      val[location] = 9.;
   }
   else if (row == 3) {
      bindx[k] = 0;   val[k++] = 1;
      bindx[k] = 1;   val[k++] = 2;
      bindx[k] = 2;   val[k++] = 1;
                      val[location] = 6.;
   }
   else perror("Invalid row number");
   bindx[location+1] = k;
}

Problem 3: 
=========
After line 98 add:

 98a   options[AZ_solver] = AZ_tfqmr;
 98b   params[AZ_tol]     = 1.0e-5;

Problem 4: 
=========
After line 68 (corresponding to original program listing) add:

 68a   if (proc_config[AZ_node] == 0) {
 68b      printf("Enter grid size (n, for n x n grid)\n");
 68c      scanf("%d",&n);
 68d   }
 68e   AZ_broadcast((char *) &n, sizeof(int), proc_config, AZ_PACK);
 68f   AZ_broadcast(NULL       ,           0, proc_config, AZ_SEND);

and enter the value 7 when prompted for the grid size.
   
Problem 5: 
=========
Change line 100 (corresponding to original program listing) to 

 100   /* Set rhs (delta function at grid corners) and initialize guess */

After line 109 (corresponding to original program listing) add:

 109a    if (update[i] == n-1  ) b[update_index[i]] = 1.0;
 109b    if (update[i] == n*n-n) b[update_index[i]] = 1.0;
 109c    if (update[i] == n*n-1) b[update_index[i]] = 1.0;

Alternatively, b can be assigned via 

        if (update[i] == 0    ) b[i] = 1.0;
        if (update[i] == n-1  ) b[i] = 1.0;
        if (update[i] == n*n-n) b[i] = 1.0;
        if (update[i] == n*n-1) b[i] = 1.0;
and AZ_reorder_vec() could be invoked after the loop to
reorder b so that it corresponds to the transformed system.

Problem 6:
=========
In our solution, we move the declaration of b (line 102 of the 
original program) to just after line 79 (of original program).

Just after the 'create_matrix_row' line (line 87 in original program), we add
      b[i] = 0.0;

Just before the 'AZ_transform' line (line 92) , we add

  for (i = 0; i < N_update; i++) {
        if (update[i] == 0    ) b[i] = 1.0;
        if (update[i] == n-1  ) b[i] = 1.0;
        if (update[i] == n*n-n) b[i] = 1.0;
        if (update[i] == n*n-1) b[i] = 1.0;
  }
Just after the 'AZ_tranform' line, we add
    AZ_reorder_vec(b, data_org, update_index, NULL);

and then we remove all references to 'b[update_index]'.


Problem 7:
=========
After line 116 (corresponding to original program listing), add:

 116a
 116b   i = AZ_find_index(  n-1, update, N_update);
 116c   if (i != -1) b[update_index[i]] = 0.0;
 116d   i = AZ_find_index(n*n-1, update, N_update);
 116e   if (i != -1) b[update_index[i]] = 0.0;
 116f
 116g   AZ_solve(x,b, options, params, NULL,bindx,NULL,NULL,NULL,val,
 116h          data_org, status, proc_config);

Note: when solving the same system mulitple times with different right 
hand sides, it is possible to reuse factorization information produced 
by the matrix scaling options and some of the preconditioners (AZ_ilu, 
AZ_lu, AZ_bilu). To do this, set options[AZ_pre_calc] to AZ_reuse
after the first solve AZ_solve().
   
Problem 8:
=========
We declared a new external function 

 21a   extern void post_create_matrix_row_5pt(int row, int location, double val[], 
 21b                     int bindx[], int update[], int update_index[], 
 21c                     int external[], int extern_index[], int data_org[]);

and added the following lines before the second call to AZ_solve().

 116f1  for (i = 0 ; i < N_update ; i++ ) {
 116f2     if (update[i]%n == n-1) 
 116f3        post_create_matrix_row_5pt(update[i],update_index[i],val,bindx, update, 
 116f4                       update_index, external, extern_index, data_org);
 116f5  }
 116f6

The uncommented new subroutine follows:
    
    void post_create_matrix_row_5pt(int row, int location, double val[], int bindx[],
                        int update[], int update_index[], int external[],
                        int extern_index[], int data_org[])
    {
       int i, k, neighbor = -1;
 
       if (row%n != 0) {
          i = AZ_find_index(row-1, update, data_org[AZ_N_internal] +
                                       data_org[AZ_N_border]);
          if (i == -1) {
             i = AZ_find_index(row-1, external, data_org[AZ_N_external]);
             if (i != -1) neighbor = extern_index[i];
          }
          else neighbor = update_index[i];
       }
 
       for (k = bindx[location] ; k < bindx[location+1] ; k++ )
          if (bindx[k] != neighbor) val[k] = 0.;
 
       val[location]     = 1.; /* matrix diagonal */
    }


Problem 9:
=========
In addition to the functions left() and right() given in the problem
it is necessary to make the following changes.

Replace 

 21    extern void create_matrix_row(int row,int i,double val[],int bindx[]);

by

 21    extern void create_matrix_row(int row,int location,double val[],int bindx[],
 21a                      int yline_head, int yval,int n);
 21b   extern int left(int yval); 
 21c   extern int right(int yval);
 21d   int    yval, yline_head, next;