az_solve.c

并行解法器,功能强大

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


/*
 * All of these routines form the interface between the code drivers and the
 * user's choice of algorithm and PDE problem.
 */

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

#define AZ_oldprecond         0
#define AZ_oldsubdomain_solve 1
#define AZ_oldreorder         2
#define AZ_oldpre_calc        3
#define AZ_SAVE_SIZE          5

void AZ_iterate(double x[], double b[], int options[], double params[],
                double status[], int proc_config[], AZ_MATRIX *Amat,
                AZ_PRECOND *precond, struct AZ_SCALING *scaling)


/*******************************************************************************
This is the new Aztec interface. This routine calls AZ_oldsolve() passing
in for example Amat->indx for the indx[] parameter in AZ_oldsolve().

NOTE: User's can still invoke AZ_solve() in the old Aztec way. AZ_solve
      also calls AZ_oldsolve(). However, matrix-free and coarse grid
       capabilities are not available via AZ_solve().


  Author:          Ray Tuminaro, SNL, 1422
  =======

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

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


  x:               Contains the result vector x.

  b:               Contains the vector b.

  options:         Determines specific solution method and other parameters.

  params:          Drop tolerance and convergence tolerance info.

  status:          On output, indicates termination status:
                    0:  terminated normally.
                   -1:  maximum number of iterations taken without achieving
                        convergence.
                   -2:  Breakdown. The algorithm can not proceed due to
                        numerical difficulties (usually a divide by zero).
                   -3:  Internal residual differs from the computed residual due
                        to a significant loss of precision.

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

  Amat:            Structure used to represent the matrix (see az_aztec.h
                   and 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.

  precond:         Structure used to represent the preconditionner
                   (see az_aztec.h ad Aztec User's Guide).

  scaling:         Structure used to represent the  scaling
                   (see az_aztec.h ad Aztec User's Guide).

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

{
  double  total_time, start_t;
  int     prec_allocated = 0;
  struct AZ_SCALING *scale2;

  if (scaling == NULL) 
     scale2 = AZ_scaling_create();
  else scale2 = scaling;

  AZ__MPI_comm_space_ok();
  if (Amat->mat_create_called != 1) {
    if (proc_config[AZ_node] == 0) {
       printf("AZ_iterate: AZ_matrix_create(int) should be called to\n");
       printf("            create matrix object (Amat) to be solved!\n");
    } 
    exit(1);
  }
  if (precond == NULL) {
     if (options[AZ_precond] == AZ_user_precond) {
       if (proc_config[AZ_node] == 0) {
          printf("AZ_iterate: Can not use NULL for precond argument when\n");
          printf("            options[AZ_precond] == AZ_user_precond.\n");
       }
       exit(1);
     }
     precond = AZ_precond_create(Amat, AZ_precondition, NULL);
     prec_allocated = 1;
  }
  if (precond->prec_create_called != 1) {
    if (proc_config[AZ_node] == 0) {
       printf("AZ_iterate: AZ_precond_create should be called to\n   ");
       printf("       create preconditioning object!\n");
    } 
    exit(1);
  }
  if (precond->Pmat->mat_create_called != 1) {
    if (proc_config[AZ_node] == 0) {
       printf("AZ_iterate: AZ_matrix_create(int) should be called to\n   ");
       printf("       create preconditioning matrix object (precond->Pmat)!\n");
    } 
    exit(1);
  }
  if (Amat->matvec == NULL) {
    if (proc_config[AZ_node] == 0) {
       printf("AZ_iterate: Matrix vector product needs to be set via ");
       printf("AZ_set_MSR(...),\n             AZ_set_VBR(...), or ");
       printf("AZ_set_MATFREE(...).\n");
    } 
    exit(1);
  }
  AZ_iterate_setup(options, params, proc_config, Amat, precond);
  AZ_sync(proc_config);
  start_t = AZ_second();
  AZ_oldsolve(x, b, options, params, status, proc_config, Amat, precond,scale2);
  total_time = AZ_gmax_double(AZ_second() - start_t, proc_config);

  status[AZ_solve_time] = total_time;
  if ((options[AZ_output] != AZ_none) && (options[AZ_output] != AZ_warnings)) {
    if (proc_config[AZ_node] == 0) {
      (void) printf("\n\n\t\tSolution time: %f (sec.)\n", total_time);
      (void) printf("\t\ttotal iterations: %d\n", (int) status[AZ_its]);
    }
  }

  AZ_flop_rates(Amat->data_org,Amat->indx,Amat->bpntr, Amat->bindx,
	      options, status, total_time, proc_config);

  AZ_iterate_finish(options, Amat, precond);
  if (prec_allocated)  AZ_precond_destroy(&precond);
  if (scaling == NULL) AZ_scaling_destroy(&scale2);
}

/* AZ_iterate*/


