az_dd_overlap.c

并行解法器,功能强大

   
    off_set = 0;
    for (i = 0 ; i < num_neigh; i++) {
       jj = proc_neigh[i];
       actual_send_leng[i] = 0;
       start_send[i]       = off_set;
       for (k = 0 ; k < send_rowcnt[jj] ; k++ ) {
          actual_send_leng[i] += (bindx[send_rownum[jj][k]+1] - 
                                     bindx[send_rownum[jj][k]]);
          for (j= bindx[send_rownum[jj][k]];j < bindx[send_rownum[jj][k]+1];j++)
             iptr[off_set++] = bindx[j];
       }
    }
    AZ_splitup_big_msg(num_neigh,(char *)iptr,(char *)&bindx[bindx[enlarged_N]],
		       sizeof(int), start_send, actual_send_leng, 
		       actual_recv_leng, proc_neigh,type,&ncnt, proc_config);

    /* Obtain the off-diagonal values */

    type            =AZ_sys_msg_type;
    AZ_sys_msg_type =(AZ_sys_msg_type+1-AZ_MSG_TYPE) % AZ_NUM_MSGS +AZ_MSG_TYPE;
   
    off_set = 0;
    for (i = 0 ; i < num_neigh; i++) {
       jj = proc_neigh[i];
       for (k = 0 ; k < send_rowcnt[jj] ; k++ ) {
          for (j =bindx[send_rownum[jj][k]];j < bindx[send_rownum[jj][k]+1];j++)
             dptr[off_set++] = val[j];
       }
    }
    AZ_splitup_big_msg(num_neigh,(char *) dptr,(char *) &val[bindx[enlarged_N]],
		       sizeof(double), start_send, actual_send_leng, 
		       actual_recv_leng, proc_neigh,type,&ncnt, proc_config);



    enlarged_N += net_gain_N;

    /* Sort the new row number index array for further processing
     *   - allocate index array for new sorted row indices (enlarged matrix)
     *   - call AZ_sort to order the update index array 
     *   (enlarged_Rownum) ==> sorted_New_Rownum
     */ 

    sorted_New_Rownum = (int*) BV_ALLOC((enlarged_N+1) * sizeof(int));
    if (sorted_New_Rownum == NULL) {
      printf("Error in allocating memory for sorted_New_Rownum.\n");
      exit(-1);
    }
    for (i=0; i<enlarged_N; i++) sorted_New_Rownum[i] = enlarged_Rownum[i];
    AZ_sort(sorted_New_Rownum, enlarged_N, NULL, NULL);

    /* Prune the links that are not in the current enlarged node list
     *  - do this only at the last iteration
     *  - use PAZ_sorted_search to identify the outcasts
     */

    if (ndist == olap_size-1) {
      k   = bindx[0];
      ind = bindx[0];
      for (i=0; i<enlarged_N; i++) {
        m = bindx[i+1];
        for (j=k; j<m; j++) {
          ret_index = PAZ_sorted_search(bindx[j],enlarged_N,
                                        sorted_New_Rownum);
          if (ret_index >= 0) {
            bindx[ind] = bindx[j];
            val[ind++] = val[j];
          }   
        }
        k = m;
        bindx[i+1] = ind;
      }
    }

    /* clean up */

    for (i=0; i<nprocs; i++) { 
      if (send_rownum[i] != NULL) {
         BV_FREE(send_rownum[i]);  send_rownum[i] = NULL;
      }
    }
  }
  BV_FREE(comm_buffer);
  BV_FREE(sorted_New_Rownum);
  BV_FREE(send_rownum);

  nzeros = bindx[enlarged_N];

  /* extract a list of external nodes -> ext_nodelist, and
   * re-order the input matrix according to the local ordering 
   *  - at exit, bindx has been reordered with the local
   *    indices 0->N_rows-1 and the external indices 
   *    N_rows->N_ext_node+N_rows-1.  The global indices for the
   *    external list are recorded in ext_nodelist  
   */

  N_ext_node = enlarged_N - N_rows;
  sorted_exts = (int *) BV_ALLOC((N_ext_node+1)*sizeof(int));
  for (i = 0 ; i < N_ext_node ; i++ ) 
     sorted_exts[i] = enlarged_Rownum[i + N_rows];

  sprintf(str,"over_map %s",context->tag); 
  *map3 = (int *) AZ_manage_memory((N_ext_node+1)*sizeof(int),AZ_ALLOC,name,
				   str,&i);

  for (i = 0 ; i < N_ext_node; i++) (*map3)[i] = i + N_rows;
  AZ_sort(sorted_exts, N_ext_node, *map3, NULL);

  PAZ_find_local_indices(N_rows,bindx,Rownum,sorted_exts,
                        N_ext_node,*map3);

  /* use all current information to construct the data_org structure */

  PAZ_set_message_info(N_ext_node,N_rows,sorted_exts, Rownum,proc_config,
                      NULL, data_org, AZ_MSR_MATRIX,name, max_per_proc,context);

  /* put the enlarged matrix in the output parameters */

  (*New_N_rows) = enlarged_N;

  /* clean up */

  BV_FREE(sorted_exts);
  BV_FREE(enlarged_Rownum);
  BV_FREE(Rownum);

  (*data_org)[AZ_matrix_type] = AZ_MSR_MATRIX;
  (*data_org)[AZ_N_internal]  = N_rows;
  (*data_org)[AZ_N_border]    = 0;
  (*data_org)[AZ_N_external]  = N_ext_node;
  (*data_org)[AZ_N_int_blk]   = N_rows;
  (*data_org)[AZ_N_bord_blk]  = 0;
  (*data_org)[AZ_N_ext_blk]   = N_ext_node;
  
} /* AZ_dd_overlap */
 
