az_subdomain_solver.c

并行解法器,功能强大

AZ_PRECOND *sub_precond;
struct AZ_SCALING *sub_scaling;
#ifdef AZ_MPI
MPI_AZComm  *tptr;
#endif
double *y;
char label[80];
int  t1, t2, t3, i, t4, t5 = 0;

/* Begin Aztec 2.1 mheroux mod */
#ifdef IFPACK
  int ione = 1;
  void *precon;
#endif
/* End Aztec 2.1 mheroux mod */

   val2   = context->A_overlapped->val;
   bindx2 = context->A_overlapped->bindx;

   switch(context->aztec_choices->options[AZ_subdomain_solve]) {

/* Begin Aztec 2.1 mheroux mod */

   case AZ_bilu_ifp:
#ifdef IFPACK
     y = (double *) malloc (N * sizeof(double));
     dcopy_(&N, x, &ione, y, &ione);
     precon = context->precon;
     ifp_apply(precon, N, 1, y, N, x, N);
     free((void *) y);
#endif
     break;

/* End Aztec 2.1 mheroux mod */

   case AZ_bilu:
      N_blk_rows = context->N_blk_rows;

      AZ_lower_triang_vbr_solve(N_blk_rows, context->A_overlapped->cpntr, 
                                context->A_overlapped->bpntr, 
				context->A_overlapped->indx,
                                bindx2, val2, x);

      AZ_upper_triang_vbr_solve(N_blk_rows, context->A_overlapped->cpntr,
                                context->A_overlapped->bpntr, 
				context->A_overlapped->indx, bindx2,
                                val2, x, context->ipvt, context->dblock);
      break;
   case AZ_ilut:
   case AZ_rilu:
   case AZ_ilu:
      AZ_lower_tsolve(x,N, val2, bindx2, context->iu, x ); 
      AZ_upper_tsolve( x, N, val2, bindx2, context->iu);
      break;
   case AZ_icc:
      AZ_lower_icc(bindx2,val2,N,x);
      AZ_upper_icc(bindx2,val2,N,x);
      break;
   case AZ_lu:
      if (N == 0) return;
      else if (N== 1) {
         x[0] *= val2[0];
         ifail = 0;
      }
      else AZ_backsolve(val2, context->pivot,x, bindx2, 
	              context->ha, context->iflag, 
                      &ifail, &(context->N_nz_factors),
		      &N, &N);
      break;
   default: 
      if (context->aztec_choices->options[AZ_subdomain_solve]
                  >= AZ_SOLVER_PARAMS) {
         printf("ERROR: Unknown subdomain solver %d\n",
                context->aztec_choices->options[AZ_subdomain_solve]);
         exit(1);
       }
       else {
          /* better to put most of this in the factorization */

          AZ_recover_sol_params(context->aztec_choices->options[
			        AZ_subdomain_solve], &sub_options, 
				&sub_params, &sub_status, &sub_matrix, 
			        &sub_precond, &sub_scaling);
          t1 = sub_options[AZ_recursion_level];
          sub_options[AZ_recursion_level]++;

          t2 = sub_options[AZ_output];
          if (context->proc_config[AZ_node] != 0 ) 
             sub_options[AZ_output] = AZ_none;

          t3 = AZ_sys_msg_type;

          /* fix data_org */

          hold_data_org = context->A_overlapped->data_org;
          new_data_org = (int *) AZ_allocate( sizeof(int) * AZ_send_list );
          if (new_data_org == NULL) {
             printf("Error: Not enough space for subdomain matrix\n");
             exit(1);
          }
          context->A_overlapped->data_org = new_data_org;
          context->A_overlapped->matvec = AZ_MSR_matvec_mult;
          new_data_org[AZ_matrix_type] = AZ_MSR_MATRIX;
          new_data_org[AZ_N_internal]  = N;
          new_data_org[AZ_N_border  ]  = 0;
          new_data_org[AZ_N_external]  = 0;
          new_data_org[AZ_N_int_blk ]  = N;
          new_data_org[AZ_N_bord_blk]  = 0;
          new_data_org[AZ_N_ext_blk ]  = 0;
          new_data_org[AZ_N_neigh   ]  = 0;
          new_data_org[AZ_total_send]  = 0;
          new_data_org[AZ_name      ]  = hold_data_org[AZ_name];
          new_data_org[AZ_internal_use]= 0;
          new_data_org[AZ_N_rows      ]= N;
          sub_precond->Pmat = context->A_overlapped;
          sub_precond->prec_function = AZ_precondition;
       
          sub_proc_config[AZ_node] = 0;
          sub_proc_config[AZ_N_procs] = 1;
#ifdef AZ_MPI
          tptr = AZ_get_comm(context->proc_config);
          AZ_set_comm(sub_proc_config, *tptr);
#endif

          sprintf(label,"y in ssolve%d", sub_options[AZ_recursion_level]);
          y = AZ_manage_memory((N+1)*sizeof(double),
                             AZ_ALLOC, AZ_SYS, label, &i);

          for (i = 0 ; i < N ; i++ ) y[i] = x[i];
          for (i = 0 ; i < N ; i++ ) x[i] = 0.0;

          t4 = sub_options[AZ_keep_info];
          sub_options[AZ_keep_info] = 1;

          if (context->aztec_choices->options[AZ_pre_calc] >= AZ_reuse) {
             t5 = sub_options[AZ_pre_calc];
             sub_options[AZ_pre_calc] = AZ_sys_reuse;
          }

          AZ_oldsolve(x, y,sub_options,sub_params, sub_status, sub_proc_config,
                       context->A_overlapped, sub_precond, sub_scaling);

          sub_options[AZ_keep_info] = t4;
          if (context->aztec_choices->options[AZ_pre_calc] == AZ_sys_reuse) 
             sub_options[AZ_pre_calc]  = t5;

          sub_options[AZ_recursion_level] = t1;
          sub_options[AZ_output] = t2;
          context->A_overlapped->data_org = hold_data_org;
          AZ_free(new_data_org);
          AZ_sys_msg_type = t3;
       }
   }
      
}


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

