az_mat_free_main.c

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

/*====================================================================
 * ------------------------
 * | CVS File Information |
 * ------------------------
 *
 * $RCSfile: az_mat_free_main.c,v $
 *
 * $Author: tuminaro $
 *
 * $Date: 2000/06/02 16:46:55 $
 *
 * $Revision: 1.13 $
 *
 * $Name:  $
 *====================================================================*/
#ifndef lint
static char rcsid[] = "$Id: az_mat_free_main.c,v 1.13 2000/06/02 16:46:55 tuminaro Exp $"
;
#endif

/*******************************************************************************
 * Copyright 1998, 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.                                                                 *
 ******************************************************************************/

/*******************************************************************************
 * Sample matrix-free driver for AZTEC package.  The software is tested by 
 * solving a simple system (Poisson equation) via Aztec's iterative solvers.
 ******************************************************************************/

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

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

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

/*
 * Set up and solve a linear system of equations using Aztec's 
 * matrix-free interface. The matrix system corresponds to 
 * Poisson's equation on a regular 2-dimensional mesh of size
 * n x n. The associated matrix has n^2 equations (one for each 
 * mesh grid point). The kth equation corresponds to grid point 
 * (i,j) where    
 *            k = n*j + i       and  i = 0, ... , n-1
 *                                   j = 0, ... , n-1
 * and the equation is of the form
 *
 *    4*u(i,j) - u(i-1,j) - u(i+1,j) - u(i,j-1) - u(i,j+1) = f(i,j)
 *
 * where u(i,j) corresponds to the unknown located at the (i,j)th 
 * grid point and f(i,j) is the right-hand side vector corresponding
 * to the (i,j)th grid point.
 * 
 * For example a 4x4 grid would correspond to a 16 equation matrix
 * with the following grid point numbering:
 *
 *                     12 13 14 15 
 *                      8  9 10 11 
 *                      4  5  6  7 
 *                      0  1  2  3 
 * 
 * The equations are solved using p^2 processors where the processors
 * are distributed as a 2-D array. For example with 4 processors, the
 * above equations would be distributed as follows:
 *      processor 0 would contain equations  0, 1, 4, 5
 *      processor 1 would contain equations  2, 3, 6, 7
 *      processor 2 would contain equations  8, 9,12,13
 *      processor 3 would contain equations 10,11,14,15
 */

