az_ilu_util.c

并行解法器,功能强大

/*====================================================================
 * ------------------------
 * | CVS File Information |
 * ------------------------
 *
 * $RCSfile: az_ilu_util.c,v $
 *
 * $Author: tuminaro $
 *
 * $Date: 1999/09/30 17:11:09 $
 *
 * $Revision: 1.15 $
 *
 * $Name:  $
 *====================================================================*/
#ifndef lint
static char rcsid[] = "$Id: az_ilu_util.c,v 1.15 1999/09/30 17:11:09 tuminaro Exp $";
#endif


/*******************************************************************************
 * Copyright 1995, Sandia Corporation.  The United States Government retains a *
 * nonexclusive license in this software as prescribed in AL 88-1 and AL 91-7. *
 * Export of this program may require a license from the United States         *
 * Government.                                                                 *
 ******************************************************************************/


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "az_aztec.h"

  /* externals */

extern void sort_blk_col_indx(int num_blks_row, int *bindx_start_row,
                       int *ordered_index);

extern void sort2(int, int *, int*);

extern void order_parallel(int M, double *val_old, double *val_new, 
                   int *bindx_old, int *bindx_new, int *indx_old, int *indx_new,
                   int *bpntr_old, int *bpntr_new, int *diag_block);
extern void get_diag(int M, int *bindx, int *bpntr, int *diag_block);

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

void sort_blk_col_indx(int num_blks_row, int *bindx_start_row,
                       int *ordered_index)

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

  Routine to sort the block entires for a given block row into increasing order.

  Author:
  =======

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

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

  num_blks_row:    Number of blocks in in this row.

  bindx_start_row: On Input:  starting address for the array of block column
                              indices for this row.
                   On Output: ordered block column indicies for this row.
                              (size >= num_blks_row).

  ordered_index:   On Input:  integer array.
                   On Output: an array of indices which describes the
                              reordering for the array "bindx_start_row"
                              (size >= num_blks_row).

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

{

  /* local variables */

  int i;


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

  /* Initialize the ordering index vector */

  for (i = 0; i < num_blks_row; i++) ordered_index[i] = i;

  /* Sort block column index array and produce the ordering index */

  sort2(num_blks_row, bindx_start_row-1, ordered_index-1);

} /* sort_blk_col_indx */

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

void sort2(int n, int ra[], int rb[])

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

  Numerical Recipes C source code modified to have first argument an integer
  array.

  Sorts the array ra[1,..,n] in ascending numerical order using heapsort
  algorithm, while making the corresponding rearrangement of the array
  rb[1,..,n].

  NOTE: The arrays start at 1 instead of 0, therefore you must pass call from C
  for a zero based array as:

                   sort(n, ra-1, rb-1);


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

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

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

  n:               Length of arrays ra and rb.

  ra:              Array to be sorted.

  rb               Second array order acording to the sorted ra array.

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

{

  /* local variables */

  int l, j, ir, i;
  int rra;
  int rrb;

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

  l  = (n >> 1) + 1;
  ir = n;

  for (;;) {
    if (l > 1) {
      rra = ra[--l];
      rrb = rb[l];
    }
    else {
      rra    = ra[ir];
      rrb    = rb[ir];
      ra[ir] = ra[1];
      rb[ir] = rb[1];

      if (--ir <= 1) {
        ra[1] = rra;
        rb[1] = rrb;
        return;
      }
    }

    i = l;
    j = l << 1;

    while (j <= ir) {
      if (j < ir && ra[j] < ra[j+1]) ++j;
      if (rra < ra[j]) {
        ra[i] = ra[j];
        rb[i] = rb[j];
        j    += (i = j);
      }
      else j = ir + 1;
    }

    ra[i] = rra;
    rb[i] = rrb;
  }

} /* sort2 */

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

void AZ_order(int M, double *val_old, double *val_new, int *bindx,
              int *indx_old, int *indx_new, int *bpntr, int *diag_block)

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

  For each row, reorders the blocks of the matrix (indices and values) in
  increasing order and constructs the array of pointers: diag_block to the
  diagonal blocks.

  Returns diag_block[i] = -1 if no diagonal block has been found at the i'th
  row.

  Author:          Lydie Prevost, SNL, 1422
  =======

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

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

  M:               The number of (block) rows in the matrix.

  val_old:         Array containing the entries of the matrix before reordering.
                   The matrix is stored block-row-by-block-row. Each block entry
                   is dense and stored by columns (VBR).

  val_new:         On output, array containing the entries of the matrix after
                   reordering. The matrix is stored block-row-by-block-row. Each
                   block entry is dense and stored by columns (VBR).

  bindx:           On input, contains the block column indices of the non-zero
                   block entries of the matrix before reordering.
                   On output, contains the block column indices of the non-zero
                   block entries of the matrix after reordering.

  indx_old:        The ith element of indx_old points to the location in val_old
                   of the (0,0) entry of the ith block entry. The last element
                   is the number of nonzero entries of matrix plus one.

  indx_new:        The ith element of indx_new points to the location in val_new
                   of the (0,0) entry of the ith block entry. The last element
                   is the number of nonzero entries of matrix plus one.

  bpntr:           The ith element of bpntr points to the first block entry of
                   the ith row in bindx. The last element is the number of
                   nonzero blocks of matrix plus one.

  diag_block:      On output, array of size M points on each diagonal block.

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

