az_solve.c

并行解法器,功能强大

  }
  else {
  switch (options[AZ_precond]) {
  case AZ_none:
    (void) printf("No preconditioning"); break;
  case AZ_Neumann:
    (void) printf("Order %d Neumann series polynomial", options[AZ_poly_ord]);
  break;
  case AZ_ls:
    (void) printf("Order %d least-squares polynomial", options[AZ_poly_ord]);
  break;
  case AZ_Jacobi:
    (void) printf("%d step block Jacobi", options[AZ_poly_ord]);
  break;
  case AZ_smoother:
    (void) printf("%d step loc avg smoother", options[AZ_poly_ord]);
  break;
  case AZ_sym_GS:
    (void) printf("%d step symmetric Gauss-Seidel", options[AZ_poly_ord]);
  break;
  case AZ_dom_decomp:
    if (options[AZ_subdomain_solve] == AZ_bilu)
       printf("BILU(%d) domain decomp. with", options[AZ_graph_fill]);
/* Begin Aztec 2.1 mheroux mod */
    else if (options[AZ_subdomain_solve] == AZ_bilu_ifp) {
       printf("IFPACK BILU(%d) ( ATHRESH = %.3e, RTHRESH = %.3e)\n ",
              options[AZ_graph_fill],params[AZ_athresh], params[AZ_rthresh]);
       printf(prefix); printf("with");
    }
/* End Aztec 2.1 mheroux mod */
    else if (options[AZ_subdomain_solve] == AZ_ilut) {
       printf("ILUT( fill-in = %.3e, drop = %.3e)\n ",
              params[AZ_ilut_fill], params[AZ_drop]);
       printf(prefix); printf("with");
    }
    else if (options[AZ_subdomain_solve] == AZ_ilu)
       printf("ILU(%d) domain decomp. with", options[AZ_graph_fill]);
    else if (options[AZ_subdomain_solve] == AZ_rilu)
       printf("RILU(%d,%.2f) domain decomp. with",options[AZ_graph_fill],
              params[AZ_omega]);
    else if (options[AZ_subdomain_solve] == AZ_lu)
       printf("LU domain decomp. with");
    else if (options[AZ_subdomain_solve] == AZ_icc)
       printf("icc(%d) domain decomp. with",options[AZ_graph_fill]);
    else if (options[AZ_subdomain_solve] < AZ_SOLVER_PARAMS)
    (void) printf("iterative subdomain solve with");
    else {
       (void) printf("Unknown subdomain solver (%d)\n",
                  options[AZ_subdomain_solve]);
       exit(1);
     }
     if      (options[AZ_overlap] == AZ_none) printf("out overlap");
     else if (options[AZ_overlap] == AZ_diag) printf(" diagonal overlap");
     else if (options[AZ_type_overlap] == AZ_symmetric) printf(" symmetric");
     if (options[AZ_overlap] >= 1) printf(" overlap = %d ",options[AZ_overlap]);
     super_special = 1;
  break;
  case AZ_icc:
    (void) printf("incomplete Choleski decomposition");
  super_special = 1;
  break;
  case AZ_user_precond:
    (void) printf("user ");
  default:
    if (options[AZ_precond] < AZ_SOLVER_PARAMS)
       (void) printf("iterative preconditioner");
  }
  }

  (void) printf("\n");

  (void) printf(prefix);

  /* lastly, print out the scaling information */

  switch (options[AZ_scaling]) {

  case AZ_none:
    (void) printf("No"); break;
  case AZ_Jacobi:
    (void) printf("block Jacobi"); break;
  case AZ_BJacobi:
    (void) printf("block Jacobi"); break;
  case AZ_row_sum:
    (void) printf("left row-sum"); break;
  case AZ_sym_diag:
    (void) printf("symmetric diagonal"); break;
  case AZ_sym_row_sum:
    (void) printf("symmetric row sum");
  }

  (void) printf(" scaling\n");

  if (super_special==1) {
    (void) printf("%sNOTE: convergence VARIES when the total number "
                  "of\n",prefix);
    (void) printf("%s      processors is changed.\n",prefix);
  }

  if ( recur == 0 ) 
     (void) printf("\t\t****************************************"
                   "***************\n");

} /* AZ_print_call_iter_solver */

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

void AZ_output_matrix(double val[], int indx[], int bindx[], int rpntr[],
                      int cpntr[], int bpntr[], int proc_config[],
                      int data_org[])

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

  Routine to perform full matrix dump on each processor.

  Author:          Scott A. Hutchinson, SNL, 1421
  =======

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

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

  a:               Array containing the nonzero entries of the matrix.

  indx,
  bindx,
  rpntr,
  cpntr,
  bpntr:           Arrays used for DMSR and DVBR sparse matrix storage.

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

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

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

