13#include "Trilinos_version.h"
14#include "ml_include.h"
25 ML_Aggregate *agg_object;
31static struct ml_info MLInfo[MAX_MI];
42static int FlgDirectSolveCoarsest = 1;
48enum direct_solver {KLU, MUMPS};
49static enum direct_solver DirectSolver = KLU;
56enum smoother_type {Jacobi, SymBlockGaussSeidel, Cheby};
57static enum smoother_type SmootherType = Cheby;
61static int FlgUseHECMWSmoother = 1;
66static int MGType = ML_MGW;
70static int MaxLevels = 10;
77enum coarsen_scheme {UncoupledMIS, METIS, ParMETIS, Zoltan, DD};
78static enum coarsen_scheme CoarsenScheme = UncoupledMIS;
82static int NumSweeps = 2;
91 int N_grids, N_levels;
92 int nlocal, nlocal_allcolumns;
94 ML_Aggregate *agg_object;
96 if (*
id <= 0 && MAX_MI < *
id) {
102 ML_Set_PrintLevel(loglevel);
108 int opt1, opt2, opt3, opt4, opt5, opt6;
114 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML using default coarse solver\n");
116 FlgDirectSolveCoarsest = 0;
117 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML coarse solver is smoother\n");
120 FlgDirectSolveCoarsest = 1;
122 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML coarse solver is KLU\n");
125 FlgDirectSolveCoarsest = 1;
126 DirectSolver = MUMPS;
127 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML coarse solver is MUMPS\n");
130 if (myrank == 0) fprintf(stderr,
"WARNING: invalid solver_opt1=%d (ignored)\n", opt1);
137 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML using default smoother\n");
139 SmootherType = Cheby;
140 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML smoother is Cheby\n");
143 SmootherType = SymBlockGaussSeidel;
144 FlgUseHECMWSmoother = 1;
145 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML smoother is SymBlockGaussSeidel\n");
148 SmootherType = Jacobi;
149 FlgUseHECMWSmoother = 1;
150 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML smoother is Jacobi\n");
153 if (myrank == 0) fprintf(stderr,
"WARNING: invalid solver_opt2=%d (ignored)\n", opt2);
160 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML using default multigrid type\n");
163 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML multigrid type is V-cycle\n");
167 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML multigrid type is W-cycle\n");
171 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML multigrid type is Full-V-cycle\n");
174 if (myrank == 0) fprintf(stderr,
"WARNING: invalid solver_opt3=%d (ignored)\n", opt3);
183 if (myrank == 0) fprintf(stderr,
"WARNING: invalid solver_opt4=%d (ignored)\n", opt4);
185 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML using default MaxLevels\n");
188 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML num of max levels is %d\n", MaxLevels);
194 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML using default coarsening scheme\n");
196 CoarsenScheme = UncoupledMIS;
197 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML coarsening scheme is UncoupledMIS\n");
200 CoarsenScheme = METIS;
201 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML coarsening scheme is METIS\n");
204 CoarsenScheme = ParMETIS;
205 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML coarsening scheme is ParMETIS\n");
208 CoarsenScheme = Zoltan;
209 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML coarsening scheme is Zoltan\n");
213 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML coarsening scheme is DD\n");
216 if (myrank == 0) fprintf(stderr,
"WARNING: invalid solver_opt5=%d (ignored)\n", opt5);
225 if (myrank == 0) fprintf(stderr,
"WARNING: invalid solver_opt6=%d (ignored)\n", opt6);
227 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML using default num of sweeps\n");
232 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML num of smoother sweeps is %d\n", NumSweeps);
238 ML_Create(&ml_object, N_grids);
241 ML_Init_Amatrix(ml_object, 0, nlocal, nlocal,
id);
243 ML_Set_Amatrix_Getrow(ml_object, 0, hecmw_ML_getrow_33, hecmw_ML_comm_33,
245 ML_Set_Amatrix_Matvec(ml_object, 0, hecmw_ML_matvec_33);
247 ML_Set_Amatrix_Getrow(ml_object, 0, hecmw_ML_getrow_nn, hecmw_ML_comm_nn,
249 ML_Set_Amatrix_Matvec(ml_object, 0, hecmw_ML_matvec_nn);
255 ML_Aggregate_Create(&agg_object);
259 int num_PDE_eqns = *Ndof;
265 }
else if (*Ndof == 2) {
270 null_vect = (
double *)
HECMW_malloc(
sizeof(
double) * null_dim * leng);
277 ML_Aggregate_Set_NullSpace(agg_object, num_PDE_eqns, null_dim, null_vect,
286 if (CoarsenScheme == UncoupledMIS) {
287 ML_Aggregate_Set_CoarsenScheme_UncoupledMIS(agg_object);
288 }
else if (CoarsenScheme == METIS) {
289 ML_Aggregate_Set_CoarsenScheme_METIS(agg_object);
290 }
else if (CoarsenScheme == ParMETIS) {
291 ML_Aggregate_Set_CoarsenScheme_ParMETIS(agg_object);
292 }
else if (CoarsenScheme == Zoltan) {
293 ML_Aggregate_Set_CoarsenScheme_Zoltan(agg_object);
294 }
else if (CoarsenScheme == DD) {
295 ML_Aggregate_Set_CoarsenScheme_DD(agg_object);
312 ML_Set_SpectralNormScheme_Calc(ml_object);
324 ML_Aggregate_Set_Dimensions(agg_object, *Ndof);
329 N_levels = ML_Gen_MultiLevelHierarchy_UsingAggregation(
330 ml_object, 0, ML_INCREASING, agg_object);
331 if (loglevel > 0 && myrank == 0) fprintf(stderr,
"INFO: ML generated num of levels is %d\n", N_levels);
342 int coarsest_level = N_levels - 1;
346 if (SmootherType == Jacobi) {
348 if (FlgUseHECMWSmoother) {
351 hecmw_ml_smoother_diag_setup_33_(
id, ierr);
353 ML_Set_Smoother(ml_object, 0, ML_BOTH,
id, hecmw_ML_smoother_diag_apply_33,
"HEC-MW");
355 hecmw_ml_smoother_diag_setup_nn_(
id, ierr);
357 ML_Set_Smoother(ml_object, 0, ML_BOTH,
id, hecmw_ML_smoother_diag_apply_nn,
"HEC-MW");
363 ML_Gen_Smoother_Jacobi(ml_object,
level,
364 ML_BOTH, NumSweeps, ML_DEFAULT);
366 }
else if (SmootherType == SymBlockGaussSeidel) {
368 if (FlgUseHECMWSmoother) {
371 hecmw_ml_smoother_ssor_setup_33_(
id, ierr);
373 ML_Set_Smoother(ml_object, 0, ML_BOTH,
id, hecmw_ML_smoother_ssor_apply_33,
"HEC-MW");
375 hecmw_ml_smoother_ssor_setup_nn_(
id, ierr);
377 ML_Set_Smoother(ml_object, 0, ML_BOTH,
id, hecmw_ML_smoother_ssor_apply_nn,
"HEC-MW");
383 ML_Gen_Smoother_SymBlockGaussSeidel(ml_object,
level,
384 ML_BOTH, NumSweeps, ML_DEFAULT, *Ndof);
388 ML_Gen_Smoother_Cheby(ml_object,
level, ML_BOTH, 20.0, NumSweeps);
394 if (FlgDirectSolveCoarsest) {
395 if (DirectSolver == MUMPS) {
396#if TRILINOS_MAJOR_VERSION < 13
397 ML_Gen_Smoother_Amesos(ml_object, coarsest_level, ML_AMESOS_MUMPS, 1, 0.0);
399 ML_Gen_Smoother_Amesos(ml_object, coarsest_level, ML_AMESOS_MUMPS, 1, 0.0, 1);
402#if TRILINOS_MAJOR_VERSION < 13
403 ML_Gen_Smoother_Amesos(ml_object, coarsest_level, ML_AMESOS_KLU, 1, 0.0);
405 ML_Gen_Smoother_Amesos(ml_object, coarsest_level, ML_AMESOS_KLU, 1, 0.0, 1);
409 if (SmootherType == Jacobi) {
410 ML_Gen_Smoother_Jacobi(ml_object, coarsest_level,
411 ML_BOTH, 3, ML_DEFAULT);
412 }
else if (SmootherType == SymBlockGaussSeidel) {
413 ML_Gen_Smoother_SymBlockGaussSeidel(ml_object, coarsest_level,
414 ML_BOTH, 3, ML_DEFAULT, *Ndof);
416 ML_Gen_Smoother_Cheby(ml_object, coarsest_level, ML_BOTH, 20.0, 2);
422 ML_Gen_Solver(ml_object, MGType, 0, N_levels - 1);
425 MLInfo[*
id - 1].ml_object = ml_object;
426 MLInfo[*
id - 1].agg_object = agg_object;
427 MLInfo[*
id - 1].ndof = *Ndof;
431 int nlocal, nlocal_allcolumns;
435 if (*
id <= 0 && MAX_MI < *
id) {
439 ml_object = MLInfo[*
id - 1].ml_object;
442 sol = (
double *)
HECMW_malloc(
sizeof(
double) * nlocal_allcolumns);
448 ML_Solve_MGV(ml_object, rhs, sol);
449 for (i = 0; i < nlocal; i++) {
456 if (*
id <= 0 && MAX_MI < *
id) {
460 ML_Aggregate_Destroy(&(MLInfo[*
id - 1].agg_object));
461 ML_Destroy(&(MLInfo[*
id - 1].ml_object));
463 if (FlgUseHECMWSmoother) {
464 if (SmootherType == Jacobi) {
465 if (MLInfo[*
id - 1].ndof == 3) {
466 hecmw_ml_smoother_diag_clear_33_(
id, ierr);
468 hecmw_ml_smoother_diag_clear_nn_(
id, ierr);
470 }
else if (SmootherType == SymBlockGaussSeidel) {
471 if (MLInfo[*
id - 1].ndof == 3) {
472 hecmw_ml_smoother_ssor_clear_33_(
id, ierr);
474 hecmw_ml_smoother_ssor_clear_nn_(
id, ierr);
483 fprintf(stderr,
"ERROR: ML not enabled\n");
487 fprintf(stderr,
"ERROR: ML not enabled\n");
491 fprintf(stderr,
"ERROR: ML not enabled\n");
void hecmw_ml_get_loglevel_(int *id, int *level)
void hecmw_ml_get_opt4_(int *id, int *opt3, int *ierr)
void hecmw_ml_get_nlocal_(int *id, int *nlocal, int *nlocal_allcolumns, int *ierr)
void hecmw_ml_get_opt6_(int *id, int *opt3, int *ierr)
void hecmw_ml_get_opt3_(int *id, int *opt3, int *ierr)
void hecmw_ml_set_opt6_(int *id, int *opt3, int *ierr)
void hecmw_ml_get_opt2_(int *id, int *opt2, int *ierr)
void hecmw_ml_get_opt5_(int *id, int *opt3, int *ierr)
void hecmw_ml_get_rbm_(int *id, double rbm[], int *ierr)
void hecmw_ml_get_opt1_(int *id, int *opt1, int *ierr)
void hecmw_ml_wrapper_clear_(int *id, int *ierr)
void hecmw_ML_wrapper_apply(int *id, double rhs[], int *ierr)
void hecmw_ml_wrapper_apply__(int *id, double rhs[], int *ierr)
void hecmw_ml_wrapper_setup__(int *id, int *sym, int *ndof, int *ierr)
void hecmw_ml_wrapper_setup_(int *id, int *sym, int *ndof, int *ierr)
void hecmw_ML_wrapper_setup(int *id, int *sym, int *Ndof, int *ierr)
void hecmw_ML_wrapper_clear(int *id, int *ierr)
void HECMW_ML_WRAPPER_CLEAR(int *id, int *ierr)
void HECMW_ML_WRAPPER_SETUP(int *id, int *sym, int *ndof, int *ierr)
void HECMW_ML_WRAPPER_APPLY(int *id, double rhs[], int *ierr)
void hecmw_ml_wrapper_apply_(int *id, double rhs[], int *ierr)
void hecmw_ml_wrapper_clear__(int *id, int *ierr)
int HECMW_Comm_rank(HECMW_Comm comm, int *rank)
HECMW_Comm HECMW_comm_get_comm(void)
int HECMW_set_error(int errorno, const char *fmt,...)
#define HECMW_malloc(size)
integer(kind=kint) myrank
PARALLEL EXECUTION.