az_tools.c

并行解法器,功能强大

     */

    AZ_sortqlists((char *) cnptr, 0, update_index, N_update, sizeof(int),
                  N_update);

    /* 3) order the external blocks */

    temp = (int *) AZ_allocate((unsigned)(N_external + 1)*sizeof(int));
    if (temp == NULL) {
      (void) fprintf(stderr,
                     "%sERROR: not enough memory to malloc temporary space\n",
                     yo);
      exit(-1);
    }

    for (i = 0; i < N_external; i++)
      temp[extern_index[i] - N_update] = cnptr[i + N_update];

    for (i = 0; i < N_external; i++) cnptr[i + N_update] = temp[i];
    AZ_free((char *) temp);

    /* 4) reconvert cnptr to give pointer information */

    AZ_convert_values_to_ptrs(cnptr, N_update + N_external, 0);
  }
  else {
    (void) fprintf(stderr, "%sERROR: matrix is not MSR or VBR\n", yo);
    exit(-1);
  }

  /*
   * Change column indices (bindx) to reflect new ordering depending depending
   * on whether or not a point is internal or external.
   */

  for (i = start; i < end; i++) {
    if (bindx[i] < N_update) bindx[i] = update_index[bindx[i]];
    else                     bindx[i] = extern_index[bindx[i] - N_update];
  }

  if (option == AZ_EXTERNS) return;

  /* reorder rows */

  if (mat_type == AZ_MSR_MATRIX) {

    /* We move the rows in four steps:
     *  1) sort the first N_update values of 'val'.
     *  2) sort the first N_update values of 'bindx'.
     *     We do this by first converting the ptrs to values
     *     representing the number of nonzero off diagonals.
     *  3) sort the off diagonal column indices.
     *  4) sort the off diagonal matrix nozeros.
     */

    j = bindx[0];
    AZ_convert_ptrs_to_values(bindx, N_update);

    AZ_sortqlists((char *) &(bindx[N_update + 1]), bindx, update_index,
                  end - N_update - 1, sizeof(int), N_update);
    AZ_sortqlists((char *) &(val[N_update + 1]), bindx, update_index,
                  end - N_update - 1, sizeof(double), N_update);
    AZ_sortqlists((char *) val, 0, update_index, N_update, sizeof(double),
                  N_update);
    AZ_sortqlists((char *) bindx, 0, update_index, N_update,
                  sizeof(int), N_update);
    AZ_convert_values_to_ptrs(bindx, N_update, j);
  }
  else {
    val_length  = indx[bnptr[N_update]];
    indx_length = bnptr[N_update];
    AZ_convert_ptrs_to_values(indx, indx_length);

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

    /* move val */

    for (i = 0; i < N_update; i++) {
      temp[i] = 0;
      for (j = bnptr[i]; j < bnptr[i + 1]; j++) temp[i] += indx[j];
    }

    AZ_sortqlists((char *) val, temp, update_index, val_length, sizeof(double),
                  N_update);
    AZ_free((char *) temp);

    AZ_convert_ptrs_to_values(bnptr, N_update);
    AZ_convert_ptrs_to_values(rnptr, N_update);

    /* move indx */

    AZ_sortqlists((char *) indx, bnptr, update_index, indx_length, sizeof(int),
                  N_update);

    /* move bindx */

    AZ_sortqlists((char *) bindx, bnptr, update_index, indx_length, sizeof(int),
                  N_update);

    /* move bnptr */

    AZ_sortqlists((char *) bnptr, 0, update_index, N_update, sizeof(int),
                  N_update);

    /* move rnptr */

    AZ_sortqlists((char *) rnptr, 0, update_index, N_update, sizeof(int),
                  N_update);

    AZ_convert_values_to_ptrs(rnptr, N_update, 0);
    AZ_convert_values_to_ptrs(bnptr, N_update, 0);
    AZ_convert_values_to_ptrs(indx, indx_length, 0);
  }

} /* AZ_reorder_matrix */

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