{

  /* local variables */

  int  iblk_row, i, j, ib1, ib2, n1, iblk, jblk, m1, ipoint, jpoint;
  int  ival = 0;
  int  k,num_nonzeros;
  int  num_total_nodes, N_external_nodes;
  int  Proc, Num_Proc;
  char str[5];
  char nstr[40];

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

  Proc               = proc_config[AZ_node];
  Num_Proc           = proc_config[AZ_N_procs];
  N_external_nodes = data_org[AZ_N_external];

  if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) {
    num_total_nodes    = data_org[AZ_N_int_blk]+data_org[AZ_N_bord_blk];
    N_external_nodes = data_org[AZ_N_ext_blk];

    /* Print out the VBR indexing information for the matrix */

    AZ_print_sync_start(Proc, AZ_TRUE, proc_config);

    (void) printf("\n----- Proc: %d indx -----\n\n", Proc);
    for (iblk = 0; iblk < num_total_nodes; iblk++) {
      for (i = *(bpntr+iblk); i < *(bpntr+iblk+1); i++)
        (void) printf("%d ", *(indx+i));

      if (iblk == num_total_nodes - 1)
        (void) printf("%d\n",*(indx+i));
      else
        (void) printf("\n");
    }

    (void) printf("\n----- Proc: %d bindx -----\n\n", Proc);
    for (iblk = 0; iblk < num_total_nodes; iblk++) {
      for (i = *(bpntr+iblk); i < *(bpntr+iblk+1); i++)
        (void) printf("%d ", *(bindx+i));
      (void) printf("\n");
    }

    (void) printf("\n----- Proc: %d rpntr -----\n\n", Proc);
    for (i = 0; i < num_total_nodes + 1; i++)
      (void) printf("%d ", *(rpntr+i));
    (void) printf("\n");

    (void) printf("\n----- Proc: %d cpntr -----\n\n", Proc);
    for (i = 0; i < num_total_nodes + N_external_nodes + 1; i++)
      (void) printf("%d ", *(cpntr+i));
    (void) printf("\n");

    (void) printf("\n----- Proc: %d bpntr -----\n\n", Proc);
    for (i = 0; i < num_total_nodes + 1; i++)
      (void) printf("%d ", *(bpntr+i));
    (void) printf("\n");

    AZ_print_sync_end(proc_config, AZ_TRUE);

    AZ_print_sync_start(Proc, AZ_TRUE, proc_config);
    (void) printf("AZ_solve debug output - full matrix dump: Processor %d\n",
                  Proc);

    /* loop over block rows */

    for (iblk_row = 0; iblk_row < num_total_nodes; iblk_row++) {

      /* number of rows in the current row block */

      m1 = rpntr[iblk_row+1] - rpntr[iblk_row];

      /* starting index of current row block */

      ival = indx[bpntr[iblk_row]];

      /* loop over all the blocks in the current block-row */

      for (j = bpntr[iblk_row]; j < bpntr[iblk_row+1]; j++) {
        jblk = bindx[j];

        /* the starting point column index of the current block */

        ib1 = cpntr[jblk];

        /* ending point column index of the current block */

        ib2 = cpntr[jblk+1];

        /* number of columns in the current block */

        n1 = ib2 - ib1;

        (void) printf("\nProc: %d Block Row: %d Block Column: %d "
                      "Row Pointer: %d Column Pointer: %d\n", Proc, iblk_row,
                      jblk, rpntr[iblk_row], rpntr[jblk]);
        (void) printf("----------------------------------------"
                      "----------------------------------------\n");

                      for (ipoint = 0; ipoint < m1; ipoint++) {
                        for (jpoint = 0; jpoint < n1; jpoint++)
                          (void) printf("a[%d]: %e ", ival+jpoint*m1+ipoint,
                                        val[ival+jpoint*m1+ipoint]);
                        (void) printf("\n");
                      }

                      ival += m1*n1;
      }
    }

    AZ_print_sync_end(proc_config, AZ_TRUE);
  }

  if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {
    num_total_nodes = data_org[AZ_N_internal]+data_org[AZ_N_border];

    N_external_nodes = data_org[AZ_N_external];
    i                  = num_total_nodes + N_external_nodes;

    if      (i < 10    ) sprintf(str, "%%1d");
    else if (i < 100   ) sprintf(str, "%%2d");
    else if (i < 1000  ) sprintf(str, "%%3d");
    else if (i < 10000 ) sprintf(str, "%%4d");
    else if (i < 100000) sprintf(str," %%5d");
    else sprintf(str, "%%d");
    sprintf(nstr, "a(,%s)=%%8.1e ", str);

    AZ_print_sync_start(Proc, AZ_TRUE, proc_config);

    (void) printf("\n----- Proc: %d -----\n\n", Proc);

    num_nonzeros = bindx[num_total_nodes];
    (void) printf("val:  ");
    for (i = 0; i < num_nonzeros; i++) {
      (void) printf("%9.1e", val[i]);
      if ((i%8) == 7) (void) printf("\n    ");
    }

    (void) printf("\nbindx:");
    for (i = 0; i < num_nonzeros; i++) {
      (void) printf("%9d", bindx[i]);
      if ((i%8) == 7) (void) printf("\n    ");
    }
    (void) printf("\n");

    for (i = 0; i < num_total_nodes; i++ ) {
      (void) printf("\nrow");
      (void) printf(str, i);
      (void) printf(":");
      (void) printf(nstr, i, val[i]);
      k = 0;
      for (j = bindx[i]; j < bindx[i+1]; j++ ) {
        (void) printf(nstr, bindx[j], val[j]);
        k++;
        if (((k%4) == 3) && (j != bindx[i+1]-1))
          (void) printf("\n      ");
      }
    }

    (void) printf("\n");
    AZ_print_sync_end( proc_config, AZ_TRUE);
  }

} /* AZ_output_matrix */

