az_tools.c

并行解法器,功能强大

  if (option == AZ_GLOBAL)
    total_ele = AZ_gsum_int(N_update, proc_config);

  /* first check that the number of offdiagonal nonzeros is positive */
  /* also compute the largest number of nonzeros in a row */

  largest = -1;
  for (i = 0; i < N_update; i++) {
    num = bnptr[i + 1] - bnptr[i];

    if (num > largest) largest = num;

    if (num < 0) {
      (void) fprintf(stderr, "%sERROR on proc %d: Number of nonzeros ", yo,
                     proc);
      (void) fprintf(stderr, "blocks in row %d = (%d - %d) = %d\n", i,
                     bnptr[i + 1], bnptr[i], bnptr[i + 1] - bnptr[i]);
      (void) fprintf(stderr, "are negative inside AZ_vbr_check()?\n");
    }
  }

  if (option == AZ_LOCAL) {
    if (largest > N_update + N_external) {
      (void) fprintf(stderr, "ERROR on proc %d: Number of blocks ", proc);
      (void) fprintf(stderr, "in a row (%d) exceeds the number of ", largest);
      (void) fprintf(stderr, "blocks on the processor %d\n",
                     N_update + N_external);
    }
  }
  else {
    if (largest > total_ele) {
      (void) fprintf(stderr, "ERROR on proc %d: Number of blocks ", proc);
      (void) fprintf(stderr, "in row %d exceeds the total number ", largest);
      (void) fprintf(stderr, "of blocks in simulation %d\n", total_ele);
    }
  }

  largest = AZ_gmax_int(largest, proc_config);
  if (proc == 0) {
    (void) fprintf(stderr, "The max number of nonzero blocks in a row = %d\n",
                   largest);
  }

  /* compute the largest column index */

  largest = -1;
  for (i = 0; i < bnptr[N_update]; i++) {
    if (bindx[i] < 0) {
      (void) fprintf(stderr, "Warning on proc %d: Negative ", proc);
      (void) fprintf(stderr, "column (%d)= %d\n", i, bindx[i]);
    }

    if (bindx[i] > largest) largest = bindx[i];
  }

  if (option == AZ_LOCAL) {
    if (largest > N_update + N_external) {
      (void) fprintf(stderr, "Warning on proc %d: Column referenced ", proc);
      (void) fprintf(stderr, "(%d) that does not exist\n", largest);
      (void) fprintf(stderr, "    # of blocks update on proc = %d\n",
                     N_update);
      (void) fprintf(stderr, "    # of external blocks = %d\n", N_external);
    }
  }
  else {
    if (largest > total_ele) {
      (void) fprintf(stderr, "Warning on proc %d: Column referenced ", proc);
      (void) fprintf(stderr, "(%d) that is larger than the total ", largest);
      (void) fprintf(stderr, "number of blocks in simulation %d\n", total_ele);
    }
  }

  largest = AZ_gmax_int(largest, proc_config);
  if (proc == 0) {
    (void) fprintf(stderr, "The largest block column index is = %d\n", largest);
  }

  /* check rnptr */

  for (i = 0; i <= N_update; i++) {
    if (rnptr[i] != cnptr[i]) {
      (void) fprintf(stderr,
                     "ERROR on proc %d: rnptr(%d) != cnptr(%d) (%d vs %d)\n",
                     proc, i, i, rnptr[i], cnptr[i]);
    }
  }

  /* check cnptr */

  largest = -1;
  for (i = 0; i < N_update + N_external; i++) {
    num = cnptr[i + 1] - cnptr[i];
    if (num > largest) largest = num;

    if (num <= 0) {
      (void) fprintf(stderr, "ERROR on proc %d: Block Size for ", proc);
      (void) fprintf(stderr, "column %d = (%d - %d) = %d\n", i, cnptr[i + 1],
                     cnptr[i], cnptr[i + 1] - cnptr[i]);
    }
  }

  largest = AZ_gmax_int(largest, proc_config);
  if (proc == 0) {
    (void) fprintf(stderr, "The largest column block size is = %d\n", largest);
  }

} /* AZ_check_vbr */

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

int AZ_find_closest_not_larger(int key, int list[], int length)

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

  Find the closest number to 'key' in 'list' which is not greater than key and
  return the index number.

  On exit, AZ_find_index() returns: i => list[i] = key or list[i] is closest
  number smaller than key.

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

  Return code:     int, see explanation above.
  ============

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

  key:             Element to be search for in list.

  list:            List (assumed to be in ascending order) to be searched.

  length:          Length of list.

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

{

  /* local variables */

  int mid, start, end;

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

  if (length == 0) return -1;

  start = 0;
  end   = length - 1;

  while (end - start > 1) {
    mid = (start + end) / 2;
    if (list[mid] > key) end = mid;
    else start = mid;
  }

  if (list[end] > key) return start;
  else return end;

} /* AZ_find_closest_not_larger */

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

void AZ_find_procs_for_externs(int N_update, int update[],
                               int external[], int N_external,
                               int proc_config[], int **extern_proc)

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

  Determine which processors update our external elements. This is done in the
  following way:

  1) the external elements are split into subpieces (to avoid using too much
     storage or overflowing message buffers) on each processor.
  2) the current subpieces for all the processors are concatentated to
     form one large list of external elements with the name of the
     processor who contributed each subpiece.
  3) The subpieces list is searched and a neighbor list is created.This list
     corresponds to the processors who contributed subpieces containing
     an external point that is updated by this processor.
  4) A message is sent out to the neighbors who in turn send all of
     their external elements back to this processor.
  5) We recieve the external elements from the neighbors, determine which
     ones we update, and send them back to the neighbors.
  6) Finally, 'extern_proc' is set for those external elements which were found.
  7) Steps 1-6 are repeated until all the external elements are discovered.

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

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

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

  N_update:        Number of elements updated on this processor.

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

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

  N_external:      Number of external elements on this processor.

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

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

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

