az_solve.c

并行解法器,功能强大

      AZ_abs_matvec_mult(x,&(newparams[AZ_weights]), Amat, proc_config);
      if ((options[AZ_pre_calc] == AZ_reuse)& (options[AZ_scaling] != AZ_none)){
         AZ_scale_f(AZ_INVSCALE_SOL, Amat, options, b, x, proc_config, scaling);
         AZ_scale_f(AZ_INVSCALE_RHS, Amat, options, &(newparams[AZ_weights]), 
                    x, proc_config, scaling);
      }
      largest = 0.;
      for (i = 0; i < NNN; i++) 
         if (newparams[i+AZ_weights] > largest) 
            largest = newparams[AZ_weights+i];
      largest *= 100.;
      for (i = 0; i < NNN; i++) 
         if (newparams[i+AZ_weights] == 0.) newparams[AZ_weights+i]=largest;
      for (i = 0; i < AZ_weights; i++) newparams[i] = params[i];
   }


  /* grab the version string and rip out the ending dots and put it into status */
  AZ_version(tstr);
  t1 = 0; j =0;
  for (i = 0; i < (int) strlen(tstr); i++)  {
     if (tstr[i] == '.') { t1++; if (t1 == 1) tstr[j++] = tstr[i];}
     else tstr[j++] = tstr[i];
  }
  tstr[j] = '\0';
  sscanf(&(tstr[6]),"%lf", &(status[AZ_Aztec_version]));
  if (!AZ_oldsolve_setup(x, b, options, newparams, status, proc_config, 
		    Amat, precond, save_old_values,scaling)) return;

  /* solve the system */

  fflush(stdout);
  switch (options[AZ_solver]) {

  case AZ_cg:

    /* conjugate gradient */

    AZ_pcg_f(b,x,&(newparams[AZ_weights]),options,params,proc_config,status,Amat,precond,conv_info);

    break;

  case AZ_gmres:

    /* GMRES */

    AZ_pgmres(b, x, &(newparams[AZ_weights]), options, params, proc_config,
              status, Amat, precond, conv_info);
 
    break;

  case AZ_fixed_pt:
    AZ_fix_pt(b, x, &(newparams[AZ_weights]), options, params, proc_config,
              status, Amat, precond, conv_info);
 
    break;

  case AZ_analyze:


    for (i=0;i < Amat->data_org[AZ_N_internal]+Amat->data_org[AZ_N_border]; i++)
       b[i] = 0.0;
    AZ_random_vector(x, data_org, proc_config);

    dtemp=AZ_gdot(Amat->data_org[AZ_N_internal]+Amat->data_org[AZ_N_border],
                     x, x, proc_config);
    dtemp = sqrt(dtemp);
    for (j=0;j<Amat->data_org[AZ_N_internal]+Amat->data_org[AZ_N_border];j++)
       x[j] = x[j]/dtemp;

    params[AZ_temp] = 1.;
    params[AZ_tol] = 1.e-40;
    i = options[AZ_max_iter];
    while (i > 0) {
       if (options[AZ_max_iter] > 10) options[AZ_max_iter] = 10;
       AZ_fix_pt(b, x, &(newparams[AZ_weights]), options, params, proc_config,
              status, Amat, precond, conv_info);
       
       dtemp=AZ_gdot(Amat->data_org[AZ_N_internal]+Amat->data_org[AZ_N_border],
                     x, x, proc_config);
       dtemp = sqrt(dtemp);
       for (j=0;j<Amat->data_org[AZ_N_internal]+Amat->data_org[AZ_N_border];j++)
          x[j] = x[j]/dtemp;
       if (options[AZ_extreme] == AZ_high) {
          if (dtemp < 2.0) params[AZ_temp] *= 100.;
       }
       else {
          if (dtemp > 1.0) { params[AZ_temp] /= (1.1*pow(dtemp,.1) ); changed++; }
          else if ( changed == 0 ) params[AZ_temp] *= 2.;
          else if (changed < 2) { params[AZ_temp] *= .7; changed += 3; }
       }
       i -= options[AZ_max_iter];
       options[AZ_max_iter] = i;
    }

       /*  find the largest point */ 
       size1 = Amat->data_org[AZ_N_internal]+Amat->data_org[AZ_N_border];
       largest = -1.; largest_index = -1;
       for (i=0;i < size1; i++) {
          if ( fabs(x[i]) > largest ) { 
            largest = fabs(x[i]); largest_index = i;
          }
       }
       global_largest = AZ_gmax_double(largest, proc_config);


#ifdef ML

       size2 = size1 + Amat->data_org[AZ_N_external];
       tempv = (double *) AZ_allocate(size2*sizeof(double));
       tempy = (double *) AZ_allocate(size2*sizeof(double));
       ibuf =  (int    *) AZ_allocate(size2*sizeof(int   ));
       dbuf =  (double *) AZ_allocate(size2*sizeof(double));
       AZ_exchange_bdry(x, Amat->data_org, proc_config);

       /* Print out the interpolation stencil at that point */

       ml = (ML *) AZ_get_precond_data(precond);
       row_length = 0;
       if (global_largest == largest) {
          allocated = size2/2;
          printf("processor %d (with %d rows) has largest value (%e,%d) in the lambda-vec\n",
                 proc_config[AZ_node], size1, largest,largest_index);
          ML_get_matrix_row(&(ml->Pmat[0]), 1, &largest_index, &allocated, &ibuf, &dbuf,
                            &row_length, 0);
          printf("interpolation operator at that point\n");
          for (i = 0; i < row_length; i++) printf("\t(%d, %e)\n",ibuf[i],dbuf[i]);
          printf("Interpolation columns corresponding to above c-points\n");
       }
       fflush(stdout);

       /* Print out interpolation column for each coarse grid    */
       /* point that appears in the above interpolation stencil. */                 

       row_length = AZ_gsum_int(row_length, proc_config);
       for (i = 0; i < row_length; i++) {
          for (j = 0; j < ml->Pmat[0].invec_leng; j++) tempv[j] = 0.0;
          for (j = 0; j < size1; j++) tempy[j] = 0.0;
          if (global_largest == largest) tempv[ibuf[i]] = 1.;
          else ibuf[i] = 0;
          ibuf[i] = AZ_gsum_int(ibuf[i], proc_config);

          ML_Operator_Apply(&(ml->Pmat[0]),ml->Pmat[0].invec_leng,tempv,size1,tempy);
          for (j = 0; j < size1; j++)
             if (tempy[j] != 0.0) printf("%d: \t (%d,%d %e)\n",proc_config[AZ_node],ibuf[i],j,tempy[j]);
       }
       fflush(stdout);

       if (global_largest == largest) {

          /* Print out the matrix stencil at the largest point */

          already_printed = (int *) tempy;
          for (i = 0; i < size2; i++) already_printed[i] = 0;

          printf("neighbors (A, lvec)\n");
          ML_get_matrix_row(&(ml->Amat[1]), 1, &largest_index, &allocated, &ibuf, &dbuf,
                            &row_length, 0);
          for (i = 0; i < row_length; i++) {
              boundary = ' ';
              if ((ibuf[i] < size1) && (Amat->val[ibuf[i]] == 1.0)) boundary = 'b';
              printf("%d: \t%d (%e , %e) %c\n",proc_config[AZ_node],ibuf[i], dbuf[i],
                     x[ibuf[i]], boundary);
              already_printed[ibuf[i]] = 1;
          }

          /* Print out the distance two values */

          printf("\ndistance 2 neighbors\n");
          ibuf2 = (int *) tempv;
          for (i = 0; i < row_length; i++) {
             ML_get_matrix_row(&(ml->Amat[1]), 1, &(ibuf[i]), &allocated, &ibuf2, &dbuf,
                               &row2_length, 0);
             for (j = 0; j < row2_length; j++) {
                boundary = ' ';
                if ((ibuf2[i] < size1) && (Amat->val[ibuf2[i]] == 1.0)) boundary = 'b';
                if ( already_printed[ ibuf2[j] ] == 0) {
                  printf("%d: \t(%d, %e) %c via %d\n",proc_config[AZ_node],ibuf2[j],
                         x[ibuf2[j]],boundary,i);
                  already_printed[ibuf2[i]] = 1;
                }
             }
          }
       }
       fflush(stdout);
#endif


              
    break;

  case AZ_GMRESR:

    /* GMRESR */

    AZ_pgmresr(b, x, &(newparams[AZ_weights]), options, params, proc_config,
              status, Amat, precond, conv_info);
 
    break;

  case AZ_cgs:

    /* conjugate gradient squared */

    AZ_pcgs(b, x, &(newparams[AZ_weights]), options, params, proc_config, status,
            Amat, precond, conv_info);
 
    break;

  case AZ_tfqmr:

    /* transpose-free quasi minimum residual */

    AZ_pqmrs(b, x, &(newparams[AZ_weights]), options, params, proc_config,status,
            Amat, precond, conv_info);
    break;

  case AZ_bicgstab:

    /* stabilized bi-conjugate gradient */
    AZ_pbicgstab(b, x,&(newparams[AZ_weights]), options, params,proc_config,status,
            Amat, precond, conv_info);
    break;

  case AZ_symmlq:

    /* conjugate gradient squared */

#ifdef eigen
    AZ_psymmlq(b, x, &(newparams[AZ_weights]), options, params, proc_config, status,
               Amat, precond, conv_info);
#else
    printf("symmlq not implemented in this version\n");
#endif
    break;

  case AZ_lu:

    data_org = Amat->data_org;
    N        = data_org[AZ_N_internal] + data_org[AZ_N_border];
    N_Blk    = data_org[AZ_N_int_blk] + data_org[AZ_N_bord_blk];

    /* direct sparse LU */

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

    if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) nonzeros = (Amat->bindx)[N];
    else if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) 
          nonzeros = (Amat->indx)[(Amat->bpntr)[N_Blk]];
    else {
       AZ_matfree_Nnzs(Amat);
       nonzeros = Amat->N_nz;
    }
    

    /* save certain options so that we can restore the user's data */

    itemp2 = options[AZ_precond];
    itemp3 = options[AZ_overlap];
    t1     = options[AZ_subdomain_solve];
    options[AZ_overlap] = AZ_none;
    options[AZ_precond] = AZ_lu;
    options[AZ_subdomain_solve] = AZ_lu;

    /* Set up the subdomain data structure */

    context.A_overlapped           = &Amat2;
    context.aztec_choices          = &aztec_choices;
    context.aztec_choices->options = options;
    context.aztec_choices->params  = newparams;
    context.tag                    = tag;
    sprintf(tag,"lu solve %d",data_org[AZ_name]);

    /* Compute the matrix bandwidth and the number of additional */
    /* nonzeros that might be required by the factorization      */

    i = AZ_compute_max_nz_per_row(Amat, N, N_Blk, &bandwidth);

    AZ_space_for_factors(0.0, nonzeros, N, &extra_factor_nonzeros,options,
                         bandwidth,0);

    /* Adjust the space for the factors based */
    /* on how much memory we can allocate.    */

    N_nz_factors = nonzeros + extra_factor_nonzeros;
    AZ_hold_space(&context,N);
    N_nz_factors = AZ_adjust_N_nz_to_fit_memory(N_nz_factors,2,1);
    AZ_free_space_holder(&context);

    /* Allocate and fill the factor arrays with MSR matrix info */

    bindx2 = (int *) AZ_allocate(N_nz_factors*sizeof(int));
    val2   = (double *) AZ_allocate(N_nz_factors*sizeof(double));
    if (val2 == NULL) AZ_perror("Not enough space for factorization\n");
    if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {
       for (i = 0 ; i < nonzeros ; i++ ) {
          bindx2[i] = (Amat->bindx)[i];
          val2[i] = (Amat->val)[i];
       }
    }
    else if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) {
       AZ_vb2msr(N_Blk,Amat->val,Amat->indx,Amat->bindx,Amat->rpntr,Amat->cpntr,
                 Amat->bpntr, val2, bindx2);
    }
    else AZ_matfree_2_msr(Amat,val2,bindx2,nonzeros);

    Amat2.bindx = bindx2;
    Amat2.val   = val2;
    Amat2.data_org = data_org;

    AZ_factor_subdomain(&context, N, N_nz_factors, &nz_used);
    AZ_solve_subdomain(x, N, &context);
    AZ_free(val2);
    AZ_free(bindx2);
    options[AZ_precond] = itemp2;
    options[AZ_overlap] = itemp3;
    options[AZ_subdomain_solve] = t1;

    status[AZ_its]      = (double ) 1.0;
    status[AZ_why]      = (double ) AZ_normal;
    status[AZ_r]        = (double ) 0.0;
    status[AZ_rec_r]    = (double ) 0.0;
    status[AZ_scaled_r] = (double ) 0.0;
    break;

  default:
    (void) fprintf(stderr,"ERROR: options[AZ_solver] has improper value(%d)\n",
                   options[AZ_solver]);
  exit(-1);
  }
  fflush(stdout);
  AZ_converge_destroy(&conv_info);

  AZ_oldsolve_finish(x, b, options, proc_config, Amat, save_old_values,scaling);
  options[AZ_ignore_scaling] = save_old_values[6];
} /* AZ_oldsolve */

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