void AZ_set_message_info(int N_external, int extern_index[], int N_update,
                         int external[], int extern_proc[], int update[],
                         int update_index[], int proc_config[], int cnptr[],
                         int *data_org[], int mat_type)

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

  Perform initializations so that local communication can occur to update the
  external elements. This initialization includes:

   1) determine the number of neighbors to which we send information as well as
      the number of neighbors from which we receive information.  These two
      should be the same. If they are not, we will set up things so that we send
      0 length messages to processors from which we receive (but have no
      information to send them).
   2) determine the total number of unknowns that we must send out.  Using this
      information we allocate space for data_org[].
   3) Initialize data_org[] (see User's Guide) so that it contains the number of
      messages to be sent/received, the node numbers of the processors to which
      we must send and receive, the length of the messages that we must send and
      that we expect to receive from each of the neighbors, and finally, a list
      of the indices of the elements that will be send to other processors (in
      the order that they will be sent).

      NOTE: Implicitly the neighbors are numbered using the ordering of the
      external elements. In particular, the external elements are ordered such
      that those elements updated by the same processor are contiguous. In this
      way, the ordering of the external elements defines an ordering for the
      neighbors.

 The information stored in 'data_org[]' is as follows:

     data_org[AZ_neighbors]         = node number of first neighbor
     data_org[AZ_neighbors+1]       = node number of second neighbor
     data_org[AZ_neighbors+2]       = node number of third  neighbor
         . . .
     data_org[AZ_rec_length]        = # of values to receive from 1st neighbor
     data_org[AZ_rec_length+1]      = # of values to receive from 2nd neighbor
     data_org[AZ_rec_length+2]      = # of values to receive from 3rd neighbor
     . . .
     data_org[AZ_send_length]       = # of values to send to 1st neighbor
     data_org[AZ_send_length+1]     = # of values to send to 2nd neighbor
     data_org[AZ_send_length+2]     = # of values to send to 3rd neighbor
     . . .
     data_org[AZ_send_list]         = elements to be sent to neighbors (starting
                                      with elements to the first neighbor [in
                                      the order that they appear on that
                                      neighbor], followed by elements to the
                                      second neighbor, etc).

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

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

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

  N_external:      Number of external elements on this processor.

  extern_index:    On output, extern_index[i] gives the local numbering of
                   global point 'external[i]'. See User's Guide.

  N_update:        Number of elements updated on this processor.

  external:        List (global indices) of external elements on this node.

  extern_proc:     extern_proc[i] is updating processor of external[i].

  update:          List (global indices) of elements updated on this node.

  update_index:    On output, update_index[i] gives the local numbering of
                   global point 'update[i]'. See User's Guide.

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

  cnptr:           VBR column pointer array (see User's Guide).

  data_org:        Array use to specifiy communication information. See User's
                   Guide.
                   On output, this array contains neighbor and message.

  mat_type:        Indicates whether this is an MSR (= AZ_MSR_MATRIX) or a
                   VBR (= AZ_VBR_MATRIX).

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

{

  /* local variables */

  int   i, j, ii, type, start, cflag, partner, newlength, new_i;
  int   total_to_be_sent, found;
  int  *new_external, *new_extern_proc;
  int  *neighbors, *tempneigh;
  int  *recv_list, *send_list;
  int   tj, end, tt, oldii;
  int   num_recv_neighbors;
  int   num_send_neighbors;
  int  *bins, shift;
  int  *send_ptr, *lens;
  int   proc,nprocs;
  int   firstone,current;
  MPI_AZRequest request[AZ_MAX_NEIGHBORS];  /* Message handle */
  unsigned int length;

  char *yo = "AZ_set_message_info: ";

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

  AZ__MPI_comm_space_ok();
  proc   = proc_config[AZ_node];
  nprocs = proc_config[AZ_N_procs];

  neighbors = (int *) AZ_allocate((unsigned) nprocs*sizeof(int));
  tempneigh = (int *) AZ_allocate((unsigned) nprocs*sizeof(int));

  /* Produce a list of the external updating processors corresponding */
  /* to each external point in the order given by 'extern_index[]'    */

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

  for (i = 0; i < N_external; i++)
    new_extern_proc[extern_index[i] - N_update] = extern_proc[i];

  /*
   * Count the number of neighbors from which we receive information to update
   * our external elements. Additionally, fill the array neighbors in the
   * following way:
   *      neighbors[i] = 0   ==>  No external elements are updated by
   *                              processor i.
   *      neighbors[i] = x   ==>  (x-1)/nprocs elements are updated from
   *                              processor i.
   */

  num_recv_neighbors = 0;
  length             = 1;

  for (i = 0; i < N_external; i++) {
    if (neighbors[new_extern_proc[i]] == 0) {
      num_recv_neighbors++;
      neighbors[new_extern_proc[i]] = 1;
    }

    if (mat_type != AZ_VBR_MATRIX)
      neighbors[new_extern_proc[i]] += nprocs;
    else {
      neighbors[new_extern_proc[i]] += (nprocs * (cnptr[i + 1 + N_update]
                                                  - cnptr[i + N_update]));

    }
  }

  /*
   * Make a list of the neighbors that will send information to update our
   * external elements (in the order that we will receive this information).
   */

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

  j = 0;
  recv_list[j++] = new_extern_proc[0];
  for (i = 1; i < N_external; i++) {
    if (new_extern_proc[i - 1] != new_extern_proc[i]) {
      recv_list[j++] = new_extern_proc[i];
    }
  }

  /* sum over all processors all the neighbors arrays */

  AZ_gsum_vec_int(neighbors, tempneigh, proc_config[AZ_N_procs], proc_config);

  /* decode the combined 'neighbors' array from all the processors */

  num_send_neighbors = neighbors[proc] % nprocs;

  /* decode 'neighbors[proc] to deduce total number of elements we must send */

  total_to_be_sent = (neighbors[proc] - num_send_neighbors) / nprocs;

  AZ_free((char *) neighbors);
  AZ_free((char *) tempneigh);

  /* send a 0 length message to each of our recv neighbors */

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


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


  /* first post receives */

  length = 0;
  for (i = 0; i < num_send_neighbors; i++) {
    partner = -1;
    mdwrap_iread((void *) external, length, &partner, &type, request+i);
  }

  /* send messages */

  for (i = 0; i < num_recv_neighbors; i++) {
    mdwrap_write((void *) external, 0, recv_list[i], type, &cflag);
  }

  /*
   * Receive message from each send neighbor to construct 'send_list'.
   */

  length = 0;
  for (i = 0; i < num_send_neighbors; i++) {
    partner = -1;
    mdwrap_wait((void *) external, length, &partner, &type, &cflag, request+i);
    if (partner != -1) send_list[i] = partner;
  }



  /*
   * Compare the two lists. In most cases they should be the same.  However, if
   * they are not add new entries to the recv list that are in the send list
   * (but not already in the recv list).
   */

  for (j = 0; j < num_send_neighbors; j++) {
    found = 0;
    for (i = 0; i < num_recv_neighbors; i++) {
      if (recv_list[i] == send_list[j])
        found = 1;
    }

    if (found == 0) {
      recv_list[num_recv_neighbors] = send_list[j];
      (num_recv_neighbors)++;
    }
  }

  AZ_free((char *) send_list);
  num_send_neighbors = num_recv_neighbors;

  /* create data_org array */

  if (!AZ_using_fortran) {
    *data_org = (int *) AZ_allocate(((unsigned) total_to_be_sent + AZ_send_list)
                               *sizeof(int));
  }
  if (*data_org == NULL) {
    (void) fprintf(stderr, "%sERROR: Not enough dynamic space.\n", yo);
    exit(-1);
  }
  if (!AZ_using_fortran) {
     for (i = 0 ; i < total_to_be_sent + AZ_send_list; i++ ) 
        (*data_org)[i] = 0;
  }

  (*data_org)[AZ_total_send] = total_to_be_sent;
  send_ptr = &((*data_org)[AZ_send_list]);

  /*