az_precond.c

并行解法器,功能强大

        /* symmetric Gauss-Seidel preconditioner only available on 1 proc */

        if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) AZ_sym_gauss_seidel();
        else if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX)
           AZ_sym_gauss_seidel_sl(val, bindx, current_rhs, data_org, options,
				  precond->context, proc_config);
     break;

     case AZ_Neumann:
     case AZ_ls:
        if (!options[AZ_poly_ord]) return;
        AZ_polynomial_expansion(current_rhs, options, proc_config, precond);
     break;

     case AZ_dom_decomp:
     case AZ_rilu:
        AZ_domain_decomp(current_rhs, precond->Pmat, options, proc_config, 
                         params, precond->context);
     break;

     case AZ_icc:
        /* incomplete Cholesky factorization */

        (void) printf("Incomplete Cholesky not available (use ilu).\n");
     break;

     case AZ_user_precond:
        precond->prec_function(current_rhs, options, proc_config, 
                               params, Amat, precond);
     break;
     case AZ_smoother:
        sprintf(label,"istatus %s",precond->context->tag);
        istatus = AZ_manage_memory(AZ_STATUS_SIZE*sizeof(double),AZ_ALLOC,
				   AZ_SYS, label,&i);
        for (i = 0 ; i < AZ_STATUS_SIZE ; i++ ) istatus[i] = 0.0;

        sprintf(label,"y %s",precond->context->tag);
        y = AZ_manage_memory((N+max_externals)*sizeof(double), AZ_ALLOC, 
			     AZ_SYS, label, &i);
        sprintf(label,"tttemp %s",precond->context->tag);
        tttemp = AZ_manage_memory((N+max_externals)*sizeof(double),AZ_ALLOC,
				  AZ_SYS, label, &i);

        for (i = 0 ; i < N ; i++ ) tttemp[i] = current_rhs[i];

        N_fixed = 0; fixed_pts = NULL;
        if (Amat->aux_ival != NULL) {
           N_fixed   = Amat->aux_ival[0][0];
           fixed_pts = Amat->aux_ival[1];
        }
        else if (options[AZ_pre_calc] != AZ_sys_reuse)
           printf("Warning: Not fixed points set for local smoothing!!\n");

        for (j = 0; j < options[AZ_poly_ord]; j++) {
           AZ_loc_avg(Amat, tttemp, y, N_fixed, fixed_pts, proc_config);
           norm1 = sqrt(AZ_gdot(N, y, y, proc_config));
           if (proc_config[AZ_node] == 0) {
              if ((j==0) && (options[AZ_output] != AZ_none) &&
                  (options[AZ_output] != AZ_last) &&
                  (options[AZ_output] != AZ_warnings))
                  printf("   %d  %e\n",j, norm1);
              else if ((j==options[AZ_poly_ord]-1) && 
		  (options[AZ_output] != AZ_none) && 
                  (options[AZ_output] != AZ_warnings))
                  printf("   %d  %e\n",j, norm1);
              else if ((options[AZ_output] > 0) && (j%options[AZ_output] == 0))
                  printf("   %d  %e\n",j, norm1);
           }
           for (i = 0 ; i < N ; i++ ) tttemp[i] = y[i];
        }
        for (i = 0 ; i < N ; i++ ) y[i] = current_rhs[i] - y[i];
        for (i = 0 ; i < N ; i++ ) current_rhs[i] = 0.0;

        opt_save1 = options[AZ_output];
        opt_save2 = options[AZ_solver];
        opt_save3 = options[AZ_precond];
        opt_save4 = options[AZ_max_iter];
        opt_save5 = options[AZ_aux_vec];

        options[AZ_output]  = AZ_warnings;
        options[AZ_solver]  = AZ_tfqmr;
        options[AZ_precond] = AZ_dom_decomp;
        options[AZ_max_iter]= 1000;
        options[AZ_aux_vec] = AZ_rand;

        options[AZ_recursion_level]++;
        AZ_oldsolve(current_rhs, y,options, params, istatus, proc_config, 
                    Amat, precond, NULL);
        options[AZ_recursion_level]--;
        options[AZ_output]  = opt_save1;
        options[AZ_solver]  = opt_save2;
        options[AZ_precond] = opt_save3;
        options[AZ_max_iter]= opt_save4;
        options[AZ_aux_vec] = opt_save5;
     break;
     default:
        if (options[AZ_precond] < AZ_SOLVER_PARAMS) {
           AZ_recover_sol_params(options[AZ_precond], &ioptions, &iparams,
                                 &istatus, &Aptr, &Pptr, &Sptr);
           sprintf(label,"y %s",precond->context->tag);
           y = AZ_manage_memory((N+max_externals)*sizeof(double),
                                AZ_ALLOC, AZ_SYS, label, &i);
           for (i = 0 ; i < N ; i++ ) y[i] = current_rhs[i];
           for (i = 0 ; i < N ; i++ ) current_rhs[i] = 0.0;

           ioptions[AZ_recursion_level] = options[AZ_recursion_level] + 1;
           if ((options[AZ_pre_calc] == AZ_sys_reuse) &&
               (ioptions[AZ_keep_info] == 1)) 
              ioptions[AZ_pre_calc] = AZ_reuse;
           AZ_oldsolve(current_rhs, y,ioptions,iparams, istatus, proc_config, 
                       Aptr, Pptr, Sptr);
        }
        else {
           (void) fprintf(stderr, "%sERROR: invalid preconditioning flag.\n"
                   "       options[AZ_precond] improperly set (%d).\n", yo,
			   options[AZ_precond]);
           exit(-1);
        }

     }
     options[AZ_pre_calc] = AZ_sys_reuse;
     precond->context->Pmat_computed = 1;

     if (multilevel_flag) {
        if (precond->next_prec == NULL) {
           multilevel_flag = 0;
           for (i = 0; i < N; i++) current_rhs[i] += x_precond[i];
        }
        else {
           for (i = 0; i < N; i++) x_precond[i] += current_rhs[i];
           AZ_compute_residual(orig_rhs, x_precond, current_rhs, 
                               proc_config, Amat);
           precond = precond->next_prec;
           options = precond->options;
           params  = precond->params;
        }
     }

  } while (multilevel_flag);

  AZ_sys_msg_type = AZ_MSG_TYPE;           /* reset all the message types.   */
                                           /* This is to make sure that all  */
                                           /* processors (even those without */
                                           /* any preconditioning work) have */
                                           /* the same message types for the */
                                           /* next message.                  */
#ifdef TIMING
  ttt = AZ_second() - ttt;
  if (input_options[AZ_recursion_level] == 0) input_precond->timing[0] += ttt;
#endif

} /* precond */

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

void AZ_calc_blk_diag_inv(double *val, int *indx, int *bindx, int *rpntr,
                          int *cpntr, int *bpntr, double *d_inv, int *d_indx,
                          int *d_bindx, int *d_rpntr, int *d_bpntr,
                          int *data_org)

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

  Routine to calculate the inverse of the block-diagonal portion of the sparse
  matrix in 'val' and the associated integer pointer vectors. This is used for
  scaling and/or preconditioning.

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

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

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

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

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

  d_inv:           Vector containing the inverses of the diagonal blocks.

  d_indx:          The 'indx' array corresponding to the inverse-block
                   diagonals.

  d_bindx:         The 'bindx' array corresponding to the inverse-block
                   diagonals.

  d_rpntr:         The 'rpntr' array corresponding to the inverse-block
                   diagonals.

  d_bpntr:         The 'bpntr' array corresponding to the inverse-block
                   diagonals.

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

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

{

  /* local variables */

  register int i, j, iblk_row, jblk, icount = 0, iblk_count = 0, ival;
  int          m1, n1, itemp;
  int          m;
  int          bpoff, idoff;
  int         *ipiv, info;
  double      *work;
  char        *yo = "AZ_calc_blk_diag_inv: ";

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

  m = data_org[AZ_N_int_blk] + data_org[AZ_N_bord_blk];

  if (m == 0) return;

  /* allocate vectors for lapack routines */

  ipiv = (int *)    AZ_allocate(rpntr[m]*sizeof(int));
  work = (double *) AZ_allocate(rpntr[m]*sizeof(double));
  if (work == NULL) AZ_perror("Not enough space for Block Jacobi\n");

  /* offset of the first block */

  bpoff = *bpntr;
  idoff = *indx;

  /* loop over block rows */

  for (iblk_row = 0; iblk_row < m; iblk_row++) {

    /* number of rows in the current row block */

    m1 = rpntr[iblk_row+1] - rpntr[iblk_row];

    /* starting index of current row block */

    ival = indx[bpntr[iblk_row] - bpoff] - idoff;

    /* loop over column block numbers, looking for the diagonal block */

    for (j = bpntr[iblk_row] - bpoff; j < bpntr[iblk_row+1] - bpoff; j++) {
      jblk = bindx[j];

      /* determine the number of columns in this block */

      n1 = cpntr[jblk+1] - cpntr[jblk];

      itemp = m1*n1;

      if (jblk == iblk_row) {   /* diagonal block */

        /* error check */

        if (n1 != m1) {
          (void) fprintf(stderr, "%sERROR: diagonal blocks are not square\n.",
                         yo);
          exit(-1);
        }
        else {

          /* fill the vectors */

          d_indx[iblk_count]  = icount;
          d_rpntr[iblk_count] = rpntr[iblk_row];
          d_bpntr[iblk_count] = iblk_row;
          d_bindx[iblk_count] = iblk_row;

          for (i = 0; i < itemp; i++) d_inv[icount++] = val[ival + i];

          /* invert the dense matrix */

          dgetrf_(&m1, &m1, &d_inv[d_indx[iblk_count]], &m1, ipiv, &info);

          if (info < 0) {
            (void) fprintf(stderr, "%sERROR: argument %d is illegal.\n", yo,
                           -info);
            exit(-1);
          }

          else if (info > 0) {
            (void) fprintf(stderr, "%sERROR: the factorization has produced a "
                           "singular U with U[%d][%d] being exactly zero.\n",
                           yo, info, info);
            exit(-1);
          }

          dgetri_(&m1, &d_inv[d_indx[iblk_count]], &m1, ipiv, work, &m1, &info);

          if (info < 0) {
            (void) fprintf(stderr, "%sERROR: argument %d is illegal.\n", yo,
                           -info);
            exit(-1);
          }

          else if (info > 0) {
            (void) fprintf(stderr, "%sERROR: U[%d][%d] is exactly zero;\n", yo,
                           info, info);
            (void) fprintf(stderr, "the matrix is singular and its inverse "
                           "could not be computed.\n");
            exit(-1);
          }
          iblk_count++;
        }
        break;
      }
      else
        ival += itemp;
    }
  }

  d_indx[iblk_count]  = icount;
  d_rpntr[iblk_count] = rpntr[iblk_row];
  d_bpntr[iblk_count] = iblk_row;

  AZ_free((void *) ipiv);
  AZ_free((void *) work);

} /* AZ_calc_blk_diag_inv */

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

void AZ_calc_blk_diag_LU(double *val, int *indx, int *bindx, int *rpntr,
                          int *cpntr, int *bpntr, double *d_inv, int *d_indx,
                          int *d_bindx, int *d_rpntr, int *d_bpntr,
                          int *data_org, int *ipvt)

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

  Routine to calculate the LU factors of the block-diagonal portion of sparse
  matrix in 'val' and the associated integer pointer vectors. This is used for
  scaling.

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

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

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

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

  indx,
  bindx,
  rpntr,
  cpntr,
  bpntr:           Arrays used for DMSR and DVBR sparse matrix storage (see