void AZ_flop_rates(int data_org[],int indx[],int bpntr[], int bindx[],
	      int options[], double status[], double total_time,
	      int proc_config[])
{
   int N, N_Blk, gn;
   double MFlops, gnnz;
  /*
   * calculate the solver MFlop/s
   */

#ifdef AZ_FLOP_CNTS
#endif

  N_Blk        = data_org[AZ_N_int_blk] + data_org[AZ_N_bord_blk];
  N            = data_org[AZ_N_internal] + data_org[AZ_N_border];
  if ( (options[AZ_output] != AZ_none) && (options[AZ_output] != AZ_warnings) &&
       (data_org[AZ_matrix_type] != AZ_USER_MATRIX)) {
    gn = AZ_gsum_int(N, proc_config);
    if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX)
      gnnz = AZ_gsum_double((double) (indx[bpntr[N_Blk]]-indx[0]),
                            proc_config);
    else if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX)
      gnnz = AZ_gsum_double((double) (bindx[N]), proc_config);
    else return;

    if (proc_config[AZ_node] == 0) {
      total_time += 1.0e-10;
      MFlops = AZ_calc_solve_flops(options, (int) status[AZ_its],total_time,gn,
                                   gnnz, data_org, proc_config);
      if (MFlops > 0.0) {
        (void) printf("\t\tSolver MFlop rate: %f MFlops/sec.\n\t\t", MFlops);
        (void) printf("Solver processor MFlop rate: %f MFlops/sec.\n\n",
                      MFlops/proc_config[AZ_N_procs]);
      }
    }
  }

#ifdef TIME_VB

  /* calculate the individual kernel times and Mflops */

  AZ_time_kernals(gn, gnnz, val, indx, bindx, rpntr,
                  rpntr, bpntr, x, b, (int) status[AZ_NUMBER_OF_ITS],
                  options, data_org, proc_config);

#endif
}
/**************************************************************/
/**************************************************************/
/**************************************************************/
void AZ_mk_identifier(double *params, int *options,
	int *data_org, char *tag) {
/*
 * Make a unique identifier that is linked to preconditioner/
 * solver/scaling/problem size options that the user has set.
 */

  double dtemp;
  int    itemp;
  char   code1,code2,code3;


  dtemp = (params[AZ_ilut_fill] + 3.1415)*(params[AZ_drop] + 2.712)*
	  (1.1 + ((double) options[AZ_graph_fill]));
  if (dtemp > 0) dtemp = pow(dtemp,.01);
  

  itemp = 4*(options[AZ_overlap] + 1) + 2*options[AZ_reorder] +
          options[AZ_type_overlap];

  if (itemp < 94) code1 = '!' + itemp; 
  else code1 = itemp%255;

  if ( options[AZ_precond] == AZ_dom_decomp )
     itemp = options[AZ_subdomain_solve];
  else itemp = options[AZ_precond];

  if (itemp < 94) code2 = '!' + itemp; 
  else code2 = itemp%255;


  itemp = options[AZ_scaling];
  if (itemp < 94) code3 = '!' + itemp; 
  else code3 = itemp%255;

  sprintf(tag,"P%d %c%c%c %.4f", data_org[AZ_N_internal]+data_org[AZ_N_border],
           code1, code2,code3, dtemp);
}

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

void AZ_mk_context(int options[], double params[], int data_org[], 
                   AZ_PRECOND *precond, int proc_config[])
{
/* Make a context to be associated with the preconditioning data structure
 * and stick it in the field: precond->context. This context should be 
 * uniquely based on the name of the matrix which the preconditioner operates 
 * on as well as the type of preconditioning options that have been chosen. 
 *
 * Note: if this context already exists, this routine will simply set
 * precond->context to point to it.
 *
 *
 */
 char tag[80];
 int  istatus;


  AZ_mk_identifier(params,options,data_org, tag);

  precond->context = (struct context *) AZ_manage_memory(sizeof(struct 
					       context), AZ_ALLOC,
                                               data_org[AZ_name],tag,&istatus);
  if (istatus == AZ_NEW_ADDRESS) {
     AZ_zero_out_context(precond->context);