void AZ_solve(double x[], double b[], int options[], double params[],
              int indx[], int bindx[], int rpntr[], int cpntr[], int bpntr[],
              double val[], int data_org[], double status[], int proc_config[])

/*******************************************************************************
  In order to assure backward compatibility  Aztec's previous AZ_solve()
  has been  renamed to AZ_oldsolve(), and this routine defines the 3 new
  parameters used by  AZ_oldsolve and called AZ_oldsolve.

  Author:          Ray Tuminaro, SNL, 1422
  =======

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

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

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

  b:               Right hand side of linear system.

  options:         Determines specific solution method and other parameters.

  params:          Drop tolerance and convergence tolerance info.

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

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

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

  status:          On output, indicates termination status:
                    0:  terminated normally.
                   -1:  maximum number of iterations taken without achieving
                        convergence.
                   -2:  Breakdown. The algorithm can not proceed due to
                        numerical difficulties (usually a divide by zero).
                   -3:  Internal residual differs from the computed residual due
                        to a significant loss of precision.

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

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


{
  AZ_MATRIX *Amat;
  AZ_PRECOND *precond;
  double  total_time, start_t;
  struct AZ_SCALING *scale2;

   scale2 = AZ_scaling_create();
   AZ__MPI_comm_space_ok();
   Amat    = AZ_matrix_create(data_org[AZ_N_internal]+data_org[AZ_N_border]);
   precond = AZ_precond_create(Amat, AZ_precondition, NULL);

   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);
   else {
      fprintf(stderr,"Unknown matrix type (%d)\n",data_org[AZ_matrix_type]);
      fprintf(stderr,"Matrix-free is now available via AZ_iterate()\n");
      exit(1);
   }

   if (options[AZ_precond] == AZ_user_precond) {
      fprintf(stderr,"Unknown preconditioning options[AZ_precond] =  (%d)\n",options[AZ_precond]);
      fprintf(stderr,"User preconditioning is now available via AZ_iterate()\n");
      exit(1);
   }

   options[AZ_recursion_level] = 0;

  if (options[AZ_pre_calc] != AZ_reuse) {
     (void) AZ_manage_memory(0,AZ_EVERYBODY_BUT_CLEAR,(Amat->data_org)[AZ_name],"kvecs",(int *) 0);
  }
  (void) AZ_manage_memory(0, AZ_CLEAR, AZ_SYS, (char *) 0, (int *) 0);

  /* output solver, scaling, and preconditioning options */

  AZ_print_call_iter_solve(options, params, proc_config[AZ_node], 0, precond);


  AZ_sync(proc_config);
  start_t = AZ_second();

  AZ_oldsolve(x, b, options, params, status, proc_config, Amat, precond,scale2);

  total_time = AZ_gmax_double(AZ_second() - start_t, proc_config);

  status[AZ_solve_time] = total_time;
  if ((options[AZ_output] != AZ_none) && (options[AZ_output] != AZ_warnings)) {
    if (proc_config[AZ_node] == 0) {
      (void) printf("\n\n\t\tSolution time: %f (sec.)\n", total_time);
      (void) printf("\t\ttotal iterations: %d\n", (int) status[AZ_its]);
    }
  }

  AZ_flop_rates(data_org,indx,bpntr, bindx, options, status, total_time, 
                proc_config);

  if (options[AZ_keep_info] == 0)
      (void) AZ_manage_memory(0,AZ_CLEAR,(Amat->data_org)[AZ_name],(char *) 0,
                              (int *) 0);

  (void) AZ_manage_memory(0, AZ_CLEAR, AZ_SYS, (char *) 0, (int *) 0);

  AZ_precond_destroy(&precond);
  AZ_matrix_destroy(&Amat);
  AZ_scaling_destroy(&scale2);

}
/* AZ_solve*/