void AZ_print_call_iter_solve(int options[], double params[], int az_proc,
			int recur, AZ_PRECOND *precond)

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

  Print out the type of solver called, scaling and preconditioning information.

  Author:          SNL
  =======

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

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

  options:         Determines specific solution method and other parameters.

  params:          Drop tolerance and convergence tolerance info.

  az_proc:         Current processor.

  recur:           Indicates the level of recursion that this call corresponds to.

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

{
  char prefix[40];
  int i, super_special = 0, str_leng = 0;

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

  if ( (options[AZ_output] == AZ_last) || (options[AZ_output] == AZ_none) ||
       (options[AZ_output] == AZ_warnings) ||  (az_proc != 0) ) return;

  prefix[str_leng++] = '\t'; prefix[str_leng++] = '\t'; prefix[str_leng++] = '*';
  prefix[str_leng++] = '*'; prefix[str_leng++] = '*'; prefix[str_leng++] = '*';
  prefix[str_leng++] = '*'; prefix[str_leng++] = ' '; 
  for (i = 0 ; i < recur ; i++ ) {
     prefix[str_leng++] = ' '; 
     prefix[str_leng++] = ' '; 
  }
  prefix[str_leng++] = '\0';

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

  (void) printf(prefix);

  /* first, print out chosen solver */

  switch (options[AZ_solver]) {

  case AZ_cg:
    (void) printf("Preconditioned CG"); break;
  case AZ_gmres:
    (void) printf("Preconditioned GMRES"); break;
  case AZ_analyze:
    (void) printf("Preconditioned analysis"); break;
  case AZ_GMRESR:
    (void) printf("Preconditioned GMRESR"); break;
  case AZ_fixed_pt:
    (void) printf("Preconditioned fixed-point iter."); break;
  case AZ_cgs:
    (void) printf("Preconditioned CGS"); break;
  case AZ_tfqmr:
    (void) printf("Preconditioned TFQMR"); break;
  case AZ_bicgstab:
    (void) printf("Preconditioned BICGSTAB"); break;
  case AZ_symmlq:
    (void) printf("Preconditioned SYMMLQ-like"); break;
  case AZ_lu:
    (void) printf("LU");
  }

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

  /* next output preconditioning options */

  (void) printf(prefix);

  if ((precond != NULL) && (precond->prec_function != AZ_precondition)) {
    if (precond->print_string == NULL) (void) printf("user ");
    else printf("%s",precond->print_string);