az_app_utils.c

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

/*====================================================================
 * ------------------------
 * | CVS File Information |
 * ------------------------
 *
 * $RCSfile: az_app_utils.c,v $
 *
 * $Author: tuminaro $
 *
 * $Date: 1999/10/22 16:08:02 $
 *
 * $Revision: 1.20 $
 *
 * $Name:  $
 *====================================================================*/
#ifndef lint
static char rcsid[] = "$Id: az_app_utils.c,v 1.20 1999/10/22 16:08:02 tuminaro Exp $";
#endif


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


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

/* ----------------- external function declarations ------------------- */

extern void create_msr_matrix(int *, double **, int **, int );
extern void AZ_transform(int *, int **, int *, double *, int *, int **,int **,
                         int **, int, int *, int *, int *, int **, int );
extern void create_fe_matrix(int *, int , int **, double **, int );
extern void create_vbr_matrix(int *, double **, int **, int, 
                              int **, int **, int **);
extern void init_options(int *, double *);
extern void init_guess_and_rhs(int *, int *, double **,double **,int *, 
                               double *, int *, int *, int *, int *, 
                               int *, int *);
extern void init_matrix_vector_structures(int *, int **, int **, int  **, 
                               int **, int **, int , double **, int **, 
                               int **, int **, int **, int **);
extern void print_global_element(int ,int *,int *,int *, int *, double *,int *);

extern int num_PDE_eqns;   /* number of PDEs being solved         */
extern int application;    /* problem being solved                */
extern int N_grid_pts;     /* N_grid_pts is the total number of   */
                           /* grid points used in the simulation. */

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

void init_options(int options[], double params[])

{

  /*
   * Choose among AZTEC options (see User's Guide).
   */

  AZ_defaults(options, params);

  options[AZ_solver]   = AZ_cgs;
  options[AZ_scaling]  = AZ_none;
  options[AZ_precond]  = AZ_ls;
  options[AZ_conv]     = AZ_r0;
  options[AZ_output]   = 1;
  options[AZ_pre_calc] = AZ_calc;
  options[AZ_max_iter] = 1550;
  options[AZ_poly_ord] = 5;
  options[AZ_overlap]  = AZ_none;
  options[AZ_kspace]   = 60;
  options[AZ_aux_vec]  = AZ_resid;

  params[AZ_tol]       = 4.00e-9;
  params[AZ_drop]      = 0.0;
  params[AZ_ilut_fill] = 1.5;
  params[AZ_omega]     = 1.;

} /* init_options */

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

void init_guess_and_rhs(int update_index[], int update[], double *x[],double
                        *ax[],int data_org[], double val[], int indx[], int
                        bindx[], int rpntr[], int cpntr[], int bpntr[], int
                        proc_config[])

/*
 * Set the initial guess and the right hand side where the right hand side
 * is obtained by doing a matrix-vector multiplication.
 *
 * Author: Ray Tuminaro, Div 1422, SNL
 * Date :  3/15/95
 *
 * Parameters
 *
 *    update_index   ==      On input, ordering of update and external
 *                           locally on this processor. For example
 *                           'update_index[i]' gives the index location
 *                           of the block which has the global index
 *                           'update[i]'.
 *    update         ==      On input, list of pts to be updated on this node
 *    data_org       ==      On input, indicates how data is set on this node.
 *                           For example, data_org[] contains information on
 *                           how many unknowns are internal, external and
 *                           border unknowns as well as which points need
 *                           to be communicated. See User's Guide for more
 *                           details.
 *    val, indx,     ==      On input, holds matrix nonzeros. See User's Guide
 *    bindx, rpntr,          for more details.
 *    cpntr, bpntr
 *    x              ==      On output, 'x' is allocated and set to all zeros.
 *    ax             ==      On output, 'ax' is allocated and is set to the
 *                           result of a matrix-vector product.
 */

{

  int i,j;
  int temp,num;
  double sum = 0.0;
  AZ_MATRIX *Amat;

  temp = data_org[AZ_N_int_blk]  + data_org[AZ_N_bord_blk];
  num  = data_org[AZ_N_internal] + data_org[AZ_N_border];

  /* allocate vectors */

  i       = num + data_org[AZ_N_external];
  *x      = (double *) AZ_allocate((i+1)*sizeof(double));
  *ax     = (double *) AZ_allocate((i+1)*sizeof(double));
  if (*ax == NULL) {
    (void) fprintf(stderr, "Not enough space in init_guess_and_rhs() for ax\n");
    exit(1);
  }
  for (j = 0 ; j < i ; j++ ) (*x)[j] = 0.0;
  for (j = 0 ; j < i ; j++ ) (*ax)[j] = 0.0;

  /* initialize 'x' to a function which will be used in matrix-vector product */

  if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) {
    for (i = 0; i < temp; i++) {
      for (j = rpntr[i]; j < rpntr[i+1]; j++) {
        (*x)[j] = (double) (update[i]) + (double)(j-rpntr[i]) /
          (double)(num_PDE_eqns);
      }
    }
  }
  else {
    for (i = 0; i < temp; i++) {
      (*x)[i] = (double) (update[i]) / (double) (num_PDE_eqns);
    }
  }

  /* Reorder 'x' so that it conforms to the transformed matrix */
 
  AZ_reorder_vec(*x,data_org,update_index,rpntr);

  if (application == 2) {

    /* take out the constant vector. Used for the */
    /* finite element problem because it is singular */

    sum = AZ_gsum_double(sum, proc_config);
    i   = AZ_gsum_int(num, proc_config);
    if (i != 0) sum = sum / ((double) i);
    for (i = 0; i < num; i++) (*x)[i] -= sum;
  }
  Amat = AZ_matrix_create(num);

  if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX)
     AZ_set_MSR(Amat, bindx, val, data_org, 0, NULL, AZ_LOCAL);
  else if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX)
     AZ_set_VBR(Amat, rpntr,cpntr, bpntr, indx,bindx, val, data_org, 0, NULL,AZ_LOCAL);


  Amat->matvec(*x, *ax, Amat, proc_config);
  AZ_matrix_destroy( &Amat );

  for (i = 0; i < num; i++) (*x)[i] = 0.0;

} /* init_guess_and_rhs */

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

void init_matrix_vector_structures(int proc_config[], int *update_index[], int
                                   *update[], int  *data_org[], int *external[],
                                   int *extern_index[], int input_option, double
                                   *val[], int *bindx[], int *indx[], int
                                   *bpntr[], int *rpntr[], int *cpntr[])