az_solve.c

并行解法器,功能强大

     if ((options[AZ_pre_calc] == AZ_reuse) && (proc_config[AZ_node] == 0)){
        fprintf(stderr,"Error:\tDid not find previous factorization (");
        fprintf(stderr, "requested \n\tby setting options[AZ_pre_calc] to ");
        fprintf(stderr, "AZ_reuse).\n\tTo find this factorization, the ");
        fprintf(stderr, "following\n\tparameters must match the previous");
        fprintf(stderr," factorization:");
        fprintf(stderr, "\n\t\t 1) Total number of unknowns.");
        fprintf(stderr, "\n\t\t 2) options[AZ_overlap]");
        fprintf(stderr, "\n\t\t 3) options[AZ_scaling]");
        fprintf(stderr, "\n\t\t 4) options[AZ_precond]");
        fprintf(stderr, "\n\t\t 5) options[AZ_reorder]");
        fprintf(stderr, "\n\t\t 6) options[AZ_type_overlap]");
        fprintf(stderr, "\n\t\t 7) options[AZ_subdomain_solve]");
        fprintf(stderr, "\n\t\t 8) options[AZ_graph_fill]");
        fprintf(stderr, "\n\t\t 9) params[AZ_ilut_fill]");
        fprintf(stderr, "\n\t\t10) params[AZ_drop]");
        fprintf(stderr, "\n\t\t11) data_org[AZ_name]\n");
printf("XXX%sXXX %d %d\n",tag,data_org[AZ_name],(int) sizeof(struct context));
(void) AZ_manage_memory(0, -43, AZ_SYS, (char *) 0, (int *) 0);
     }
     if (options[AZ_pre_calc] == AZ_reuse) exit(1);

     tag[0] = 'T';
     precond->context->tag = (char *) AZ_manage_memory(sizeof(char)*80,
					      AZ_ALLOC, data_org[AZ_name],
                                              tag,&istatus);
     tag[0] = 'P';
     sprintf(precond->context->tag,"%s",tag);
  }
}
/************************************************************************/
/************************************************************************/
/************************************************************************/

void AZ_rm_context(int options[], double params[], int data_org[])
{
/* Remove context associated with the preconditioning data structure.
 * See 'AZ_mk_context' for more details.
 *
 *
 */
 char tag[80];
 int  istatus;


  AZ_mk_identifier(params,options,data_org, tag);

  tag[0] = 'T';
  AZ_manage_memory(sizeof(char)*80, AZ_SELECTIVE_CLEAR,
                            data_org[AZ_name], tag,&istatus);
  tag[0] = 'P';
  AZ_manage_memory(sizeof(struct context), AZ_SELECTIVE_CLEAR,
                            data_org[AZ_name],tag,&istatus);
}
     
int AZ_oldsolve_setup(double x[], double b[], int options[], double params[], 
	double status[], int proc_config[], AZ_MATRIX *Amat, 
	AZ_PRECOND *precond, int save_old_values[],
	struct AZ_SCALING *scaling)


/*******************************************************************************
  Aztec's previous AZ_solve() is renamed to AZ_oldsolve() with 3 new
  parameters appended to it: Amat, precond, scaling.
  This routine is never called directly by an application. It is only
  used internally by Aztec.

  Solve the system of equations given in the VBR format using an iterative
  method specified by 'options[AZ_solver]'. Store the result in 'x'.

  Author:          Ray Tuminaro, SNL, 1422
  =======

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

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

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

  b:               Right hand side of linear system.

  options:         Determines specific solution method and other parameters.

  params:          Drop tolerance and convergence tolerance info.

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

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

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

  status:          On output, indicates termination status:
                    0:  terminated normally.
                   -1:  maximum number of iterations taken without achieving
                        convergence.
                   -2:  Breakdown. The algorithm can not proceed due to
                        numerical difficulties (usually a divide by zero).
                   -3:  Internal residual differs from the computed residual due
                        to a significant loss of precision.

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

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

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

  scaling:         Structure used to represent the  scaling
                   (see az_aztec.h ad Aztec User's Guide).


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

{

  /* local variables */

  int     i;
  int     Proc, Num_Proc;
  int     *data_org;
  int     prec_context, *toptions;
  int *tdata_org;
  AZ_PRECOND *tprecond;
  double *tparams;
  char tag[80];
  

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

  if (options[AZ_recursion_level] == 0) precond->timing[0] = 0.0;