/* *********************************************************************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* Given a local matrix with N_local number of rows (N_local, bindx) and
 * a list of the local nodes, construct the external node list 
 *
 * Input :
 *
 *   N_local       : the number of local nodes
 *   sorted_Rownum : a list of row indices (sorted) for the local rows
 *   bindx         : the connectivity of the matrix in MSR
 *
 * Output :
 *
 *   N_ext_node : the number of external nodes found
 *   ext_nodelist : a list of node numbers for the external nodes
 * ----------------------------------------------------------------------
 */

void PAZ_compose_external(int N_local, int *sorted_Rownum, int *bindx, 
                          int *N_ext_node, int **ext_nodelist)
{  
  int  i, ret_index, enode_leng, *enode_list;
  int k,kk;

  enode_leng = 0;
  for (i=N_local+1; i<bindx[N_local]; i++) {
    ret_index = PAZ_sorted_search(bindx[i],N_local,sorted_Rownum);
    if (ret_index < 0) enode_leng++;
  }

  /* allocate initially certain amount of memory for holding the list */

  enode_list = (int*) BV_ALLOC((1+enode_leng) * sizeof(int));
  if (enode_list == NULL) 
    AZ_perror("Error in allocating memory for enode_list.\n");

  /* fetch the column indices in bindx to find out if they are in the
   * local list (sorted_Rownum).  If not, record them as external nodes.
   * In order not to duplicate the external nodes on the external node
   * list, we call PAZ_insert_list to do the job.  We also keep track
   * that the list does not overflow, and if an overflow is detected,
   * re-allocate more memory to hold the additional nodes.
   */    

  enode_leng = 0;
  for (i=N_local+1; i<bindx[N_local]; i++) {
    ret_index = PAZ_sorted_search(bindx[i],N_local,sorted_Rownum);
    if (ret_index < 0) 
       enode_list[enode_leng++] = bindx[i];
  }
  AZ_sort(enode_list, enode_leng, NULL, NULL);

  /* remove duplicates */
 
  kk = 0;
  for (k = 1; k < enode_leng ; k++) {
    if (enode_list[kk] != enode_list[k]) {
      kk++;
      enode_list[kk] = enode_list[k];
    }
  }
 
  if (enode_leng != 0) kk++;
  enode_leng = kk;

  

  /* clean up */

  (*N_ext_node)   = enode_leng;
  (*ext_nodelist) = enode_list; 
}
/* *********************************************************************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* Given a sorted list of indices and the key, find the position of the
 * key in the list.  If not found, return where it would have been stored.
 *
 * author : Charles Tong (8117)
 *------------------------------------------------------------------------
 */

int PAZ_sorted_search(int key, int nlist, int *list)
{
  int  nfirst, nlast, nmid, found, index = 0;

  if (nlist <= 0) return -1;
  nfirst = 0;  
  nlast  = nlist-1; 
  if (key > list[nlast])  return -(nlast+1);
  if (key < list[nfirst]) return -(nfirst+1);
  found = 0;
  while ((found == 0) && ((nlast-nfirst)>1)) {
    nmid = (nfirst + nlast) / 2; 
    if (key == list[nmid])     {index  = nmid; found = 1;}
    else if (key > list[nmid])  nfirst = nmid;
    else                        nlast  = nmid;
  }
  if (found == 1)               return index;
  else if (key == list[nfirst]) return nfirst;
  else if (key == list[nlast])  return nlast;
  else                          return -(nfirst+1);
}

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