{
/*****************************************************************************/
/* Generic Aztec variables common to most matrix-free applications           */
/*****************************************************************************/
                     /*                                                      */
int                  /*                                                      */
   proc_config[      /* Processor information.                               */
     AZ_PROC_SIZE],  /*                                                      */
   options[          /* Used to pass parameters to Aztec.                    */
     AZ_OPTIONS_SIZE];/*                                                     */
double               /*                                                      */
   params[           /* Used to pass in parameters to Aztec                  */
     AZ_PARAMS_SIZE],/*                                                      */
   status[           /* Returned from AZ_iterate() indicating success or     */
     AZ_STATUS_SIZE];/* failure of the iterative solvers                     */ 
AZ_MATRIX *Amat=NULL,/* Structure representing matrix to be solved.          */
          *Pmat=NULL;/* Structure representing preconditioning matrix.       */
                     /*                                                      */
AZ_PRECOND *Prec;    /* Structure representing preconditioner.               */
                     /*                                                      */
double *x, *rhs;     /* Initial guess and right-hand side to linear system.  */

/*****************************************************************************/
/* Example specific variables                                                */
/*****************************************************************************/
int N_equations;     /* Number of equations in the linear system residing on */
                     /* this processor. In our example, this is 'nx*ny'.     */
                     /*                                                      */
int N_coord;         /* Total number of points in each coordinate direction. */
                     /* That is, the matrix dimension is of size N_coord^2   */
                     /*                                                      */
int  Nproc_coord;    /* Number of processor in each coordinate direction     */
                     /* That is, the linear system is solved by a processor  */
                     /* array of size Nproc_coord x Nproc_coord.             */
                     /*                                                      */
int nx, ny;          /* Number of grid points in each coordinate direction   */
                     /* within the local grid owned by this processor.       */
                     /*                                                      */
struct pass_data     /* Used to pass data through Aztec into the USER's      */
       pass_data;    /* matrix-vector product routine.                       */
struct exchange      /* Used to exchange data within the USER's matrix-vector*/
      comm_structure;/* product routine.                                     */

int i, N_ghost, precond_choice;
/* ----------------------- execution begins --------------------------------*/

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

#ifdef AZ_MPI
  MPI_Init(&argc,&argv);
   AZ_set_proc_config(proc_config, MPI_COMM_WORLD);
#else
   AZ_set_proc_config(proc_config, AZ_NOT_MPI);
#endif

  /* initialize AZTEC options */

   AZ_defaults(options, params);


   precond_choice = example_specific_read_input(options, proc_config,
                                                &N_coord,&Nproc_coord);

  /* Set up the example specific data structure 'comm_structure' which  */
  /* handles communication. In addition, set (nx, ny) to the number of  */
  /* grid points in the x and y direction for this processor.           */
  
   example_specific_comm_setup(N_coord, &nx, &ny, proc_config[AZ_node], 
                               Nproc_coord, &comm_structure);
   N_equations = nx*ny;

  /* Set up an EXAMPLE SPECIFIC data structure to pass information      */
  /* through Aztec to the USER's routines. NOTE: while specific to      */
  /* this example, most codes will do something similar.                */

   pass_data.nx    = nx;                      
   pass_data.ny    = ny;                      
   pass_data.nproc = (int) sqrt( ((double) proc_config[AZ_N_procs]) + .01);
   pass_data.py    = proc_config[AZ_node]/pass_data.nproc;
   pass_data.px    = proc_config[AZ_node]%pass_data.nproc;
   pass_data.comm_structure  = &comm_structure;

  /* initialize Aztec matrix to be solved */

   Amat = AZ_matrix_create(N_equations);
   AZ_set_MATFREE(Amat, &pass_data, example_specific_matvec);

   Prec = NULL;    /* When preconditioning structure is NULL, AZ_iterate()  */
                   /* applies Aztec's preconditioners to application matrix */
                   /* (i.e. user does not supply a preconditioning routine  */
                   /* or an additional matrix for preconditioning.          */

   switch (options[AZ_precond]) {

   
   case AZ_none:           /***** NO preconditioning. ***********************/
     break;

   case AZ_ls:             /***** Polynomial preconditioning (least-squares */
   case AZ_Neumann:        /***** or Neumann). ******************************/
     AZ_set_MATFREE_matrix_norm(Amat, 8.0); /* Aztec needs upper bound for  */
     break;                                 /* matrix norm                  */

   case AZ_Jacobi:         /***** Jacobi preconditioning. *******************/
      /*                                                                    */
      /* Matrix diagonal is needed and can be supplied via:                 */
      /*     1) a getrow() function associated with Amat                    */
      /*     2) a new matrix (MSR, VBR, or USER with getrow function)       */
      /*        associated with preconditioner.                             */
      /* We choose the second option in this example.                       */


      Pmat = AZ_matrix_create(N_equations);
      AZ_set_MATFREE_getrow(Pmat, NULL, example_specific_diagonal_getrow, 
                            NULL, 0, proc_config);
      Prec = AZ_precond_create(Pmat, AZ_precondition, NULL);
      break;

   case AZ_dom_decomp:     /***** Domain Decomposition preconditioning. *****/
      /*                                                                    */
      /* Local matrix must be supplied via a getrow associated with Amat or */
      /* via a new matrix associated with preconditioner.                   */

      if (precond_choice == NONOVERLAPDD_PRECOND) {

         /* Since no overlapping is done, communication information and     */
         /* matrix data corresponding to ghost unknowns is not needed.      */

         Pmat = AZ_matrix_create(N_equations);
         AZ_set_MATFREE_getrow(Pmat, &pass_data,simple_example_specific_getrow,