az_scaling.c

并行解法器,功能强大

/*====================================================================
 * ------------------------
 * | CVS File Information |
 * ------------------------
 *
 * $RCSfile: az_scaling.c,v $
 *
 * $Author: tuminaro $
 *
 * $Date: 2000/06/02 16:48:32 $
 *
 * $Revision: 1.31 $
 *
 * $Name:  $
 *====================================================================*/
#ifndef lint
static char rcsid[] = "$Id: az_scaling.c,v 1.31 2000/06/02 16:48:32 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 <string.h>
#include "az_aztec.h"

/* static functions */

static void calc_blk_diag_Chol(double *val, int *indx, int *bindx, int *rpntr,
                               int *cpntr, int *bpntr, double *L, int *d_indx,
                               int *d_bindx, int *d_rpntr, int *d_bpntr,
                               int *data_org);
extern void sym_row_sum_scaling(int proc_config[]);

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

void AZ_scale_f(int action, AZ_MATRIX *Amat, int options[], double b[], double x[], 
	        int proc_config[], struct AZ_SCALING *scaling)

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

  Scale the matrix, rhs and initial guess as specified by the user.

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

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

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

  val:             Array containing the nonzero entries of the matrix (see
                    Aztec Usert).

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

  b:               Right hand side of linear system.

  x:               On input, contains the initial guess. On output contains the
                   solution to the linear system.

  options:         Determines specific solution method and other parameters.

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

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

  action:          Flag which determines whether to scale the matrix/rhs/sol or 
                   just the rhs or just the solution or to inverse scale the
                   rhs or solution.

  Amat:            Structure used to represent the matrix (see az_aztec.h
                   and Aztec User's Guide).


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

{

  /* local variables */

  char *yo = "AZ_scale: ";
  int  temp;

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

  temp = options[AZ_scaling];
  if ( (options[AZ_scaling] != AZ_none) && 
       (Amat->data_org[AZ_matrix_type] != AZ_MSR_MATRIX) &&
       (Amat->data_org[AZ_matrix_type] != AZ_VBR_MATRIX)  ) {
    options[AZ_scaling] = AZ_none;
  }

  if ( (options[AZ_ignore_scaling] == AZ_TRUE) && (options[AZ_pre_calc] <= AZ_recalc) 
      && (action == AZ_SCALE_MAT_RHS_SOL)) 
    scaling->A_norm = AZ_gmax_matrix_norm(Amat->val, Amat->indx, Amat->bindx, 
                                          Amat->rpntr, Amat->cpntr, Amat->bpntr, 
                                          proc_config, Amat->data_org);

  switch (options[AZ_scaling]) {

  case AZ_none:
    break;

  case AZ_Jacobi:
    AZ_block_diagonal_scaling(action, Amat, Amat->val, Amat->indx, Amat->bindx, 
                              Amat->rpntr, Amat->cpntr, Amat->bpntr, Amat->data_org, 
			      b, options, proc_config, scaling);
    break;

  case AZ_BJacobi:
    AZ_block_diagonal_scaling(action, Amat, Amat->val, Amat->indx, Amat->bindx, 
                              Amat->rpntr, Amat->cpntr, Amat->bpntr, Amat->data_org, 
			      b, options, proc_config, scaling);
    break;

  case AZ_row_sum:
    AZ_row_sum_scaling(action, Amat, b, options, scaling);
    break;

  case AZ_sym_diag:
    AZ_sym_diagonal_scaling(action,Amat,b,x,options,proc_config,scaling);
    break;

  case AZ_sym_row_sum:
    AZ_sym_row_sum_scaling_sl(action, Amat, 
                              b, x, options, proc_config, scaling);
    break;

  case AZ_sym_BJacobi:

    /* symmetric block Jacobi scaling */
#ifdef next_version
    AZ_sym_block_diagonal_scaling(val, indx, bindx, rpntr, cpntr, bpntr, b,
                                    options, data_org, proc_config);
    else if (action == AZ_RESCALE_SOL)
      AZ_sym_rescale_vbr(x, data_org);
#endif
    break;

  default:
    (void) fprintf(stderr, "%sERROR: invalid scaling option.\n"
                   "          options[AZ_scaling] = %d\n", yo,
                   options[AZ_scaling]);
  exit(-1);
  }
  options[AZ_scaling] = temp;

} /* AZ_scale */

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

void AZ_block_diagonal_scaling(int action, AZ_MATRIX *Amat, double val[], int indx[], 
			       int bindx[], int rpntr[], int cpntr[], int bpntr[], 
			       int data_org[], double b[], int options[], int proc_config[],
			       struct AZ_SCALING *scaling)

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

  Routine to block Jacobi scale the sparse matrix problem.  Note: this scales
  the matrix and the right-hand side and the resulting matrix is non-symmetric.

  If the matrix is in MSR format, it is treated as point entry (block size 1)
  and standard Jacobi scaling is performed.  Else, if the matrix is in VBR
  format, block Jacobi scaling is performed.

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

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

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

  val:             Array containing the nonzero entries of the matrix (see
                    Aztec Usert).

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

  b:               Right hand side of linear system.

  options:         Determines specific solution method and other parameters.

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

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

  action:          Flag which determines whether to scale the matrix or to
                   rescale the solution.

  Amat:            Structure used to represent the matrix (see az_aztec.h
                   and Aztec User's Guide).


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

{

  /* local variables */

  static AZ_MATRIX *Amat_inv;
  register int   iblk_row, i, j, k, ival, d_ival, jblk, ib;
  int            bpoff, idoff;
  int            d_bpoff, d_idoff;
  int            m1, n1, ib1, ib2;
  int            ione = 1, itemp;
  int            m, Proc, N;
  int            tsize;
  int            j_last, bindx_row;
  int            max_blk;
  static int    *d3_indx, *d3_bindx, *d3_rpntr, *d3_bpntr;
  static double *d3_inv, *sc_vec;
  double         *c, *work;
  char           None[2], Left[2], Lower[2], Unit[2], Upper[2];
  char          *yo = "AZ_block_diagonal_scaling: ";
  char          label[80];
int *ipiv,info, iminus_one = -1;
double one = 1.0;

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

  if ( (action == AZ_INVSCALE_SOL) || (action == AZ_SCALE_SOL)) return;

  /* initialize */

  N    = data_org[AZ_N_internal] + data_org[AZ_N_border];
  m    = data_org[AZ_N_int_blk]  + data_org[AZ_N_bord_blk];
  Proc = proc_config[AZ_node];

  scaling->action = AZ_left_scaling;

  if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {

    /***** MSR Matrix => point Jacobi scaling *****/

    sprintf(label,"sc_vec%d",options[AZ_recursion_level]);
    sc_vec = (double *) AZ_manage_memory((N+data_org[AZ_N_external])*sizeof(double), AZ_ALLOC,
                                         data_org[AZ_name], label,
                                         &itemp);

    if (((action == AZ_SCALE_RHS) || (action == AZ_INVSCALE_RHS) || 
         (options[AZ_pre_calc] >= AZ_reuse)) && (itemp == AZ_NEW_ADDRESS)) {
      (void) fprintf(stderr, "%sERROR: Previous scaling not found for matrix "
                     "with\ndata_org[AZ_name] = %d. Check\n"
                     "options[AZ_pre_calc]\n", yo,
                     data_org[AZ_name]);
      exit(-1);
    }

    if ( (options[AZ_pre_calc] <= AZ_recalc) && (action == AZ_SCALE_MAT_RHS_SOL)) {
      for (ib = 0; ib < N; ib++) {
        if (fabs(val[ib]) > DBL_MIN) sc_vec[ib] = 1.0 / val[ib];
        else                         sc_vec[ib] = 1.0;

        val[ib] = 1.0;
        j_last  = bindx[ib+1] - bindx[ib];
        bindx_row = bindx[ib];

        for (j = 0; j < j_last; j++) {
          k       = bindx_row + j;
          val[k] *= sc_vec[ib];
        }
      }
    }
    if ( (action == AZ_SCALE_MAT_RHS_SOL) || (action == AZ_SCALE_RHS) ) {
       for (ib = 0; ib < N; ib++) b[ib] *= sc_vec[ib];
    }
    if ( action == AZ_INVSCALE_RHS)  {
       for (ib = 0; ib < N; ib++) b[ib] /= sc_vec[ib];
    }
  }

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

    /***** VBR Matrix => block Jacobi scaling *****/

    sprintf(None, "N");

    /* First, compute the block-diagonal inverse (if not already computed) */

    tsize = 0;
    for (i = 0; i < m; i++)
      tsize += (rpntr[i+1] - rpntr[i]) * (cpntr[i+1] - cpntr[i]);

    sprintf(label,"d3_indx%d",options[AZ_recursion_level]);
    d3_indx  = (int *)    AZ_manage_memory((m+1)*sizeof(int), AZ_ALLOC,
                                           data_org[AZ_name], label, &itemp);
    sprintf(label,"d3_bindx%d",options[AZ_recursion_level]);
    d3_bindx = (int *)    AZ_manage_memory(m*sizeof(int), AZ_ALLOC,
                                           data_org[AZ_name], label, &itemp);
    sprintf(label,"d3_rpntr%d",options[AZ_recursion_level]);
    d3_rpntr = (int *)    AZ_manage_memory((m+1)*sizeof(int), AZ_ALLOC,
                                           data_org[AZ_name], label, &itemp);
    sprintf(label,"d3_bpntr%d",options[AZ_recursion_level]);
    d3_bpntr = (int *)    AZ_manage_memory((m+1)*sizeof(int), AZ_ALLOC,
                                           data_org[AZ_name], label, &itemp);
    sprintf(label,"d3_inv%d",options[AZ_recursion_level]);
    d3_inv   = (double *) AZ_manage_memory(tsize*sizeof(double), AZ_ALLOC,
                                           data_org[AZ_name], label, &itemp);
    sprintf(label,"Amat_inv%d",options[AZ_recursion_level]);
    Amat_inv = (AZ_MATRIX *) AZ_manage_memory(sizeof(AZ_MATRIX), AZ_ALLOC,
                                           data_org[AZ_name], label, &itemp);

    sprintf(label,"ipv %d",options[AZ_recursion_level]);
    ipiv = (int *) AZ_manage_memory((N+1)*sizeof(int), AZ_ALLOC,
                                           data_org[AZ_name], label, &itemp);

    Amat_inv->rpntr       = d3_rpntr;   Amat_inv->cpntr    = d3_rpntr;
    Amat_inv->bpntr       = d3_bpntr;   Amat_inv->bindx    = d3_bindx;
    Amat_inv->indx        = d3_indx;    Amat_inv->val      = d3_inv;
    Amat_inv->data_org    = data_org;
    Amat_inv->matvec      = Amat->matvec;
    Amat_inv->matrix_type = Amat->matrix_type;


    if (((action == AZ_SCALE_RHS) || (action == AZ_INVSCALE_RHS) || 
         (options[AZ_pre_calc] >= AZ_reuse)) && (itemp == AZ_NEW_ADDRESS)) {
      (void) fprintf(stderr, "%sERROR: Previous scaling not found for matrix "
                     "with\ndata_org[AZ_name] = %d. Check\n"
                     "options[AZ_pre_calc]\n", yo,
                     data_org[AZ_name]);
      exit(-1);
    }



    if ( (options[AZ_pre_calc] <= AZ_recalc) && (action == AZ_SCALE_MAT_RHS_SOL)) {

      AZ_calc_blk_diag_LU(val, indx, bindx, rpntr, cpntr, bpntr, d3_inv,
                           d3_indx, d3_bindx, d3_rpntr, d3_bpntr, data_org, ipiv);

      /* offset of the first block */

      bpoff = *bpntr;
      idoff = *indx;
      d_bpoff = *d3_bpntr;
      d_idoff = *d3_indx;

      /* scale the matrix 'A' */

      max_blk = 0;
      for (i = 0; i < m + data_org[AZ_N_ext_blk] ; i++) {
        if ( cpntr[i+1]-cpntr[i] > max_blk )
          max_blk = cpntr[i+1] - cpntr[i];
      }

      work = (double *) AZ_allocate(max_blk*max_blk*sizeof(double));
      if (work == NULL) {
        (void) fprintf(stderr, "%sERROR: not enough memory for diagonal\n"
                       "      scaling. Not able to allocate work\n"
                       "      array. Must run a smaller problem\n", yo);
        exit(-1);
      }

      /* loop over the block rows */

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

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

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

        /* starting index of current row block */

        ival = indx[bpntr[iblk_row] - bpoff] - idoff;

        /* starting index of current block row for diagonal scaling blocks */

        d_ival  = d3_indx[d3_bpntr[iblk_row] - d_bpoff] - d_idoff;

        /* loop over the (block) columns in the current (block) row */

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

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

          ib1 = cpntr[jblk];

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

          ib2 = cpntr[jblk+1];

          /* number of columns in the current block */

          n1    = ib2 - ib1;
          itemp = m1*n1;

          if (jblk == iblk_row) {

            /* diagonal block => set to identity */

            if (m1 != n1) {
              if (Proc == 0) {
                (void) fprintf(stderr, "%sERROR: diagonal blocks are not\n"
                               "       square inside scaling\n", yo);
              }
              exit(-1);
            }

            for (i = 0; i < m1; i++)
              for (k = 0; k < n1; k++)
                if (i == k)
                  val[ival + i + m1*k] = 1.0;
                else
                  val[ival + i + m1*k] = 0.0;
          }

          else {
            if (itemp > max_blk*max_blk) {
              (void) fprintf(stderr, "%sERROR: block size (%d) is too big =>\n",
                             yo, itemp);
              exit(-1);
            }

            /* Matrix solve */

            dgetrs_(None, &m1, &n1, &d3_inv[d_ival], &m1, 
                    &(ipiv[rpntr[iblk_row]]), &val[ival], &m1, &info,1);
          }
          ival += itemp;
        }
      }

      AZ_free((void *) work);



    }
    d_bpoff = *d3_bpntr;
    d_idoff = *d3_indx;

    /* lastly, scale the rhs */

    c = (double *) AZ_allocate((unsigned) N * sizeof(double));