solver.c

cfd求解器使用与gmsh网格的求解

      w  = (double*) Malloc((M->N+1) * sizeof(double));
      jw = (int*) Malloc(2 * (M->N+1) * sizeof(int));    
      ilut_(&M->N, a, ja, ia, &p->Nb_Fill, &p->Dropping_Tolerance,
	    M->S.alu, M->S.jlu, M->S.ju, &nnz_ilu, w, jw, &ierr);    
      Free(w); Free(jw); break;
      
    case ILUTP :
      w  = (double*) Malloc((M->N+1) * sizeof(double));
      jw = (int*) Malloc(2 * (M->N+1) * sizeof(int));    
      ilutp_(&M->N, a, ja, ia, &p->Nb_Fill, &p->Dropping_Tolerance, 
	     &p->Permutation_Tolerance, &M->N, M->S.alu, M->S.jlu, 
	     M->S.ju, &nnz_ilu, w, jw, M->S.permp, &ierr);    
      Free(jw); Free(w); break;
      
    case ILUD :
      w  = (double*) Malloc((M->N+1) * sizeof(double));
      jw = (int*) Malloc(2 * (M->N+1) * sizeof(int));    
      ilud_(&M->N, a, ja, ia, &p->Diagonal_Compensation, 
	    &p->Dropping_Tolerance, M->S.alu, M->S.jlu, 
	    M->S.ju, &nnz_ilu, w, jw, &ierr);    
      Free(w); Free(jw); break;
      
    case ILUDP :
      w     = (double*) Malloc((M->N+1) * sizeof(double));
      jw    = (int*) Malloc(2 * (M->N+1) * sizeof(int));    
      iludp_(&M->N, a, ja, ia, &p->Diagonal_Compensation, 
	     &p->Dropping_Tolerance, &p->Permutation_Tolerance, 
	     &M->N, M->S.alu, M->S.jlu, M->S.ju, &nnz_ilu, 
	     w, jw, M->S.permp, &ierr);    
      Free(jw); Free(w); break;
      
    case ILUK :    
      levels = (int*) Malloc(nnz_ilu * sizeof(int));
      w      = (double*) Malloc((M->N+1) * sizeof(double));
      jw     = (int*) Malloc(3 * (M->N+1) * sizeof(int));
      iluk_(&M->N, a, ja, ia, &p->Nb_Fill, 
	    M->S.alu, M->S.jlu, M->S.ju, 
	    levels, &nnz_ilu, w, jw, &ierr);    
      Free(levels); Free(w); Free(jw); break;
      
    case ILU0 :
      jw = (int*) Malloc((M->N+1) * sizeof(int));    
      ilu0_(&M->N, a, ja, ia, M->S.alu, M->S.jlu, M->S.ju, jw, &ierr);    
      Free(jw); break;
      
    case MILU0 :
      jw = (int*) Malloc((M->N+1) * sizeof(int));    
      milu0_(&M->N, a, ja, ia, M->S.alu, M->S.jlu, M->S.ju, jw, &ierr);    
      Free(jw); break;      
      
    }
    
    switch (ierr){
    case  0 : 
      break;
    case -1 :
      Msg(GERROR, "Input matrix may be wrong");
      break;
    case -2 : /* Matrix L in ILU overflows work array 'al' */
    case -3 : /* Matrix U in ILU overflows work array 'alu' */
      nnz_ilu += nnz_ilu/2 ;
      Msg(INFO, "Reallocating ILU (NZ: %d)", nnz_ilu);
      Free(M->S.alu) ;
      M->S.alu = (scalar*) Malloc(nnz_ilu * sizeof(scalar));
      Free(M->S.jlu) ;
      M->S.jlu = (int*) Malloc(nnz_ilu * sizeof(int));
      goto reallocate ;
    case -4 :
      Msg(GERROR, "Illegal value of nb_fill in ILU");
      break;
    case -5 :
      Msg(GERROR, "Zero row encountered in ILU");
      break;
    default :
      Msg(GERROR, "Zero pivot on line %d in ILU",ierr);
      break;
    }
    
    if(p->Preconditioner != NONE)
      print_matrix_info_MSR(M->N, M->S.alu, M->S.jlu);
    
    if(p->Matrix_Printing == 2 || p->Matrix_Printing == 3){
      Msg(SPARSKIT, "ILU printing\n");
      psplot_(&M->N, M->S.jlu, M->S.jlu, &trente_et_un, &deux);      
    }      
  
    Msg(RESOURCES, "");
  
  }
  
  /* RHS reordering */    
  
  for(i=0;i<M->N;i++){
    rhs[i] = b[M->S.rpermr[i] - 1];
  }
    
  /* Iterations */
  
  ipar[1] = 0;
  ipar[2] = (p->Preconditioner == NONE) ? 0 : p->Preconditioner_Position;
  ipar[3] = 1;
  ipar[4] = 0;
  ipar[5] = p->Krylov_Size;
  ipar[6] = p->Nb_Iter_Max;
  
  fpar[1] = p->Stopping_Test;
  fpar[2] = 0.0;
  fpar[11] = 0.0;
  
  switch (p->Algorithm){      
  case CG      : Msg(SPARSKIT, "Conjugate Gradient (CG)\n"); 
                 ipar[4] = 5 * M->N; break;
  case CGNR    : Msg(SPARSKIT, "CG Normal Residual equation (CGNR)\n");
                 ipar[4] = 5 * M->N; break;
  case BCG     : Msg(SPARSKIT, "Bi-Conjugate Gradient (BCG)\n"); 
                 ipar[4] = 7 * M->N; break;
  case DBCG    : Msg(SPARSKIT, "BCG with partial pivoting (DBCG)\n"); 
                 ipar[4] = 11 * M->N; break;
  case BCGSTAB : Msg(SPARSKIT, "BCG stabilized (BCGSTAB)\n");
                 ipar[4] = 8 * M->N; break;
  case TFQMR   : Msg(SPARSKIT, "Transpose-Free Quasi-Minimum Residual (TFQMR)\n"); 
                 ipar[4] = 11 * M->N; break;
  case FOM     : Msg(SPARSKIT, "Full Orthogonalization Method (FOM)\n"); 
                 ipar[4] = (M->N+3) * (ipar[5]+2) + (ipar[5]+1) * ipar[5]/2; break;
  case GMRES   : Msg(SPARSKIT, "Generalized Minimum RESidual (GMRES)\n");  
                 ipar[4] = (M->N+3) * (ipar[5]+2) + (ipar[5]+1) * ipar[5]/2; break;
  case FGMRES  : Msg(SPARSKIT, "Flexible version of Generalized Minimum RESidual (FGMRES)\n"); 
                 ipar[4] = 2*M->N * (ipar[5]+1) + (ipar[5]+1)*ipar[5]/2 + 3*ipar[5] + 2; break;
  case DQGMRES : Msg(SPARSKIT, "Direct version of Quasi Generalize Minimum RESidual (DQGMRES)\n"); 
                 ipar[4] = M->N + (ipar[5]+1) * (2*M->N+4); break;
  case PGMRES  : Msg(SPARSKIT, "Alternative Generalized Minimum RESidual (GMRES)\n");
                 ipar[4] = (M->N+4) * (ipar[5]+2) + (ipar[5]+1) * ipar[5]/2; break;
  default      : Msg(GERROR, "Unknown algorithm for sparse matrix solver"); break;
  }
  
  w = (double*) Malloc(ipar[4] * sizeof(double));

  its = 0;
  end = 0;
  res = 0.0;


  if( Flag_FMM && p->Scaling != NONE ){
    Msg(SPARSKIT, "FMM Scaling\n") ;    
    FMM_Scaling(M->rowscal, M->colscal);
  }

  if( Flag_FMM && p->Renumbering_Technique == RCMK ){
    Msg(SPARSKIT, "FMM DTAx with Renumbering\n") ; 
    FMM_Renumbering( M->N, M->S.permr, M->S.permp ) ;
  }


  while(1){

    switch(p->Algorithm){
    case CG      : cg_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;      
    case CGNR    : cgnr_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
    case BCG     : bcg_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
    case DBCG    : dbcg_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
    case BCGSTAB : bcgstab_(&M->N, rhs, sol, &ipar[1], &fpar[1], w);	break;
    case TFQMR   : tfqmr_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
    case FOM     : fom_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
    case GMRES   : gmres_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;
    case FGMRES  : fgmres_(&M->N, rhs, sol, &ipar[1], &fpar[1], w); break;      
    case DQGMRES : dqgmres_(&M->N, rhs, sol, &ipar[1], &fpar[1], w);	break;
    case PGMRES  : pgmres_ (&M->N, &p->Krylov_Size, rhs, sol, w, &p->Stopping_Test,
			    &p->Nb_Iter_Max, &six, a, ja, ia, 
			    M->S.alu, M->S.jlu, M->S.ju, &ierr); 
                   end = 1; break;
    }
    
    if(!end){
      
      if(ipar[7] != its){
	if(its) Msg(ITER, " %4d  %.7e  %.7e\n", its, res, res/res1);
	its = ipar[7] ;
      }
      
      res = fpar[5];
      if(its==1) res1 = fpar[5] ;
       
      switch(ipar[1]){
      case 1 : 
	amux_(&M->N, &w[ipar[8]-1], &w[ipar[9]-1], a, ja, ia); break;
      case 2 : 
	atmux_(&M->N, &w[ipar[8]-1], &w[ipar[9]-1], a, ja, ia); break;
      case 3 : case 5 : 
	lusol_(&M->N, &w[ipar[8]-1], &w[ipar[9]-1], M->S.alu, M->S.jlu, M->S.ju); break;
      case 4 : case 6 : 
	lutsol_(&M->N, &w[ipar[8]-1], &w[ipar[9]-1], M->S.alu, M->S.jlu, M->S.ju); break;
      case 0 : 
	end = 1; break;
      case -1 : 
	Msg(WARNING, "Iterative solver has iterated too many times"); end = 1; break;
      case -2 : 
	Msg(WARNING, "Iterative solver was not given enough work space");
	Msg(WARNING, "The work space should at least have %d elements", ipar[4]);
	end = 1; break;
      case -3 : 
	Msg(WARNING, "Iterative solver is facing a break-down"); end = 1; break;
      default : 
	Msg(WARNING, "Iterative solver terminated (code = %d)", ipar[1]); end = 1; break;
      }
      
      if (Flag_FMM && !Flag_DTA){ 
	FMM_MatVectorProd(&w[ipar[8]-1], &w[ipar[9]-1]) ;
      }
      
    }    
    if(end) break;
    
  }

  
  /* Convergence results monitoring */
  
  Msg(ITER, " %4d  %.7e  %.7e\n", ipar[7], fpar[6], fpar[6]/res1);
  
  amux_(&M->N, sol, w, a, ja, ia);
  
 
  for(i=0 ; i<M->N ; i++){
    w[M->N+i] = sol[i] - 1.0 ;
    w[i] -= rhs[i] ;
  }      
  
 
  Msg(SPARSKIT, "%d Iterations / Residual: %g\n", ipar[7], dnrm2_(&M->N,w,&un));
  /*
  Msg(ITER, "Conv. Rate: %g, |Res|: %g, |Err|: %g\n", 
      fpar[7], dnrm2_(&M->N,w,&un), dnrm2_(&M->N,&w[M->N],&un));
  */
  Free(w);

  /* Inverse renumbering */
  
  for (i=0;i<M->N;i++) {
    j = M->S.permr[i] - 1;
    k = M->S.permp[j+M->N] - 1;        
    x[i] = sol[k];
  }

  
  if( Flag_FMM && p->Renumbering_Technique == RCMK ){
    Msg(SPARSKIT, "FMM InverseRenumbering\n") ; 
    FMM_InverseRenumbering(M->S.rpermr) ;
    Flag_FMM = 3;
  } 

  if( Flag_FMM && p->Scaling != NONE ){  
    Msg(SPARSKIT, "FMM UnScaling\n") ; 
    FMM_UnScaling(M->rowscal, M->colscal);
  } 
  

  
  /* Free memory */
  
  Free(rhs);
  Free(sol);
  
  if (!M->ILU_Exists){
    if(p->Preconditioner != NONE) {
      Free(M->S.alu);
      Free(M->S.jlu);
      Free(M->S.ju);
    }
    if (p->Renumbering_Technique == RCMK) {
      Free(M->S.rpermr);
      Free(M->S.permr);
      Free(M->S.permp);
      Free(M->S.a_rcmk);
      Free(M->S.ia_rcmk);
      Free(M->S.ja_rcmk);
    }          
  }

  if(M->T == DENSE){
    List_Delete(M->S.a);
    List_Delete(M->S.jptr);
    List_Delete(M->S.ai);
  }

  if(p->Scaling)
    scale_vector (COLUMN, M, x) ;

}



