az_scaling.c

并行解法器,功能强大

        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]);
        }

        row_sum = row_sum * sgn(val[irow]);

        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 / 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];
        }
      }
    }

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

      /* loop over the block rows */

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

        /* find out how many rows are in this block row */

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

        /* find out how many block columns are in this block row */

        num_col_blks = bpntr[iblk_row+1] - bpntr[iblk_row];

        /* loop over all the rows in this block row */

        for (irow = 0; irow < num_blk_rows; irow++) {
          I = rpntr[iblk_row] + irow;     /* true matrix row number */

          /* loop over all the column blocks in this block row */

          for (jblk_col = 0; jblk_col < num_col_blks; jblk_col++) {

            /* find out which column block this is */

            icol_blk = bindx[iblk];
            indx_ptr = indx[iblk++];

            /* find out how many columns are in this block */

            num_blk_cols = cpntr[icol_blk+1] - cpntr[icol_blk];

            /* loop over all the columns in this block */

            for (jcol = 0; jcol < num_blk_cols; jcol++) {
              J     = cpntr[icol_blk] + jcol;     /* true matrix column */
              t_index = indx_ptr + jcol*num_blk_rows + irow;

              /* diagonal entry => get sign */

              if (I == J) sign = sgn(val[t_index]);

              row_sum += fabs(val[t_index]);
            }
          }

          /* reset the block counter */

          iblk -= num_col_blks;

          if ( fabs(sign) < (1.0 - sqrt(DBL_EPSILON)) ) {
            (void) fprintf(stderr, "%sERROR: sign not set => no diagonal "
                           "entry.\n         inside row_sum.\n", yo);
            exit(-1);
          }

          if (fabs(row_sum) == 0.0) {
            (void) fprintf(stderr,"%sERROR: row %d is all 0's.\n", yo, I);
            exit(-1);
          }

          inv_row_sum = sign / row_sum;
          sc_vec[I]   = inv_row_sum;
          row_sum     = sign = 0.0;

          /* scale the matrix */

          for (jblk_col = 0; jblk_col < num_col_blks; jblk_col++) {

            /* find out which column block this is */

            icol_blk = bindx[iblk];
            indx_ptr = indx[iblk++];

            /* find out how many columns are in this block */

            num_blk_cols = cpntr[icol_blk+1] - cpntr[icol_blk];

            /* loop over all the columns in this block */

            for (jcol = 0; jcol < num_blk_cols; jcol++) {
              J           = cpntr[icol_blk] + jcol;     /* true matrix column */
              t_index       = indx_ptr + jcol * num_blk_rows + irow;
              val[t_index] *= inv_row_sum;
            }
          }

          /* reset the block counter */

          iblk -= num_col_blks;
        }

        /* last row in this block row => offset the correction above */

        iblk += num_col_blks;
      }
    }
  }
  if ( (action == AZ_SCALE_MAT_RHS_SOL) || (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];
  }

} /* AZ_row_sum_scaling */

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

void sym_row_sum_scaling(int proc_config[])

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

  Routine to symmetricaly row-sum scale sparse matrix problem; Note: this scales
  the entire matrix problem Ax = b, the routine sym_rescale must be used to
  transform solution back to recover soution to original problem.

  Author:          , SNL,
  =======

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

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

  proc_config:     Machine configuration.  proc_config[AZ_node] is the node
                   number.  proc_config[AZ_N_procs] is the number of processors.

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

{

  /* local variables */

  char *yo = "sym_row_sum_scaling: ";

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

  if (proc_config[AZ_node]  == 0) {
    (void) fprintf(stderr, "%sWARNING: sym_row_sum_scaling not implemented"
                   "           for VBR matrices\n"
                   "No preconditioning performed\n", yo);
  }

} /* sym_row_sum_scaling */

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

void AZ_sym_diagonal_scaling(int action, AZ_MATRIX *Amat, 
		        double b[], double x[], 
                        int options[], 
			int proc_config[], struct AZ_SCALING *scaling)

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

  Routine to symmetrically diagonally scale sparse matrix problem; Note: this
  scales the entire matrix problem Ax = b, the routine sym_rescale must be used
  to transform solution back to recover soution to original problem.

  Author:          John N. Shadid, SNL, 1421 (MSR format)
  =======          Lydie Prevost, SNL 1422 (VBR format) 

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

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

  val:             Array containing the nonzero entries of the matrix (see
                    Aztec Usert).

  bindx:           Arrays used for DMSR and DVBR sparse matrix storage (see
                   file  Aztec Usert).

  b:               Right hand side of linear system.

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

  options:         Determines specific solution method and other parameters.

  proc_config:     Machine configuration.  proc_config[AZ_node] is the node
                   number.  proc_config[AZ_N_procs] is the number of processors.

  x:               Current solution vector.

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

