az_scaling.c

并行解法器,功能强大

       (void) fprintf(stderr, "%sWARNING: sym_row_sum_scaling not implemented"
                   "\n           for VBR matrices\n"
                   "No scaling performed\n", yo);
     }
     return;
  }
  N = data_org[AZ_N_internal] + data_org[AZ_N_border];

  if (action == AZ_INVSCALE_SOL) {
     AZ_sym_reinvscale_sl(x, data_org, options, proc_config, scaling); return;
  }
  if (action == AZ_SCALE_SOL) {
     AZ_sym_rescale_sl(x, data_org, options, proc_config, scaling); return;
  }

  sprintf(label,"sc_vec%d",options[AZ_recursion_level]);
  sc_vec = (double *) AZ_manage_memory((N + data_org[AZ_N_external]) *
                                       sizeof(double), AZ_ALLOC,
                                       data_org[AZ_name], label, &i);
  scaling->action = AZ_left_and_right_scaling;

  if ((options[AZ_pre_calc] >= AZ_reuse) && (i == AZ_NEW_ADDRESS)) {
    (void) fprintf(stderr, "%sERROR: Previous scaling not found for matrix "
                   "with\ndata_org[AZ_name] = %d. Check\n"
                   "options[AZ_pre_calc]\n", yo,
                   data_org[AZ_name]);
    exit(-1);
  }

  if ( (options[AZ_pre_calc] <= AZ_recalc) && (action == AZ_SCALE_MAT_RHS_SOL)) {
    for(irow = 0; irow < N; irow++) {

      /* get row sum */

      j_last  = bindx[irow+1] - bindx[irow];
      bindx_row = bindx[irow];
      row_sum = fabs(val[irow]);

      for (jcol = 0; jcol < j_last; jcol++) {
        k        = bindx_row + jcol;
        row_sum += fabs(val[k]);
      }

      if (fabs(row_sum) < DBL_MIN) {
        (void) fprintf(stderr, "%sERROR: Row %d is all zero's\n", yo, irow);
        exit(-1);
      }

      sc_vec[irow] = 1.0 / sqrt(fabs(row_sum));

      /* scale matrix row */

      val[irow] *= sc_vec[irow];

      for (jcol = 0; jcol < j_last; jcol++) {
        k       = bindx_row + jcol;
        val[k] *= sc_vec[irow];
      }
    }

    AZ_exchange_bdry(sc_vec, data_org, proc_config);

    /* do right diagonal scaling */
    /* index through rows of matrix */

    for (irow = 0; irow < N; irow++) {
      val[irow] *= sc_vec[irow];
      j_last     = bindx[irow+1] - bindx[irow];
      bindx_row    = bindx[irow];

      for (j = 0; j < j_last; j++) {
        k       = bindx_row + j;
        val[k] *= sc_vec[bindx[k]];
      }
    }
  }

  if ( (action == AZ_SCALE_RHS) ) {
       for (irow = 0; irow < N; irow++) b[irow] *= sc_vec[irow];
  }
  if ( action == AZ_INVSCALE_RHS)  {
       for (irow = 0; irow < N; irow++) b[irow] /= sc_vec[irow];
  }
  if (action == AZ_SCALE_MAT_RHS_SOL) {
       for (irow = 0; irow < N; irow++) b[irow] *= sc_vec[irow];
       for (irow = 0; irow < N; irow++) x[irow] /= sc_vec[irow];
  }

} /* sym_row_sum_scaling */

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

void AZ_sym_rescale_sl(double x[], int data_org[], int options[],
	int proc_config[], struct AZ_SCALING *scaling)

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

  Routine to symmetrically diagonally rescale the sparse matrix problem;
  Note: this rescales the entire matrix problem Ax = b.

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

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

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

  x:               Current solution vector.

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

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

{

  /* local variables */

  register int i;
  int          N, k;
  double      *sc_vec;
  char        *yo = "AZ_sym_rescale_sl: ";
  char        label[80];

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

  if ((data_org[AZ_matrix_type] != AZ_MSR_MATRIX) &&
      (data_org[AZ_matrix_type] != AZ_VBR_MATRIX) ) {
     (void) fprintf(stderr,"%sWARNING: Matrix type is neither MSR nor VBR\n",
                    yo);
     return;
  }



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

  sprintf(label,"sc_vec%d",options[AZ_recursion_level]);
  sc_vec = (double *) AZ_manage_memory((N + data_org[AZ_N_external]) *
                                       sizeof(double), AZ_ALLOC,
                                       data_org[AZ_name], label, &k);
  scaling->action = AZ_left_and_right_scaling;

  if (k == AZ_NEW_ADDRESS) {
    (void) fprintf(stderr, "%sWARNING: Scaling vector not found: "
                   "Not rescaling solution\n", yo);
    return;
  }

  for (i = 0; i < N; i++) x[i] = x[i] / sc_vec[i];

  AZ_exchange_bdry(x, data_org, proc_config);

} /* AZ_sym_rescale_sl */
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/

