az_tools.c

并行解法器,功能强大

   * Create 'new_external' which explicitly put the external elements in the
   * order given by 'extern_index'
   */

  new_external = (int *) AZ_allocate((unsigned) (N_external+1)*sizeof(int));
  if (new_external == NULL) {
    (void) fprintf(stderr, "%sERROR: Not enough dynamic space.\n", yo);
    exit(-1);
  }

  for (i = 0; i < N_external; i++) {
    new_external[extern_index[i] - N_update] = external[i];
  }

  /*
   * Send each processor the global index list of the external elements in the
   * order that I will want to receive them when updating my external elements
   */



  type            = AZ_sys_msg_type;
  AZ_sys_msg_type = (AZ_sys_msg_type+1-AZ_MSG_TYPE) % AZ_NUM_MSGS + AZ_MSG_TYPE;

  lens = (int *) AZ_allocate((num_recv_neighbors+1)*sizeof(int));
  if (lens == NULL) {
    (void) fprintf(stderr, "%sERROR: Not enough dynamic space.\n", yo);
    exit(-1);
  }

  for (i = 0; i < num_recv_neighbors; i++) {
    length     = sizeof(int);
    partner    = recv_list[i];
    mdwrap_iread((void *) &(lens[i]), length, &partner, &type, request+i);
  }


  j = 0;
  for (i = 0; i < num_recv_neighbors; i++) {
    start  = j;
    newlength = 0;

    /* go through list of external elements until updating processor changes */

    while ((j < N_external) && (new_extern_proc[j] == recv_list[i])) {

      if (mat_type != AZ_VBR_MATRIX) newlength++;
      else newlength += (cnptr[j + 1 + N_update] - cnptr[j + N_update]);

      j++;
      if (j == N_external) break;
    }

    (*data_org)[AZ_rec_length+i] = newlength;
    (*data_org)[AZ_neighbors+i]  = recv_list[i];

    length = j - start;
    mdwrap_write((void *) &length, sizeof(int), recv_list[i], type, &cflag);
  }

  /* receive from each neighbor the global index list of external ele */

  for (i = 0; i < num_recv_neighbors; i++) {
    length     = sizeof(int);
    partner    = (*data_org)[AZ_neighbors+i];
    mdwrap_wait((void *) &(lens[i]), length, &partner, &type, &cflag,request+i);
    (*data_org)[AZ_send_length + i] = lens[i];
  }
  AZ_free(lens);







  type            = AZ_sys_msg_type;
  AZ_sys_msg_type = (AZ_sys_msg_type+1-AZ_MSG_TYPE) % AZ_NUM_MSGS + AZ_MSG_TYPE;


  j = 0;
  for (i = 0; i < num_recv_neighbors; i++) {
    length     = (  (*data_org)[AZ_send_length + i] ) * sizeof(int);
    partner    = (*data_org)[AZ_neighbors+i];
    mdwrap_iread((void *) &(send_ptr[j]), length, &partner, &type, request+i);
    j += (*data_org)[AZ_send_length + i];
  }

  j = 0;
  for (i = 0; i < num_recv_neighbors; i++) {
    start  = j;
    newlength = 0;

    /* go through list of external elements until updating processor changes */

    while ((j < N_external) && (new_extern_proc[j] == recv_list[i])) {

      if (mat_type != AZ_VBR_MATRIX) newlength++;
      else newlength += (cnptr[j + 1 + N_update] - cnptr[j + N_update]);

      j++;
      if (j == N_external) break;
    }
    mdwrap_write((void *) &(new_external[start]), (j-start)* sizeof(int),
             recv_list[i], type, &cflag);
  }

  /* receive from each neighbor the global index list of external ele */

  j = 0;
  for (i = 0; i < num_recv_neighbors; i++) {
    length     = (  (*data_org)[AZ_send_length + i] ) * sizeof(int);
    partner    = (*data_org)[AZ_neighbors+i];
    mdwrap_wait((void *) &(send_ptr[j]), length, &partner, &type, &cflag,
		 request+i);
    j += (*data_org)[AZ_send_length + i];
  }



  tj = 0;
  if (mat_type == AZ_VBR_MATRIX) {
    (*data_org)[AZ_matrix_type] = AZ_VBR_MATRIX;

    /* move the information to the back of the send_ptr */

    for (i = j - 1; i >= 0; i--) send_ptr[total_to_be_sent+i-j] = send_ptr[i];
    tj = total_to_be_sent - j;
  }
  else (*data_org)[AZ_matrix_type] = AZ_MSR_MATRIX;

  total_to_be_sent = j;

  /* replace global indices by local indices */

  oldii     = ii = 0;
  tt        = tj;
  firstone  = AZ_send_length;
  current   = firstone;

  while (((*data_org)[current] == 0) && (current-firstone < num_recv_neighbors))
    current++;

  bins = (int *) AZ_allocate((unsigned) (N_update / 4 + 10)*sizeof(int));
  if (bins == NULL) {
    (void) fprintf(stderr, "%sERROR: Not enough dynamic space.\n", yo);
    exit(-1);
  }
  for (i = 0 ; i < N_update / 4 + 10 ; i++ ) bins[i] = 0;

  AZ_init_quick_find(update, N_update, &shift, bins);
  for (i = tj; i < total_to_be_sent + tj; i++) {

    /*
     * j = AZ_find_index(send_ptr[i], update, N_update);
     */

    j = AZ_quick_find(send_ptr[i], update, N_update, shift, bins);
    if (j == -1) {
      (void) fprintf(stderr, "%sERROR: A point (%d) requested from processor "
                     "%d\nwas not found among the update ele of proc %d\n", yo,
                     send_ptr[i], (*data_org)[AZ_neighbors+current-firstone],
                     proc);
      exit(-1);
    }

    if (mat_type != AZ_VBR_MATRIX) {
      send_ptr[ii] = update_index[j];
      ii +=  1;
    }
    else {
      start = cnptr[update_index[j]];
      end   = cnptr[update_index[j] + 1] - 1;
      for (new_i = start; new_i <= end; new_i++) send_ptr[ii++] = new_i;

      if (i - tt + 1 == (*data_org)[current]) {
        (*data_org)[current] = ii - oldii;
        current++;

        while ( ((*data_org)[current] == 0) &&
                (current-firstone < num_recv_neighbors)) current++;
        oldii = ii;
        tt    = i + 1;
      }
    }
  }

  AZ_free((char *) bins);
  AZ_free((char *) recv_list);
  AZ_free((char *) new_external);
  AZ_free((char *) new_extern_proc);

  (*data_org)[AZ_N_neigh] = num_send_neighbors;

} /* AZ_set_message_info */

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

