az_domain_decomp.c

来自「并行解法器,功能强大」· C语言 代码 · 共 726 行 · 第 1/2 页

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

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


void AZ_domain_decomp(double x[], AZ_MATRIX *Amat, int options[],
                   int proc_config[], double params[],
		   struct context *context)


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

  Precondition 'x' using an overlapping domain decomposition method where a 
  solver specified by options[AZ_subdomain_solve] is used on the subdomains. 
  Note: if a factorization needs to be computed on the first iteration, this
  will be done and stored for future iterations.

  Author:          Lydie Prevost, SNL, 9222
  =======          Revised by R. Tuminaro (4/97), SNL, 9222

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

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

  N_unpadded:      On input, number of rows in linear system (unpadded matrix) 
                   that will be factored (after adding values for overlapping).

  Nb_unpadded:     On input, number of block rows in linear system (unpadded) 
                   that will be factored (after adding values for overlapping).

  N_nz_unpadded:   On input, number of nonzeros in linear system (unpadded)
                   that will be factored (after adding values for overlapping).
             
  x:               On output, x[] is preconditioned by performing the subdomain
                   solve indicated by options[AZ_subdomain_solve].

  val    indx       
  bindx  rpntr:    On input, arrays containing matrix nonzeros (see manual). 
  cpntr  bpntr            

  options:         Determines specific solution method and other parameters.  In
                   this routine, we are concerned with options[AZ_overlap]:

                      == AZ_none: nonoverlapping domain decomposition
                      == AZ_diag: use rows corresponding to external variables 
                                  but only keep the diagonal for these rows.
                      == k      : Obtain rows that are a distance k away from
                                  rows owned by this processor.
                                  
  data_org:        Contains information on matrix data distribution and 
                   communication parameters (see manual).