void AZ_sym_reinvscale_sl(double x[], int data_org[], int options[],
	int proc_config[], struct AZ_SCALING *scaling)

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

  Routine to symmetrically diagonally rescale the sparse matrix problem;
  Note: this rescales the entire matrix problem Ax = b.

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

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

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

  x:               Current solution vector.

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

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

{

  /* local variables */

  register int i;
  int          N, k;
  double      *sc_vec;
  char        *yo = "AZ_sym_rescale_sl: ";
  char        label[80];

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

  if ((data_org[AZ_matrix_type] != AZ_MSR_MATRIX) &&
      (data_org[AZ_matrix_type] != AZ_VBR_MATRIX) ) {
     (void) fprintf(stderr,"%sWARNING: Matrix type is neither MSR nor VBR\n",
                    yo);
     return;
  }



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

  sprintf(label,"sc_vec%d",options[AZ_recursion_level]);
  sc_vec = (double *) AZ_manage_memory((N + data_org[AZ_N_external]) *
                                       sizeof(double), AZ_ALLOC,
                                       data_org[AZ_name], label, &k);
  scaling->action = AZ_left_and_right_scaling;

  if (k == AZ_NEW_ADDRESS) {
    (void) fprintf(stderr, "%sWARNING: Scaling vector not found: "
                   "Not rescaling solution\n", yo);
    return;
  }

  for (i = 0; i < N; i++) x[i] = x[i] * sc_vec[i];

  AZ_exchange_bdry(x, data_org, proc_config);

} /* AZ_sym_rescale_sl */

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

static void calc_blk_diag_Chol(double *val, int *indx, int *bindx, int *rpntr,
                               int *cpntr, int *bpntr, double *L, int *d_indx,
                               int *d_bindx, int *d_rpntr, int *d_bpntr,
                               int *data_org)

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

  Routine to calculate the Cholesky factors 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 Usert).

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

  L:               Vector containing the upper Cholesky factors 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 Usert).

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

{

  /* local variables */

  register int i, j, iblk_row, jblk, icount = 0, iblk_count = 0, ival;
  int          m1, n1, itemp, iL;
  int          m;
  int          bpoff, idoff;
  int          info;
  char        *yo = "calc_blk_diag_Chol: ", *uplo = "L";

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

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

  if (m == 0) return;

  /* 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++) L[icount++] = val[ival + i];

          /* Compute the Cholesky factors for this block */

          iL = d_indx[d_bpntr[iblk_row] - *d_bpntr] - *d_indx;
          dpotrf_(uplo, &m1, L+iL, &m1, &info, strlen(uplo));

          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 L with L[%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;

} /* calc_blk_diag_Chol */

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

#ifdef next_version
void AZ_sym_rescale_vbr(double x[], int data_org[], int options[])

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

  Routine to symmetrically block-diagonally rescale the sparse matrix problem;
  Note: this rescales the entire matrix problem Ax = b.

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

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

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

  x:               Current solution vector.

  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 */

  int          N, k;
  double      *sc_vec;
  char        *yo = "AZ_sym_rescale_vbr: ";

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

  if (data_org[AZ_matrix_type] != AZ_VBR_MATRIX) return;

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

  sprintf(label,"sc_vec%d",options[AZ_recursion_level]);
  sc_vec = (double *) AZ_manage_memory((N + data_org[AZ_N_external]) *
                                       sizeof(double), AZ_ALLOC,
                                       data_org[AZ_name], label, &k);
  scaling->action = AZ_left_and_right_scaling;

  if (k == AZ_NEW_ADDRESS) {
    (void) fprintf(stderr, "%sWARNING: Scaling vector not found - "
                   "not rescaling solution\n", yo);
    return;
  }
  /*
    for (i = 0; i < N; i++) x[i] = x[i] * sc_vec[i];

    AZ_exchange_bdry(x, data_org, proc_config);
    */
} /* AZ_sym_rescale_vbr */
#endif

struct AZ_SCALING *AZ_scaling_create()
{
   struct AZ_SCALING *temp;

   temp = (struct AZ_SCALING *) AZ_allocate(sizeof(struct AZ_SCALING));
   if (temp == NULL) {
      printf("AZ_scaling_create: Not enough space\n");
      exit(1);
   }
   temp->A_norm = 0.0;
   temp->action = AZ_none;
   temp->mat_name = 0;
   temp->scaling_opt = AZ_none;
   return temp;
}
void AZ_scaling_destroy(struct AZ_SCALING **temp)
{
   if (*temp != NULL) AZ_free(*temp);
   *temp = NULL;
}