az_matrix_util.c

并行解法器,功能强大

/*====================================================================
 * ------------------------
 * | CVS File Information |
 * ------------------------
 *
 * $RCSfile: az_matrix_util.c,v $
 *
 * $Author: tuminaro $
 *
 * $Date: 2000/06/02 16:46:55 $
 *
 * $Revision: 1.27 $
 *
 * $Name:  $
 *====================================================================*/
#ifndef lint
static char rcsid[] = "$Id: az_matrix_util.c,v 1.27 2000/06/02 16:46:55 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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>

#include "az_aztec.h"

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

double AZ_gmax_matrix_norm(double val[], int indx[],int bindx[], int rpntr[],
                           int cpntr[], int bpntr[], int proc_config[],
                           int data_org[])

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

  This routine returns maximum matrix norm for VBR matrix A: this is a parallel
  version for a distributed matrix.

  Author:          Scott A. Hutchinson, SNL, 1421
  =======

  Return code:     double, maximum matrix norm.
  ============

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

  val:             Array containing the nonzero entries of the matrix (see file
                   Aztec User's Guide).

  indx,
  bindx,
  rpntr,
  cpntr,
  bpntr:           Arrays used for DMSR and DVBR sparse matrix storage (see
                   file 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.

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

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

{

  /* local variables */

  register int indx_ptr, irow, iblk_row, jcol, jblk_col, icol_blk, iblk = 0;
  int           num_blk_rows, num_col_blks, num_blk_cols;
  double        row_sum = 0.0, row_max = 0.0;
  int           k;
  int           j_last, bindx_row;

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

  if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {

    for (irow = 0; irow < data_org[AZ_N_internal] + data_org[AZ_N_border];
         irow++) {

      /* compute diagonal contribution */

      row_sum = fabs(val[irow]);

      /* nonzero off diagonal contibution */

      j_last  = bindx[irow+1] - bindx[irow];
      bindx_row = bindx[irow];

      for (jcol = 0; jcol < j_last; jcol++) {
        k        = bindx_row + jcol;
        row_sum += fabs(val[k]);
      }
      row_max = max(row_sum, row_max);
    }

    row_max = AZ_gmax_double(row_max, proc_config);

    return row_max;
  }

  else if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) {

    /* loop over the block rows */

    for (iblk_row = 0;
         iblk_row < data_org[AZ_N_int_blk] + data_org[AZ_N_bord_blk];
         iblk_row++) {

      /* find out how many rows are in this block row */

      num_blk_rows = rpntr[iblk_row+1] - rpntr[iblk_row];

      /* find out how many block columns are in this block row */

      num_col_blks = bpntr[iblk_row+1] - bpntr[iblk_row];

      /* loop over all the rows in this block row */

      for (irow = 0; irow < num_blk_rows; irow++) {

        /* loop over all the column blocks in this block row */

        for (jblk_col = 0; jblk_col < num_col_blks; jblk_col++) {

          /* find out which column block this is */

          icol_blk = bindx[iblk];
          indx_ptr = indx[iblk++];

          /* find out how many columns are in this block */

          num_blk_cols = cpntr[icol_blk+1] - cpntr[icol_blk];

          /* loop over all the columns in this block */

          for (jcol = 0; jcol < num_blk_cols; jcol++)
            row_sum += fabs(val[indx_ptr + jcol*num_blk_rows + irow]);
        }

        iblk   -= num_col_blks;
        row_max = max(row_sum, row_max);
        row_sum = 0.0;
      }

      iblk += num_col_blks;
    }

    row_max = AZ_gmax_double(row_max, proc_config);

    return row_max;
  }

  else {
    (void) fprintf(stderr, "ERROR: invalid matrix type %d\n",
                   data_org[AZ_matrix_type]);
    exit(1);
  }
  return(0.0);

} /* AZ_gmax_matrix_norm */

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

double AZ_gvector_norm(int n, int p, double *x, int proc_config[])

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

  Function which returns the lp norm of the vector x, i.e., if p = 2, the
  standard Euclidean norm is returned. NOTE, to get the l-infinity norm,
  set p = -1.

  Author:          Scott A. Hutchinson, SNL, 1421
  =======

  Return code:     double, requested norm value.
  ============

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

  n:               Order of vector x.

  p:               Order of the norm to perform, i.e., ||x||p.

  x:               Vector of length n (on this processor).

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

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

{

  /* local variables */

  register int    i, j;
  register double sum, power;
  int             index, one = 1, *n_ptr;
  double          norm;

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

  /* error checking */

  if (p <= 0 && p != -1) return -1.0;

  /* initialize */

  n_ptr = &n;

  switch (p) {

  case -1:                      /* infinity norm */

    index = idamax_(n_ptr, x, &one) - 1;
    if (index < 0) return -1.0;
    norm  = AZ_gmax_double(fabs(x[index]), proc_config);
    break;

  case 1:                       /* sum norm */
    sum  = dasum_(n_ptr, x, &one);
    norm = AZ_gsum_double(sum, proc_config);
    break;

  case 2:                       /* Euclidean norm */
    norm = sqrt(AZ_gdot(n, x, x, proc_config));
    break;

  default:                      /* any other p-norm */

    sum = 0.0;
    for (i = 0; i < n; i++) {
      power = *x;
      for (j = 0; j < p; j++)
        power *= *x;
      sum += fabs(power);
      x++;
    }

    norm = pow(AZ_gsum_double(sum, proc_config), 1.0 / (double) p);
  }

  return norm;

} /* vector norm */

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

void AZ_print_vbr_matrix(int matrix_flag, int Proc, int itotal_nodes,
                      int ext_nodes, double val[], int indx[], int bindx[],
                      int rpntr[], int bpntr[])

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

  Prints out the VBR matrix and pointers.

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

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

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

  matrix_flag:     = 0 no matrix output, = 1 output matrix.

  Proc:            Current processor number.

  itotal_nodes:    Number of internal + border nodes on this processor.

  ext_nodes:       Number of external nodes.

  val:             Array containing the nonzero entries of the matrix (see file
                   Aztec User's Guide).

  indx,
  bindx,
  rpntr,
  cpntr,
  bpntr:           Arrays used for DMSR and DVBR sparse matrix storage (see
                   file Aztec User's Guide).

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

{

  /* local variables */

  int iblk, i, iblk_row, m1, n1, ipoint;
  int ival, jblk, j, ib1, ib2, jpoint;

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

  /* print out the VBR indexing information for the matrix */

  (void) printf("\n----- Proc: %d indx -----\n\n", Proc);

  for (iblk = 0; iblk < itotal_nodes; iblk++) {
    for (i = *(bpntr+iblk); i < *(bpntr+iblk+1); i++)
      (void) printf("%d ", *(indx+i));

    if (iblk == itotal_nodes - 1){
      (void) printf("%d\n", *(indx+i));
    }
    else
      (void) printf("\n");
  }

  (void) printf("\n----- Proc: %d bindx -----\n\n", Proc);

  for (iblk = 0; iblk < itotal_nodes; iblk++) {
    for (i = *(bpntr+iblk); i < *(bpntr+iblk+1); i++)
      (void) printf("%d ", *(bindx+i));
    (void) printf("\n");
  }

  (void) printf("\n----- Proc: %d rpntr -----\n\n", Proc);

  for (i = 0; i < itotal_nodes + ext_nodes + 1; i++)
    (void) printf("%d ", *(rpntr+i));
  (void) printf("\n");

  (void) printf("\n----- Proc: %d bpntr -----\n\n", Proc);

  for (i = 0; i < itotal_nodes + 1; i++)
    (void) printf("%d ", *(bpntr+i));
  (void) printf("\n");

  /* dump of matrix in a block output format */

  if (matrix_flag) {

    /* loop over block rows */

    for (iblk_row = 0; iblk_row < itotal_nodes; iblk_row++) {

      /* number of rows in the current row block */

      m1 = rpntr[iblk_row+1] - rpntr[iblk_row];

      /* starting index of current row block */

      ival = indx[bpntr[iblk_row]];

      /* loop over all the blocks in the current block-row */

      for (j = bpntr[iblk_row]; j < bpntr[iblk_row+1]; j++) {
        jblk = bindx[j];

        /* the starting point column index of the current block */

        ib1 = rpntr[jblk];

        /* ending point column index of the current block */

        ib2 = rpntr[jblk+1];

        /* number of columns in the current block */

        n1 = ib2 - ib1;

        (void) printf("\nProc: %d Block Row: %d Block Column: %d "
                      "Row Pointer: %d Column Pointer: %d\n", Proc, iblk_row,
                      jblk, rpntr[iblk_row], rpntr[jblk]);

        (void) printf("---------------------------------------"
                      "---------------------------------------\n");

        for (ipoint = 0; ipoint < m1; ipoint++) {
          for (jpoint = 0; jpoint < n1; jpoint++)
            (void) printf("val[%d]: %e ", ival+jpoint*m1+ipoint,
                          val[ival+jpoint*m1+ipoint]);
          (void) printf("\n");
        }

        ival += m1*n1;
      }
    }
  }

} /* print_vbr_matrix */

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

void AZ_dtrans(int *m, int *n, double *A)

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

  Perform an in-place transpose of a general m x n matrix stored in "A".

  Author:          Scott A. Hutchinson, SNL, 1421
  =======

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

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

  m, n:            Number of rows, columns in "A", respectively.

  A:               Matrix to be transposed.

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

{

  /* local variables */

  register int  i, j, w_index = 0, A_index;
  int           itemp;
  double       *work;

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

  work = (double *) AZ_allocate(*m * *n * sizeof(double));