title: OpenACC and CUDA aware MPI
body: <p>I want to move on the device the whole while loop in the main. The problems emerges when I add <code>#pragma acc host_data use_device(err)</code> to <code>MPI_Allreduce (&amp;err, &amp;err, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);</code>.</p>
  <p>The error is that the reduction on err doesn't work so that the code exit after one step from the loop.</p>
  <p>After the <code>MPI_Allreduce()</code>, even using <code>#pragma acc update self(err)</code>, err is still equal to zero.</p>
  <p>I'm compiling with <code>mpicc -acc -ta=tesla:managed -Minfo=accel -w jacobi.c</code>
  And running with <code>mpirun -np 2 -mca pml ^ucx ./a.out</code></p>
  <p>Could you help me to find the error?</p>
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;math.h&gt;
  
  #define PARALLEL
  #define NX_GLOB     128   /* Global number of interior points */
  #define NY_GLOB     128  /* Global number of interior points */
  #define NGHOST   1
  #define NDIM     2
  
  #ifdef PARALLEL
    #include &lt;mpi.h&gt;
    MPI_Comm MPI_COMM_CART;
  #endif
  
  typedef struct MPI_Decomp_{
    int nprocs[NDIM];     /*  Number of processors in each dimension */
    int periods[NDIM];    /*  Periodicity flag in each dimension     */
    int coords[NDIM];     /*  Cartesian coordinate in the MPI topology */
    int gsize[NDIM];      /*  Global domain size (no ghosts)  */
    int lsize[NDIM];      /*  Local domain size (no ghosts)   */
    int start[NDIM];      /*  Local start index in each dimension           */
    int procL[NDIM];      /*  Rank of left-lying  process in each direction */
    int procR[NDIM];      /*  Rank of right-lying process in each direction */
    int rank;             /*  Local process rank */
    int size;             /*  Communicator size  */
  } MPI_Decomp;
  
  
  
  void BoundaryConditions(double **, double *, double *, int, int, MPI_Decomp *);
  void DomainDecomposition(MPI_Decomp *);
  void WriteSolution (double **, int, int, MPI_Decomp *);
  double **Allocate_2DdblArray(int, int);
  int    **Allocate_2DintArray(int, int);
  void   Show_2DdblArray(double **, int, int, const char *);
  void   Show_2DintArray(int **, int, int, const char *);
  
  int nx_tot, ny_tot;
  
  int main(int argc, char ** argv)
  {
    int    nx, i, ibeg, iend;
    int    ny, j, jbeg, jend;
    int    k, rank=0, size=1;
    double xbeg = 0.0, xend = 1.0;
    double ybeg = 0.0, yend = 1.0;
    double dx = (xend - xbeg)/(NX_GLOB + 1);
    double dy = (yend - ybeg)/(NY_GLOB + 1);
    double *xg, *yg, *x, *y, **phi, **phi0;
    double err, tol;
    MPI_Decomp  mpi_decomp;
    double err_glob;
    int procL[NDIM] = {-1,-1};
    int procR[NDIM] = {-1,-1};
  
  /* --------------------------------------------------------
     0. Initialize the MPI execution environment
     -------------------------------------------------------- */
  
  #ifdef PARALLEL
    MPI_Datatype row_type, col_type;
  
    MPI_Init(&amp;argc, &amp;argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &amp;rank);
    MPI_Comm_size(MPI_COMM_WORLD, &amp;size);
  
    DomainDecomposition(&amp;mpi_decomp);
    nx = mpi_decomp.lsize[0];
    ny = mpi_decomp.lsize[1];
  #else
    mpi_decomp.gsize[0] = mpi_decomp.lsize[0] = nx = NX_GLOB;
    mpi_decomp.gsize[1] = mpi_decomp.lsize[1] = ny = NY_GLOB;
    mpi_decomp.procL[0] = mpi_decomp.procL[1] = -1;
    mpi_decomp.procR[0] = mpi_decomp.procR[1] = -1;
  #endif
  
  /* --------------------------------------------------------
     1. Set local grid indices
     -------------------------------------------------------- */
  
    ibeg   = NGHOST;
    iend   = ibeg + nx - 1;
    nx     = iend - ibeg + 1;
    nx_tot = nx + 2*NGHOST;
  
    jbeg   = NGHOST;
    jend   = jbeg + ny - 1;
    ny     = jend - jbeg + 1;
    ny_tot = ny + 2*NGHOST;
  
  /* --------------------------------------------------------
     2. Generate global and local grids
     -------------------------------------------------------- */
  
    xg = (double *) malloc ( (NX_GLOB+2*NGHOST)*sizeof(double));
    yg = (double *) malloc ( (NY_GLOB+2*NGHOST)*sizeof(double));
    for (i = 0; i &lt; (NX_GLOB+2*NGHOST); i++) xg[i] = xbeg + (i-ibeg+1)*dx;
    for (j = 0; j &lt; (NY_GLOB+2*NGHOST); j++) yg[j] = ybeg + (j-jbeg+1)*dy;
    #ifdef PARALLEL
    x = xg + mpi_decomp.start[0];
    y = yg + mpi_decomp.start[1];
    #else
    x = xg;
    y = yg;
    #endif
  
  /* --------------------------------------------------------
     3. Allocate memory on local processor and
        assign initial conditions.
     -------------------------------------------------------- */
  
    phi  = Allocate_2DdblArray(ny_tot, nx_tot);
    phi0 = Allocate_2DdblArray(ny_tot, nx_tot);
  
    for (j = jbeg; j &lt;= jend; j++){
    for (i = ibeg; i &lt;= iend; i++){
      phi0[j][i] = 0.0;
    }}
  
  #ifdef PARALLEL
    MPI_Type_contiguous (nx_tot, MPI_DOUBLE, &amp;row_type);
    MPI_Type_vector     (ny_tot, 1, nx_tot, MPI_DOUBLE, &amp;col_type);
    MPI_Type_commit (&amp;row_type);
    MPI_Type_commit (&amp;col_type);
  #endif
  
  /* --------------------------------------------------------
     4. Main iteration cycle
     -------------------------------------------------------- */
  
    tol = 1.e-5;
    err = 1.0;
    k   = 0;
  
    #pragma acc enter data copyin(phi[:ny_tot][:nx_tot], phi0[:ny_tot][:nx_tot], x[:NX_GLOB+2*NGHOST], y[NX_GLOB+2*NGHOST], err, err_glob)
    while  (err &gt; tol){
  
    /* -- 4a. Set boundary conditions first -- */
  
      BoundaryConditions(phi0, x, y, nx, ny, &amp;mpi_decomp);
  
    /* -- 4b. Jacobi's method and residual (interior points) -- */
  
      err = 0.0;
  
      #pragma acc parallel loop collapse(2) reduction(+:err) present(phi[:ny_tot][:nx_tot], phi0[:ny_tot][:nx_tot])
      for (j = jbeg; j &lt;= jend; j++){
      for (i = ibeg; i &lt;= iend; i++){
        phi[j][i] = 0.25*(  phi0[j][i-1] + phi0[j][i+1]
                          + phi0[j-1][i] + phi0[j+1][i] );
  
        err += dx*dy*fabs(phi[j][i] - phi0[j][i]);
      }}
  
  
      #pragma acc parallel loop collapse(2) present(phi[:ny_tot][:nx_tot], phi0[:ny_tot][:nx_tot])
      for (j = jbeg; j &lt;= jend; j++){
      for (i = ibeg; i &lt;= iend; i++){
        phi0[j][i] = phi[j][i];
      }}
  
      #ifdef PARALLEL
      // double err_glob;
      #pragma acc host_data use_device(err, err_glob)
      {
      MPI_Allreduce (&amp;err, &amp;err_glob, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
      }
      err = err_glob;
      #endif
  
      // #pragma acc update host(err)
      if (rank == 0){
        printf (&quot;k = %d; err = %8.3e\n&quot;,k, err);
      }
      k++;
    }
    #pragma acc exit data copyout(phi[:ny_tot][:nx_tot], phi0[:ny_tot][:nx_tot], err, err_glob)
  
    WriteSolution (phi, nx, ny, &amp;mpi_decomp);
  
    #ifdef PARALLEL
    MPI_Finalize();
    #endif
    return 0;
  }
  
  #ifdef PARALLEL
  /* ********************************************************************* */
  void DomainDecomposition(MPI_Decomp *mpi_decomp)
  /*
   *
   *********************************************************************** */
  {
    int dim, i;
    int rank, size;
    int *coords  = mpi_decomp-&gt;coords;
    int *gsize   = mpi_decomp-&gt;gsize;
    int *lsize   = mpi_decomp-&gt;lsize;
    int *nprocs  = mpi_decomp-&gt;nprocs;
    int *periods = mpi_decomp-&gt;periods;
    int *procL   = mpi_decomp-&gt;procL;
    int *procR   = mpi_decomp-&gt;procR;
    int *start   = mpi_decomp-&gt;start;
    int new_coords[NDIM];
  
  /* --------------------------------------------------------
     1. Get rank &amp; size
     -------------------------------------------------------- */
  
    MPI_Comm_rank(MPI_COMM_WORLD, &amp;rank);
    MPI_Comm_size(MPI_COMM_WORLD, &amp;size);
  
    mpi_decomp-&gt;rank = rank;
    mpi_decomp-&gt;size = size;
  
  /* --------------------------------------------------------
     2. Obtain number of processor along each dimension.
        Use maximally squared decomp.
     -------------------------------------------------------- */
  
    nprocs[0] = (int)sqrt(size);
    nprocs[1] = size/nprocs[0];
    if (nprocs[0]*nprocs[1] != size){
      if (rank == 0) printf (&quot;! Cannot decompose\n&quot;);
      MPI_Finalize();
      exit(1);
    }
    if (rank == 0){
      printf (&quot;Decomposition achieved with %d X %d procs\n&quot;,nprocs[0],nprocs[1]);
    }
  
    periods[0] = 0;
    periods[1] = 0;
  
  /* --------------------------------------------------------
     3. Create Cartesian topology
     -------------------------------------------------------- */
  
    MPI_Cart_create(MPI_COMM_WORLD, NDIM, nprocs, periods,
                                          0, &amp;MPI_COMM_CART);
    MPI_Cart_get(MPI_COMM_CART, NDIM, nprocs, periods, coords);
  
  /* --------------------------------------------------------
     4. Fill structure members
     -------------------------------------------------------- */
  
    gsize[0] = NX_GLOB;
    gsize[1] = NY_GLOB;
    lsize[0] = NX_GLOB/nprocs[0];
    lsize[1] = NY_GLOB/nprocs[1];
    start[0] = coords[0]*lsize[0];
    start[1] = coords[1]*lsize[1];
  
  /* --------------------------------------------------------
     5. Determine ranks of neighbour processors
     -------------------------------------------------------- */
  
    for (dim = 0; dim &lt; NDIM; dim++) {
      for (i = 0; i &lt; NDIM; i++) new_coords[i] = coords[i];
  
      new_coords[dim] = coords[dim] + 1;
      if (new_coords[dim] &lt; nprocs[dim]) {
        MPI_Cart_rank ( MPI_COMM_CART, new_coords, &amp;(procR[dim]) );
      } else {
        procR[dim] = MPI_PROC_NULL;
      }
  
      new_coords[dim] = coords[dim] - 1;
      if (new_coords[dim] &gt;= 0) {
       MPI_Cart_rank ( MPI_COMM_CART, new_coords, &amp;(procL[dim]) );
      } else {
        procL[dim] = MPI_PROC_NULL;
      }
    }
  
  /* --------------------------------------------------------
     6. Print processor information.
        (Use MPI_Bcast() to print in sequence)
     -------------------------------------------------------- */
  
    int proc, go;
    for (proc = 0; proc &lt; size; proc++){
      go = proc;
      MPI_Bcast(&amp;go, 1, MPI_INT, 0, MPI_COMM_WORLD);
      if (rank == go) {
        printf (&quot;[Rank %d]\n&quot;,rank);
        printf (&quot;  coords = [%d, %d], lsize = [%d, %d]\n&quot;,
                   coords[0], coords[1], lsize[0], lsize[1]);
        for (dim = 0; dim &lt; NDIM; dim++){
          printf (&quot;  (procL, procR)[%d] = %d, %d\n&quot;, dim, procL[dim], procR[dim]);
        }
      }
      MPI_Barrier(MPI_COMM_WORLD);
    }
  
    return;
  }
  #endif
  
  /* ********************************************************************* */
  void BoundaryConditions(double **phi, double *x, double *y,
                          int nx, int ny, MPI_Decomp *mpi_decomp)
  /*
   *********************************************************************** */
  {
    int i,j;
    int ibeg   = NGHOST;
    int iend   = ibeg + nx - 1;
  
    int jbeg   = NGHOST;
    int jend   = jbeg + ny - 1;
  
    int *procL = mpi_decomp-&gt;procL;
    int *procR = mpi_decomp-&gt;procR;
  #ifdef PARALLEL
    int rank = mpi_decomp-&gt;rank;
    int size = mpi_decomp-&gt;size;
    double send_buf[NX_GLOB + 2*NGHOST];
    double recv_buf[NX_GLOB + 2*NGHOST];
  
  /*  Used for testing
      for (j = 0; j &lt;= jend+1; j++){
      for (i = 0; i &lt;= iend+1; i++){
        phi[j][i] = -1;
      }}
  
      for (j = jbeg; j &lt;= jend; j++){
      for (i = ibeg; i &lt;= iend; i++){
        phi[j][i] = rank;
      }}
  */
  
    #pragma acc enter data create(send_buf[:NX_GLOB+2*NGHOST], recv_buf[NX_GLOB+2*NGHOST])
  
    // Left buffer
    i = ibeg;
  
    #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], send_buf[:NX_GLOB+2*NGHOST])
    for (j = jbeg; j &lt;= jend; j++) send_buf[j] = phi[j][i];
    #pragma acc host_data use_device(send_buf, recv_buf)
    {
    MPI_Sendrecv (send_buf, jend+1, MPI_DOUBLE, procL[0], 0,
                  recv_buf, jend+1, MPI_DOUBLE, procL[0], 0,
                  MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    }
    #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], recv_buf[NX_GLOB+2*NGHOST])
    for (j = jbeg; j &lt;= jend; j++) phi[j][i-1] = recv_buf[j];
  
    // Right buffer
    i = iend;
    #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], send_buf[:NX_GLOB+2*NGHOST])
    for (j = jbeg; j &lt;= jend; j++) send_buf[j] = phi[j][i];
    #pragma acc host_data use_device(send_buf, recv_buf)
    {
    MPI_Sendrecv (send_buf, jend+1, MPI_DOUBLE, procR[0], 0,
                  recv_buf, jend+1, MPI_DOUBLE, procR[0], 0,
                  MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    }
    #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], recv_buf[NX_GLOB+2*NGHOST])
    for (j = jbeg; j &lt;= jend; j++) phi[j][i+1] = recv_buf[j];
  
  
    // Bottom buffer
    j = jbeg;
    #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], send_buf[:NX_GLOB+2*NGHOST])
    for (i = ibeg; i &lt;= iend; i++) send_buf[i] = phi[j][i];
    // #pragma acc update self(send_buf[:NX_GLOB+2*NGHOST])
    #pragma acc host_data use_device(send_buf, recv_buf)
    {
    MPI_Sendrecv (send_buf, iend+1, MPI_DOUBLE, procL[1], 0,
                    recv_buf, iend+1, MPI_DOUBLE, procL[1], 0,
                    MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    }
    #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], recv_buf[NX_GLOB+2*NGHOST])
    for (i = ibeg; i &lt;= iend; i++) phi[j-1][i] = recv_buf[i];
  
    // Top buffer
    j = jend;
    #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], send_buf[:NX_GLOB+2*NGHOST])
    for (i = ibeg; i &lt;= iend; i++) send_buf[i] = phi[j][i];
    #pragma acc host_data use_device(send_buf, recv_buf)
    {
    MPI_Sendrecv (send_buf, iend+1, MPI_DOUBLE, procR[1], 0,
                  recv_buf, iend+1, MPI_DOUBLE, procR[1], 0,
                  MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    }
    #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], recv_buf[NX_GLOB+2*NGHOST])
    for (i = ibeg; i &lt;= iend; i++) phi[j+1][i] = recv_buf[i];
  
    #pragma acc exit data copyout(send_buf[:NX_GLOB+2*NGHOST], recv_buf[NX_GLOB+2*NGHOST])
  
  #endif
  
  /* -- Left -- */
  
    if (procL[0] &lt; 0){
      i = ibeg-1;
      #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], y[:NY_GLOB+2*NGHOST])
      for (j = jbeg; j &lt;= jend; j++) phi[j][i] = 1.0-y[j];
    }
  
  /* -- Right -- */
  
    if (procR[0] &lt; 0){
      i = iend+1;
      #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], y[:NY_GLOB+2*NGHOST])
      for (j = jbeg; j &lt;= jend; j++) phi[j][i] = y[j]*y[j];
    }
  
  /* -- Bottom -- */
  
    if (procL[1] &lt; 0){
      j = jbeg-1;
      #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], x[:NX_GLOB+2*NGHOST])
      for (i = ibeg; i &lt;= iend; i++) phi[j][i] = 1.0-x[i];
    }
  
  /* -- Top -- */
  
    if (procR[1] &lt; 0){
      j = jend+1;
      #pragma acc parallel loop present(phi[:ny_tot][:nx_tot], x[:NX_GLOB+2*NGHOST])
      for (i = ibeg; i &lt;= iend; i++) phi[j][i] = x[i];
    }
  
    return;
  
  #ifdef PARALLEL
  // Print
    MPI_Barrier(MPI_COMM_WORLD);
    int go, proc;
    for (proc = 0; proc &lt; size; proc++){
      go = proc;
      MPI_Bcast(&amp;go, 1, MPI_INT, 0, MPI_COMM_WORLD);
  
      if (rank == go) {
        printf (&quot;Boundary [Rank %d]\n&quot;,rank);
        for (j = jend+1; j &gt;= 0; j--){
          for (i = 0; i &lt;= iend+1; i++){
            printf (&quot;%6.2f  &quot;, phi[j][i]);
          }
          printf (&quot;\n&quot;);
        }
      }
    }
  
  MPI_Finalize();
  exit(0);
  #endif
  }
  
  /* ********************************************************************* */
  void WriteSolution (double **phi, int nx, int ny, MPI_Decomp *md)
  /*
   *********************************************************************** */
  {
    int i,j;
    int ibeg   = NGHOST;
    int iend   = ibeg + nx - 1;
  
    int jbeg   = NGHOST;
    int jend   = jbeg + ny - 1;
  
    static int nfile = 0;
    char fname[32];
  
    sprintf (fname,&quot;laplace2D_MPIACC.txt&quot;,nfile);
  
  /*
  for (j = jbeg-1; j &lt;= jend+1; j++) for (i = ibeg-1; i &lt;= iend+1; i++) {
    phi[j][i] = -1;
  }
  for (j = jbeg; j &lt;= jend; j++) for (i = ibeg; i &lt;= iend; i++) {
    phi[j][i] = md-&gt;rank;
  }
  */
  #ifdef PARALLEL
    MPI_File fh;
    MPI_Datatype type_local, type_domain;
    int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY;
    int gsize[2], lsize[2], start[2];
  
  /* --------------------------------------------------------
     1. Create a local array type without the ghost zones
        This datatype will be passed to MPI_File_write()
     -------------------------------------------------------- */
  
    gsize[0] = md-&gt;lsize[0] + 2*NGHOST;
    gsize[1] = md-&gt;lsize[1] + 2*NGHOST;
  
    lsize[0] = md-&gt;lsize[0];
    lsize[1] = md-&gt;lsize[1];
  
    start[0] = NGHOST;
    start[1] = NGHOST;
  
    MPI_Type_create_subarray (NDIM, gsize, lsize, start,
                              MPI_ORDER_FORTRAN, MPI_DOUBLE, &amp;type_local);
    MPI_Type_commit (&amp;type_local);
  
  /* --------------------------------------------------------
     2. Create the subarry in the global domain.
        This datatype is used to set the file view.
     -------------------------------------------------------- */
  
    gsize[0] = NX_GLOB;
    gsize[1] = NY_GLOB;
  
    lsize[0] = md-&gt;lsize[0];
    lsize[1] = md-&gt;lsize[1];
  
    start[0] = lsize[0]*md-&gt;coords[0];  // equal to md-&gt;start[0]
    start[1] = lsize[1]*md-&gt;coords[1];  // equal to md-&gt;start[1]
  
    MPI_Type_create_subarray (NDIM, gsize, lsize, start,
                              MPI_ORDER_FORTRAN, MPI_DOUBLE, &amp;type_domain);
    MPI_Type_commit (&amp;type_domain);
  
  /* --------------------------------------------------------
     3. Write to disk
     -------------------------------------------------------- */
  
    MPI_File_delete(fname, MPI_INFO_NULL);
    MPI_File_open(MPI_COMM_CART, fname, amode, MPI_INFO_NULL, &amp;fh);
    MPI_File_set_view(fh, 0, MPI_DOUBLE, type_domain, &quot;native&quot;, MPI_INFO_NULL);
    MPI_File_write_all(fh, phi[0], 1, type_local, MPI_STATUS_IGNORE);
    MPI_File_close(&amp;fh);
    MPI_Type_free (&amp;type_local);
    MPI_Type_free (&amp;type_domain);
  #else
    FILE *fp;
    printf (&quot;&gt; Writing %s\n&quot;,fname);
    fp = fopen(fname, &quot;wb&quot;);
  
    for (j = jbeg; j &lt;= jend; j++){
      fwrite (phi[j] + ibeg, sizeof(double), nx, fp);
    }
  
    fclose(fp);
  #endif
  
    nfile++;
  }
  
  /* ********************************************************************* */
  double **Allocate_2DdblArray(int nx, int ny)
  /*
   * Allocate memory for a double precision array with
   * nx rows and ny columns
   *********************************************************************** */
  {
    int i,j;
    double **buf;
  
    buf    = (double **)malloc (nx*sizeof(double *));
    buf[0] = (double *) malloc (nx*ny*sizeof(double));
    for (j = 1; j &lt; nx; j++) buf[j] = buf[j-1] + ny;
  
    return buf;
  }
  /* ********************************************************************* */
  int **Allocate_2DintArray(int nx, int ny)
  /*
   * Allocate memory for an integer-type array with
   * nx rows and ny columns
   *********************************************************************** */
  {
    int i,j;
    int **buf;
  
    buf    = (int **)malloc (nx*sizeof(int *));
    buf[0] = (int *) malloc (nx*ny*sizeof(int));
    for (j = 1; j &lt; nx; j++) buf[j] = buf[j-1] + ny;
  
    return buf;
  }
  
  
  /* ********************************************************************* */
  void Show_2DdblArray(double **A, int nx, int ny, const char *string)
  /*
   *********************************************************************** */
  {
    int i, j;
  
    printf (&quot;%s\n&quot;,string);
    printf (&quot;------------------------------\n&quot;);
    for (i = 0; i &lt; nx; i++) {
      for (j = 0; j &lt; ny; j++) {
        printf (&quot;%8.2f  &quot;, A[i][j]);
      }
      printf (&quot;\n&quot;);
    }
    printf (&quot;------------------------------\n&quot;);
  }
  /* ********************************************************************* */
  void Show_2DintArray(int **A, int nx, int ny, const char *string)
  /*
   *********************************************************************** */
  {
    int i, j;
  
    printf (&quot;%s\n&quot;,string);
    for (j = 0; j &lt; ny; j++) printf (&quot;-----&quot;);
    printf (&quot;\n&quot;);
  
    for (i = 0; i &lt; nx; i++) {
      for (j = 0; j &lt; ny; j++) {
        printf (&quot;%03d  &quot;, A[i][j]);
      }
      printf (&quot;\n&quot;);
    }
  
    for (j = 0; j &lt; ny; j++) printf (&quot;-----&quot;);
    printf (&quot;\n&quot;);
  }
  </code></pre>