az_subdomain_solver.c

并行解法器,功能强大

/*====================================================================
 * ------------------------
 * | CVS File Information |
 * ------------------------
 *
 * $RCSfile: az_subdomain_solver.c,v $
 *
 * $Author: tuminaro $
 *
 * $Date: 2000/06/02 16:48:32 $
 *
 * $Revision: 1.18 $
 *
 * $Name:  $
 *====================================================================*/
#ifndef lint
static char rcsid[] = "$Id: az_subdomain_solver.c,v 1.18 2000/06/02 16:48:32 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"
/* Begin Aztec 2.1 mheroux mod */
#ifdef IFPACK
#include "az_ifpack.h"
#endif
/* End Aztec 2.1 mheroux mod */

/*
 * 
 * These routines contain all the solver specific stuff for solving on subdomains.
 *
 * Note: See az_solve.c for an example of using these routines for an lu solver.
 *
 */
/****************************************************************************/
/****************************************************************************/
/****************************************************************************/

void AZ_factor_subdomain(struct context *context, int N, int N_nz,
	int *nz_used)
{
/****************************************************************************
  Given an overlapped subdomain matrix, factor it according to the
  chosen algorithm and store the result back in subdomain. Additionally, 
  store the number of nonzeros used in the factorization in nz_used.

  Notes:
    1) Matrix comes in as an MSR matrix.
    2) context contains several fields which need to be appropriately
       set. These fields are specific to the individual solvers.
    3) The factorization overwrites the matrix. However, different
       solvers will store the factorization in different formats.

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

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

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

  context        On input, context contains the matrix to be 
                   factored in context.A_overlapped (MSR format), 
                   On output, context contains the factored matrix
                   which is stored in a format specific to the solver and 
                   any additional parameters required by the backsolver.

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

  N_nz             On input, the number of nonzero values in the matrix
                   to be factored.

  nz_used          On output, the number of nonzero values in the matrix
                   representing the factorization.

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

	int i, j, *bindx, *bpntr, *iw, ifail;
        double *fake_rhs, *aflag;
        double *cr, *unorm, *a, *val;
        int    *ind, *jnz, *ja, *rnr, ifill;
        double dtemp = (context->aztec_choices->params)[AZ_omega];
        int    N_blk_rows, name = context->A_overlapped->data_org[AZ_name], 
               N_nz_matrix;
        char   str[80];

/* Begin Aztec 2.1 mheroux mod */
#ifdef IFPACK
	void *precon, *bmat;
	double rthresh, athresh;
	int N_int_blk, N_bord_blk, graph_fill;
#endif
/* End Aztec 2.1 mheroux mod */

        bindx = context->A_overlapped->bindx;

        *nz_used = bindx[N];

        switch(context->aztec_choices->options[AZ_subdomain_solve]) {
/* Begin Aztec 2.1 mheroux mod */
        case AZ_bilu_ifp:
#ifdef IFPACK
           if (N == 0) return;

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

           /* for bilu(k) with k > 1 , figure out the new sparsity pattern */

           AZ_sort_msr(bindx, val, N);

           /* Let IFPACK handle fillin */
	      graph_fill = (context->aztec_choices->options)[AZ_graph_fill];
           (context->aztec_choices->options)[AZ_graph_fill] = 0;
           /* recover some space so that there will */
           /* be enough room to convert back to vbr */

           i = AZ_compress_msr(&(context->A_overlapped->bindx), 
                         &(context->A_overlapped->val), context->N_nz_allocated,
                         *nz_used, name, context);
           context->N_nz = *nz_used;
           context->N_nz_allocated = *nz_used;

           AZ_msr2vbr_mem_efficient(N, &(context->A_overlapped->bindx), 
                                 &(context->A_overlapped->val), 
                                 &(context->A_overlapped->cpntr), 
                                 &(context->A_overlapped->bpntr), 
                                 &(context->A_overlapped->indx), &N_blk_rows, 
                                 (context->A_overlapped->data_org)[AZ_name],
                                 context->tag,i);

	   context->A_overlapped->matrix_type = AZ_VBR_MATRIX;
   
	   /*ifp_initialize();*/
  
	   /* Create IFPACK encapsulation of Amat */

	   context->A_overlapped->rpntr = context->A_overlapped->cpntr;
	   N_int_blk = context->A_overlapped->data_org[AZ_N_int_blk];
	   N_bord_blk = context->A_overlapped->data_org[AZ_N_bord_blk];
	   context->A_overlapped->data_org[AZ_N_int_blk] = N_blk_rows;
	   context->A_overlapped->data_org[AZ_N_bord_blk] = 0;
	   (context->aztec_choices->options)[AZ_graph_fill] = graph_fill;

	   az2ifp_blockmatrix(&bmat, context->A_overlapped); 

	   context->A_overlapped->data_org[AZ_N_int_blk] = N_int_blk;
	   context->A_overlapped->data_org[AZ_N_bord_blk] = N_bord_blk;

	   rthresh =  (context->aztec_choices->params)[AZ_rthresh];
	   athresh =  (context->aztec_choices->params)[AZ_athresh];
           ifill = (context->aztec_choices->options)[AZ_graph_fill];
	   ifp_preconditioner(&precon, bmat, IFP_BILUK, (double) ifill, 0.0,
			    IFP_SVD, rthresh, athresh);
	   context->precon = precon;
           break;

/* End Aztec 2.1 mheroux mod */

#else
        AZ_perror("IFPACK not linked.  Must compile with -DIFPACK");
#endif
        case AZ_bilu:
           if (N == 0) return;

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

           /* for bilu(k) with k > 1 , figure out the new sparsity pattern */

           AZ_sort_msr(bindx, val, N);
           ifill = (context->aztec_choices->options)[AZ_graph_fill];
           if (ifill > 0) {
              *nz_used = AZ_fill_sparsity_pattern(context, ifill, 
                                                  bindx, val, N);

           }
           /* recover some space so that there will */
           /* be enough room to convert back to vbr */

           i = AZ_compress_msr(&(context->A_overlapped->bindx), 
                         &(context->A_overlapped->val), context->N_nz_allocated,
                         *nz_used, name, context);
           context->N_nz = *nz_used;
           context->N_nz_allocated = *nz_used;

           AZ_msr2vbr_mem_efficient(N, &(context->A_overlapped->bindx), 
                                 &(context->A_overlapped->val), 
                                 &(context->A_overlapped->cpntr), 
                                 &(context->A_overlapped->bpntr), 
                                 &(context->A_overlapped->indx), &N_blk_rows, 
                                 (context->A_overlapped->data_org)[AZ_name],
                                 context->tag,i);

	   context->A_overlapped->matrix_type = AZ_VBR_MATRIX;
   
           bindx = context->A_overlapped->bindx;
           bpntr = context->A_overlapped->bpntr;
           val   = context->A_overlapped->val;

	   sprintf(str,"ipvt %s",context->tag);
           context->ipvt  = (int *) AZ_manage_memory((N+1)*sizeof(int),
                                    AZ_ALLOC, name, str, &i);
           sprintf(str,"dblock %s",context->tag);
           context->dblock= (int *) AZ_manage_memory((N_blk_rows+1)*
                                                 sizeof(int), AZ_ALLOC, name,
                                                 str, &i);

           context->N_blk_rows = N_blk_rows;

           /* set dblock to point to the diagonal block in each block row */

           for (i = 0 ; i < N_blk_rows ; i++ ) {
              for (j = bpntr[i] ; j < bpntr[i+1] ; j++ ) {
                 if (bindx[j] == i) context->dblock[i] = j;
              }
           }

           AZ_fact_bilu(N_blk_rows, context->A_overlapped, context->dblock,
                        context->ipvt);
           break;

	case AZ_ilut:
           cr = (double *) AZ_allocate((2*N+3+context->max_row)*sizeof(int)+
                                     (2*N+2+context->max_row)*sizeof(double));
           if (cr == NULL) AZ_perror("Out of space in ilut.\n");
           unorm = &(cr[N+2]);
           a     = &(unorm[N]);
           ind   = (int *) &(a[context->max_row]);
           jnz   = &(ind[N+3]);
           ja    = &(jnz[N]);
           sprintf(str,"iu %s",context->tag);
           context->iu    = (int *) AZ_manage_memory((N+1)*sizeof(int),
                                             AZ_ALLOC, name, str, &i);
           AZ_fact_ilut(&N, context->A_overlapped, a, ja,
                        (context->aztec_choices->params)[AZ_drop], 
                        context->extra_fact_nz_per_row, N_nz - bindx[N],
                        context->iu,cr,unorm,ind, nz_used, jnz);
           AZ_free(cr);
           break;
	case AZ_ilu:
           dtemp = 0.0;
	case AZ_rilu:
           if (N == 0) return;
           sprintf(str,"iu %s",context->tag);
           bindx = context->A_overlapped->bindx;
           val   = context->A_overlapped->val;

           /* for ilu(k) with k > 1 , figure out the new sparsity pattern */

           AZ_sort_msr(bindx, val, N);
           ifill = (context->aztec_choices->options)[AZ_graph_fill];
           if (ifill > 0) {
              *nz_used = AZ_fill_sparsity_pattern(context, ifill, 
                                                  bindx, val, N);
           }
           context->iu= (int *) AZ_manage_memory((N+1)*sizeof(int),AZ_ALLOC,
                                                    name, str, &i);
           iw = (int *) AZ_allocate((N+1)*sizeof(int));
           if (iw == NULL) AZ_perror("Out of space in ilu.\n");
           AZ_fact_rilu(N, nz_used, context->iu, iw, context->A_overlapped, 
                        dtemp);
           AZ_free(iw);
           break;
	case AZ_icc:
           sprintf(str,"iu %s",context->tag);
           bindx = context->A_overlapped->bindx;
           val   = context->A_overlapped->val;

           /* for ilu(k) with k > 1 , figure out the new sparsity pattern */

           AZ_sort_msr(bindx, val, N);
           ifill = (context->aztec_choices->options)[AZ_graph_fill];
           if (ifill > 0)
              *nz_used = AZ_fill_sparsity_pattern(context, ifill, 
                                                  bindx, val, N);

           AZ_fact_chol(context->A_overlapped->bindx,
                        context->A_overlapped->val,N);
           break;
	case AZ_lu:
           if (N == 0) return;
           aflag = (double *) AZ_allocate(8*sizeof(double));
           rnr   = (int *) AZ_allocate(N_nz*sizeof(int));
           if (rnr == NULL) AZ_perror("Out of space in lu.\n");

           sprintf(str,"iflag %s",context->tag);
           context->iflag = (int *) AZ_manage_memory(10*sizeof(int), AZ_ALLOC,
                                                       name, str ,&i);