az_gmres.c

并行解法器,功能强大

      dble_tmp = sqrt(hh[i][i] * hh[i][i] + hh[i1][i] * hh[i1][i]);

      /* Estimate condition number of the GMRES */
      /* least-squares problem using ICE.       */

      if (i == 0) {
        big = dble_tmp;
        small = big;
        svbig[0] = doubleone;
        svsml[0] = doubleone;
      }
      else {
        for (k = 0; k < i; k++) vectmp[k] = hh[k][i];
        vectmp[i] = dble_tmp;
        ijob = 1;
        dlaic1_(&ijob, &i, svbig, &big, vectmp, &vectmp[i], &sestpr, &ss, &cc);
        big = sestpr;
        dscal_(&i, &ss, svbig, &one);
        svbig[i] = cc;
        ijob = 2;
        dlaic1_(&ijob, &i, svsml, &small, vectmp, &vectmp[i], &sestpr, &ss,
                &cc);
        small = sestpr;
        dscal_(&i, &ss, svsml, &one);
        svsml[i] = cc;
      }

      if ((small == 0.0) || (dble_tmp  < DBL_EPSILON * r_2norm) ||
          (big/small > 1.0e+11) ) {

        /* take most recent solution and get out */

        for (k = 0; k < i1; k++) tmp[k] = rs[k];
        AZ_get_x_incr(options, data_org, proc_config, params, i, hh, tmp, 
                      temp, v, Amat, precond, x, &first_time, &converged, kspace);

        AZ_scale_true_residual(x, b,
                               v[kspace], weight, &actual_residual,
                               &true_scaled_r, options, data_org, proc_config, 
                               Amat, convergence_info);

        if (dble_tmp  < DBL_EPSILON * r_2norm) i = AZ_breakdown;
        else                                   i = AZ_ill_cond;

        AZ_terminate_status_print(i, iter, status, rec_residual, params,
                                  true_scaled_r, actual_residual, options,
                                  proc_config);
        return;
      }

      dble_tmp = 1.0 / dble_tmp;
      c[i]     =  hh[i][i]  * dble_tmp;
      s[i]     =  hh[i1][i] * dble_tmp;
      rs[i1]   = -s[i] * rs[i];
      rs[i]   *=  c[i];

      if (r_avail) {
        dble_tmp  = c[i] * rs[i1];
        dble_tmp1 = s[i] * s[i];

        for (k = 0; k < N; k++)
          res[k] = dble_tmp*v[i1][k] + dble_tmp1*res[k];
      }

      /* determine residual norm & test convergence */

      hh[i][i]     = c[i] * hh[i][i] + s[i] * hh[i1][i];
      r_2norm      = fabs(rs[i1]);
      rec_residual = r_2norm;

      /*
       * Compute a few global scalars:
       *     1) ||r||                corresponding to options[AZ_conv]
       *     2) scaled ||r||         corresponding to options[AZ_conv]
       * NOTE: if r_avail = AZ_TRUE or AZ_FIRST is passed in, we perform
       * step 1), otherwise ||r|| is taken as rec_residual.
       */

      AZ_compute_global_scalars(Amat, x, b,
                                res, weight, &rec_residual, &scaled_r_norm,
                                options, data_org, proc_config, &r_avail, dummy,
                                dummy, dummy, convergence_info);

      converged = scaled_r_norm < epsilon;

      if ( (iter%print_freq == 0) && proc == 0)
        (void) fprintf(stdout, "%siter: %4d           residual = %e\n",prefix,iter,
                       scaled_r_norm);

      i++;      /* subspace dim. counter dim(K) = i - 1 */

      if ( (i == kspace) || converged || iter == options[AZ_max_iter]) {

        /* update x and set temp to delta x */

        for (k = 0; k <= i1; k++) tmp[k] = rs[k];

        AZ_get_x_incr(options, data_org, proc_config, params, i, hh, tmp, 
                temp, v, Amat, precond,  x, &first_time, &converged, kspace);

      }

      if (converged) {

        /* compute true residual using 'v[kspace]' as a temporary vector */

        AZ_scale_true_residual(x, b,
                               v[kspace], weight, &actual_residual,
                               &true_scaled_r, options, data_org, proc_config, 
                               Amat, convergence_info);

        converged = true_scaled_r < params[AZ_tol];

        if (!converged && (AZ_get_new_eps(&epsilon, scaled_r_norm,
                                          true_scaled_r,
                                          proc_config) == AZ_QUIT)) {

          /*
           * Computed residual has converged, actual residual has not converged,
           * AZ_get_new_eps() has decided that it is time to quit.
           */

          AZ_terminate_status_print(AZ_loss, iter, status, rec_residual, params,
                                    true_scaled_r, actual_residual, options,
                                    proc_config);
          return;
        }

        /* restore previous solution */

        if ( (!converged) && (i != kspace) )
          for (k = 0; k < N; k++) x[k] = x[k] - temp[k];
      }

      if ( (i == kspace) && !converged) {
        if (r_avail)
          for (k = 0; k < N; k++) v[0][k] = res[k];
        else {
          Amat->matvec(temp, v[kspace], Amat, proc_config);

          dscal_(&N, &first_2norm, v[0], &one);
          daxpy_(&N, &minusone, v[kspace], &one, v[0], &one);
        }
      }
    }
  }

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


  /* check if we exceeded maximum number of iterations */

  if (converged) {
    i = AZ_normal;
    scaled_r_norm = true_scaled_r;
  }
  else
    i = AZ_maxits;

  AZ_terminate_status_print(i, iter, status, rec_residual, params,
                            scaled_r_norm, actual_residual, options,
                            proc_config);

} /* AZ_pgmres */

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