void AZ_space_for_factors(double input_fill, int N_nz, int N, 
	int *extra_factor_nonzeros, int options[],int bandwidth,
        int  max_nz_per_row)
{
/****************************************************************************

  Compute the additional number of nonzeros required to do the factorization
  specified by options[AZ_subdomain_solve]. 
 
  Author:          Ray Tuminaro, SNL, 9222
 
  Return code:     void
  ============
 
  Parameter list:
  ===============
 
  input_fill:      On input, input_fill*N_nz is roughly the number of
                   nonzeros that will be allowed in the matrix factors
                   for ilut.

  N_nz:            On input, number of nonzeros in the padded matrix
                   that will be factored.
  
  N:               On input, the order of the padded matrix to be factored.

  extra_factor_nonzeros: 
                   On output, the additional space that will be added to
                   accommodate fill-in during the factorization.

  options:         On input, user specified input options.

 */

  int fill, i;
  double new_nz, nz_per_row, t;
  int    max_per_row, temp;
  
 
  if (options[AZ_subdomain_solve] == AZ_ilut) {
     input_fill -= 1.0;
     if (N == 0) *extra_factor_nonzeros = 0;
     else {

        new_nz = input_fill * ((double) N_nz);
        nz_per_row = new_nz/((double) N);
        fill = (int) floor( nz_per_row/2. + .5);
        t = ((double) N_nz)/ ((double) N);
        max_per_row = N - (int) ceil(t);
        if ( 2*fill > max_per_row) fill = max_per_row/2;

        *extra_factor_nonzeros =  2*N*fill + 1;
         while ( *extra_factor_nonzeros < 0) {
           fill--; *extra_factor_nonzeros =  2*N*fill + 1;
        }
     }
     temp = N*bandwidth;
     if (temp < 0) temp = *extra_factor_nonzeros;
     if (temp < *extra_factor_nonzeros) *extra_factor_nonzeros = temp;
  }
/* Begin Aztec 2.1 mheroux mod */
  else if ((options[AZ_subdomain_solve] == AZ_rilu) ||
           (options[AZ_subdomain_solve] == AZ_ilu ) ||
           (options[AZ_subdomain_solve] == AZ_icc ) ||
           (options[AZ_subdomain_solve] == AZ_bilu_ifp) ||
           (options[AZ_subdomain_solve] == AZ_bilu) ){
     fill = options[AZ_graph_fill];
/* End Aztec 2.1 mheroux mod */
     if (fill < 0) {
         printf("options[AZ_graph_fill] must be greater or equal to 0\n");
         printf("when using an incomplete factorization\n");
         exit(1);
     }
     if (fill == 0) *extra_factor_nonzeros = 3;
     else {
        temp = max_nz_per_row;
        for (i = 0 ; i < fill ; i++) { 
           temp *= max_nz_per_row;
           if (temp > bandwidth) break;
        }
        if (temp > bandwidth) temp = bandwidth;
        temp = temp*N;
        *extra_factor_nonzeros = temp - N_nz + 200;
     }
  }
  else if (options[AZ_subdomain_solve] == AZ_lu) {
     *extra_factor_nonzeros = N*bandwidth - N_nz + 200;
                             /* for small matrices y12m seems to need some    */
                            /* additional space. It might use a dense solver */
                            /* in some case???? who knows...                 */
  }
  else *extra_factor_nonzeros = 1;

  /* make sure things don't overflow */

  temp = 2*(N_nz + *extra_factor_nonzeros)*sizeof(double);
  if ( temp < 0 ) {
      temp = 2;
      while ( temp < (2*temp) ) temp = 2*temp;
      *extra_factor_nonzeros = temp/(2*sizeof(double));
  }

}

void AZ_zero_out_context(struct context *context)
{

   context->iu     = NULL;
   context->iflag  = NULL;         
   context->ha     = NULL;
   context->ipvt   = NULL;
   context->dblock = NULL; 
   context->space_holder  = NULL;
   context->pivot  = NULL;
   context->A_overlapped  = NULL;
   context->aztec_choices = NULL;
   context->x_pad         = NULL; 
   context->ext_vals      = NULL; 
   context->x_reord       = NULL;
   context->padded_data_org = NULL;
   context->map           = NULL;
   context->inv_ordering  = NULL;
   context->tag           = NULL;
   context->proc_config   = NULL;

   context->extra_fact_nz_per_row = 0;
   context->N_large_int_arrays = 0;
   context->N_large_dbl_arrays = 0;
   context->N_nz_factors = 0;
   context->N_nz_matrix = 0; 
   context->N_blk_rows = 0;
   context->max_row = 0;
   context->N = 0;
   context->N_unpadded = 0;
   context->N_nz = 0; 
   context->Pmat_computed = 0;
}