az_util.c

并行解法器,功能强大

            j = (int) ( current->name[nn] - '!');
            if (j != i) {
              printf("==> check preconditioner and subdomain solver choices\n");
              return(NULL);
            }
            i = (int) ( name[++nn] - '!');
            j = (int) ( current->name[nn] - '!');
            if (j != i) {
              printf("==> check scaling choices\n");
              return(NULL);
            }
           printf("==> check AZ_ilut_fill, AZ_drop & AZ_graph_fill choices\n");
           return(NULL);
         }
      }
      prev    = current;
      current = current->next;
    }
  }
  else {
    (void) fprintf(stderr, "Error: Invalid action(%d) in AZ_manage_memory()\n",
                   action);
    exit(-1);
  }
  return((double *) NULL);

} /* AZ_manage_memory */

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

void AZ_mysleep(int i)

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

  A dummy routine that just wastes some time.  We essentially try and mimick the
  unix sleep() function which is not available on the nCUBE.

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

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

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

  i:               On input, indicates how much time to spend in this routine.
                   The larger the value of 'i', the more time is spent in this
                   routine. On the Ncube 2, this usually corresponds to about i
                   seconds.

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

{

  /* local variables */

  int    j, k, top, inner;
  double a;

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

  a = 10.01;
  top = 700 * i;
  inner = 10 * 10 * 10;

  for (k = 0; k < top; k++)
    for (j = 0; j < inner; j++)
      a = 2.0 / (a - 1);

} /* AZ_mysleep */

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

void AZ_p_error(char *str, int proc)

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

  Print the string if we are processor 0.

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

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

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

  str:             String to be printed.

  proc:            Current processor number.

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

{
  if (proc == 0) (void) fprintf(stdout, "%s", str);
}

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

int AZ_get_new_eps(double *epsilon, double recursive, double actual,
                   int proc_config[])

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

  Routine which decides what to do when the computed residual has converged but
  not the true residual.

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

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

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

  epsilon:         Residual tolerance.

  recursive:       On input, the norm of the residual produced by an iterative
                   method via recursion (e.g. r = r_0 + delta_1 + delta_2 +
                   delta_3 ... )

  actual:          On input, the norm of the explicitly computed residual (i.e.
                   r = b - A x ). In the absence of rounding error, recursive
                   and actual should be the same value. However, due to rounding
                   errors these two might differ.

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

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

{

  /* local variables */

  double difference;

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

  difference = fabs(actual - recursive);
  if (difference > *epsilon) return AZ_QUIT;
  else {
    *epsilon = *epsilon - 1.5 * difference;
    while (*epsilon < 0.0) *epsilon += .1*difference;
   }

  if (proc_config[AZ_node] == 0)
    (void) printf("\n\t\tTrying to reduce actual residual "
                  "further\n\t\t     (recursive = %e, actual = %e)\n\n",
                  recursive, actual);

  return AZ_CONTINUE;

} /* AZ_get_new_eps */

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

void AZ_terminate_status_print(int situation, int iter, double status[],
                               double rec_residual, double params[],
                               double scaled_r_norm, double actual_residual,
                               int options[], int proc_config[])

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

  Routine to output conditions under which iterative solver was terminated if
  other than specified convergence.

  Author:          Scott A. Hutchinson, SNL, 1421
  =======

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

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

  situation:       On input, indicates why iterative method terminated.
                   = AZ_normal   : iterative solver terminates normally (i.e.
                                   convergence criteria met).
                   = AZ_breakdown: iterative solver has broken down.
                   = AZ_loss     : iterative solver has terminated due to a lack
                                   of accuracy in the recursive residual
                                   (caused by rounding errors).
                   = AZ_maxits   : iterative solver has reached the maximum
                                   number of requested iterations without
                                   convergence.
                   = AZ_ill_cond : the upper triangular factor of the
                                   Hessenberg within GMRES is ill-conditioned.
                                   There is a possibility that the matrix
                                   is singular and we have the least-squares
                                   solution.

  iter:            Number of iterations completed.

  status:          !!!!!THIS MAY NOT BE ACCURATE FOR THIS ROUTINE!!!!!
                   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.

  rec_residual:    On input, the norm of the residual produced by an iterative
                   method via recursion (e.g. r = r_0 + delta_1 + delta_2 +
                   delta_3 ... )

  params:          Drop tolerance and convergence tolerance info.

  scaled_r_norm:   Residual expression (requested by the user via
                   OPTIONS[AZ_conv] ... see user's guide and
                   AZ_compute_global_scalars()) which is used in printing and
                   compared with PARAMS[AZ_TOL] when determining convergence.

  actual_residual: On input, the norm of the explicitly computed residual (i.e.
                   r = b - A x ).  In the absence of rounding error, recursive
                   and actual should be the same value. However, due to rounding
                   errors these two might differ.

  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.

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

