az_precond.c

并行解法器,功能强大

/*====================================================================
 * ------------------------
 * | CVS File Information |
 * ------------------------
 *
 * $RCSfile: az_precond.c,v $
 *
 * $Author: tuminaro $
 *
 * $Date: 1999/11/15 18:37:29 $
 *
 * $Revision: 1.34 $
 *
 * $Name:  $
 *====================================================================*/
#ifndef lint
static char rcsid[] = "$Id: az_precond.c,v 1.34 1999/11/15 18:37:29 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 <float.h>
#include "az_aztec.h"

/*---------------- External Definitions -------------------------------------*/

#ifdef eigen
extern void AZ_do_Jacobi(double val[], int indx[], int bindx[], int rpntr[],
                     int cpntr[], int bpntr[], double x[], double b[],
                     double temp[], int options[], int data_org[],
                     int proc_config[], double params[], int flag);
#endif
extern void AZ_calc_blk_diag_inv(double *val, int *indx, int *bindx, int *rpntr,
                     int *cpntr, int *bpntr, double *d_inv, int *d_indx,
                     int *d_bindx, int *d_rpntr, int *d_bpntr,
                     int data_org[]);
extern void jacobi(double val[], double x[], int data_org[]);

extern int AZ_sys_msg_type;

/*---------------------------------------------------------------------------*/


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

