az_qmrcgs.c

并行解法器,功能强大

    if (iter==1) status[AZ_first_precond] = AZ_second() - init_time;

    for (i = 0; i < N; i++) ubar[i] = ubar[i] + beta * qbar[i];

    Amat->matvec(ubar, Aubar, Amat, proc_config);

    daxpy_(&N, &beta, v, &one, Aqbar, &one);
    for (i = 0; i < N; i++) v[i] = Aubar[i] + beta * Aqbar[i];

    sigma = AZ_gdot(N, rtilda, v, proc_config);

    if (fabs(sigma) < brkdown_tol) { /* possible problem */
      if (AZ_breakdown_f(N, rtilda, v, sigma, proc_config)) {

        /* break down */

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

        AZ_terminate_status_print(AZ_breakdown, iter, status, est_residual,
                                  params, true_scaled_r, actual_residual,
                                  options, proc_config);
        return;
      }
      else brkdown_tol = 0.1 * fabs(sigma);
    }

    alpha = rhon / sigma;

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

    dcopy_(&N, v, &one, qbar, &one);

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

    for (i = 0; i < N; i++) qbar[i] = ubar[i] - alpha * qbar[i];
    Amat->matvec(qbar, Aqbar, Amat, proc_config);

    for (i = 0; i < N; i++) r_cgs[i] = r_cgs[i] - alpha*(Aubar[i] + Aqbar[i]);

    /* QMRS scaling and iterates weights 5.11 */

    norm_r_n_cgs = sqrt(AZ_gdot(N, r_cgs, r_cgs, proc_config));

    /* m is odd in Freund's paper */

    omega = sqrt(norm_r_nm1_cgs * norm_r_n_cgs);
    nu_m  = omega / tau_mm1;
    c     = 1.0 / sqrt(1.0 + nu_m * nu_m);
    tau_m = tau_mm1 * nu_m * c;
    eta_m = c * c * alpha;

    if (brkdown_will_occur) {
      AZ_scale_true_residual(x, b, v,
                             weight, &actual_residual, &true_scaled_r, options,
                             data_org, proc_config, Amat, convergence_info);

      AZ_terminate_status_print(AZ_breakdown, iter, status, est_residual,
                                params, true_scaled_r, actual_residual, options,
                                proc_config);
      return;
    }

    dtemp = nu_mm1 *nu_mm1 * eta_mm1 / alpha;
    for (i = 0; i < N; i++) d[i] = ubar[i] + dtemp * d[i];
    daxpy_(&N, &eta_m, d, &one, x, &one); /* x = x - eta_m d  */

    if (r_avail) {
      for (i = 0; i < N; i++) Ad[i] = Aubar[i] + dtemp * Ad[i];
    }

    /* save some values */

    eta_mm1 = eta_m;  tau_mm1        = tau_m;
    nu_mm1  = nu_m;   norm_r_nm1_cgs = norm_r_n_cgs;

    /* m is even in Freund's paper */

    omega = norm_r_n_cgs;

    if (tau_mm1 == 0.0) nu_m = 0.0;
    else                nu_m  = omega / tau_mm1;

    c     = 1.0 / sqrt(1.0 + nu_m * nu_m);
    tau_m = tau_mm1 * nu_m * c;

    if (options[AZ_check_update_size]) {
       eta_m = eta_m/(c*c*alpha);
       for (i = 0; i < N; i++) ubar[i] = eta_m*d[i];
    }
    eta_m = c * c * alpha;

    dtemp = nu_mm1 * nu_mm1 * eta_mm1 / alpha;
    for (i = 0; i < N; i++) d[i] = qbar[i] + dtemp * d[i];
    daxpy_(&N, &eta_m, d, &one, x, &one); /* x = x - eta_m d  */

    if (r_avail) {
      for (i = 0; i < N; i++) Ad[i] = Aqbar[i] + dtemp * Ad[i];
    }

    /* save some values */

    eta_mm1 = eta_m;  tau_mm1        = tau_m;
    nu_mm1  = nu_m;   norm_r_nm1_cgs = norm_r_n_cgs;
    rhonm1  = rhon;

    if (r_avail) {
      for (i = 0; i < N; i++) Aubar[i] = r_cgs[i] - (eta_m - alpha) * Ad[i];

      /* Note: Aubar temporarily holds qmr residual */

    }
    else {

      /*
       * We want to estimate the 2-norm of the qmr residual. Freund gives the
       * bound ||r|| <= tau_m * sqrt(2*iter+1).  We use this bound until we get
       * close to the solution. At that point we compute the real residual norm
       * and use this to estimate the norm of ||W|| in Freund's paper.
       */

      dtemp = sqrt((double) (2 * iter + 1));

      if ((scaled_r_norm < epsilon * dtemp) && !offset) {
        AZ_scale_true_residual(x, b,
                               Aubar, weight, &actual_residual, &true_scaled_r,
                               options, data_org, proc_config, Amat,
			       convergence_info);

        if (tau_m != 0.0) W_norm = actual_residual / tau_m;
        if (W_norm < 1.0) W_norm = 1.0;

        offset       = 2 * iter + 1;
        est_residual = actual_residual;
      }
      else est_residual = sqrt((double)(2 * iter + 1 - offset) +
                               W_norm * W_norm) * tau_m;
    }

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

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

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

    /* convergence tests */

    converged = scaled_r_norm < epsilon;
    if (options[AZ_check_update_size] & converged) {
      daxpy_(&N, &doubleone , d, &one, ubar, &one); 
      converged = AZ_compare_update_vs_soln(N, -1.,eta_m, ubar, x,
                                           params[AZ_update_reduction],
                                           options[AZ_output], proc_config, &first_time);
    }

    if (converged) {
      AZ_scale_true_residual(x, b, Aubar,
                             weight, &actual_residual, &true_scaled_r, options,
                             data_org, proc_config, Amat,convergence_info);

      converged = true_scaled_r < params[AZ_tol];

      /*
       * Note: epsilon and params[AZ_tol] may not be equal due to a previous
       * call to AZ_get_new_eps().
       */

      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, est_residual, params,
                                  true_scaled_r, actual_residual, options,
                                  proc_config);
        return;
      }
    }

    beta = rhon / rhonm1;
  }

  iter--;

  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, est_residual, params,
                            scaled_r_norm, actual_residual, options,
                            proc_config);

} /* pqmrs */