#ifdef AZ_MPI
  if ( proc_config[AZ_Comm_Set] != AZ_Done_by_User) {
      printf("Error: Communicator not set. Use AZ_set_comm()\n");
      printf("       (e.g. AZ_set_comm(proc_config,MPI_COMM_WORLD)).\n");
      exit(1);
  }
#endif

  /* save a few variables that Aztec will mess with so that */
  /* we can restore them before returning to the user.      */


  save_old_values[AZ_oldprecond] = options[AZ_precond];
  save_old_values[AZ_oldsubdomain_solve] = options[AZ_subdomain_solve];
  save_old_values[AZ_oldreorder] = options[AZ_reorder];
  save_old_values[AZ_oldpre_calc] = options[AZ_pre_calc];
  i  = 1;
  if (options[AZ_precond] == AZ_lu) sprintf(tag,"AZ_lu");
  else if (options[AZ_precond] == AZ_ilu ) sprintf(tag,"AZ_ilu");
  else if (options[AZ_precond] == AZ_ilut) sprintf(tag,"AZ_ilut");
  else if (options[AZ_precond] == AZ_bilu) sprintf(tag,"AZ_bilu");
/* Begin Aztec 2.1 mheroux mod */
  else if (options[AZ_precond] == AZ_bilu_ifp) sprintf(tag,"AZ_bilu_ifp");