{

  /* local variables */

  int  *extern_mess, *neigh_list, *tempneigh;
  int  i, j, k = 0, ii, jj, kk;
  int  num_send_neighbors, num_recv_neighbors;
  int  type, cflag, partner, length;
  int  num_sublist_elements, all_sublists_length;
  int  found = 0, flag = 0;
  int  stride = 0;
  int  first_extern;
  int  buf_size;
  int  nprocs, proc;
  MPI_AZRequest request[AZ_MAX_NEIGHBORS];  /* Message handle */
  MPI_AZRequest send_request[2*AZ_MAX_NEIGHBORS];
  int num_sent, exch_neighbor;
  unsigned int exch_length = 0;
  int to_recv;
int oldk, iii, *ext_copy, **exch_buffers, exch_count, 
    *newbuffer, send_counter;


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

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

  if (!AZ_using_fortran) {
    *extern_proc = (int *) AZ_allocate((unsigned) (N_external+1)*sizeof(int));
    if ( *extern_proc == NULL) {
       (void) fprintf(stderr, "%d: Not enough memory for extern_proc\n",proc);
       exit(-1);
     }
  }

  for (i = 0; i < N_external; i++) (*extern_proc)[i] = -1;

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

  if ( (neigh_list == NULL) || (tempneigh == NULL)) {
       (void) fprintf(stderr, " Not enough memory to find externals\n");
       exit(-1);
  }

  /*
   * 'extern_mess' must hold the concatentated sublists as well as all the
   * externals on a processor.
   */

  if (N_external > AZ_MAX_MESSAGE_SIZE) i = AZ_MAX_MESSAGE_SIZE;
  else                                  i = N_external;

  buf_size = AZ_gmax_int(i + 1, proc_config);
  buf_size = max(buf_size, (AZ_TEST_ELE+1) * nprocs);

  extern_mess = (int *) AZ_allocate((unsigned) buf_size*sizeof(int));
  if (extern_mess == NULL) {
    (void) fprintf(stderr, " Not enough memory to find the external procs\n");
    exit(-1);
  }
  for (i = 0 ; i < buf_size ; i++) extern_mess[i] = 0;

  first_extern = 0;

  while (flag == 0) {

    /*
     * Initialize  extern_mess to a sublist of external (of size
     * AZ_TEST_ELE). Additionally, encode the processor number at the head of
     * this list as -(proc+1)
     */

    extern_mess[0]       = -proc - 1;
    num_sublist_elements = min(AZ_TEST_ELE, N_external - first_extern);

    if (num_sublist_elements != 0)
      stride = (N_external - first_extern) / num_sublist_elements;

    for (i = 0; i < num_sublist_elements; i++)
      extern_mess[i + 1] = external[i * stride + first_extern];

    all_sublists_length = num_sublist_elements + 1;

    /*
     * Append together all the sublists found on all processors.  After this
     * operation each processor has a list of all the external elements
     * contained in all the sublists.
     */

    AZ_gappend(extern_mess, &all_sublists_length, buf_size, proc_config);

    /*
     * Figure out which external elements we update by searching all the
     * elements in the concatinated sublists.  Store the processor numbers in
     * the front of 'extern_mess' corresponding to the external elements which
     * we update.
     */

    for (i = 0; i < nprocs; i++) neigh_list[i] = 0;
    num_send_neighbors = 0;

    for (i = 0; i < all_sublists_length; i++) {

      /*
       * Is extern_mess[i] in the current sequence or is it the head of a new
       * subsequence.
       */

      if (extern_mess[i] >= 0) {
        if (found == -1) {      /* no point in current subsequence */

          /* updated by this processor */

          found = AZ_find_index(extern_mess[i], update, N_update);
          if (found != -1) {

            /* this point updated by this processor */

            neigh_list[-k - 1]                = 1;
            extern_mess[num_send_neighbors++] = k;
          }
        }
      }
      else {
        found = -1;
        k     = extern_mess[i]; /* k = processor of current sublist */
      }
    }

    /* find the total number of neighbors from who we receive */

    AZ_gsum_vec_int(neigh_list, tempneigh, proc_config[AZ_N_procs],proc_config);
    num_recv_neighbors = neigh_list[proc];

    if ((num_sublist_elements != 0) && (num_recv_neighbors == 0)) {
      (void) fprintf(stderr, "A column referenced on one processor was not\n");
      (void) fprintf(stderr,
                     "found on any other processors. This means that\n");
      (void) fprintf(stderr, "no processor is assigned to the row equal to \n");
      (void) fprintf(stderr, "this column number.\n");
      (void) fprintf(stderr, "Note: matrices must be square.\n");
      (void) fprintf(stderr, "The following columns were not found:\n");

      for (i = 0 ; i < num_sublist_elements; i++ ) {
        (void) fprintf(stderr, "%3d ", external[i * stride + first_extern]);
      }
      (void) fprintf(stderr,"\n");
        exit(-1);
    }

    /* send a zero length message to processors in the neigh_list */

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


    /* post receives */

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

    /* send messages */

    for (i = 0; i < num_send_neighbors; i++) {
      k = 1;
      mdwrap_write((void *) &k, 0, -extern_mess[i] - 1, type, &cflag);
    }

    /* receive the messages and create a new neighbor list */

    for (i = 0; i < num_recv_neighbors; i++) {
      partner  = -1;
      length   = mdwrap_wait((void *) &k, 0, &partner, &type, &cflag,request+i);
      neigh_list[i] = partner;
    }


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

    k