az_util.c

并行解法器,功能强大

      (void) fprintf(stdout, "\t\t||r||_inf\n");
    (void) fprintf(stdout, "\t----------------------------- : %e\t%e\n",
                   scaled_r_norm, eps);
    (void) fprintf(stdout, "\t||A||_inf ||x||_1 + ||b||_inf\n");
    break;
    case AZ_expected_values:
    case AZ_weighted:
      (void) fprintf(stdout, "\t||r||_WRMS:\t\t%e\t%e\n", scaled_r_norm, eps);
    break;
    default:
      (void) fprintf(stderr, "terminate_status: ERROR: convergence test %d "
                     "not implemented\n", conv_flag);
    exit(-1);
    }

    (void) fprintf(stderr,"\n\t*********************************************"
                   "******************\n\n");
  }

  /* free memory */

  if (solver_name != NULL)
    AZ_free((void *) solver_name);

} /* AZ_terminate_status_print */

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

int AZ_breakdown_f(int N, double v[], double w[], double inner,
                   int proc_config[])

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

  Determine if the 2 vectors v and w are orthogonal. In particular,

       |<v,w>| <  100 ||v||_2 ||w||_2 DBL_EPSILON

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

  Return code:     int
  ============

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

  N:               Length of vectors v and w.

  v, w:            Two vectors of length N to be checked for othogonality.

  inner:           <v,w>

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

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

{

  double v_norm, w_norm;

  v_norm = AZ_gvector_norm(N, 2, v, proc_config);
  w_norm = AZ_gvector_norm(N, 2, w, proc_config);

  return (fabs(inner) < 100.0 * v_norm * w_norm * DBL_EPSILON);

} /* breakdown */

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