void AZ_precondition(double x[], int input_options[], int proc_config[],
                     double input_params[], AZ_MATRIX *Amat, 
		     AZ_PRECOND *input_precond)


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

  This routine calls appropriate sparse matrix preconditioner.

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

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

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


  x:               On input, contains the current solution. On output contains
                   the preconditioned solution to the linear system.

  options:         Determines specific solution method and other parameters.

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

  params:          Drop tolerance and convergence tolerance info.

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

  precond:         Structure used to represent the preconditionner
                   (see file az_aztec.h and Aztec User's Guide).

 * --------------------------------------------------------------------

 Related routines:

   scaling routines:
        AZ_block_diagonal_scaling -- block-diagonally scales sparse matrix
                                     problem.
        AZ_row_sum_scaling        -- row sum scales sparse matrix problem.
        sym_diagonal_scaling      -- diagonaly scales symm. sparse problem.
        sym_row_sum_scaling       -- row sum scales symmetric sparse problem.

   preconditioners:
        jacobi                 -- point Jacobi method.
        AZ_polynomial_expansion-- Polynomial expansion; Neumann series and
                                  least squares.
        domain decomposition   -- Block solvers (LU , ILU or ILUT) used on 
                                  each processor. The blocks are either
                                  non-overlapping or overlapping.
        icc                    -- incomplete sparse Choleski (symmetric
                                  version).

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

{

  /* local variables */

  int            ione = 1;
  double        *temp;
  int            m, N, k, length;
  int            i, step, j;
  static int    *d2_indx,*d2_bindx,*d2_rpntr,*d2_bpntr;
  static double *d2_inv;
  static AZ_MATRIX *Dmat;
  int            tsize, multilevel_flag = 0, max_externals;
  static int     previous_factors = -1;
  double        *v, *y;
  char          *yo = "precond: ";
  int          *data_org, *bindx, *indx, *cpntr, *rpntr, *bpntr;
  double       *val;
  char         label[64],suffix[32];
  char         tag[80];
  double       *current_rhs, *orig_rhs = NULL, *x_precond = NULL;
  int          *options, *ioptions, N_fixed, *fixed_pts;
  double       *params,  *iparams, *istatus;
  AZ_MATRIX    *Aptr, *Pmat;
  AZ_PRECOND   *Pptr, *precond;
  struct AZ_SCALING *Sptr;
  int          opt_save1, opt_save2, opt_save3, opt_save4, opt_save5, *itemp;
  double       *tttemp, norm1, *dtemp;
#ifdef TIMING
  double       ttt;
#endif


#ifdef eigen
  double         *tb, *tr;
#endif

  /**************************** execution begins ******************************/
#ifdef TIMING
  ttt = AZ_second();
#endif

  precond = input_precond;

  sprintf(suffix," in precond%d",input_options[AZ_recursion_level]);  
                                              /* set string that will be used */
                                              /* for manage_memory label      */

  data_org = precond->Pmat->data_org;
  options  = input_options;
  params   = input_params;

  m    = data_org[AZ_N_int_blk] + data_org[AZ_N_bord_blk];
  N    = data_org[AZ_N_internal] + data_org[AZ_N_border];
  max_externals = Amat->data_org[AZ_N_external];
  if (max_externals < data_org[AZ_N_external]) 
     max_externals = data_org[AZ_N_external];

  current_rhs = x; 
  if (options[AZ_precond] == AZ_multilevel) {

     /* make extra vectors to hold rhs and residual */

     sprintf(tag,"orig_rhs %s",precond->context->tag);
     orig_rhs = AZ_manage_memory((N+max_externals)*sizeof(double),
                               AZ_ALLOC, AZ_SYS,tag,&i);
     sprintf(tag,"x_prec %s",precond->context->tag);
     x_precond    = AZ_manage_memory((N+max_externals)*sizeof(double),
                               AZ_ALLOC, AZ_SYS, tag,&i);
     for (i = 0 ; i < N; i++) x_precond[i] = 0.0;
     for (i = 0 ; i < N; i++) orig_rhs[i] = current_rhs[i];
     multilevel_flag = 1;
     options = precond->options;
     params  = precond->params;
  }

  do {
     data_org = precond->Pmat->data_org;
     val      = precond->Pmat->val;
     bindx    = precond->Pmat->bindx;
     cpntr    = precond->Pmat->cpntr;
     indx     = precond->Pmat->indx;
     rpntr    = precond->Pmat->rpntr;
     bpntr    = precond->Pmat->bpntr;
     if (max_externals < data_org[AZ_N_external]) 
        max_externals = data_org[AZ_N_external];

     switch (options[AZ_precond]) {
     case AZ_none:
     break;

     case AZ_Jacobi:
        if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {
           for (i = 0; i < N; i++) current_rhs[i] /= val[i];

           if (options[AZ_poly_ord] > 1) {
              sprintf(tag,"v_prec %s",precond->context->tag);
              v = AZ_manage_memory((N+max_externals)*sizeof(double),
                                    AZ_ALLOC, AZ_SYS, tag, &i);
              sprintf(tag,"y_prec %s",precond->context->tag);
              y = AZ_manage_memory(N*sizeof(double), AZ_ALLOC, AZ_SYS, tag,&i);
              for (i = 0; i < N; i++) v[i] = current_rhs[i];

              for (step = 1; step < options[AZ_poly_ord]; step++) {
                 Amat->matvec(v, y, Amat, proc_config);
                 for(i = 0; i < N; i++) v[i] += current_rhs[i] - y[i] / val[i];
              }
              for (i = 0; i < N; i++) current_rhs[i] = v[i];
           }
        }
        else if (data_org[AZ_matrix_type] == AZ_USER_MATRIX) {
           if (options[AZ_pre_calc] < AZ_sys_reuse) {
              sprintf(tag,"d2_inv %s",precond->context->tag);
              d2_inv   = (double *) AZ_manage_memory(N*sizeof(double),AZ_ALLOC,
						data_org[AZ_name],tag,&i);
              Pmat = precond->Pmat;
              if ( (Pmat->N_nz < 0) || (Pmat->max_per_row < 0)) 
                 AZ_matfree_Nnzs(Pmat);

              if ( (Pmat->getrow == NULL) && (N != 0) ) {
                 printf("Error: Only matrices with getrow() defined via ");
                 printf("AZ_set_MATFREE_getrow(...) can do Jacobi preconditioning\n");
                 exit(1);
              }
              sprintf(tag,"dtemp %s",precond->context->tag);
              dtemp = (double *) AZ_manage_memory(Pmat->max_per_row*
				                sizeof(double),AZ_ALLOC,
						data_org[AZ_name],tag,&i);
              sprintf(tag,"itemp %s",precond->context->tag);
              itemp = (int    *) AZ_manage_memory(Pmat->max_per_row*
				                sizeof(int   ),AZ_ALLOC,
						data_org[AZ_name],tag,&i);
  
	      for (i = 0; i < N; i++) {
                 Pmat->getrow(itemp,dtemp,&length,Pmat,1,&i,Pmat->max_per_row);
                 for (k =0; k < length; k++) 
                    if (itemp[k] == i) break;

                 if (k == length) d2_inv[i] = 0.0; /* no diagonal */
                 else d2_inv[i] = 1./dtemp[k];
              }
           }
           for (i = 0; i < N; i++) current_rhs[i] *= d2_inv[i];

           if (options[AZ_poly_ord] > 1) {
              sprintf(tag,"v_prec %s",precond->context->tag);
              v = AZ_manage_memory((N+max_externals)*sizeof(double),
                                    AZ_ALLOC, AZ_SYS, tag, &i);
              sprintf(tag,"y_prec %s",precond->context->tag);
              y = AZ_manage_memory(N*sizeof(double), AZ_ALLOC, AZ_SYS, tag,&i);
              for (i = 0; i < N; i++) v[i] = current_rhs[i];

              for (step = 1; step < options[AZ_poly_ord]; step++) {
                 Amat->matvec(v, y, Amat, proc_config);
                 for(i = 0; i < N; i++) v[i] += current_rhs[i] - y[i]*d2_inv[i];
              }
              for (i = 0; i < N; i++) current_rhs[i] = v[i];
           }
        }
        else if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) {
           /* block Jacobi preconditioning */

           if (options[AZ_pre_calc] < AZ_sys_reuse) {
              /* First, compute block-diagonal inverse */
              /* (only if not already computed)        */

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

                 sprintf(tag,"d2_indx %s",precond->context->tag);
                 d2_indx  = (int *) AZ_manage_memory((m+1)*sizeof(int),AZ_ALLOC,
                                            data_org[AZ_name], tag, &i);
                 sprintf(tag,"d2_bindx %s",precond->context->tag);
                 d2_bindx = (int *) AZ_manage_memory(m*sizeof(int), AZ_ALLOC,
                                            data_org[AZ_name], tag, &i);
                 sprintf(tag,"d2_rpntr %s",precond->context->tag);
                 d2_rpntr = (int *) AZ_manage_memory((m+1)*sizeof(int),AZ_ALLOC,
                                            data_org[AZ_name], tag, &i);
                 sprintf(tag,"d2_bpntr %s",precond->context->tag);
                 d2_bpntr = (int *) AZ_manage_memory((m+1)*sizeof(int),AZ_ALLOC,
                                            data_org[AZ_name], tag, &i);
                 sprintf(tag,"d2_inv %s",precond->context->tag);
                 d2_inv   = (double *) AZ_manage_memory(tsize*sizeof(double),
                                            AZ_ALLOC, data_org[AZ_name],tag,&i);
                 d2_bpntr[0] = 0;
                 sprintf(tag,"dmat_calk_binv %s",precond->context->tag);
                 Dmat     = (AZ_MATRIX *) AZ_manage_memory(sizeof(AZ_MATRIX), 
                                            AZ_ALLOC,data_org[AZ_name],tag,&i);

                 Dmat->rpntr         = d2_rpntr;   Dmat->cpntr    = d2_rpntr;
                 Dmat->bpntr         = d2_bpntr;   Dmat->bindx    = d2_bindx;
                 Dmat->indx          = d2_indx;    Dmat->val      = d2_inv;
                 Dmat->data_org      = data_org;
                 Dmat->matvec        = precond->Pmat->matvec;
                 Dmat->matrix_type   = precond->Pmat->matrix_type;

                 if (options[AZ_pre_calc] != AZ_reuse) {
                    AZ_calc_blk_diag_inv(val, indx, bindx, rpntr, cpntr, bpntr,
                                         d2_inv, d2_indx, d2_bindx, d2_rpntr, 
                                         d2_bpntr, data_org);
                 }
                 else if (i == AZ_NEW_ADDRESS) {
                   fprintf(stderr, "Error: options[AZ_pre_calc]==AZ_reuse and"
                         "previous factors\n       not found. Check"
                         "data_org[AZ_name].\n");
                   exit(-1);
                 }
           }
           else if (previous_factors != data_org[AZ_name]) {
              fprintf(stderr, "Warning: Using a previous factorization as a"
                       "preconditioner\neven though matrix"
                       "(data_org[AZ_name]) has changed\n");
           }
           previous_factors = data_org[AZ_name];

           /* scale rhs */

           sprintf(tag,"v_prec %s",precond->context->tag);
           v = AZ_manage_memory((N+max_externals)*sizeof(double),
                           AZ_ALLOC, AZ_SYS, tag, &i);

           Dmat->matvec(current_rhs, v, Dmat, proc_config);

           dcopy_(&N, v, &ione, current_rhs, &ione);

           if (options[AZ_poly_ord] > 1) {
              sprintf(tag,"y_prec %s",precond->context->tag);
              y = AZ_manage_memory((N+max_externals)*sizeof(double),
                             AZ_ALLOC, AZ_SYS, tag, &i);

              sprintf(tag,"temp_prec %s",precond->context->tag);
              temp = AZ_manage_memory(N*sizeof(double), AZ_ALLOC,AZ_SYS,tag,&i);

              for (step = 1; step < options[AZ_poly_ord]; step++) {
                 Amat->matvec(v, y, Amat, proc_config);
                 Dmat->matvec(y, temp, Dmat, proc_config);

                 for (i = 0; i < N; i++) v[i] += current_rhs[i] - temp[i];
              }

              for (i = 0; i < N; i++) current_rhs[i] = v[i];
           }
        }
     break;
     case AZ_sym_GS: