az_subdomain_solver.c

并行解法器,功能强大

           sprintf(str,"ha %s",context->tag);
           context->ha = (int *) AZ_manage_memory(11*(N+1)*sizeof(int),
                                             AZ_ALLOC, name, str, &i);
           sprintf(str,"pivot %s",context->tag);
           context->pivot = (double *) AZ_manage_memory((N+1)*sizeof(double),
                                             AZ_ALLOC, name, str,&i);

           aflag[0] = 16.0;    aflag[2] = 1.0e8;   aflag[3] = 1.0e-12;   
           aflag[1] = (context->aztec_choices->params)[AZ_drop];

           /* set up flags for the sparse factorization solver */

           context->iflag[0] = 1;         context->iflag[1] = 2;
           context->iflag[2] = 1;         context->iflag[3] = 0;
           context->iflag[4] = 2;    
           /*    Note: if matrix is pos def, iflag[2] = 2 is cheaper */

           N_nz_matrix = bindx[N] - 1;

           AZ_msr2lu(N, context->A_overlapped, rnr);

           /* Mark bindx so we can see what was not used later */

           for (i =  N_nz_matrix ; i < N_nz ; i++) bindx[i] = -7;

           /* factor the matrix */ 

           if (N == 1) {
             context->A_overlapped->val[0]=1./context->A_overlapped->val[0];
           }
           else {
              context->N_nz_factors = N_nz;
              fake_rhs = (double *) AZ_allocate(N*sizeof(double));
              if (fake_rhs == NULL) {
                 printf("Not enough memory inside subdomain_solve\n");
              }
              for (i = 0 ; i < N ; i++ ) fake_rhs[i] = 0.0;
              AZ_fact_lu(fake_rhs, context->A_overlapped,aflag, 
                         context->pivot, rnr, context->ha, 
			 context->iflag, &N_nz_matrix,
                         &ifail, &(context->N_nz_factors),
                         &N, &N);

              (context->iflag)[4] = 3; 
              AZ_free(fake_rhs);

              /* find out what was not used by checking what was not touched */

              *nz_used = N_nz;
              for (i = N_nz_matrix; i < N_nz ; i++ ) {
                 if (bindx[i] != -7) *nz_used = i;
              }
              (*nz_used)++;
              context->N_nz_factors = *nz_used;
           }
           AZ_free(rnr);
           AZ_free(aflag);
           break;
        default:
           if (context->aztec_choices->options[AZ_subdomain_solve]
                  >= AZ_SOLVER_PARAMS) {
              fprintf(stderr,"Unknown subdomain solver(%d)\n",
                   context->aztec_choices->options[AZ_subdomain_solve]);
              exit(1);
           }
        }      
}

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

void AZ_free_space_holder(struct context *context)
{
/****************************************************************************
  This routine is used in conjunction with AZ_hold_space().
  Essentially, this routine deallocates memory allocated via
  AZ_hold_space(). The whole point of these two routines is to
  allocated all the space needed during the factorization process
  EXCLUDING all arrays whose size is related to the number of 
  nonzeros. Once this is done, we can determine how much space 
  there is left for the large arrays required for the factorization
  and split the remaining space amoung these large arrays. In this
  way LU routines where it is difficult to know the space requirements
  ahead of time can try to use as large an array as possible.

  Note: after factorization 'realloc' is used to reduce the array sizes.
  

  Author:          Ray Tuminaro, SNL, 9222 (3/98)

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

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

  context          On input, context->aztec_choices->
		   options[AZ_subdomain_solve]
                   contains the preconditioner choice while
                   context->space_holder holds memory previously
                   allocated via AZ_hold_space().
                   On output, context->space_holder is deallocated.

                 
*******************************************************************************/
    int which = context->aztec_choices->options[AZ_subdomain_solve];

/* Begin Aztec 2.1 mheroux mod */
    if ( (which == AZ_ilut) || (which == AZ_lu ) || (which == AZ_bilu) ||
         (which == AZ_bilu_ifp) ||
         (which == AZ_rilu )|| (which == AZ_ilu) || (which == AZ_icc) )
       AZ_free(context->space_holder);
/* End Aztec 2.1 mheroux mod */
}

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

