az_qmrcgs.c

并行解法器,功能强大

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

void AZ_pqmrs(double b[], double x[], double weight[], int options[], 
	double params[], int proc_config[], double status[], AZ_MATRIX *Amat, 
	AZ_PRECOND *precond, struct AZ_CONVERGE_STRUCT *convergence_info)

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

  Freund's transpose free QMR routine to solve the nonsymmetric matrix problem
  Ax = b. NOTE: this routine differs from Freund's paper in that we compute
  ubar (= M^-1 u ) and qbar (= M^-1 q) instead of u and q defined in Freund's
  paper.

  IMPORTANT NOTE: While an estimate of the 2-norm of the qmr residual is
  available (see comment below), the actual qmr residual is not normally
  computed as part of the qmr algorithm. Thus, if the user uses a convergence
  condition (see AZ_compute_global_scalars()) that is based on the 2-norm of the
  residual there is no need to compute the residual (i.e. r_avail = AZ_FALSE).
  However, if another norm of r is requested, AZ_compute_global_scalars() will
  set r_avail = AZ_TRUE and the algorithm will compute the residual.

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

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

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

  b:               Right hand side of linear system.

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

  weight:          Vector of weights for convergence norm #4.

  options:         Determines specific solution method and other parameters.

  params:          Drop tolerance and convergence tolerance info.

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

  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.

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

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

{

  /* local variables */

  register int    i;
  int             N, NN, converged, one = 1, iter= 1,r_avail = AZ_FALSE, j;
  int             precond_flag, print_freq, proc;
  int             brkdown_will_occur = AZ_FALSE;
  double          alpha, beta = 0.0, true_scaled_r=0.0;
  double          *ubar, *v, *r_cgs, *rtilda, *Aubar, *qbar, *Aqbar, *d, *Ad = NULL;
  double          rhonm1, rhon, est_residual, actual_residual = -1.0;
  double          scaled_r_norm, sigma, epsilon, brkdown_tol = DBL_EPSILON;
  double          omega, c, norm_r_n_cgs, norm_r_nm1_cgs;
  double          tau_m, nu_m, eta_m, init_time = 0.0;
  double          tau_mm1, nu_mm1 = 0.0, eta_mm1 = 0.0, doubleone = 1.0;
  register double dtemp;
  double          W_norm = 0.0;
  int             offset = 0;
  int          *data_org, str_leng, first_time = AZ_TRUE;
  char         label[64],suffix[32], prefix[64];

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

  sprintf(suffix," in qmrcgs%d",options[AZ_recursion_level]);
                                                /* set string that will be used */
                                                /* for manage_memory label      */
  /* set prefix for printing */

  str_leng = 0;
  for (i = 0; i < 16; i++) prefix[str_leng++] = ' ';
  for (i = 0 ; i < options[AZ_recursion_level]; i++ ) {
     prefix[str_leng++] = ' '; prefix[str_leng++] = ' '; prefix[str_leng++] = ' ';
     prefix[str_leng++] = ' '; prefix[str_leng++] = ' ';
  }
  prefix[str_leng] = '\0';

  data_org = Amat->data_org;

  /* pull needed values out of parameter arrays */

  N            = data_org[AZ_N_internal] + data_org[AZ_N_border];
  precond_flag = options[AZ_precond];
  epsilon      = params[AZ_tol];
  proc         = proc_config[AZ_node];
  print_freq   = options[AZ_print_freq];

  /* allocate memory for required vectors */

  NN     = N + data_org[AZ_N_external];
  if (NN == 0) NN++; /* make sure everyone allocates something */
  NN = NN + (NN%2);
      /* make sure things are aligned on double words for paragon */


  sprintf(label,"ubar%s",suffix);
  ubar   = (double *) AZ_manage_memory(8*NN*sizeof(double),
				      AZ_ALLOC,AZ_SYS,label,&j);
  v      = &(ubar[1*NN]);
  Aubar  = &(ubar[2*NN]);
  d      = &(ubar[3*NN]);
  qbar   = &(ubar[4*NN]);
  rtilda = &(ubar[5*NN]);
  Aqbar  = &(ubar[6*NN]);
  r_cgs  = &(ubar[7*NN]);

  AZ_compute_residual(b, x, r_cgs, proc_config, Amat);

  /* d, qbar, Aqbar, v = 0 */

  for (i = 0; i < N; i++)
    d[i] = qbar[i] = Aqbar[i] = v[i] = 0.0;

  /* set rtilda */

  if (options[AZ_aux_vec] == AZ_resid) dcopy_(&N, r_cgs, &one, rtilda, &one);
  else AZ_random_vector(rtilda, data_org, proc_config);

  /*
   * Compute a few global scalars:
   *     1) ||r_cgs||              corresponding to options[AZ_conv]
   *     2) scaled ||r_cgs||       corresponding to options[AZ_conv]
   *     3) rhon = <rtilda, r_cgs>
   * Note: step 1) is performed if r_avail = AZ_TRUE on entry or
   *       AZ_FIRST_TIME is passed in. Otherwise, ||r_cgs|| is taken as
   *       est_residual.
   */

  AZ_compute_global_scalars(Amat, x, b, r_cgs,
                            weight, &est_residual, &scaled_r_norm, options,
                            data_org, proc_config, &r_avail, r_cgs, rtilda,
                            &rhon, convergence_info);
  true_scaled_r = scaled_r_norm;

  if ((options[AZ_output] != AZ_none) && (options[AZ_output] != AZ_last) &&
      (options[AZ_output] != AZ_warnings) && (proc == 0))
    (void) fprintf(stdout, "%siter:    0           residual = %e\n",prefix,scaled_r_norm);

  norm_r_nm1_cgs = est_residual;
  tau_mm1        = norm_r_nm1_cgs;
  rhonm1         = rhon;

  /* Set up aux-vector if we need to compute the qmr residual */

  if (r_avail) {
    sprintf(label,"Ad%s",suffix);
    Ad = (double *) AZ_manage_memory(NN*sizeof(double),AZ_ALLOC,
				     AZ_SYS, label, &j);
    for (i = 0; i < N; i++) Ad[i] = 0.0;
  }

  converged = scaled_r_norm < epsilon;

  for (iter = 1; iter <= options[AZ_max_iter] && !converged; iter++) {

    if (fabs(rhon) < brkdown_tol) { /* possible breakdown problem */

      if (AZ_breakdown_f(N, r_cgs, rtilda, rhon, proc_config))
        brkdown_will_occur = AZ_TRUE;
      else
        brkdown_tol = 0.1 * fabs(rhon);
    }

    /* ubar  = M^-1 r_cgs + beta*qbar               */
    /* Aubar = A ubar                               */
    /* v     = A ubar + beta ( A qbar + beta pnm1 ) */
    /*       = Aubar  + beta ( Aqbar +  beta v)     */

    dcopy_(&N, r_cgs, &one, ubar, &one);

    if (iter==1) init_time = AZ_second();

    if (precond_flag)
    precond->prec_function(ubar,options,proc_config,params,Amat,precond);