void AZ_convert_values_to_ptrs(int array[], int length, int start)

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

  Change 'array[]' so that on exit
     1) array[0] = start
     2) array[i+1] - array[i] = value of array[i] on entry to this routine.

  Author:          Ray S. Tuminaro, SNL, 1422
  =======

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

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

  array:           On entry to this routine and array[0] = start.
                   On output, array[i+1] - array[i] = value of array[i].

  length:          Length of array[].

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

{

  /* local variables */

  int i;

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

  for (i = 1; i < length; i++) array[i] += array[i - 1];
  for (i = length; i > 0; i--)  array[i]  = array[i - 1] + start;

  array[0] = start;

} /* AZ_convert_values_to_ptrs */

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

void AZ_convert_ptrs_to_values(int array[], int length)

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

  Author:          Ray S. Tuminaro, SNL, 1422
  =======

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

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

  array:           On entry to this routine and array[0] = start.
                   On output, array[i+1] - array[i] = value of array[i].

  length:          Length of array[].

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

{
  int i;

  for (i = 0; i < length; i++) array[i] = array[i + 1] - array[i];
}

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

void AZ_msr2vbr(double val[], int indx[], int rnptr[], int cnptr[], int bnptr[],
                int bindx[], int msr_bindx[], double msr_val[],
                int total_blk_rows, int total_blk_cols, int blk_space,
                int nz_space, int blk_type)

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

  Convert the MSR matrix defined in [msr_val,msr_bindx] to a VBR matrix.

  Author:          Ray S. Tuminaro, SNL, 1422
  =======

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

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

  val,
  rnptr,
  bindx,
  indx,
  bnptr,
  cnptr:           Sparse matrix arrays. See User's Guide.
                   On input, the matrix corresponds to the initial ordering
                   (e.g. row i corresponds to global row update[i]).
                   On output, the matrix rows and columns are renumbered to
                   correspond to the ordering given by 'update_index' and
                   'extern_index'. (e.g. global row update[i] now appears
                   as row update_index[i] in the matrix).

  msr_val,
  msr_bindx:       On input, MSR matrix to be converted to VBR.
                   See User's Guide.

  total_blk_rows:  Number of block rows in resulting local VBR matrix.

  total_blk_cols:  Number of block columns in resulting local VBR matrix.

  blk_space:       Length of storage allocated for bindx[] and indx[]. An error
                   message will be printed if we try to write past these arrays.

  nz_space:        Length of storage allocated for val[]. An error message will
                   be printed if we try to write past this array.

  blk_type:        If blk_type > 0, blk_type indicates that all block rows (and
                   colunns) have the same size given by 'blk_type'. If
                   blk_type < 0, the block rows have different sizes.

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

{

  /* local variables */

  int therow, thecol;
  int i, j;

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

  for (i = 0; i < total_blk_rows; i++) rnptr[i] = cnptr[i];

  AZ_convert_values_to_ptrs(rnptr, total_blk_rows, 0);
  AZ_convert_values_to_ptrs(cnptr, total_blk_cols, 0);

  indx[0] = bnptr[0] = 0;

  /* go through each block row */

  for (i = 0; i < total_blk_rows; i++) {
    bnptr[i + 1] = bnptr[i];

    for (therow = rnptr[i]; therow < rnptr[i + 1]; therow++) {

      /* add the diagonal entry */

      thecol = therow;
      AZ_add_new_ele(cnptr, therow, i, bindx, bnptr, indx, val, therow,
                     msr_val[therow], total_blk_cols, blk_space, nz_space,
                     blk_type);

      /* add off diagonal entries */

      for (j = msr_bindx[therow]; j < msr_bindx[therow + 1]; j++) {
        thecol = msr_bindx[j];
        AZ_add_new_ele(cnptr, thecol, i, bindx, bnptr, indx, val, therow,
                       msr_val[j], total_blk_cols, blk_space, nz_space,
                       blk_type);
      }
    }
  }

} /* AZ_msr2vbr */

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

int AZ_find_block_col(int cnptr[], int column, int max_blocks, int blk_size)

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

  Return the local index of the block column witch contains the point column
  given by local index 'column'.

  Author:          Ray S. Tuminaro, SNL, 1422
  =======

  Return code:     int, local index of the block column.
  ============

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

  cnptr:           cnptr[0] = 0 and cnptr[i+1] - cnptr[i] gives the column
                   dimension of global block column
                     'update[i]' if  i <  N_update
                     'external[k]' if  i >= N_update
                   where k = i - N_update.

  column:          Local column index of the column for which we are trying to
                   find the block column containing it.

  blk_size:        blk_size > 0 ==> all the blocks are the same size so we can
                   use a shortcut in computing the block column index.
                   blk_size = 0 ==> short cut not used.

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