*******************************************************************************/
{
  int N_unpadded, Nb_unpadded, N_nz_unpadded;
  double *x_pad = NULL, *x_reord = NULL, *ext_vals = NULL;
  int N_nz, N_nz_padded, nz_used;
  int mem_orig, mem_overlapped, mem_factor;
  int name, i, bandwidth;
  int *ordering = NULL;
  double start_t;
  int estimated_requirements;
  char str[80];
int *garbage;

  int N;
  int *padded_data_org = NULL, *bindx, *data_org;
  double *val;
  int *inv_ordering = NULL;
  int *map = NULL;
  AZ_MATRIX *A_overlapped = NULL;
  struct aztec_choices aztec_choices;


  /**************************** execution begins ******************************/
  data_org = Amat->data_org;
  bindx    = Amat->bindx;
  val      = Amat->val;
  N_unpadded = data_org[AZ_N_internal] + data_org[AZ_N_border];
  Nb_unpadded = data_org[AZ_N_int_blk]+data_org[AZ_N_bord_blk];
  if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) 
     N_nz_unpadded = bindx[N_unpadded];
  else if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX)
     N_nz_unpadded = (Amat->indx)[(Amat->bpntr)[Nb_unpadded]];
  else {
     if (Amat->N_nz < 0) 
        AZ_matfree_Nnzs(Amat);
     N_nz_unpadded = Amat->N_nz;
  }

  
  aztec_choices.options  = options;
  aztec_choices.params   = params;
  context->aztec_choices = &aztec_choices;
  context->proc_config   = proc_config;
  name                   = data_org[AZ_name];


  if ((options[AZ_pre_calc] >= AZ_reuse) && (context->Pmat_computed)) {
     N               = context->N;
     N_nz            = context->N_nz;
     A_overlapped    = context->A_overlapped;
     A_overlapped->data_org  = data_org;
     A_overlapped->matvec = Amat->matvec;
  }
  else {
     sprintf(str,"A_over %s",context->tag);
     A_overlapped = (AZ_MATRIX *) AZ_manage_memory(sizeof(AZ_MATRIX), 
                                                   AZ_ALLOC, name, str, &i);
     AZ_matrix_init(A_overlapped, 0);

     context->A_overlapped     = A_overlapped;
     A_overlapped->data_org    = data_org;
     A_overlapped->matvec      = Amat->matvec;
     A_overlapped->matrix_type = AZ_MSR_MATRIX;

     AZ_init_subdomain_solver(context);

     AZ_compute_matrix_size(Amat, options, N_nz_unpadded, N_unpadded, 
			 &N_nz_padded, data_org[AZ_N_external],
		 	 &(context->max_row), &N, &N_nz, params[AZ_ilut_fill], 
                         &(context->extra_fact_nz_per_row),
                         Nb_unpadded,&bandwidth);

     
        estimated_requirements = N_nz;
        if (N_nz*2 > N_nz) N_nz = N_nz*2;	/* check for overflow */
						/* Add extra memory to N_nz. */
                                                /* This extra memory is used */
                                                /* as temporary space during */
                                                /* overlapping calculation   */

        /* Readjust N_nz by allocating auxilliary arrays and allocate */
        /* the MSR matrix to check that there is enough space.        */

        /* block off some space for map and padded_data_org in dd_overlap */

        garbage = (int *) AZ_allocate((AZ_send_list + 2*(N-N_unpadded) +10)*
                                      sizeof(int));
        AZ_hold_space(context, N);

   
        if (N_nz*((int) sizeof(double)) < N_nz) N_nz=N_nz/2; /* check for overflow */
        if (N_nz*((int) sizeof(double)) < N_nz) N_nz=N_nz/2; /* check for overflow */
        if (N_nz*((int) sizeof(double)) < N_nz) N_nz=N_nz/2; /* check for overflow */
        if (N_nz*((int) sizeof(double)) < N_nz) N_nz=N_nz/2; /* check for overflow */
        if (N_nz*((int) sizeof(double)) < N_nz) N_nz=N_nz/2; /* check for overflow */

        N_nz = AZ_adjust_N_nz_to_fit_memory(N_nz,
                                 context->N_large_int_arrays,
                                 context->N_large_dbl_arrays);
        context->N_nz_factors = N_nz;

        if (N_nz <= N_nz_unpadded) {
           printf("Error: Not enough space for domain decomposition\n");
           AZ_exit(1);
        }


        if (estimated_requirements > N_nz ) estimated_requirements = N_nz;

        /* allocate matrix via AZ_manage_memory() */

        sprintf(str,"bindx %s",context->tag);
        A_overlapped->bindx =(int    *) AZ_manage_memory(N_nz*sizeof(int),
                                                AZ_ALLOC, name, str, &i);
        sprintf(str,"val %s",context->tag);
        A_overlapped->val =(double *) AZ_manage_memory(N_nz*sizeof(double),
                                                AZ_ALLOC, name, str, &i);
        context->N_nz_allocated = N_nz;
        AZ_free_space_holder(context);
        AZ_free(garbage);

        /* convert to MSR if necessary */

        if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX)
          AZ_vb2msr(Nb_unpadded,val,Amat->indx,bindx,Amat->rpntr,Amat->cpntr,
		    Amat->bpntr, A_overlapped->val, A_overlapped->bindx);
        else if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {
          for (i = 0 ; i < N_nz_unpadded; i++ ) {
             A_overlapped->bindx[i] = bindx[i];
             A_overlapped->val[i]   = val[i];
          }
        }
        else AZ_matfree_2_msr(Amat,A_overlapped->val,A_overlapped->bindx,N_nz);

        mem_orig = AZ_gsum_int(A_overlapped->bindx[N_unpadded],proc_config);

        start_t = AZ_second();
        AZ_pad_matrix(context, proc_config, N_unpadded, &N, 
                      &(context->map), &(context->padded_data_org), &N_nz, 
                      estimated_requirements);

/*
        if (proc_config[AZ_node] == 0)
           printf("matrix padding took %e seconds\n",AZ_second()-start_t);
*/


        mem_overlapped = AZ_gsum_int(A_overlapped->bindx[N],proc_config);
  
        if (options[AZ_reorder]) {
           start_t = AZ_second();
           AZ_find_MSR_ordering(A_overlapped->bindx,&ordering,N,
                                &(context->inv_ordering),name,context);
/*
           if (proc_config[AZ_node] == 0) 
              printf("took %e seconds to find ordering\n", AZ_second()-start_t);
*/
           start_t = AZ_second();
           AZ_mat_reorder(N,A_overlapped->bindx,A_overlapped->val,&ordering,
                          context->inv_ordering);
/*
           if (proc_config[AZ_node] == 0) 
              printf("took %e seconds to reorder\n", AZ_second()-start_t);
*/
                /* NOTE: ordering is freed inside AZ_mat_reorder */
        }

        /* Do a factorization if needed.  */

        start_t = AZ_second();
        AZ_factor_subdomain(context, N, N_nz, &nz_used);