void AZ_random_vector(double u[], int data_org[], int proc_config[])

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

  Set the the vector u to a random vector.

  Author:          John N. Shadid, SNL, 1421
  =======

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

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

  u:               Vector to be initialized.

  data_org:        Array containing information on the distribution of the
                   matrix to this processor as well as communication parameters
                   (see Aztec User's Guide).

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

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

{

  /* local variables */

  static int seed = 493217272;
  int        i, N;

  /*********************** BEGIN EXECUTION *********************************/

  N    = (proc_config[AZ_node]+7) * (proc_config[AZ_node]+13) *
    (19 + proc_config[AZ_node]);

  seed = (int) (AZ_srandom1(&N)* (double) seed);
  N    = data_org[AZ_N_internal] + data_org[AZ_N_border];

  for (i = 0; i < N; i++) u[i] = AZ_srandom1(&seed);

} /* AZ_random_vector */

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

double AZ_srandom1(int *seed)

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

  Random number generator.

  Author:
  =======

  Return code:     double, random number.
  ============

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

  seed:            Random number seed.

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

{

  /* local variables */

  int    a = 16807;
  int    m = 2147483647;
  int    q = 127773;
  int    r = 2836;

  int    lo, hi, test;
  double rand_num;

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

  hi   = *seed / q;
  lo   = *seed % q;
  test = a * lo - r * hi;

  if (test > 0) *seed = test;
  *seed = test + m;

  rand_num = (double) *seed / (double) m;

  return rand_num;

} /* AZ_srandom1 */



int allo_count = 0, free_count = 0;

#ifdef MEM_CHECK
/* sophisticated wrappers for allocating memory */

struct widget {                  /* widget is used to maintain a linked   */
   int order;                    /* list of all memory that was allocated */
   int size;
   char *address;
   struct widget *next;
};
struct widget *widget_head =  NULL;  /* points to first element of allocated */
                                     /* memory linked list.                  */

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

void AZ_print_it() {
/*
 * Print out the allocated memory linked list
 */

   struct widget *current;

   current = widget_head;
   while (current != NULL) {
      printf("(%d,%d,%u)\n",current->order, current->size, current->address);
      current = current->next;
   }
}
extern void AZ_print_it();

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

char *AZ_allocate(unsigned int isize) {

/* 
 *  Allocate memory and record the event by placing an entry in the 
 *  widget_head list. Also recored the size of this entry as well.
 *
 *  Note: we actually allocate more memory than is requested (7 doubles more).
 *  This additional memory is used to record the 'size' and to mark the
 *  memory with a header and trailer which we can later check to see if
 *  they were overwritten.
 *
 */

    char *ptr, *header_start, *header_end;
    struct widget *widget;
    int *size_ptr, i, size;
    double *dptr;

    size = (int) isize;

    size = size + 7*sizeof(double);
    widget = (struct widget *) malloc(sizeof(struct widget));
    if (widget == NULL) return(NULL);
    ptr = (char *) malloc(size);
    if (ptr == NULL) {
       free(widget);
       return(NULL);
    }
    allo_count++;

    /* put trash in the space to make sure nobody is expecting zeros */
    for (i = 0 ; i < size/sizeof(char) ; i++ ) 
       ptr[i] = 'f';


    /* record the entry */

    widget->order = allo_count;
if (size == 7*sizeof(double) ) {
printf("allocating 0 space %u (%d)\n",ptr,size);
 i = 0;
 size = 1/i;
 widget = NULL;
}
    widget->size  = size - 7*sizeof(double);
    widget->next  = widget_head;
    widget_head   = widget;
    widget->address = ptr;

    size_ptr = (int *) ptr;
    size_ptr[0] = size - 7*sizeof(double);
    dptr     = (double *) ptr;

    /* mark the header */

    header_start = (char *) &(dptr[1]);

    for (i = 0 ; i < 3*sizeof(double)/sizeof(char) ; i++ )
       header_start[i] = 'x';

    /* mark the trailer */

    header_end = &(ptr[ (size/sizeof(char)) - 1]);
    header_start = (char *) &(dptr[4]);
    header_start = & (header_start[(size-7*sizeof(double))/sizeof(char)]);

    while (header_start <= header_end) {
       *header_start = 'x';
       header_start++;
    }

    return( (char *) &(dptr[4]) );
}

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

void AZ_free(void *vptr) {
/*
 * Free memory and remove the corresponding entry from the widget_head
 * list. Additionally, check that the size stored in the header is correct
 * and that there has not been any memory overwritten in the header or
 * the trailer.
 *
 */

   struct widget *current, *prev;
   double *dptr;
   int *iptr, size, i;
   char *header_start, *header_end, *ptr;

    ptr = (char *) vptr;
    free_count++;
    if (ptr == NULL) {
       printf("Trying to free a NULL ptr\n");
i = 0;
size = 1/i;
widget_head = NULL;
    }
    else {
       current = widget_head;
       prev    = NULL;
       dptr = (double *) ptr;
       --dptr;
       --dptr;
       --dptr;
       --dptr;
       ptr = (char *) dptr;
       while (current != NULL) {
          if (current->address == ptr) break;
          else { prev = current; current = current->next; }
       }
       if (current == NULL) {
          printf("the pointer %u was not found and thus can not be freed.\n",
                  ptr);
          exit(1);
       }
       else {
           /* check to see if the header is corrupted */
           iptr = (int *) ptr;
           header_start = (char *) &(dptr[1]);

           for (i = 0 ; i < 3*sizeof(double)/sizeof(char) ; i++ ) {
              if (header_start[i] != 'x') {
                 printf("header is corrupted for %u (%d,%d)\n",ptr,
                         current->size,current->order);
                 size =  0;
                 size = 1/size;
              }
           }
           size = iptr[0];

           /* check to see if the sizes are different */

           if (current->size != size) {
              printf("Freeing %u whose size has changed (%d,%d)\n",
                     current->address,current->size,size);
              exit(1);
           }

           /* check to see if the trailer is corrupted */

           header_end = &(ptr[ ((size+7*sizeof(double))/sizeof(char)) - 1]);
           header_start = (char *) &(dptr[4]);
           header_start = &(header_start[size/sizeof(char)]);

           while (header_start <= header_end) {
              if ( *header_start != 'x') {
                 printf("trailer is corrupted for %u (%d,%d)\n",
                         ptr, size,
                         current->order);
                 size =  0;
                 size = 1/size;
              }
              header_start++;
           }

           /* free the space and the widget */

           free(ptr);
           if (widget_head == current) widget_head = current->next;
           else prev->next = current->next;
           free(current);

       }
   }

}

char *AZ_realloc(void *vptr, unsigned int new_size) {

   struct widget *current, *prev;
   int i, *iptr, size, *new_size_ptr;
   char *header_start, *header_end, *ptr;
   char *data1, *data2, *new_ptr, *new_header_start, *new_header_end;
   int newmsize, smaller;
   double *dptr, *new_dptr;

    ptr = (char *) vptr;
    if (ptr == NULL) {
       printf("Trying to realloc a NULL ptr\n");
       exit(1);
    }
    else {
       current = widget_head;
       prev    = NULL;
data1 = ptr;
       dptr = (double *) ptr;
       --dptr;
       --dptr;
       --dptr;
       --dptr;
       ptr = (char *) dptr;
       while (current !=