{

  /* local variables */

  register int j, k, irow, icol;
  int          N, m;
  int          j_last, bindx_row, i;
  double       *sc_vec;
  int count;
  char         label[80];

  char        *yo = "AZ_sym_diagonal_scaling: ";
  int         *indx, *bindx, *rpntr, *cpntr, *bpntr, *data_org;
  double      *val;


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

  val  = Amat->val;
  indx = Amat->indx;
  bindx = Amat->bindx;
  rpntr = Amat->rpntr;
  cpntr = Amat->cpntr;
  bpntr = Amat->bpntr;
  data_org = Amat->data_org;
  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;
  }

  N = data_org[AZ_N_internal] + data_org[AZ_N_border];
  m    = data_org[AZ_N_int_blk]  + data_org[AZ_N_bord_blk];



  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 (((action == AZ_SCALE_RHS) || (action == AZ_INVSCALE_RHS) || 
       (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\n", yo,
                   data_org[AZ_name]);
    exit(-1);
  }



  if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {

    /***** MSR Matrix *****/


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

        /* scale matrix */

        j_last  = bindx[irow+1] - bindx[irow];
        bindx_row = bindx[irow];

        if (fabs(val[irow]) < DBL_MIN) {
          (void) fprintf(stderr, "%sERROR: diagonal of row %d is zero\n", yo,
                         irow);
          exit(-1);
        }

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

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

      /* do right diagonal scaling */

      AZ_exchange_bdry(sc_vec, data_org, proc_config);

      /* 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]];
        }
      }
    }

  }

  else {

    /***** VBR Matrix *****/
   
    /* find the diagonal terms */

    if ( (options[AZ_pre_calc] <= AZ_recalc) && (action == AZ_SCALE_MAT_RHS_SOL)) {

     /* index through block rows of matrix */
   
         for (irow = 0; irow < m; irow++) {
   
         /* for all the blocks in a row */
   
         for (k = bpntr[irow]; k < bpntr[irow+1]; k++) {
            icol = bindx[k];
   
            count = 0;
   
            for (j = cpntr[icol]; j < cpntr[icol+1]; j++) {
               for (i = rpntr[irow]; i < rpntr[irow+1]; i++) {
                  if ( icol == irow && i ==j ){
                     sc_vec[i] =  1.0 / sqrt(fabs(val[indx[k]+count]));
                  }
                  count++;
               }
            }
   
         }
      }


     /* do left and right diagonal scaling */
   
      AZ_exchange_bdry(sc_vec, data_org, proc_config);
   
     /* index through rows of matrix */
   
      for (irow = 0; irow < m; irow++) {
   
         /* for all the blocks in a row */
   
         for (k = bpntr[irow]; k < bpntr[irow+1]; k++) {
            icol = bindx[k];
   
            count = 0;
   
            for (j = cpntr[icol]; j < cpntr[icol+1]; j++) {
               for (i = rpntr[irow]; i < rpntr[irow+1]; i++) {
                  val[indx[k]+count] *=  sc_vec[i] * sc_vec[j] ;
                  count++;
               }
            }
   
         }
      }
   
    }
  }

  /* rescale right hand side and solution */

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

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

void AZ_sym_row_sum_scaling_sl(int action, AZ_MATRIX *Amat, 
                              double b[],
                              double x[], int options[], 
			      int proc_config[], struct AZ_SCALING *scaling)

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

  Routine to symmetrically diagonally scale sparse matrix problem; Note: this
  scales the entire matrix problem Ax = b, the routine sym_rescale must be used
  to transform solution back to recover soution to original problem.

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

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

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

  val:             Array containing the nonzero entries of the matrix (see
                    Aztec Usert).

  bindx:           Arrays used for DMSR and DVBR sparse matrix storage (see
                   file  Aztec Usert).

  b:               Right hand side of linear system.

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

  options:         Determines specific solution method and other parameters.

  proc_config:     Machine configuration.  proc_config[AZ_node] is the node
                   number.  proc_config[AZ_N_procs] is the number of processors.

  x:               Current solution vector.

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

{

  /* local variables */

  register int j, k, irow;
  int          N;
  int          i, jcol, j_last, bindx_row;
  double       row_sum, *sc_vec, *val;
  char        *yo = "AZ_sym_row_sum_scaling_sl: ";
  char         label[80];
  int         *bindx, *data_org;



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

  val = Amat->val;
  bindx = Amat->bindx;
  data_org = Amat->data_org;
  if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) {
     if ( (options[AZ_pre_calc] >= AZ_reuse) && (proc_config[AZ_node]  == 0)) {