{

  /* local variables */

  static int iterations = 0;
  char      *solver_name;
  int        solver_flag, conv_flag;
  double     eps;

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

  eps = params[AZ_tol];

  /* set status */

  if (scaled_r_norm < eps) situation = AZ_normal;

  status[AZ_its] = (double ) iter;
  status[AZ_why] = (double ) situation;
  status[AZ_r]   = actual_residual;

  if (actual_residual == -1.0) status[AZ_r] = rec_residual;

  status[AZ_rec_r]    = rec_residual;
  status[AZ_scaled_r] = scaled_r_norm;

  if (situation == AZ_normal) return; /* nothing else to do */

  if (options[AZ_output] ==  AZ_none) return;

  /* initialize */

  solver_flag = options[AZ_solver];
  conv_flag   = options[AZ_conv];
  if (!iterations) iterations = iter;

  switch (solver_flag) {
  case AZ_cg:
    solver_name = (char *) AZ_allocate(3*sizeof(char));
    (void) sprintf(solver_name, "cg");
    break;
  case AZ_fixed_pt:
    solver_name = (char *) AZ_allocate(9*sizeof(char));
    (void) sprintf(solver_name, "fixed-pt");
    break;
  case AZ_GMRESR:
    solver_name = (char *) AZ_allocate(7*sizeof(char));
    (void) sprintf(solver_name, "gmresr");
    break;
  case AZ_analyze:
    solver_name = (char *) AZ_allocate(8*sizeof(char));
    (void) sprintf(solver_name, "analyze");
    break;
  case AZ_gmres:
    solver_name = (char *) AZ_allocate(6*sizeof(char));
    (void) sprintf(solver_name, "gmres");
    break;
  case AZ_cgs:
    solver_name = (char *) AZ_allocate(4*sizeof(char));
    (void) sprintf(solver_name, "cgs");
    break;
  case AZ_tfqmr:
    solver_name = (char *) AZ_allocate(7*sizeof(char));
    (void) sprintf(solver_name, "tfqmr");
    break;
  case AZ_bicgstab:
    solver_name = (char *) AZ_allocate(9*sizeof(char));
    (void) sprintf(solver_name, "bicgstab");
    break;
  case AZ_symmlq:
    solver_name = (char *) AZ_allocate(7*sizeof(char));
    (void) sprintf(solver_name, "symmlq");
    break;
  case AZ_lu:
    solver_name = (char *) AZ_allocate(4*sizeof(char));
    (void) sprintf(solver_name, "lu");
    break;
  default:
    (void) fprintf(stderr,
                   "Error: invalid solver flag %d passed to terminate_status\n",
                   solver_flag);
  exit(-1);
  }

  if (proc_config[AZ_node] == 0) {
    (void) fprintf(stderr, "\n\n");
    (void) fprintf(stderr,"\t*************************************************"
                   "**************\n\n");

    switch (situation) {
    case AZ_ill_cond:
      (void) fprintf(stderr, "\tWarning: the GMRES Hessenberg matrix is "
                     "ill-conditioned.  This may \n\tindicate that the "
                     "application matrix is singular. In this case, GMRES\n"
                     "\tmay have a least-squares solution.\n");
    break;
    case AZ_breakdown:
      if (!solver_flag) {
        (void) fprintf(stderr, "\tWarning: direction vector is no longer "
                       "A-conjugate \n\tto previous vector but solution has "
                       "not converged.\n");
      }
      else {
        (void) fprintf(stderr, "\tWarning: a breakdown in this "
                       "method\n\thas occurred and solution has not "
                       "converged.\n");
      }
      break;

    case AZ_loss:
      (void) fprintf(stderr, "\tWarning: recursive residual indicates "
                     "convergence\n\tthough the true residual is too large.\n");
      (void) fprintf(stderr, "\n\tSometimes this occurs when storage is ");
      (void) fprintf(stderr, "overwritten (e.g. the\n\tsolution vector was not ");
      (void) fprintf(stderr, "dimensioned large enough to hold\n\texternal ");
      (void) fprintf(stderr, "variables). Other times, this is due to roundoff. ");
      (void) fprintf(stderr, "In\n\tthis case, the solution has either ");
      (void) fprintf(stderr, "converged to the accuracy\n\tof the machine or ");
      (void) fprintf(stderr, "intermediate roundoff errors ");
      (void) fprintf(stderr, "occurred\n\tpreventing full convergence. In the ");
      (void) fprintf(stderr, "latter case, try solving\n\tagain using the new ");
      (void) fprintf(stderr, "solution as an initial guess.\n");
      break;

    case AZ_maxits:
      (void) fprintf(stderr, "\tWarning: maximum number of iterations "
                     "exceeded\n\twithout convergence\n");
    break;

    default:
      (void) fprintf(stderr, "\tError: improper code passed from solver %s\n\n",
                     solver_name);
    (void) fprintf(stderr,"\t***********************************************%s",
                   "**********\n\n");
    exit(-1);
    }

    (void) fprintf(stdout,"\n\tSolver:\t\t\t%s\n", solver_name);
    (void) fprintf(stdout,"\tnumber of iterations:\t%d\n\n", iter);

    if (actual_residual != -1.0)
      (void) fprintf(stdout,"\tActual residual = %11.4e",actual_residual);

    (void) fprintf(stdout,"\tRecursive residual = %11.4e\n\n",rec_residual);
    (void) fprintf(stdout,"\tCalculated Norms\t\t\t\tRequested Norm\n");
    (void) fprintf(stdout,"\t--------------------------------------------");
    (void) fprintf(stdout,"\t--------------\n\n");

    switch (conv_flag) {
    case AZ_noscaled:
      (void) fprintf(stdout, "\t||r||_2 :\t\t%e\t%e\n", scaled_r_norm,
                     eps);
    break;
    case AZ_r0:
      (void) fprintf(stdout, "\t||r||_2 / ||r0||_2:\t\t%e\t%e\n", scaled_r_norm,
                     eps);
    break;
    case AZ_rhs:
      (void) fprintf(stdout, "\t||r||_2 / ||b||_2:\t\t%e\t%e\n", scaled_r_norm,
                     eps);
    break;
    case AZ_Anorm:
      (void) fprintf(stdout, "\t||r||_2 / ||A||_inf:\t\t%e\t%e\n",
                     scaled_r_norm, eps);
    break;
    case AZ_sol: