az_ilu_util.c

并行解法器,功能强大

  bpntr_old:       The i'th element of bpntr points to the first block entry of
                   the i'th row in bindx_old. The last element is the number of
                   nonzero blocks of matrix plus one.

  bpntr_new:       On, output, the i'th element of bpntr points to the first
                   block entry of the i'th row in bindx_new. The last element
                   is the number of nonzero blocks of matrix plus one.

  diag_block:      On output, array of size M points on each diagonal block.

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

{

  /* local variables */

  int  i, kk, j;
  int  num_blks_row_old, num_blks_row_new;
  int *sort, ptr, compt, temp;

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

  /* Allocate work space */

  sort = (int *) AZ_allocate(sizeof(int) * (M));
  if (sort == NULL) {
    (void) fprintf(stderr, "Error: not enough memory inside order_parallel\n"
                   "       must run a smaller problem\n");
    exit(-1);
  }

  /* Initialize */

  for (i = 0; i < M; i++) diag_block[i] = -1;
  bpntr_new[0] = bindx_new[0] = 0;

  for (i = 0; i < M; i++) {     /* loop over the rows */

    /* constructs the new array of block column indices in this row */

    num_blks_row_old = bpntr_old[i+1] - bpntr_old[i];

    /* copy old block column index array and then sort it this defines the new
       block column index array */

    for (j = 0; j < num_blks_row_old; j++)
      bindx_new[bpntr_new[i] + j] = bindx_old[bpntr_old[i] + j];

    sort_blk_col_indx(num_blks_row_old, &bindx_new[bpntr_new[i]], sort);

    /* Count the blocks that multiply internal and border unknowns */

    num_blks_row_new = 0;
    for (j = 0; j < num_blks_row_old; j++) {
      if (bindx_new[bpntr_new[i] + j] >= M) break;

      num_blks_row_new++;
    }

    bpntr_new[i+1] = bpntr_new[i] + num_blks_row_new;

    /* for each block of this row compute the new index vector and construct the
       pointers to the diagonal block i */

    for (kk = bpntr_new[i]; kk < bpntr_new[i+1]; kk++) {

      /* Define new indx vector */

      if (kk - bpntr_new[i] == 0) {
        indx_new[0] = indx_old[0];
      }
      else {
        temp         = sort[kk-1-bpntr_old[i]] + bpntr_old[i];
        indx_new[kk] = indx_new[kk-1] + (indx_old[temp+1] - indx_old[temp]);
      }

      /* Get diagonal block pointers */

      if (bindx_new[kk] == i) diag_block[i] = kk;
    }

    /* constructs the new array containing the entries of the matrix */

    for (kk = bpntr_new[i]; kk < bpntr_new[i+1]; kk++) {
      ptr   = indx_old[sort[kk-bpntr_old[i]] + bpntr_old[i]];
      compt = -1;

      for (j = indx_new[kk]; j < indx_new[kk+1]; j++) {
        compt++;
        val_new[j] = val_old[ptr + compt];
      }
    }
  }

  AZ_free((void *) sort);

} /* order_parallel */

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

void get_diag(int M, int *bindx, int *bpntr, int *diag_block)

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

  Constructs the array of pointers: diag_block to the diagonal blocks.  Returns
  diag_block[i] = -1 if no diagonal block has been found at the ith row.

  Author:          Lydie Prevost, SNL, 1421
  =======

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

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

  M:               Number of (block) rows in the matrix and L.

  bindx:           Contains the block column indices of the non-zero block
                   entries.

  bpntr:           The i'th element of bpntr points to the first block entry of
                   the i'th row in bindx. The last element is the number of
                   nonzero blocks of matrix plus one.

  diag_block:      On output, array of size M points on each diagonal block.

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

{
  int i, kk;

  for (i = 0; i < M; i++) diag_block[i] = -1;

  for (i = 0; i < M; i++) {
    for (kk = bpntr[i]; kk < bpntr[i+1]; kk++) {
      if (bindx[kk] == i)
        diag_block[i] = kk;
    }
  }

} /* get_diag */




#define add_scaled_row(row)                             \
    accum_col[Ncols++] = row;                           \
    for (kk = bindx[row] ; kk <= last[row]; kk++ ) {    \
       col = bindx[kk];                                 \
       if (col >= 0) accum_col[Ncols++] = col;          \
    }

/***************************************************************************/
/***************************************************************************/
/***************************************************************************/
void AZ_MSR_mult_patterns(int bindx[], int N, int last[], int bindx_length,
	int *accum_col)
{
/*
 * Multiply two MSR matrices (actually only the sparsity patterns)
 * and store the resulting sparse matrix. That is,
 *                 C = A * B
 *
 * Note: A, B and C are all stored in the same array (bindx). Specifically,
 *       all the elements in bindx[] correspond to A. However, some columns
 *       are encoded as negative numbers. To obtain the real column number
 *       one must do the following:  real_column = -2 - bindx[k]. The matrix,
 *       B, corresponds only to the positive column indices in bindx[].
 *
 * Parameters
 * ======
 *    bindx[]       On input, bindx[] holds the two input matrices as 
 *                  described above. On output, bindx[] holds the matrix
 *                  corresponding to the product of the two input matrices.
 *                  However, any matrix entry which was not in the matrix B
 *                  is encoded as a negative column number (see above).
 *
 *    N             On input, size of the matrices to multiply.
 *
 *    last          On input, uninitialized workspace of size N.
 *
 *    bindx_length  On input, the size of the array allocated for bindx.
 *                  In general, the number of elements in bindx[] will grow and
 *                  so bindx[] should contain additional room for this growth.
 *
 *    accum_col     On input, uninitialized workspace of size 2*N.
 */

   int    i, k, kk, next_nz, Ncols;
   int    row, col, *signs;
   int    start_row, end_row, first_one, orig_col;
   int largest_col, smallest_col;


   /* move off-diagonals to the back of the array */
   /* and  initialize start/end ptrs              */

   kk = bindx_length-1;
   end_row   = bindx[N]-1; 
   for (i = N-1 ; i >= 0 ; i--) {
      start_row = bindx[i];
      last[i] = kk;
      for (k = end_row; k >= start_row; k-- ) bindx[kk--] = bindx[k];
      end_row   = start_row - 1; 
      bindx[i] = kk+1;
   }

   /* initialize the arrays */

   for (i = 0 ; i < 2*N; i++) accum_col[i] = 0;
   signs = &(accum_col[N]);

   next_nz   = N+1;

   largest_col = 0;
   for (i = 0 ; i < N ; i++) {
      if (largest_col < i) largest_col = i;
      Ncols = 0;
      add_scaled_row(i);

      for (k = bindx[i] ; k <= last[i]; k++ ) {
         if (Ncols >= N) {
            AZ_sort(accum_col, Ncols, NULL, NULL);
            AZ_rm_duplicates(accum_col, &Ncols);
         }
         row        = bindx[k];
         if (row < 0) row = -2 - row;
         add_scaled_row(row);
      }
      AZ_sort(accum_col, Ncols, NULL, NULL);
      AZ_rm_duplicates(accum_col, &Ncols);

      for (k = 0 ; k < Ncols; k++) signs[accum_col[k]] = -1;

      /* compute the smallest and largest column */
      /* that could appear in a factorization    */

      smallest_col = i;
      first_one    = next_nz;
      if (bindx[i] <= last[i]) {
         orig_col = bindx[bindx[i]];
         if (orig_col < 0)  orig_col = -2 - orig_col; 
         if (smallest_col > orig_col) smallest_col = orig_col;

         orig_col = bindx[last[i]];
         if (orig_col < 0)  orig_col = -2 - orig_col; 
         if (largest_col < orig_col) largest_col = orig_col;
      }

      /* record sign of column to be stored  */

      for (k = bindx[i]; k <= last[i]; k++) {
         orig_col = bindx[k];
         if (orig_col >= 0)  signs[orig_col] = 1;    
      }

      if (next_nz+Ncols-2 > last[i]) {
         printf("Not enough room for the larger sparsity pattern\n");
         exit(1);
      }
      for (k = 0 ; k < Ncols ; k++) {
         col = accum_col[k];
         if (col != i) {
            if (signs[col] == -1) col = -2 - col;

            if ((accum_col[k] <= largest_col) && 
                (accum_col[k] >= smallest_col ) ) 
                   bindx[next_nz++] = col;
         }
      }
      bindx[i] = first_one;
      last[i]  = next_nz - 1;
   }
   bindx[N] = last[N-1]+1;

}
       
/***************************************************************************/
/***************************************************************************/
/***************************************************************************/
void AZ_rm_duplicates(int array[], int *N)
{
/*
 * remove any duplicates that might appear in the SORTED
 * array 'array'.
 *
 */
  int k, kk;

  kk = 0;
  for (k = 1; k < *N; k++) {
    if (array[kk] != array[k]) {
      kk++;
      array[kk] = array[k];
    }
  }
  if (*N != 0) kk++;

  *N= kk;
}

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

int AZ_fill_sparsity_pattern(struct context *context, int ifill, int bindx[], 
			     double val[], int N)
{
/*
 * Expand the MSR matrix (bindx,val) by adding zeros such that the new
 * matrix has the same sparsity pattern as when ILU(ifill) is performed
 * on the original matrix.
 *
 */
int   length, last_one, flag, i;
int   *work1, *work2;
double temp;

   length = context->N_nz_allocated;
   last_one = bindx[N]-1;

   /* allocate the work space arrays. If there is enough space in */
   /* val[] use that for one of the two work arrays               */

   if ( (length - last_one - 2)*sizeof(double) > (N+1)*sizeof(int) ) {
      flag = 0;
      work1 = (int *) &(val[last_one+1]);
   }
   else {
      work1 = (int *) AZ_allocate((N+1)*sizeof(int));
      flag = 1;
   }
   work2 = (int *) AZ_allocate(2*(N+1)*sizeof(int));
   if (work2 == NULL) AZ_perror("Out of space in ilu.\n");

   /* Take the power of the matrix. That is, A = A^ifill. */

   for (i = 0 ; i < ifill; i++)
      AZ_MSR_mult_patterns(bindx, N, work1, length,work2);

   AZ_free(work2);
   if (flag) AZ_free(work1);

   /* Move the nonzero values into their proper location in the  */
   /* new expanded matrix. Also, decode any columns which appear */
   /* as negative numbers (see AZ_MSR_mult_patterns()).          */

   for (i = bindx[N]-1; i >=  bindx[0]; i--) {
      if (bindx[i] >= 0) {
         temp            = val[last_one];
         val[last_one--] = 0.0;
         val[i]          = temp;
      }
      else {
         bindx[i] = -2 - bindx[i];
         val[i] = 0.0;
      }
   }

   return(bindx[N]);
}

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

void AZ_sort_msr(int bindx[], double val[], int N)
{
/*
 * Sort the column numbers within each MSR row.
 */

   int i, start, last;

   for (i = 0 ; i < N ; i++ ) {
      start = bindx[i];
      last  = bindx[i+1];
      AZ_sort( &(bindx[start]), last - start , NULL, &(val[start]));
   }
}