/* ------------------------------------------------------------------------ */
/*  U t i l s                                                               */
/* ------------------------------------------------------------------------ */

void print_parametres (Solver_Params *p){
  printf(" Matrix_Format           : %d\n", p->Matrix_Format);
  printf(" Matrix_Printing         : %d\n", p->Matrix_Printing);
  printf(" Renumbering_Technique   : %d\n", p->Renumbering_Technique);
  printf(" Preconditioner          : %d\n", p->Preconditioner);
  printf(" Preconditioner_Position : %d\n", p->Preconditioner_Position);
  printf(" Nb_Fill                 : %d\n", p->Nb_Fill);
  printf(" Dropping_Tolerance      : %g\n", p->Dropping_Tolerance);
  printf(" Permutation_Tolerance   : %g\n", p->Permutation_Tolerance);
  printf(" Diagonal_Compensation   : %g\n", p->Diagonal_Compensation);
  printf(" Algorithm               : %d\n", p->Algorithm);
  printf(" Krylov_Size             : %d\n", p->Krylov_Size);
  printf(" IC_Acceleration         : %g\n", p->IC_Acceleration);
  printf(" Iterative_Improvement   : %d\n", p->Iterative_Improvement);
  printf(" Nb_Iter_Max             : %d\n", p->Nb_Iter_Max);
  printf(" Stopping_Test           : %g\n", p->Stopping_Test);
}