az_scaling.c

并行解法器,功能强大

    if (c == NULL) {
       (void) fprintf(stderr, "%sERROR: not enough memory for block diagonal\n"
                     "       scaling. Not able to allocate c\n"
                     "       array. Must run a smaller problem\n", yo);
       exit(-1);
    }
    if ( action == AZ_INVSCALE_RHS)  {
       sprintf(Left,  "L");
       sprintf(Upper, "U");
       sprintf(Lower, "L");
       sprintf(Unit, "U");
       for (iblk_row = 0; iblk_row < m; iblk_row++) {
          m1 = rpntr[iblk_row+1] - rpntr[iblk_row];
          d_ival  = d3_indx[d3_bpntr[iblk_row] - d_bpoff] - d_idoff;
          dtrmm_( Left, Upper, None, None, &m1,
                  &ione, &one, &d3_inv[d_ival], &m1, &(b[rpntr[iblk_row]]), &m1,
                  strlen(Left), strlen(Upper), strlen(None), strlen(None));
          dtrmm_( Left, Lower, None, Unit, &m1, &ione,
                  &one, &d3_inv[d_ival], &m1, &(b[rpntr[iblk_row]]), &m1,
                  strlen(Left), strlen(Lower), strlen(None), strlen(Unit));
          dlaswp_( &ione, &(b[rpntr[iblk_row]]), &m1, &ione, &m1, &(ipiv[rpntr[iblk_row]]), &iminus_one );

       }
    }
    if ( (action == AZ_SCALE_MAT_RHS_SOL) || (action == AZ_SCALE_RHS) ) {
       for (iblk_row = 0; iblk_row < m; iblk_row++) {
          m1 = rpntr[iblk_row+1] - rpntr[iblk_row];
          d_ival  = d3_indx[d3_bpntr[iblk_row] - d_bpoff] - d_idoff;
          dgetrs_(None, &m1, &ione, &d3_inv[d_ival], &m1, &(ipiv[rpntr[iblk_row]]), 
                  &(b[rpntr[iblk_row]]), &m1, &info,1);
       }
    }
    AZ_free((void *) c);

  }

} /* AZ_block_diagonal_scaling */

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

void AZ_sym_block_diagonal_scaling(double val[], int indx[], int bindx[],
                                   int rpntr[], int cpntr[], int bpntr[],
                                   double b[], int options[], 
                                   int data_org[], int proc_config[]
                                   /* , struct AZ_SCALING *scaling*/)

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

  Routine to symmetric block Jacobi scale the sparse matrix problem.  Note: this
  scales the matrix, the solution and the right-hand side.

  If the matrix is in MSR format, it is treated as point entry (block size 1)
  and standard Jacobi scaling is performed.  Else, if the matrix is in VBR
  format, block symmetric Jacobi scaling is performed.

  For the DVBR format, the following system is created from Ax = b:

       (trans(inv(L)) A inv(L) (L x) = trans(inv(L)) b

  where L is the Cholesky factor of the block diagonal portion of A.

  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).

  b:               Right hand side of linear system.

  options:         Determines specific solution method and other parameters.

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

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

  action:          Flag which determines whether to scale the matrix or to
                   rescale the solution.

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

{

  /* local variables */

  register int   iblk_row, i, j, k, ival, jblk;
  register int   it, jt, icount, iwork, ilow, ib = 0;
  int            bpoff, idoff;
  int            d_bpoff, d_idoff;
  int            m1, n1, ib1, ib2, idL;
  int            ione = 1, itemp, itemp2;
  int            m, Proc;
  int            tsize;
  int            max_blk;
  static int    *d3_indx, *d3_bindx, *d3_rpntr, *d3_bpntr;
  static double *L;
  double         done = 1.0;
  double        *work;
  char           None[2];
  char          *side = "L", *uplo = "L", *transa = "N", *diag = "N";
  char          *yo = "sym_AZ_block_diagonal_scaling: ";
  char           label[80];

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

  printf("Error: Symmetric block scaling not implemented\n");
  exit(-1);

  /* initialize */

  m    = data_org[AZ_N_int_blk]  + data_org[AZ_N_bord_blk];
  Proc = proc_config[AZ_node];

  sprintf(None, "N");

  /*
   * First, compute the block-diagonal Cholesky factorization, its inverse and
   * its inverse's transpose (if not already computed).
   */

  tsize = 0;
  for (i = 0; i < m; i++)
    tsize += (rpntr[i+1] - rpntr[i]) * (cpntr[i+1] - cpntr[i]);

  sprintf(label,"d3_indx%d",options[AZ_recursion_level]);
  d3_indx  = (int *)    AZ_manage_memory((m+1)*sizeof(int), AZ_ALLOC,
                                         data_org[AZ_name], label, &itemp);
  sprintf(label,"d3_bindx%d",options[AZ_recursion_level]);
  d3_bindx = (int *)    AZ_manage_memory(m*sizeof(int), AZ_ALLOC,
                                         data_org[AZ_name], label, &itemp);
  sprintf(label,"d3_rpntr%d",options[AZ_recursion_level]);
  d3_rpntr = (int *)    AZ_manage_memory((m+1)*sizeof(int), AZ_ALLOC,
                                         data_org[AZ_name], label, &itemp);
  sprintf(label,"d3_bpntr%d",options[AZ_recursion_level]);
  d3_bpntr = (int *)    AZ_manage_memory((m+1)*sizeof(int), AZ_ALLOC,
                                         data_org[AZ_name],
                                         label, &itemp);
  sprintf(label,"L%d",options[AZ_recursion_level]);
  L        = (double *) AZ_manage_memory(tsize*sizeof(double), AZ_ALLOC,
                                         data_org[AZ_name], label, &itemp);

  if ((options[AZ_pre_calc] >= AZ_reuse) && (itemp == 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 necessary, calculate the Cholesky factors (L) of the diagonal blocks and
   * store in L and the d3_ pointer vectors.
   */

  if (options[AZ_pre_calc] <= AZ_recalc) {
    calc_blk_diag_Chol(val, indx, bindx, rpntr, cpntr, bpntr, L, d3_indx,
                       d3_bindx, d3_rpntr, d3_bpntr, data_org);

    /* offset of the first block */

    bpoff = *bpntr;
    idoff = *indx;

    d_bpoff = *d3_bpntr;
    d_idoff = *d3_indx;

    /* symmetrically scale the matrix 'A' using the Cholesky factors */

    max_blk = 0;
    for (i = 0; i < m + data_org[AZ_N_ext_blk] ; i++) {
      if ( cpntr[i+1]-cpntr[i] > max_blk )
        max_blk = cpntr[i+1] - cpntr[i];
    }
    work = (double *) AZ_allocate(max_blk*max_blk*sizeof(double));
    if (work == NULL) {
      (void) fprintf(stderr, "%sERROR: not enough memory for diagonal\n"
                     "      scaling. Not able to allocate work\n"
                     "      array. Must run a smaller problem\n", yo);
      exit(-1);
    }

    /* loop over the block rows */

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

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

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

      /* starting index of current row block */

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

      /* loop over the (block) columns in the current (block) row */

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

        /* the starting point column index of the current block */

        ib1 = cpntr[jblk];

        /* ending point column index of the current block */

        ib2 = cpntr[jblk+1];

        /* number of columns in the current block */

        n1    = ib2 - ib1;
        itemp = m1*n1;

        if (jblk == iblk_row) {

          /* diagonal block => set to identity */

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

          for (i = 0; i < m1; i++) {
            itemp2 = ival + i;
            for (k = 0; k < n1; k++)
              if (i == k)
                val[itemp2 + m1*k] = 1.0;
              else
                val[itemp2 + m1*k] = 0.0;
          }
        }

        else if (jblk < iblk_row) {

          /* lower off-diagonal block */

          if (itemp > max_blk*max_blk) {
            (void) fprintf(stderr, "%sERROR: block size (%d) is too big =>\n",
                           yo, itemp);
            exit(-1);
          }

          /*
           * Fill the work array with the proper block in "A" and take its
           * transpose.
           */

          icount = 0;
          for (it = 0; it < m1; it++) {
            for (jt = 0; jt < n1; jt++) {
              work[icount] = val[ival+icount];
              icount++;
            }
          }

          AZ_dtrans(&m1, &n1, work);

          /* starting index for the diagonal L block for this first operation */

          idL = d3_indx[d3_bpntr[jblk] - d_bpoff] - d_idoff;

          /* perform a backsolve on L*work' = A' to get 'work' array */

          dtrsm_(side, uplo, transa, diag, &m1, &n1, &done, L+idL, &m1, work,
                 &m1, strlen(side), strlen(uplo), strlen(transa), strlen(diag));

          /* need the transpose of the work array */

          AZ_dtrans(&m1, &n1, work);

          /* starting index for the diagonal 'L' block for the next operation */

          idL = d3_indx[d3_bpntr[iblk_row] - d_bpoff] - d_idoff;

          /* perform a backsolve on L*work2 = work */

          dtrsm_(side, uplo, transa, diag, &m1, &n1, &done, L+idL, &m1, work,
                 &m1, strlen(side), strlen(uplo), strlen(transa), strlen(diag));

          /* copy the transpose of this result into the proper block in 'val' */

          iwork  = AZ_get_sym_indx(iblk_row, jblk, indx, bindx, bpntr);
          icount = 0;
          for (i = 0; i < n1; i++)
            for (k = 0; k < m1; k++)
              *(val + iwork + i + k*n1) = *(work + icount++);
        }

        ival += itemp;
      }

      /* scale the RHS */

      idL = d3_indx[d3_bpntr[iblk_row] - d_bpoff] - d_idoff;
      dtrsm_(side, uplo, transa, diag, &m1, &ione, &done, L+idL, &m1, b+ib,
             &m1, strlen(side), strlen(uplo), strlen(transa), strlen(diag));
      ib += m1;
    }

    AZ_free((void *) work);

    /*
     * Copy the lower block to their corresponding upper blocks for symmetry.
     */

    /* loop over the block rows */

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

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

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

      /* starting index of current row block */

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

      /* loop over the (block) columns in the current (block) row */

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

        /* the starting point column index of the current block */

        ib1 = cpntr[jblk];

        /* ending point column index of the current block */

        ib2 = cpntr[jblk+1];

        /* number of columns in the current block */

        n1    = ib2 - ib1;
        itemp = m1*n1;

        if (jblk > iblk_row) {

          /*
           * Lower off-diagonal block => copy to corresponding upper block.
           * NOTE: we do not pass in cpntr as the matrix should be symmetric.
           */

          ilow = AZ_get_sym_indx(iblk_row, jblk, indx, bindx, bpntr);

          icount = 0;
          for (i = 0; i < n1; i++)
            for (k = 0; k < m1; k++)
              val[ilow + i + k*n1] = val[ival + icount++];
        }

        ival += itemp;
      }
    }
  }

} /* sym_AZ_block_diagonal_scaling */

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

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

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

  Routine to left row-sum scale the sparse matrix problem; Note: this scales
  the entire matrix problem Ax = b

  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).

  b:               Right hand side of linear system.

  options:         Determines specific solution method and other parameters.

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

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

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

{

  /* local variables */

  register int indx_ptr, irow, iblk_row, jcol, jblk_col, icol_blk, iblk = 0;
  int          num_blk_rows, num_col_blks, num_blk_cols, I, J, t_index;
  int          m, k, N;
  int          j_last, bindx_row;
  double       row_sum = 0.0, sign = 0.0, inv_row_sum, *sc_vec;
  char        *yo = "AZ_row_sum_scaling: ";
  char        label[80];
  int         *indx, *bindx, *rpntr, *cpntr, *bpntr, *data_org;
  double      *val;

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

  if ( (action == AZ_INVSCALE_SOL) || (action == AZ_SCALE_SOL)) return;

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

    if (((action == AZ_SCALE_RHS) || (action == AZ_INVSCALE_RHS) || 
         (options[AZ_pre_calc] >= AZ_reuse)) && (k == 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)) {
    if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {
      for (irow = 0; irow < N; irow++) {

        /* get row sum */

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