#include <string.h>
#ifdef USING_ML
#include "ml_include.h"
#endif

void AZ_oldsolve(double x[], double b[], int options[], double params[], 
	double status[], int proc_config[], AZ_MATRIX *Amat, 
	AZ_PRECOND *precond, struct AZ_SCALING *scaling)


/*******************************************************************************
  Aztec's previous AZ_solve() is renamed to AZ_oldsolve() with 3 new
  parameters appended to it: Amat, precond, scaling.
  This routine is never called directly by an application. It is only
  used internally by Aztec.

  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:          Ray Tuminaro, SNL, 1422
  =======

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

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

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

  b:               Right hand side of linear system.

  options:         Determines specific solution method and other parameters.

  params:          Drop tolerance and convergence tolerance info.

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

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

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

  status:          On output, indicates termination status:
                    0:  terminated normally.
                   -1:  maximum number of iterations taken without achieving
                        convergence.
                   -2:  Breakdown. The algorithm can not proceed due to
                        numerical difficulties (usually a divide by zero).
                   -3:  Internal residual differs from the computed residual due
                        to a significant loss of precision.

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

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

  precond:         Structure used to represent the preconditionner
                   (see az_aztec.h ad Aztec User's Guide).

  scaling:         Structure used to represent the  scaling
                   (see az_aztec.h ad Aztec User's Guide).


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

{

  /* local variables */

  int     i, itemp2, itemp3, j;
  int     nonzeros;
  int     N, N_Blk;
  int     *data_org;
  struct  context context;
  struct aztec_choices aztec_choices;
  int    nz_used;
  int bandwidth, extra_factor_nonzeros;
  AZ_MATRIX Amat2;
  int *bindx2, N_nz_factors, t1, NNN;
  double *val2, *newparams, largest, dtemp;
  char tag[80];
  int   save_old_values[7], changed = 0;
  char  tstr[15];
  struct AZ_CONVERGE_STRUCT *conv_info;
  int size1, largest_index;
  double global_largest;
#ifdef ML
  int size2, *ibuf, allocated, row_length, *already_printed, *ibuf2, row2_length;
  double *tempv, *tempy, *tempz, *dbuf;
  char boundary;
  ML *ml;
#endif


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

  newparams = params;
  conv_info = AZ_converge_create();
  data_org  = Amat->data_org;
  conv_info->scaling = scaling;
  AZ__MPI_comm_space_ok();
  status[AZ_Aztec_version] = -1.;
  save_old_values[6] = options[AZ_ignore_scaling];
  if ( options[AZ_conv] == AZ_expected_values) {
      options[AZ_ignore_scaling] = AZ_TRUE;
      NNN = data_org[AZ_N_internal] + data_org[AZ_N_border];
      sprintf(tag, "some weights %d %d %d", data_org[AZ_name], 
              options[AZ_recursion_level],NNN);
      newparams = (double *) AZ_manage_memory((AZ_PARAMS_SIZE+NNN)*
					      sizeof(double), AZ_ALLOC,
                                               data_org[AZ_name],tag,&i);
      if ((options[AZ_pre_calc] == AZ_reuse)& (options[AZ_scaling] != AZ_none))
         AZ_scale_f(AZ_SCALE_SOL, Amat, options, b, x, proc_config, scaling);