/* End Aztec 2.1 mheroux mod */
  else if (options[AZ_precond] == AZ_rilu) sprintf(tag,"AZ_rilu");
  else i = 0;
  if (i == 1) {
     printf("To use this preconditioner, you must now set\n");
     printf("      options[AZ_precond] = AZ_dom_decomp;\n");
     printf("      options[AZ_subdomain_solve] = %s;\n",tag);
     exit(1);
  }
     

  data_org = Amat->data_org;

  Proc         = proc_config[AZ_node];
  Num_Proc     = proc_config[AZ_N_procs];

  data_org[AZ_internal_use] = 1;

  /* check for inconsistent user input options */

  if (!AZ_check_options(options, Proc, data_org, Num_Proc, params,
      Amat, precond)) {
 
    status[AZ_its]      = (double )  0.0;
    status[AZ_why]      = (double )  AZ_param;
    status[AZ_r]        = (double ) -1.0;
    status[AZ_rec_r]    = (double ) -1.0;
    status[AZ_scaled_r] = (double ) -1.0;
    return(0);
  }

  /* Test to see if we should capture matrix, rhs and partitioning info
     in an ASCII data file.  If the file "AZ_write_matrix_now" exists in
     the current working directory, then the files 
     - AZ_capture_matrix.dat
     - AZ_capture_rhs.dat
     - AZ_capture_partition.dat (VBR only)

     will be appended with the current matrix in (i,j,val) format, the
     current RHS and the current partition information.  The existence
     of "AZ_write_matrix_now" is check each time.  Thus, capturing can
     be turned on and off at will during the run of a simulation.
  */

  AZ_capture_matrix(Amat->val, Amat->indx, Amat->bindx, Amat->rpntr, 
		    Amat->cpntr, Amat->bpntr, proc_config,
		    data_org, b);


  /* If desired, print out the matrix and indexing arrays */

  if ((data_org[AZ_matrix_type] == AZ_MSR_MATRIX) ||
      (data_org[AZ_matrix_type] == AZ_VBR_MATRIX)){
     if (options[AZ_output] == AZ_all)
        AZ_print_out((int *) NULL, (int *) NULL, (int *) NULL, (int *) NULL,
                     Amat->val, Amat->indx, Amat->bindx, Amat->rpntr, 
                     Amat->cpntr, Amat->bpntr, proc_config, AZ_input_form, 
                     data_org[AZ_matrix_type], data_org[AZ_N_int_blk] + 
                     data_org[AZ_N_bord_blk],data_org[AZ_N_ext_blk],0);
  }

  toptions  = options;
  tparams   = params;
  tprecond  = precond;
  tdata_org = precond->Pmat->data_org;

  prec_context = 1;
  if (options[AZ_precond] == AZ_multilevel) 
     prec_context = AZ_multilevel;

  while ( prec_context ) {

     /* adjust print frequency according to user input */

     i = toptions[AZ_output];

     if  (i == AZ_warnings) toptions[AZ_print_freq] =toptions[AZ_max_iter] + 10;
     else if (i == AZ_none) toptions[AZ_print_freq] =toptions[AZ_max_iter] + 10;
     else if (i == AZ_all ) toptions[AZ_print_freq] = 1;
     else if (i == AZ_last) toptions[AZ_print_freq] =toptions[AZ_max_iter] + 1;
     else                   toptions[AZ_print_freq] =toptions[AZ_output];
     if ((i != AZ_none) && (i != AZ_warnings) &&
         (proc_config[AZ_node] == 0)) printf("\n");

     /*  Setup data structure to record domain decomposition factorization */
     /*  information that can be used in future solves.                    */
     /*  Note: If this information already exists, we recover the previous */
     /*        information.                                               */

     AZ_mk_context(toptions, tparams, tdata_org, tprecond, proc_config);
     
     /* check if multiple preconditioners are being passed in */
     /* if this is the case we must afix a subdomain context  */
     /* to each one.                                          */

     if ( prec_context == AZ_multilevel) {
       if (tprecond->next_prec == NULL) prec_context = 0;
       else {
          tprecond = tprecond->next_prec; 
          toptions = tprecond->options;
          tparams  = tprecond->params;
          tdata_org= tprecond->Pmat->data_org;
       }
     } 
     else prec_context = 0;
  }

  /* scale matrix, rhs and initial guess if required */

   AZ_scale_f(AZ_SCALE_MAT_RHS_SOL, Amat, options, b, x, proc_config, scaling);
   return(1);
}



void AZ_oldsolve_finish(double x[], double b[], int options[], 
	int proc_config[], AZ_MATRIX *Amat, int save_old_values[],
	struct AZ_SCALING *scaling)
{


  AZ_scale_f(AZ_INVSCALE_SOL, Amat, options, b, x, proc_config, scaling);
  AZ_scale_f(AZ_INVSCALE_RHS, Amat, options, b, x, proc_config, scaling);

  options[AZ_precond]         = save_old_values[AZ_oldprecond];
  options[AZ_subdomain_solve] = save_old_values[AZ_oldsubdomain_solve];
  options[AZ_reorder]         = save_old_values[AZ_oldreorder];
  options[AZ_pre_calc]        = save_old_values[AZ_oldpre_calc];

  fflush(stdout);

} /* AZ_oldsolve_finish  */

void AZ_iterate_setup(int options[], double params[], int proc_config[], 
	AZ_MATRIX *Amat, AZ_PRECOND *precond)