{

  /* local variables */

  int     i, kk, j, ii;
  int     num_blks_row;
  int    *sort, old_blk_index, counter;
  int     new_blk;
  int    *temp_ind, size_temp_ind = 10, size_temp_val = 40;
  double *temp_val;
  int     total_vals;
  int     start, end;

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

  temp_ind = (int    *) AZ_allocate(size_temp_ind*sizeof(int));
  temp_val = (double *) AZ_allocate(size_temp_val*sizeof(double));
  sort     = (int    *) AZ_allocate(sizeof(int) * (M));

  if ( (temp_val == NULL) || (sort == NULL)) 
     AZ_perror("Out of space inside AZ_sort()\n");

  for (i = 0; i < M; i++) diag_block[i] = -1;

  for (i = 0; i < M; i++) {     /* loop over the rows */

    /* constructs the new array of block column indices in this row */

    num_blks_row = bpntr[i+1] - bpntr[i];

    if (num_blks_row+1 > size_temp_ind) {
      size_temp_ind = num_blks_row + 1;
      AZ_free(temp_ind);
      temp_ind = (int *) AZ_allocate(size_temp_ind * sizeof(int));
    }

    for (ii = bpntr[i]; ii <= bpntr[i+1]; ii++)
      temp_ind[ii - bpntr[i]] = indx_old[ii];

    total_vals = indx_old[bpntr[i+1]] - indx_old[bpntr[i]];

    sort_blk_col_indx(num_blks_row, bindx+bpntr[i], sort);

    /* for each block of this row computes the new indices and constructs the
       pointers on the diagonal block i */

    indx_new[0] = indx_old[0];
    for (kk = 0; kk < num_blks_row; kk++) {
      new_blk = kk + bpntr[i];

      /* index into old VBR matrix and get the block size */

      indx_new[new_blk+1] = indx_new[new_blk] +
        (temp_ind[sort[kk]+1] - temp_ind[sort[kk]]);

      if (bindx[new_blk] == i) diag_block[i] = new_blk;
    }

    /* constructs the new array containing the entries of the matrix */

    if (total_vals > size_temp_val) {
      size_temp_val = total_vals;
      AZ_free(temp_val);
      temp_val = (double *) AZ_allocate(size_temp_val * sizeof(double));
    }

    start   = indx_old[bpntr[i]];
    end     = indx_old[bpntr[i+1]];
    counter = 0;

    for (ii = start; ii < end; ii++)
      temp_val[counter++] = val_old[ii];

    for (kk = 0; kk < num_blks_row; kk++) {

      /*
       * Get old block index into the coefficient array and new block location.
       */

      old_blk_index = temp_ind[sort[kk]] - temp_ind[0];
      new_blk       = kk + bpntr[i];

      counter = 0;
      for (j = indx_new[new_blk]; j < indx_new[new_blk+1]; j++) {
        val_new[j] = temp_val[old_blk_index + counter++];
      }
    }
  }

  AZ_free((void *) sort);
  AZ_free((void *) temp_ind);
  AZ_free((void *) temp_val);

} /* order */

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

void order_parallel(int M, double *val_old, double *val_new, int *bindx_old,
                    int *bindx_new, int *indx_old, int *indx_new,
                    int *bpntr_old, int *bpntr_new, int *diag_block)

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

  For each row, reorders the blocks of the matrix (indices and values) in
  increasing order and constructs the array of pointers: diag_block to the
  diagonal blocks.

  Returns diag_block[i]=-1  if no diagonal block has been found at the i'th row.

  Author:          Lydie Prevost, SNL, 1422
  =======

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

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

  M:               The number of (block) rows in the matrix.

  val_old:         Array containing the entries of the matrix before reordering.
                   The matrix is stored block-row-by-block-row. Each block entry
                   is dense and stored by columns (VBR).

  val_new:         On output, array containing the entries of the matrix after
                   reordering. The matrix is stored block-row-by-block-row. Each
                   block entry is dense and stored by columns (VBR).

  bindx_old:       Contains the block column indices of the non-zero block
                   entries of the matrix before reordering.

  bindx_new:       On output, contains the block column indices of the reordered
                   non-zero block entries of the matrix before reordering.

  indx_old:        The ith element of indx_old points to the location in val_old
                   of the (0,0) entry of the ith block entry. The last element
                   is the number of nonzero entries of matrix plus one.

  indx_new:        The ith element of indx_new points to the location in val_new
                   of the (0,0) entry of the ith block entry. The last element
                   is the number of nonzero entries of matrix plus one.