az_util.c

并行解法器,功能强大

/*====================================================================
 * ------------------------
 * | CVS File Information |
 * ------------------------
 *
 * $RCSfile: az_util.c,v $
 *
 * $Author: tuminaro $
 *
 * $Date: 2000/06/02 16:48:32 $
 *
 * $Revision: 1.65 $
 *
 * $Name:  $
 *====================================================================*/
#ifndef lint
static char rcsid[] = "$Id: az_util.c,v 1.65 2000/06/02 16:48:32 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 <string.h>
#include <math.h>
#include <float.h>
#include "az_aztec.h"

/*
 * File containing utility functions for solvers.  Note: Some of the fem high
 * level solvers such as the nonlinear solver call these routines as well.
 */

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

void AZ_compute_residual(double b[], double x[], double r[], 
			 int proc_config[], AZ_MATRIX *Amat)

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

  Compute the residual r = b - Ax.

  Author:          John N. Shadid, SNL, 1421
  =======

  Return code:     void
  ============

  Parameter list:
  ===============

  val:             Array containing the nonzero entries of the matrix (see
                   Aztec User's Guide).

  indx,
  bindx,
  rpntr,
  cpntr,
  bpntr:           Arrays used for DMSR and DVBR sparse matrix storage (see
                   file Aztec User's Guide).

  b:               Right hand side of linear system.

  x:               On input, contains the initial guess. On output contains the
                   solution to the linear system.

  r:               On output, residual vector.

  options:         Determines specific solution method and other parameters.

  proc_config:     Machine configuration.  proc_config[AZ_node] is the node
                   number.  proc_config[AZ_N_procs] is the number of processors.

  params:          Drop tolerance and convergence tolerance info.

  Amat:            Structure used to represent the matrix (see az_aztec.h
                   and Aztec User's Guide).



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

{

  /* local variables */

  register int i;
  int N;

  /**************************** execution begins ******************************/

  N = Amat->data_org[AZ_N_internal] + Amat->data_org[AZ_N_border];

  Amat->matvec(x, r, Amat, proc_config);


  for(i = 0; i < N; i++) r[i] = b[i] - r[i];

} /* AZ_compute_residual */

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

double AZ_gmax_vec(int N, double vec[], int proc_config[])

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

  Routine to return the maximum element of a distributed vector "vec".

  Author:          John N. Shadid, SNL, 1421
  =======

  Return code:     double, maximum value in vector 'vec'
  ============

  Parameter list:
  ===============

  N:               Length of vector 'vec'.

  vec:             Vector of length 'N'.

  proc_config:     Machine configuration.  proc_config[AZ_node] is the node
                   number.  proc_config[AZ_N_procs] is the number of processors.

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

{

  /* local variables */

  register int i;
  double       rmax = 0.0;

  /**************************** execution begins ******************************/

  for (i = 0; i < N; i++) rmax = max(rmax, vec[i]);
  rmax = AZ_gmax_double(rmax, proc_config);

  return rmax;

} /* AZ_gmax_vec */

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

double AZ_gdot(int N, double r[], double z[], int proc_config[])

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

  Routine to perform dot product of r and z with unit stride. This routine call
  the BLAS routine ddot to do the local vector dot product and then uses the
  global summation routine AZ_gsum_double to obtain the reguired global result.

  Author:          John N. Shadid, SNL, 1421
  =======

  Return code:     double, dot product of vectors 'r' and 'z'
  ============

  Parameter list:
  ===============

  N:               Length of vector 'vec'.

  r, z:            Vectors of length 'N'.

  proc_config:     Machine configuration.  proc_config[AZ_node] is the node
                   number.  proc_config[AZ_N_Procs] is the number of processors.

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

{

  static int one = 1;
  int        add_N;

  add_N = N;

  return AZ_gsum_double(ddot_(&add_N, r, &one, z, &one), proc_config);

} /* dot */

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

#ifdef TIME_VB

void AZ_time_kernals(int gN, double gnzeros, double val[], int indx[],
                     int bindx[], int rpntr[], int cpntr[], int bpntr[],
                     double x[], double y[], int ntimes, int options[],
                     int data_org[], int proc_config[], double params[],
                     AZ_MATRIX *Amat)

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

  Solve the system of equations given in the VBR format using an iterative
  method specified by 'options[AZ_solver]'. Store the result in 'x'.

  Author:          John N. Shadid, SNL, 1421
  =======

  Return code:     void
  ============

  Parameter list:
  ===============

  gN:              Global order of the linear system of equations.

  gnzeros:         Global number of nonzeros in val.

  val:             Array containing the nonzero entries of the matrix (see
                   Aztec User's Guide).

  indx,
  bindx,
  rpntr,
  cpntr,
  bpntr:           Arrays used for DMSR and DVBR sparse matrix storage (see
                   file Aztec User's Guide).

  x, y:            Vectors of length gN.

  ntimes:          Number of times to perform each operation.

  options:         Determines specific solution method and other parameters.

  data_org:        Array containing information on the distribution of the
                   matrix to this processor as well as communication parameters
                   (see Aztec User's Guide).

  proc_config:     Machine configuration.  proc_config[AZ_node] is the node
                   number.  proc_config[AZ_N_procs] is the number of processors.

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

{

  /* local variables */

  register int i, j;
  double       start_t, total_t;
  double       sprx_comm_t, sprx_comp_t, dot_comm_t, dot_comp_t;
  double       sprx_t, sprx_overlap_cop_t;
  double       daxpy_comp_t, dzpax_comp_t;
  double       Djunk;
  double      *ptr_vec1, *ptr_vec2, *ptr_vec3, *z1, *z2, *z3, alpha = 1.0;
  int          one = 1, bpntr_index;
  double       Mflop, Mflops_comp, Mflops_node_comp;
  double       Mflops_comm, Mflops_node_comm;
  double       sparax_overlap_border_comp_t, sparax_overlap_internal_comp_t;
  double       read_nonoverlap_t, read_nonoverlap_t_max, read_t_max;
  double       time_overlap_max;
  double       gather_t, write_t, read_t, time;
  double       max, min, avg;
  double       overall_t, start_overall_t;

  int          Num_Proc, Proc;
  int          Num_Internal_Blks;
  int          iout;

  char        *message_recv_add[AZ_MAX_NEIGHBORS];
  char        *message_send_add[AZ_MAX_NEIGHBORS];
  int          message_recv_length[AZ_MAX_NEIGHBORS];
  int          message_send_length[AZ_MAX_NEIGHBORS];
  AZ_MATRIX   Amat2;

  /*
    message_send_add:
    message_send_add[i] points to the beginning of the list of
    values to be sent to the ith neighbor (i.e. data_org[AZ_neighbors+i] or
    sometimes locally defined as proc_num_neighbor[i]). That is,
    *(message_send_add[i] + j) is the jth value to be sent to the
    ith neighbor.

    message_send_length:
    message_send_length[i] is the number of bytes to be sent to
    the ith neighbor (i.e. data_org[AZ_neighbors+i] or sometimes
    locally defined as proc_num_neighbor[i]).

    message_recv_add:
    message_recv_add[i] points to the beginning of the list of
    locations which are to receive values sent by the ith neighbor (i.e.
    data_org[AZ_neighbors+i] or sometimes locally defined as
    proc_num_neighbor[i]). That is, *(message_recv_add[i] + j) is the
    location where the jth value sent from the ith neighbor will be stored.
    message_recv_length:
    message_recv_length[i] is the number of bytes to be sent to
    the ith neighbor (i.e. data_org[AZ_neighbors+i] or sometimes
    locally defined as proc_num_neighbor[i]).
    */

  int          temp1, temp2;

  /**************************** execution begins ******************************/

  Proc              = proc_config[AZ_node];
  Num_Proc          = proc_config[AZ_N_procs];
  Num_Internal_Blks = data_org[AZ_N_int_blk];
  iout              = options[AZ_output];

  if (data_org[AZ_matrix_type] == AZ_USER_MATRIX){
      if (az_proc == 0) {
        (void) fprintf(stderr, "ERROR: matrix-vector timing not available for"
                       "matrix-free.\n"
                       "      (data_org[AZ_matrix_type] = %d)\n\n",
                       data_org[AZ_matrix_type]);
      }
      exit (-1);
  }

  if (data_org[AZ_matrix_type] != AZ_VBR_MATRIX) {
    (void) fprintf(stderr, "I'm not sure if the timing stuff works for \n"
                   "nonVBR matrices. For example I don't think that\n"
                   "we have overlapped communication/computations.\n"
                   "However, probably a lot of the stuff works?\n");
    exit(-1);
  }

  if (ntimes < 0)   ntimes = -ntimes;
  if (ntimes < 200) ntimes = 200;

  if (Proc == 0) {
    (void) printf("++++++++++ timing run ++++++++++\n");
    (void) printf("\nnprocs: %d ntimes: %d\n", Num_Proc, ntimes);
    (void) printf("N: %d\tNzeros: %e\t\n\n", gN, gnzeros);
  }

  /* sparax */

  if (iout > 0 && Proc == 0) {
    (void) printf("timing nonoverlapped matrix vector multiply\n");
  }

  start_t = AZ_second();
  for (i = 1; i <= ntimes; i++)
    Amat->matvec(x, r, Amat, proc_config);


  sprx_t = AZ_second() - start_t;

  /* exchange boundary */

  if (iout > 0 && Proc ==0){
    (void) printf("timing unstructured communication\n");
  }

  start_t = AZ_second();
  for (i = 1; i <= ntimes; i++) {
    AZ_exchange_bdry(x, data_org, proc_config);
    AZ_exchange_bdry(y, data_org, proc_config);
  }

  sprx_comm_t = 0.5*(AZ_second() - start_t);
  sprx_comp_t = sprx_t - sprx_comm_t;

  /* overlapped communication */

  if (iout > 0 && Proc == 0) {
    (void) printf("timing overlapped sparse matrix vector multiply\n");
  }
  start_t = AZ_second();
  for (i = 1; i <= ntimes; i++)
    Amat->matvec(x, r, Amat, proc_config);

  time = AZ_second() - start_t;

  gather_t                       = (double) 0.0;
  write_t                        = (double) 0.0;
  read_t                         = (double) 0.0;
  read_nonoverlap_t              = (double) 0.0;
  sparax_overlap_internal_comp_t = (double) 0.0;
  sparax_overlap_border_comp_t   = (double) 0.0;

  if (iout > 0 && Proc == 0) {
    (void) printf("time the individual routines for sparax_overlap\n");
  }

  start_overall_t = AZ_second();

  for (i = 1; i <= ntimes; i++) {

    /* time the individual routines for sparax_overlap */

    start_t = AZ_second();
    AZ_gather_mesg_info(x, data_org, message_recv_add, message_send_add,
                        message_recv_length, message_send_length);
    gather_t += AZ_second() - start_t;

    start_t = AZ_second();
    AZ_write_local_info(data_org, message_recv_add, message_send_add,