void AZ_hold_space(struct context *context, int N) 
{
/****************************************************************************
  This routine is used in conjunction with AZ_free_space_holder().
  Essentially, this routine allocates memory while AZ_free_space_holder()
  deallocates.  The whole point of these two routines is to
  allocated all the space needed during the factorization process
  EXCLUDING all arrays whose size is related to the number of 
  nonzeros. Once this is done, we can determine how much space 
  there is left for the large arrays required for the factorization
  and split the remaining space amoung these large arrays. In this
  way LU routines where it is difficult to know the space requirements
  ahead of time can try to use as large an array as possible.

  Note: after factorization 'realloc' is used to reduce the array sizes.

  Author:          Ray Tuminaro, SNL, 9222 (3/98)

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

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

  context          On input, context->aztec_choices->
	           options[AZ_subdomain_solve]
                   indicates the solver being used.
                   On output, context->space_holder points to a
                   block of memory which can hold all the 'nonlarge' arrays
                   required by this solver.
     
  N                On input, the size of the matrix to be allocated.

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


   switch(context->aztec_choices->options[AZ_subdomain_solve]) {
   case AZ_ilut:
      context->space_holder = (int *) AZ_allocate((2*N+2+context->max_row)*
                                                   sizeof(double) + sizeof(int)*
                                                   (3*N+8+context->max_row));
      if (context->space_holder  == NULL) AZ_perror("No space in ilut.\n");

      /*   Space for cr (N+2 doubles), unorm (N doubles), a (max_row doubles),*/
      /*   ind (N+3 ints), jnz (N ints), iu (N+1 ints), ja(max_row ints),     */
      /*   + 4 ints for manage memory header                                  */
      break;
   case AZ_lu:
        context->space_holder = (int *) AZ_allocate((2*N+9)* sizeof(double) + 
                                           (11*(N+1)+22)*sizeof(int) );

      /*   Space for aflag (8 doubles), ifail (10 ints), ha (11*(N+1) ints), */
      /*   pivot (N+1 doubles), fake_rhs (N doubles)+12 ints for manage      */
      /*   memory header                                                     */
      /*   Note: Arrays of size N_nz (e.g. rnr) are not allocated by this    */
      /*   routine. Instead the subdomain field N_int_arrays is set.         */

        if (context->space_holder == NULL) AZ_perror("Out of space in lu.\n");
      break;
   case AZ_bilu:
/* Begin Aztec 2.1 mheroux mod */
   case AZ_bilu_ifp:
/* End Aztec 2.1 mheroux mod */
        context->space_holder= (int *) AZ_allocate((N+1)*sizeof(double));
        if (context->space_holder == NULL) AZ_perror("No space for bilu.\n");

        /* BILU is a little funny in that the maximum amount of memory */
        /* required does not occur during the factorization. Instead   */
        /* it occurs when converting MSR to VBR. At this point, an     */
        /* additional array of length N is needed.                     */
      break;

   case AZ_icc:
       context->space_holder= (int *) AZ_allocate((2*(N+1))*sizeof(int)+
                                                     (N+1)*sizeof(double));
       if (context->space_holder == NULL) AZ_perror("Out of space in ilu.\n");
       break;
   case AZ_ilu:
   case AZ_rilu:
       context->space_holder= (int *) AZ_allocate((2*(N+1)+4)*sizeof(int));

       /*   Space for iu (N+1 ints), iw (N+1 ints) + 4 ints for manage_memory */

       if (context->space_holder == NULL) AZ_perror("Out of space in ilu.\n");
       break;
    default:
      ;
   }
}

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

void AZ_init_subdomain_solver(struct context *context) 
{
/****************************************************************************
  Initialize the data structure 'context' to whatever values will
  be need by the particular routine.

  NOTE: The fields 'N_large_int_arrays' and 'N_large_dbl_arrays' must be 
  set to the number of arrays whose lengthes are equal to the number
  of nonzeros in the factorization. Most solvers will set these equal to
  one indicating that an msr bindx array is needed as well as an msr val
  array. However, the 'lu' solver for example requires an additional integer
  array to hold nonzeros and so N_large_int_arrays is set to 2. These values
  will be used when we try to allocate as much space as is needed to hold
  the factorization.

  Author:          Ray Tuminaro, SNL, 9222 (3/98)

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

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

  context          On output, the various fields are set to solver specific
                   information that is needed in an initialization phase. 

  options          On input, is the user specified algorithm options given to 
                   AZ_solve() or AZ_iterate().

  params           On input, is the user specified algorithm options given to 
                   AZ_solve() or AZ_iterate().

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

   context->N_large_int_arrays = 1;
   context->N_large_dbl_arrays = 1;
   if (context->aztec_choices->options[AZ_subdomain_solve] == AZ_lu)
      context->N_large_int_arrays= 2;
}

extern int AZ_sys_msg_type;

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

void AZ_solve_subdomain(double x[],int N, struct context *context)
{
/****************************************************************************
  Given a vector 'x' representing the right hand side, solve the system
  using whatever subdomain solver is indicated by 'context->which'
  and whatever factorization information has already been computed.

  Author:          Ray Tuminaro, SNL, 9222 (3/98)

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

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

  x                On input, the right hand side of the subdomain system that
                   is to be solved. 
                   On output, the solution of the subdomain system.

  N                On input, the size of the linear system to be solved.

  bindx2,val2      On input, matrix or factorization information to be used 
                   by the solver. For most schemes, this information is in
                   MSR format. However, the lu and bilu scheme would have
                   this information in another format.
                   Note: additional array information can be passed through
                   context.

  context          On input, the various fields are set to solver specific
                   information corresponding to algorithm parameters as
                   well as a previously done factorization.

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

double *val2;
int    *bindx2;
int N_blk_rows;
int ifail;
int *sub_options, sub_proc_config[AZ_PROC_SIZE], *hold_data_org, *new_data_org;
double *sub_params, *sub_status;
AZ_MATRIX *sub_matrix;