/*
        start_t        = AZ_second()-start_t;
        max_time = AZ_gmax_double(start_t,proc_config);
        min_time = AZ_gmin_double(start_t,proc_config);
        if (proc_config[AZ_node] == 0) 
           printf("time for subdomain solvers ranges from %e to %e\n",
                  min_time,max_time);
*/
  
        if ( A_overlapped->matrix_type == AZ_MSR_MATRIX)
           AZ_compress_msr(&(A_overlapped->bindx), &(A_overlapped->val),
                     context->N_nz_allocated, nz_used, name, context);


        context->N_nz = nz_used;
        context->N    = N;
        context->N_nz_allocated = nz_used;

        mem_factor = AZ_gsum_int(nz_used,proc_config);

        if (proc_config[AZ_node] == 0)
           AZ_print_header(options,mem_overlapped,mem_orig,mem_factor);

        if (options[AZ_overlap] >= 1) {
           sprintf(str,"x_pad %s",context->tag);
           context->x_pad  = (double *) AZ_manage_memory(N*sizeof(double),
                                                   AZ_ALLOC, name, str, &i);
           sprintf(str,"ext_vals %s",context->tag);
           context->ext_vals = (double *) AZ_manage_memory((N-N_unpadded)*
                                             sizeof(double), AZ_ALLOC, name, 
                                             str, &i);
        }
        if (options[AZ_reorder]) {
           sprintf(str,"x_reord %s",context->tag);
           context->x_reord = (double *) AZ_manage_memory(N*sizeof(double),
                                             AZ_ALLOC, name, str, &i);
        }

     }

  /* Solve L u = x where the solution u overwrites x */

    x_reord         = context->x_reord;
    inv_ordering    = context->inv_ordering;
    x_pad           = context->x_pad;
    ext_vals        = context->ext_vals;
    padded_data_org = context->padded_data_org;
    map             = context->map;

   if (x_pad == NULL) x_pad = x;

   if (options[AZ_overlap] >= 1) {

      for (i = 0 ; i < N_unpadded ; i++) x_pad[i] = x[i];
      AZ_exchange_bdry(x_pad,padded_data_org, proc_config);
      for (i = 0 ; i < N-N_unpadded ; i++ ) 
         ext_vals[map[i]-N_unpadded] = x_pad[i+N_unpadded];
      for (i = 0 ; i < N-N_unpadded ; i++ ) x_pad[i + N_unpadded] = ext_vals[i];
   }
   else if (options[AZ_overlap] == AZ_diag) 
	AZ_exchange_bdry(x_pad,data_org, proc_config);

   if (x_reord == NULL) x_reord = x_pad;
   if (options[AZ_reorder]) {
      for (i = 0 ; i < N ; i++ ) x_reord[inv_ordering[i]] = x_pad[i];
   }

   AZ_solve_subdomain(x_reord,N, context);

   if (options[AZ_reorder])
      for (i = 0; i < N; i++) x_pad[i] = x_reord[inv_ordering[i]];

   AZ_combine_overlapped_values(options[AZ_type_overlap],padded_data_org, 
                             options, x_pad, map,ext_vals,name,proc_config);

   if (x_pad != x) 
     for (i = 0 ; i < N_unpadded ; i++ ) x[i] = x_pad[i];

} /* subdomain driver*/

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

void AZ_print_header(int options[], int mem_overlapped,
                          int mem_orig, int mem_factor)
{
if ((options[AZ_overlap] < 1) && 
    (options[AZ_subdomain_solve] != AZ_ilut)) return;
   if ((options[AZ_output] != AZ_none ) && (options[AZ_output] != AZ_warnings)){
      printf("\n\t\t*******************************************************\n");
      if (options[AZ_overlap] > 0) {
         printf("\t\t*****       Subdomain overlapping requires %.3e times\n", 
                ((double) mem_overlapped)/ ((double) mem_orig));
         printf("\t\t*****       the memory used for the nonoverlapped\n");
         printf("\t\t*****       subdomain matrix.\n");
      }
      if (options[AZ_subdomain_solve] == AZ_ilut) {
         printf("\t\t***** ilut: The ilut factors require %.3e times \n\t\t", 
                 ((double) mem_factor)/((double) mem_overlapped));
         printf("*****       the memory of the overlapped subdomain matrix.");
      }
      printf("\n\t\t*******************************************************\n");