az_precond.c

并行解法器,功能强大

                   file Aztec User's Guide).

  d_inv:           Vector containing the LU of the diagonal blocks.

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

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

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

  d_bpntr:         The 'bpntr' array corresponding to the LU-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         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 */

  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, &(ipvt[rpntr[iblk_row]]), &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);
          }
          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 *) work);

} /* AZ_calc_blk_diag_inv */

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

void jacobi(double val[], double x[], int data_org[])

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

  Simple Jacobi iteration (undamped).  Not yet implemented for DVBR formatted
  matrices.

  Author:          John N. Shadid, 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).

  x:               On input, contains the current solution to the linear system.
                   On output contains the Jacobi preconditioned solution.

  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;
  int          N;

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

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

  for (i = 0; i < N; i++)
    *x++ /=  *val++;

} /* jacobi */

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

extern void AZ_sym_gauss_seidel(void)

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

  Symmetric Gauss-Siedel preconditioner

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

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

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

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

{

  /* local variables */

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

  (void) fprintf(stderr, "WARNING: sym Gauss-Seidel preconditioning not\n"
                 "         implemented for VBR matrices\n");
  exit(-1);

} /* AZ_sym_gauss_seidel */

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

void AZ_sym_gauss_seidel_sl(double val[],int bindx[],double x[],int data_org[],
			    int options[], struct context *context,
			    int proc_config[])

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

  Symmetric Gauss-Siedel preconditioner.

  Author:          John N. Shadid, 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).

  x:               On input, contains the current solution to the linear system.
                   On output contains the Jacobi preconditioned solution.

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

  options:         Determines specific solution method and other parameters.

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

{

  /* local variables */

  register int    *bindx_ptr;
  register double sum, *ptr_val;
  int             i, bindx_row, j_last, N, step, ione = 1, j;
  double          *b, *ptr_b;
  char            tag[80];

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

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

  sprintf(tag,"b/sGS %s",context->tag);
  b = AZ_manage_memory(N*sizeof(double), AZ_ALLOC, AZ_SYS, tag, &i);

  dcopy_(&N, x, &ione, b, &ione);
  ptr_val = val;

  for (i = 0; i < N; i++) {
    (*ptr_val) = 1.0 / (*ptr_val);
    x[i]     = 0.0;
    ptr_val++;
  }

  for (step = 0; step < options[AZ_poly_ord]; step++) {
    AZ_exchange_bdry(x, data_org, proc_config);

    bindx_row = bindx[0];
    bindx_ptr = &bindx[bindx_row];
    ptr_val   = &val[bindx_row];
    ptr_b   = b;

    for (i = 0; i < N; i++) {
      sum    = *ptr_b++;
      j_last = bindx[i+1] - bindx[i];

      for (j = 0; j < j_last; j++) {
        sum -= *ptr_val++ * x[*bindx_ptr++];
      }
      x[i] = sum * val[i];
    }

    bindx_row = bindx[N];
    bindx_ptr = &bindx[bindx_row-1];
    ptr_val   = &val[bindx_row-1];

    for (i = N - 1; i >= 0; i--) {
      sum = b[i];
      j_last  = bindx[i+1] - bindx[i];

      for (j = 0; j < j_last; j++) {
        sum -= *ptr_val-- * x[*bindx_ptr--];
      }
      x[i] = sum * val[i];
    }
  }

  for (i = 0; i < N; i++)
    val[i] = 1.0 / val[i];

} /* AZ_sym_gauss_seidel_sl */

#ifdef eigen
void AZ_do_Jacobi(double val[], int indx[], int bindx[], int rpntr[],
                     int cpntr[], int bpntr[], double x[], double b[],
                     double temp[], int options[], int data_org[],
                     int proc_config[], double params[], int flag)
{
double *v;
int i,step;
int N;
 
  N    = data_org[AZ_N_internal] + data_org[AZ_N_border];
 
    if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {
 
      if ( (options[AZ_poly_ord] != 0) && (flag == 1) )
         for (i = data_org[AZ_N_internal]; i < N; i++) x[i] = b[i]/val[i];
 
      if (options[AZ_poly_ord] > flag) {
        v = AZ_manage_memory((N+data_org[AZ_N_external])*sizeof(double),
                             AZ_ALLOC, AZ_SYS, "v in do_jacobi", &i);
 
        for (i = 0; i < N; i++) v[i] = x[i];
 
        for (step = flag; step < options[AZ_poly_ord]; step++) {
          Amat->matvec(v, temp, Amat, proc_config);

          for(i = 0; i < N; i++) v[i] += (b[i] - temp[i]) / val[i];
        }
        for (i = 0; i < N; i++) x[i] = v[i];
      }
    }
    else {
       (void) fprintf(stderr,"AZ_slu with option[AZ_poly_ord] > 0 only \n");
       (void) fprintf(stderr,"implemented for MSR matrices.\n");
       exit(-1);
    }

}
#endif