/*******************************************************************************
This is the new Aztec interface. This routine calls AZ_oldsolve() passing
in for example Amat->indx for the indx[] parameter in AZ_oldsolve().

NOTE: User's can still invoke AZ_solve() in the old Aztec way. AZ_solve
      also calls AZ_oldsolve(). However, matrix-free and coarse grid
       capabilities are not available via AZ_solve().


  Author:          Ray Tuminaro, SNL, 1422
  =======

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

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


  x:               Contains the result vector x.

  b:               Contains the vector b.

  options:         Determines specific solution method and other parameters.

  params:          Drop tolerance and convergence tolerance info.

  status:          On output, indicates termination status:
                    0:  terminated normally.
                   -1:  maximum number of iterations taken without achieving
                        convergence.
                   -2:  Breakdown. The algorithm can not proceed due to
                        numerical difficulties (usually a divide by zero).
                   -3:  Internal residual differs from the computed residual due
                        to a significant loss of precision.

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

  Amat:            Structure used to represent the matrix (see az_aztec.h
                   and 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.

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

  scaling:         Structure used to represent the  scaling
                   (see az_aztec.h ad Aztec User's Guide).

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

{

  if (Amat->matrix_type == AZ_MSR_MATRIX) {
      Amat->matvec = AZ_MSR_matvec_mult;
      Amat->data_org[AZ_N_int_blk] = Amat->data_org[AZ_N_internal];
      Amat->data_org[AZ_N_bord_blk] = Amat->data_org[AZ_N_border];
  }

  else if (Amat->matrix_type == AZ_VBR_MATRIX) 
      Amat->matvec = AZ_VBR_matvec_mult;

  else if (Amat->matrix_type == AZ_USER_MATRIX){
      Amat->data_org[AZ_N_int_blk] = Amat->data_org[AZ_N_internal];
      Amat->data_org[AZ_N_bord_blk] = Amat->data_org[AZ_N_border];
  }


  Amat->data_org[AZ_matrix_type] = Amat->matrix_type;

/*
  if (options[AZ_precond] != AZ_user_precond)
      precond->prec_function = AZ_precondition;
*/

  options[AZ_recursion_level] = 0;

  if (options[AZ_pre_calc] != AZ_reuse) {
     (void) AZ_manage_memory(0,AZ_EVERYBODY_BUT_CLEAR,(Amat->data_org)[AZ_name],"kvecs",
                             (int *) 0);
     (void) AZ_manage_memory(0,AZ_EVERYBODY_BUT_CLEAR,(precond->Pmat->data_org)[AZ_name],
                             "kvecs", (int *) 0);
  }
  (void) AZ_manage_memory(0, AZ_CLEAR, AZ_SYS, (char *) 0, (int *) 0);

  /* output solver, scaling, and preconditioning options */

  AZ_print_call_iter_solve(options, params, proc_config[AZ_node], 0, precond);

}

void AZ_iterate_finish(int options[], AZ_MATRIX *Amat, AZ_PRECOND *precond)
{

  if (options[AZ_keep_info] == 0) {
      (void) AZ_manage_memory(0,AZ_CLEAR,(Amat->data_org)[AZ_name],
                              (char *) 0,(int *) 0);
      (void) AZ_manage_memory(0,AZ_CLEAR,(precond->Pmat->data_org)[AZ_name],
                              (char *) 0,(int *) 0);
  }

  (void) AZ_manage_memory(0, AZ_CLEAR, AZ_SYS, (char *) 0, (int *) 0);
}

int AZ_initialize(double x[], double b[], int options[], 
	double params[], double status[], int proc_config[], AZ_MATRIX *Amat, 
	AZ_PRECOND *precond, int save_old_values[], struct AZ_SCALING *scaling)
{
   AZ_iterate_setup(options, params, proc_config, Amat, precond);
   return(AZ_oldsolve_setup(x, b, options, params, status, proc_config, Amat, 
	             precond, save_old_values, scaling));
   
}

void AZ_finalize(double x[], double b[], int options[], int 
   proc_config[], AZ_MATRIX *Amat, AZ_PRECOND *precond, int save_old_values[],
   struct AZ_SCALING *scaling)
{
   AZ_oldsolve_finish(x, b, options, proc_config, Amat, save_old_values,scaling);
   AZ_iterate_finish(options, Amat, precond);
}