void PAZ_order_ele(int extern_index[], int N_update, int externs[], 
                   int N_external, int *m1, int *m3, int max_per_proc)

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

  Find an ordering for the elements in 'update' and 'external' (if option ==
  AZ_EXTERNS, we only order 'external'). We first order 'update' by numbering
  the internal elements first followed by the border elements. On output
  update_index[] contains the numbers 0->N_update-1 ordered such that
  update_index[k] < update_index[j] if update[k] is internal and update[j] is
  border. Next, we order the external elements such that elements that are
  updated by the same processor are contiguous.

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

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

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

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

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

  internal:        On output, equals the number of internal elements on this
                   processor.

  border:          On output, equals the number of border elements on this
                   processor.

  N_update:        Number of elements updated on this processor.

  bpntr, bindx:    Nonzero matrix information (see User's Guide).
                   Note: for MSR arrays the user invokes this routine by calling
                   AZ_order_ele(...,bindx,bindx, ....). That is both 'bpntr' and
                   'bindx' point to 'bindx'.

  externs:     

  N_external:      Number of external elements on this processor.

  option:          option = AZ_ALL ==> order internal, border and extern ele.
                   option = AZ_EXTERNS ==> order only external elements.
  m_type:          On input, m_type = AZ_MSR_MATRIX
                     or      m_type = AZ_VBR_MATRIX

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

{

  /* local variables */

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

  int  i, j, count;

  /*
   * Sift through the external elements. For each newly encountered external
   * point assign it the next index in the sequence. Then look for other
   * external elements who are update by the same node and assign them the next
   * set of index numbers in the sequence (ie. elements updated by the same node
   * have consecutive indices).
   */

  count = N_update;
  for (i = 0; i < N_external; i++) extern_index[i] = -1;

  for (i = 0; i < N_external; i++) {
    if (extern_index[i] == -1) {
      extern_index[i] = count++;
      m3[ extern_index[i] - N_update ] = m1[i];

      for (j = i + 1; j < N_external; j++) {
        if ((externs[j]/max_per_proc) == (externs[i]/max_per_proc)) {
           extern_index[j] = count++;
           m3[ extern_index[j] - N_update ] = m1[j];
        }
      }
    }
  }

} /* PAZ_order_ele */

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

void PAZ_find_local_indices(int N_update, int bindx[], int update[],
                           int *sorted_ext, int N_external, int map[])

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

  Given the global column indices for the matrix and a list of elements updated
  on this processor, compute the external elements and store them in the list
  'external' and change the global column indices to local column indices. In
  particular,

  On input, the column index bindx[k] is converted to j on output where

          update[j] = bindx[k]
   or
          external[j - N_update] = bindx[k]

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

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

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

  N_update:        Number of elements updated on this processor.

  bindx:           MSR or VBR column indices. On input, they refer to global
                   column indices. On output, they refer to local column indices
                   as described above. See User's Guide for more information.

  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.

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

  bnptr:           'bpntr[N_update]' indicates the location of the
                   last VBR nonzero.

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

{

  /* local variables */

  int  j, k;
  int *bins,shift;
  int  start,end;

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

  /* set up some bins so that we will be able to use AZ_quick_find() */

  bins = (int *) BV_ALLOC((N_update / 4 + 10)*sizeof(int));
  if  (bins == NULL) {
    (void) fprintf(stderr, "ERROR: Not enough temp space\n");
    exit(-1);
  }

  AZ_init_quick_find(update, N_update, &shift, bins);

  /*
   * Compute the location of the first and last column index that is stored in
   * the bindx[].
   */

  start = bindx[0]; end = bindx[bindx[0]-1]; 
  

  /*
   * Estimate the amount of space we will need by counting the number of
   * references to columns not found among 'update[]'.  At the same time replace
   * column indices found in update[] by the appropriate index into update[].
   * Add N_update to columns not found in 'update[]' (this effectively marks
   * them as external elements).
   *
   * Note the space estimate will be an over-estimate as we do not take into
   * account that there will be duplicates in the external list.
   */

  for (j = start; j < end; j++) {
    k = AZ_quick_find(bindx[j], update, N_update,shift,bins);

    if (k != -1) bindx[j] = k;
    else {
       k = AZ_find_index(bindx[j], sorted_ext,N_external);