void AZ_get_x_incr(int options[], int data_org[], int
                   proc_config[], double params[], int i, double **hh, double
                   *rs, double *temp, double **v, AZ_MATRIX *Amat, AZ_PRECOND 
                   *precond, double x[], int *first_time, int *converged,
		   int kspace)

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

  This routine is normally invoked from GMRES and is used to compute the
  increment to the solution (as well as the new solution) that should be applied as a 
  result of solving the upper hessenberg system produces by GMRES.

  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.

  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 file 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.

  params:          Drop tolerance and convergence tolerance info.

  i:               The size of the Krylov subspace that has be generated.

  hh:              The upper triangular matrix produced by gmres (after reducing
                   the upper hessenberg matrix) after sweeping out i steps.

  rs:              The projected right hand side produced by gmres.

  temp:            vector of length data_org[AZ_num_int_unk] +
                   data_org[AZ_num_bord_unk].  On output to AZ_get_x_incr(),
                   temp contains the increment that must be added the current
                   approximation to get the new gmres approximate solution.

  v:               Orthogonal vectors produced by Gram-Schmidt process that
                   span the Krylov space swept out.


  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).


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

{

  /* local variables */

  int precond_flag;
  int    ii, k, k1, j, N;
  double t, doubleone = 1.0, update_norm = 1.0;
  int    one = 1;

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

  if (i == 0) return;

  N = data_org[AZ_N_internal] + data_org[AZ_N_border];

  /* solve upper triangular system, compute solution */

  i--;  /* set i = Krylov subspace dimension */

  rs[i] /= hh[i][i];
  for (ii = 1; ii <= i; ii++) {
    k  = i - ii;
    k1 = k + 1;
    t  = rs[k];
    for (j = k1; j <= i; j++) t -= hh[k][j] * rs[j];
    rs[k] = t / hh[k][k];
  }

  /*
   * done with back substitution; form linear combination to get solution
   */

  precond_flag = options[AZ_precond];
  if (options[AZ_check_update_size] & *converged) {
     for (j = 0; j < N; j++) temp[j] = v[i][j];
     if (precond_flag) precond->prec_function(temp,options,proc_config,params, Amat,precond);
     update_norm = fabs(rs[i]*sqrt(AZ_gdot(N, temp, temp, proc_config)));
  }

  for (j = 0; j < N; j++) temp[j] = 0.0;

  for (j = 0; j <= i; j++) {
    t = rs[j];
    daxpy_(&N, &t, v[j], &one, temp, &one);
  }

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

  daxpy_(&N, &doubleone, temp, &one, x, &one);

  if (options[AZ_check_update_size] & *converged) {
     *converged =AZ_compare_update_vs_soln(N,update_norm, rs[i],temp,x,
                       params[AZ_update_reduction],options[AZ_output],proc_config,first_time);

     /* restore previous solution */

     if ( (!(*converged)) && (i != kspace) ) {
        doubleone = -1.;
        daxpy_(&N, &doubleone, temp, &one, x, &one);
     }
  }

} /* AZ_get_x_incr */