diff --git a/src/parcsr_ls/HYPRE_parcsr_amg.c b/src/parcsr_ls/HYPRE_parcsr_amg.c index 5c8dfd456f..f5c9232306 100644 --- a/src/parcsr_ls/HYPRE_parcsr_amg.c +++ b/src/parcsr_ls/HYPRE_parcsr_amg.c @@ -2241,6 +2241,18 @@ HYPRE_BoomerAMGSetIsolatedFPoints(HYPRE_Solver solver, isolated_fpt_index) ); } +/*-------------------------------------------------------------------------- + * HYPRE_BoomerAMGSetUseAuxStrengthMatrix + *--------------------------------------------------------------------------*/ + +HYPRE_Int +HYPRE_BoomerAMGSetUseAuxStrengthMatrix(HYPRE_Solver solver, + HYPRE_Int use_aux_strength_mat) +{ + return (hypre_BoomerAMGSetUseAuxStrengthMatrix( (void *) solver, + use_aux_strength_mat) ); +} + /*-------------------------------------------------------------------------- * HYPRE_BoomerAMGSetCumNnzAP *--------------------------------------------------------------------------*/ @@ -2262,3 +2274,4 @@ HYPRE_BoomerAMGGetCumNnzAP( HYPRE_Solver solver, { return ( hypre_BoomerAMGGetCumNnzAP( (void *) solver, cum_nnz_AP ) ); } + diff --git a/src/parcsr_ls/HYPRE_parcsr_ls.h b/src/parcsr_ls/HYPRE_parcsr_ls.h index 5f0d0620cd..ac14ea1243 100644 --- a/src/parcsr_ls/HYPRE_parcsr_ls.h +++ b/src/parcsr_ls/HYPRE_parcsr_ls.h @@ -1473,6 +1473,14 @@ HYPRE_Int HYPRE_BoomerAMGSetIsolatedFPoints(HYPRE_Solver solver, HYPRE_Int num_isolated_fpt, HYPRE_BigInt *isolated_fpt_index); +/** + * (Optional) if use_aux_strength_mat equals 1, the strength matrix is created from + * an auxilliary M-matrix that is generated from the original system matrix A. + **/ +HYPRE_Int +HYPRE_BoomerAMGSetUseAuxStrengthMatrix(HYPRE_Solver solver, + HYPRE_Int use_aux_strength_mat); + /** * (Optional) if Sabs equals 1, the strength of connection test is based * on the absolute value of the matrix coefficients diff --git a/src/parcsr_ls/Makefile b/src/parcsr_ls/Makefile index 9eaa82ce43..a8c0dc2cba 100644 --- a/src/parcsr_ls/Makefile +++ b/src/parcsr_ls/Makefile @@ -143,6 +143,7 @@ FILES =\ par_schwarz.c\ par_stats.c\ par_strength.c\ + par_aux_strength.c\ par_sv_interp.c\ par_sv_interp_ln.c\ par_vardifconv.c\ diff --git a/src/parcsr_ls/_hypre_parcsr_ls.h b/src/parcsr_ls/_hypre_parcsr_ls.h index 9cef339a33..9c5cff4cba 100644 --- a/src/parcsr_ls/_hypre_parcsr_ls.h +++ b/src/parcsr_ls/_hypre_parcsr_ls.h @@ -310,6 +310,14 @@ typedef struct HYPRE_Solver dslu_solver; #endif + /* use auxilliary strength matrix */ + HYPRE_Int use_aux_strength_matrix; + + /* store complexity info */ + HYPRE_Real grid_complexity; + HYPRE_Real operator_complexity; + HYPRE_Real memory_complexity; + } hypre_ParAMGData; /*-------------------------------------------------------------------------- @@ -581,6 +589,14 @@ typedef struct #define hypre_ParAMGDataDSLUSolver(amg_data) ((amg_data)->dslu_solver) #endif +/* use auxilliary matrix for defining strength */ +#define hypre_ParAMGDataUseAuxStrengthMatrix(amg_data) ((amg_data)->use_aux_strength_matrix) + +/* store complexity info */ +#define hypre_ParAMGDataGridComplexity(amg_data) ((amg_data)->grid_complexity) +#define hypre_ParAMGDataOperatorComplexity(amg_data) ((amg_data)->operator_complexity) +#define hypre_ParAMGDataMemoryComplexity(amg_data) ((amg_data)->memory_complexity) + #endif /****************************************************************************** * Copyright (c) 1998 Lawrence Livermore National Security, LLC and other @@ -2058,6 +2074,8 @@ HYPRE_Int HYPRE_BoomerAMGSetIsolatedFPoints( HYPRE_Solver solver, HYPRE_Int num_ HYPRE_BigInt *isolated_fpt_index ); HYPRE_Int HYPRE_BoomerAMGSetFPoints( HYPRE_Solver solver, HYPRE_Int num_fpt, HYPRE_BigInt *fpt_index ); +HYPRE_Int HYPRE_BoomerAMGSetUseAuxStrengthMatrix(HYPRE_Solver solver, + HYPRE_Int use_aux_strength_mat); HYPRE_Int HYPRE_BoomerAMGSetCumNnzAP ( HYPRE_Solver solver, HYPRE_Real cum_nnz_AP ); HYPRE_Int HYPRE_BoomerAMGGetCumNnzAP ( HYPRE_Solver solver, HYPRE_Real *cum_nnz_AP ); @@ -2684,6 +2702,8 @@ HYPRE_Int hypre_BoomerAMGSetFPoints( void *data, HYPRE_Int isolated, HYPRE_Int n HYPRE_Int hypre_BoomerAMGSetCumNnzAP ( void *data, HYPRE_Real cum_nnz_AP ); HYPRE_Int hypre_BoomerAMGGetCumNnzAP ( void *data, HYPRE_Real *cum_nnz_AP ); +HYPRE_Int hypre_BoomerAMGSetUseAuxStrengthMatrix( void *data, HYPRE_Int use_aux_strength_mat); + /* par_amg_setup.c */ HYPRE_Int hypre_BoomerAMGSetup ( void *amg_vdata, hypre_ParCSRMatrix *A, hypre_ParVector *f, hypre_ParVector *u ); @@ -3098,6 +3118,9 @@ HYPRE_Int hypre_IntersectTwoArrays ( HYPRE_Int *x, HYPRE_Real *x_data, HYPRE_Int HYPRE_Int hypre_IntersectTwoBigArrays ( HYPRE_BigInt *x, HYPRE_Real *x_data, HYPRE_Int x_length, HYPRE_BigInt *y, HYPRE_Int y_length, HYPRE_BigInt *z, HYPRE_Real *output_x_data, HYPRE_Int *intersect_length ); +HYPRE_Int hypre_IntersectTwoBigIntegerArrays ( HYPRE_BigInt *x, HYPRE_Int x_length, + HYPRE_BigInt *y, HYPRE_Int y_length, HYPRE_BigInt *z, + HYPRE_Int *intersect_length ); HYPRE_Int hypre_SortedCopyParCSRData ( hypre_ParCSRMatrix *A, hypre_ParCSRMatrix *B ); HYPRE_Int hypre_BoomerAMG_MyCreateS ( hypre_ParCSRMatrix *A, HYPRE_Real strength_threshold, HYPRE_Real max_row_sum, HYPRE_Int num_functions, HYPRE_Int *dof_func, hypre_ParCSRMatrix **S_ptr ); @@ -3365,7 +3388,8 @@ HYPRE_Int hypre_BoomerAMGCreateSabsHost ( hypre_ParCSRMatrix *A, HYPRE_Real stre HYPRE_Int hypre_BoomerAMGCreate2ndSDevice( hypre_ParCSRMatrix *S, HYPRE_Int *CF_marker, HYPRE_Int num_paths, HYPRE_BigInt *coarse_row_starts, hypre_ParCSRMatrix **C_ptr); HYPRE_Int hypre_BoomerAMGMakeSocFromSDevice( hypre_ParCSRMatrix *A, hypre_ParCSRMatrix *S); - +HYPRE_Int hypre_BoomerAMGCreateAuxS(hypre_ParCSRMatrix *A, hypre_ParCSRMatrix *S, + hypre_ParCSRMatrix **S_aux_ptr, HYPRE_Int method); /* par_sv_interp.c */ HYPRE_Int hypre_BoomerAMGSmoothInterpVectors ( hypre_ParCSRMatrix *A, HYPRE_Int num_smooth_vecs, diff --git a/src/parcsr_ls/par_amg.c b/src/parcsr_ls/par_amg.c index b10876e1d4..60c7680f95 100644 --- a/src/parcsr_ls/par_amg.c +++ b/src/parcsr_ls/par_amg.c @@ -66,6 +66,7 @@ hypre_BoomerAMGCreate( void ) HYPRE_Int rap2; HYPRE_Int keepT; HYPRE_Int modu_rap; + HYPRE_Int use_aux_strength_mat; /* solve params */ HYPRE_Int min_iter; @@ -305,6 +306,8 @@ hypre_BoomerAMGCreate( void ) keepT = 0; modu_rap = 0; + use_aux_strength_mat = -1; + if (hypre_GetExecPolicy1(memory_location) == HYPRE_EXEC_DEVICE) { keepT = 1; @@ -546,6 +549,8 @@ hypre_BoomerAMGCreate( void ) hypre_ParAMGDataKeepTranspose(amg_data) = keepT; hypre_ParAMGDataModularizedMatMat(amg_data) = modu_rap; + hypre_ParAMGDataUseAuxStrengthMatrix(amg_data) = use_aux_strength_mat; + /* information for preserving indices as coarse grid points */ hypre_ParAMGDataCPointsMarker(amg_data) = NULL; hypre_ParAMGDataCPointsLocalMarker(amg_data) = NULL; @@ -563,6 +568,11 @@ hypre_BoomerAMGCreate( void ) hypre_ParAMGDataDSLUSolver(amg_data) = NULL; #endif + /* init complexity info */ + hypre_ParAMGDataGridComplexity(amg_data) = -1.0; + hypre_ParAMGDataOperatorComplexity(amg_data) = -1.0; + hypre_ParAMGDataMemoryComplexity(amg_data) = -1.0; + HYPRE_ANNOTATE_FUNC_END; return (void *) amg_data; @@ -5212,6 +5222,17 @@ hypre_BoomerAMGSetFPoints(void *data, return hypre_error_flag; } +HYPRE_Int +hypre_BoomerAMGSetUseAuxStrengthMatrix( void *data, + HYPRE_Int use_aux_strength_mat) +{ + hypre_ParAMGData *amg_data = (hypre_ParAMGData*) data; + + hypre_ParAMGDataUseAuxStrengthMatrix(amg_data) = use_aux_strength_mat; + + return hypre_error_flag; +} + HYPRE_Int hypre_BoomerAMGSetCumNnzAP( void *data, HYPRE_Real cum_nnz_AP ) @@ -5228,7 +5249,6 @@ hypre_BoomerAMGSetCumNnzAP( void *data, return hypre_error_flag; } - HYPRE_Int hypre_BoomerAMGGetCumNnzAP( void *data, HYPRE_Real *cum_nnz_AP ) @@ -5244,3 +5264,4 @@ hypre_BoomerAMGGetCumNnzAP( void *data, return hypre_error_flag; } + diff --git a/src/parcsr_ls/par_amg.h b/src/parcsr_ls/par_amg.h index 715cf32cdf..540fce5da5 100644 --- a/src/parcsr_ls/par_amg.h +++ b/src/parcsr_ls/par_amg.h @@ -292,6 +292,14 @@ typedef struct HYPRE_Solver dslu_solver; #endif + /* use auxilliary strength matrix */ + HYPRE_Int use_aux_strength_matrix; + + /* store complexity info */ + HYPRE_Real grid_complexity; + HYPRE_Real operator_complexity; + HYPRE_Real memory_complexity; + } hypre_ParAMGData; /*-------------------------------------------------------------------------- @@ -563,4 +571,12 @@ typedef struct #define hypre_ParAMGDataDSLUSolver(amg_data) ((amg_data)->dslu_solver) #endif +/* use auxilliary matrix for defining strength */ +#define hypre_ParAMGDataUseAuxStrengthMatrix(amg_data) ((amg_data)->use_aux_strength_matrix) + +/* store complexity info */ +#define hypre_ParAMGDataGridComplexity(amg_data) ((amg_data)->grid_complexity) +#define hypre_ParAMGDataOperatorComplexity(amg_data) ((amg_data)->operator_complexity) +#define hypre_ParAMGDataMemoryComplexity(amg_data) ((amg_data)->memory_complexity) + #endif diff --git a/src/parcsr_ls/par_amg_setup.c b/src/parcsr_ls/par_amg_setup.c index 544bdd66d9..ddf78a8dd9 100644 --- a/src/parcsr_ls/par_amg_setup.c +++ b/src/parcsr_ls/par_amg_setup.c @@ -50,7 +50,7 @@ hypre_BoomerAMGSetup( void *amg_vdata, hypre_IntArray **CF_marker_array; hypre_IntArray **dof_func_array; hypre_IntArray *dof_func; - HYPRE_Int *dof_func_data; + HYPRE_Int *dof_func_data = NULL; HYPRE_Real *relax_weight; HYPRE_Real *omega; HYPRE_Real schwarz_relax_wt = 1; @@ -231,6 +231,10 @@ hypre_BoomerAMGSetup( void *amg_vdata, HYPRE_Real wall_time = 0.0; /* for debugging instrumentation */ HYPRE_Int add_end; + HYPRE_Int use_aux_strength_mat = hypre_ParAMGDataUseAuxStrengthMatrix(amg_data); + hypre_ParCSRMatrix *S_aux = NULL; + hypre_ParCSRMatrix *S_H = NULL; + #ifdef HYPRE_USING_DSUPERLU HYPRE_Int dslu_threshold = hypre_ParAMGDataDSLUThreshold(amg_data); #endif @@ -983,12 +987,68 @@ hypre_BoomerAMGSetup( void *amg_vdata, hypre_ParAMGDataSmoother(amg_data) = smoother; } + /*------------------------------------------------------------------------------------------------ + * Initialize Auxilliary matrix for computing strength matrix if needed + *-----------------------------------------------------------------------------------------------*/ + if (use_aux_strength_mat >= 0) + { + HYPRE_Int method = 0; + switch (use_aux_strength_mat) + { + case 1: + method = 1; + break; + case 2: + method = 2; + break; + case 10: + method = 10; + break; + case 11: + method = 11; + break; + case 12: + method = 12; + break; + default: + method = 0; + } + if (my_id == 0 && amg_print_level > 1) + { + hypre_printf("\n\nUsing Auxilliary matrix for computing strength matrix: method %d \n", method); + } + use_aux_strength_mat = method; + } + if (use_aux_strength_mat == 10 || use_aux_strength_mat == 12) + { + /* Build auxilliary matrix for strength */ + hypre_BoomerAMGCreateAuxS(A, NULL, &S_aux, use_aux_strength_mat); + } + else if (use_aux_strength_mat == 11) + { + /* create iniitial strenth matrix */ + hypre_BoomerAMGCreateS(A, strong_threshold, max_row_sum, + num_functions, dof_func_data, &S); + /* Build auxilliary matrix from strength matrix */ + hypre_BoomerAMGCreateAuxS(A, S, &S_aux, 1); + /* create new strength matrix from auxilliary matrix */ + hypre_ParCSRMatrixDestroy(S); + S = NULL; + } + /*----------------------------------------------------- * Enter Coarsening Loop *-----------------------------------------------------*/ - while (not_finished_coarsening) { + /* WM: debug - hack to get rid of zeros in A, which can screw with the sparseity pattern of S_aux */ + hypre_CSRMatrix *delete_zeros = hypre_CSRMatrixDeleteZeros(hypre_ParCSRMatrixDiag(A_array[level]), + 1.0e-12); + if (delete_zeros) + { + hypre_CSRMatrixDestroy(hypre_ParCSRMatrixDiag(A_array[level])); + hypre_ParCSRMatrixDiag(A_array[level]) = delete_zeros; + } /* only do nodal coarsening on a fixed number of levels */ if (level >= nodal_levels) { @@ -1128,8 +1188,70 @@ hypre_BoomerAMGSetup( void *amg_vdata, { if (!useSabs) { - hypre_BoomerAMGCreateS(A_array[level], strong_threshold, max_row_sum, - num_functions, dof_func_data, &S); + if (use_aux_strength_mat == 0 || use_aux_strength_mat == 2) + { + /* Build auxilliary matrix for strength */ + hypre_BoomerAMGCreateAuxS(A_array[level], NULL, &S_aux, use_aux_strength_mat); + /* build strength matrix */ + hypre_BoomerAMGCreateS(S_aux, strong_threshold, max_row_sum, + num_functions, dof_func_data, &S); + /* WM: debug - print out S and S_aux */ + if (amg_print_level == 3) + { + char file[256]; + hypre_sprintf(file, "matrices/S_aux_str%.2f_auxs%d_%d", strong_threshold, use_aux_strength_mat, + level); + hypre_ParCSRMatrixPrintIJ(S_aux, 0, 0, file); + hypre_sprintf(file, "matrices/S_str%.2f_auxs%d_%d", strong_threshold, use_aux_strength_mat, level); + hypre_ParCSRMatrixPrintIJ(S, 0, 0, file); + } + /* destroy auxilliary matrix */ + hypre_ParCSRMatrixDestroy(S_aux); + S_aux = NULL; + } + else if (use_aux_strength_mat == 1) + { + /* create iniitial strenth matrix */ + hypre_BoomerAMGCreateS(A_array[level], strong_threshold, max_row_sum, + num_functions, dof_func_data, &S); + /* Build auxilliary matrix from strength matrix */ + hypre_BoomerAMGCreateAuxS(A_array[level], S, &S_aux, 1); + /* create new strength matrix from auxilliary matrix */ + hypre_ParCSRMatrixDestroy(S); + S = NULL; + hypre_BoomerAMGCreateS(S_aux, strong_threshold, max_row_sum, + num_functions, dof_func_data, &S); + /* destroy auxilliary matrix */ + hypre_ParCSRMatrixDestroy(S_aux); + S_aux = NULL; + } + else if (use_aux_strength_mat == 10 || use_aux_strength_mat == 11 || use_aux_strength_mat == 12) + { + hypre_BoomerAMGCreateS(S_aux, strong_threshold, max_row_sum, + num_functions, dof_func_data, &S); + /* WM: debug - print out S and S_aux */ + if (amg_print_level == 3) + { + char file[256]; + hypre_sprintf(file, "matrices/S_aux_str%.2f_auxs%d_%d", strong_threshold, use_aux_strength_mat, + level); + hypre_ParCSRMatrixPrintIJ(S_aux, 0, 0, file); + hypre_sprintf(file, "matrices/S_str%.2f_auxs%d_%d", strong_threshold, use_aux_strength_mat, level); + hypre_ParCSRMatrixPrintIJ(S, 0, 0, file); + } + } + else + { + hypre_BoomerAMGCreateS(A_array[level], strong_threshold, max_row_sum, + num_functions, dof_func_data, &S); + /* WM: debug - print out S and S_aux */ + if (amg_print_level == 3) + { + char file[256]; + hypre_sprintf(file, "matrices/S_str%.2f_auxs%d_%d", strong_threshold, use_aux_strength_mat, level); + hypre_ParCSRMatrixPrintIJ(S, 0, 0, file); + } + } } else { @@ -2855,7 +2977,7 @@ hypre_BoomerAMGSetup( void *amg_vdata, { hypre_MatvecCommPkgCreate(A_H); } - + hypre_printf("Building RAP here ... nongal \n"); /* Delete AP */ hypre_ParCSRMatrixDestroy(Q); } @@ -3078,6 +3200,17 @@ hypre_BoomerAMGSetup( void *amg_vdata, { hypre_BoomerAMGBuildCoarseOperatorKT(P_array[level], A_array[level], P_array[level], keepTranspose, &A_H); + + if (use_aux_strength_mat == 10 || use_aux_strength_mat == 11 || use_aux_strength_mat == 12) + { + hypre_BoomerAMGBuildCoarseOperatorKT(P_array[level], S_aux, + P_array[level], keepTranspose, &S_H); + /* Do we need to drop small entries in S_H ?? */ + /* reset S_aux */ + hypre_ParCSRMatrixDestroy(S_aux); + S_aux = S_H; + } + } if (Pnew && ns == 1) @@ -3088,6 +3221,22 @@ hypre_BoomerAMGSetup( void *amg_vdata, } } + /* WM: debug - use my own naming and just get P out once... */ + if (amg_print_level == 3) + { + char file[256]; + if (level == 0) + { + hypre_sprintf(file, "matrices/A_str%.2f_auxs%d_%d", strong_threshold, use_aux_strength_mat, 0); + hypre_ParCSRMatrixPrintIJ(A_array[0], 0, 0, file); + } + hypre_sprintf(file, "matrices/A_str%.2f_auxs%d_%d", strong_threshold, use_aux_strength_mat, + level + 1); + hypre_ParCSRMatrixPrintIJ(A_H, 0, 0, file); + + hypre_sprintf(file, "matrices/P_str%.2f_auxs%d_%d", strong_threshold, use_aux_strength_mat, level); + hypre_ParCSRMatrixPrintIJ(P_array[level], 0, 0, file); + } #if DEBUG_SAVE_ALL_OPS if (level == 0) { @@ -3197,6 +3346,14 @@ hypre_BoomerAMGSetup( void *amg_vdata, } HYPRE_ANNOTATE_REGION_END("%s", "Coarse solve"); HYPRE_ANNOTATE_MGLEVEL_END(level); + + /* destroy S_aux */ + if (use_aux_strength_mat == 0 || use_aux_strength_mat == 11) + { + hypre_ParCSRMatrixDestroy(S_aux); + S_aux = NULL; + } + hypre_GpuProfilingPopRange(); if (level > 0) diff --git a/src/parcsr_ls/par_aux_strength.c b/src/parcsr_ls/par_aux_strength.c new file mode 100644 index 0000000000..855c9f0198 --- /dev/null +++ b/src/parcsr_ls/par_aux_strength.c @@ -0,0 +1,954 @@ +/****************************************************************************** + * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other + * HYPRE Project Developers. See the top-level COPYRIGHT file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + ******************************************************************************/ + +/****************************************************************************** + * + *****************************************************************************/ + +#include "_hypre_parcsr_ls.h" + +/*==========================================================================*/ +/*==========================================================================*/ +/** + Generates an auxiilliary matrix, S_aux, with M-matrix properties such that + S_aux = A - B, where B is obtained by distributing positive off-diagonal + contributions to 'strong' neighbors ("positive" here means having the same + sign as the diagonal). + + NOTE: Current implementation assumes A is a square matrix with a regular + partitioning of rows and columns. + + {\bf Input files:} + _hypre_parcsr_ls.h + + @return Error code. + + @param A [IN] + coefficient matrix + @param S_ptr [OUT] + strength matrix + + @see */ +/*--------------------------------------------------------------------------*/ + +HYPRE_Int +hypre_BoomerAMGCreateAuxMMatrix(hypre_ParCSRMatrix *A, + hypre_ParCSRMatrix **S_aux_ptr, + HYPRE_Int method) +{ +#ifdef HYPRE_PROFILE + hypre_profile_times[HYPRE_TIMER_ID_CREATES] -= hypre_MPI_Wtime(); +#endif + + MPI_Comm comm = hypre_ParCSRMatrixComm(A); + HYPRE_MemoryLocation memory_location_P = hypre_ParCSRMatrixMemoryLocation(A); + hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A); + hypre_ParCSRCommHandle *comm_handle; + + /* diag part of A */ + hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A); + HYPRE_Complex *A_diag_data = hypre_CSRMatrixData(A_diag); + HYPRE_Int *A_diag_i = hypre_CSRMatrixI(A_diag); + HYPRE_Int *A_diag_j = hypre_CSRMatrixJ(A_diag); + /* off-diag part of A */ + hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A); + HYPRE_Complex *A_offd_data = hypre_CSRMatrixData(A_offd); + HYPRE_Int *A_offd_i = hypre_CSRMatrixI(A_offd); + HYPRE_Int *A_offd_j = hypre_CSRMatrixJ(A_offd); + + /* diag part of S */ + hypre_CSRMatrix *S_diag = NULL; + HYPRE_Int *S_diag_i = NULL; + HYPRE_Int *S_diag_j = NULL; + // HYPRE_Int skip_diag = S ? 0 : 1; + /* off-diag part of S */ + hypre_CSRMatrix *S_offd = NULL; + HYPRE_Int *S_offd_i = NULL; + HYPRE_Int *S_offd_j = NULL; + + HYPRE_BigInt *row_starts = hypre_ParCSRMatrixRowStarts(A); + HYPRE_Int num_variables = hypre_CSRMatrixNumRows(A_diag); + HYPRE_BigInt global_num_vars = hypre_ParCSRMatrixGlobalNumRows(A); + HYPRE_Int num_cols_A_diag = hypre_CSRMatrixNumCols(A_diag); + HYPRE_Int num_cols_A_offd = hypre_CSRMatrixNumCols(A_offd); + + HYPRE_BigInt *A_col_starts = hypre_ParCSRMatrixColStarts(A); + HYPRE_BigInt *A_offd_colmap = hypre_ParCSRMatrixColMapOffd(A); + + HYPRE_Int A_diag_nnz = A_diag_i[num_variables]; + HYPRE_Int A_offd_nnz = A_offd_i[num_variables]; + + HYPRE_Int num_nonzeros_diag; + HYPRE_Int num_nonzeros_offd = 0; + HYPRE_Int num_cols_offd = 0; + + hypre_ParCSRMatrix *S_aux; + hypre_CSRMatrix *S_aux_diag; + HYPRE_Int *S_aux_diag_i = NULL; + HYPRE_Int *S_aux_diag_j = NULL; + HYPRE_Complex *S_aux_diag_data = NULL; + /* HYPRE_Real *S_diag_data; */ + hypre_CSRMatrix *S_aux_offd; + HYPRE_Int *S_aux_offd_i = NULL; + HYPRE_Int *S_aux_offd_j = NULL; + HYPRE_Complex *S_aux_offd_data = NULL; + HYPRE_BigInt *S_aux_offd_colmap = NULL; + /* HYPRE_Real *S_offd_data; */ + /* off processor portions of S */ + hypre_CSRMatrix *A_ext = NULL; + HYPRE_Int *A_ext_i = NULL; + HYPRE_Real *A_ext_data = NULL; + HYPRE_BigInt *A_ext_j = NULL; + HYPRE_Int A_ext_nnz = 0; + + HYPRE_Real diag, row_scale, row_sum; + HYPRE_Int i, j, k, jj, jA, jS, kp, kn, m, n, P_n_max, N_n_max, cnt_n, cnt_p; + + HYPRE_Int ierr = 0; + + HYPRE_Int *dof_func_offd; + HYPRE_Int num_sends; + HYPRE_Int *int_buf_data; + HYPRE_Int index, start; + + HYPRE_Int *P_i = NULL; + HYPRE_Int *P_n = NULL; + HYPRE_BigInt *N_i = NULL; + HYPRE_Int *N_n = NULL; + HYPRE_BigInt *B_i = NULL; + HYPRE_Int *B_n = NULL; + HYPRE_Complex *B_a = NULL; + HYPRE_Complex *B_a_sum = NULL; + + HYPRE_Complex pval, bval; + + HYPRE_Int *prefix_sum_workspace; + + HYPRE_MemoryLocation memory_location = hypre_ParCSRMatrixMemoryLocation(A); + + HYPRE_Int num_procs, my_id; + hypre_MPI_Comm_size(comm, &num_procs); + hypre_MPI_Comm_rank(comm, &my_id); + + // Begin: + /* Allocate memory for auxilliary strength matrix */ + S_aux = hypre_ParCSRMatrixCreate(comm, global_num_vars, global_num_vars, + row_starts, row_starts, + num_cols_A_offd, A_diag_nnz, A_offd_nnz); + + S_aux_diag = hypre_ParCSRMatrixDiag(S_aux); + hypre_CSRMatrixI(S_aux_diag) = hypre_CTAlloc(HYPRE_Int, num_variables + 1, memory_location); + hypre_CSRMatrixJ(S_aux_diag) = hypre_CTAlloc(HYPRE_Int, A_diag_nnz, memory_location); + hypre_CSRMatrixData(S_aux_diag) = hypre_CTAlloc(HYPRE_Complex, A_diag_nnz, memory_location); + + S_aux_diag_i = hypre_CSRMatrixI(S_aux_diag); + S_aux_diag_j = hypre_CSRMatrixJ(S_aux_diag); + S_aux_diag_data = hypre_CSRMatrixData(S_aux_diag); + + S_aux_offd = hypre_ParCSRMatrixOffd(S_aux); + hypre_CSRMatrixI(S_aux_offd) = hypre_CTAlloc(HYPRE_Int, num_variables + 1, memory_location); + + S_aux_offd_i = hypre_CSRMatrixI(S_aux_offd); + + + if (num_cols_A_offd) + { + hypre_CSRMatrixJ(S_aux_offd) = hypre_CTAlloc(HYPRE_Int, A_offd_nnz, memory_location); + hypre_CSRMatrixData(S_aux_offd) = hypre_CTAlloc(HYPRE_Complex, A_offd_nnz, memory_location); + + S_aux_offd_j = hypre_CSRMatrixJ(S_aux_offd); + S_aux_offd_data = hypre_CSRMatrixData(S_aux_offd); + + S_aux_offd_colmap = hypre_TAlloc(HYPRE_BigInt, num_cols_A_offd, memory_location); + hypre_ParCSRMatrixColMapOffd(S_aux) = S_aux_offd_colmap; + +#ifdef HYPRE_USING_OPENMP + #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE +#endif + for (i = 0; i < num_cols_A_offd; i++) + { + S_aux_offd_colmap[i] = A_offd_colmap[i]; + } + } + + // extract off-processor portion of A + if (num_procs > 1) + { + A_ext = hypre_ParCSRMatrixExtractBExt(A, A, 1); + A_ext_data = hypre_CSRMatrixData(A_ext); + A_ext_i = hypre_CSRMatrixI(A_ext); + A_ext_j = hypre_CSRMatrixBigJ(A_ext); + A_ext_nnz = A_ext_i[num_cols_A_offd]; + + } + + /* allocate work arrays */ + HYPRE_Int P_nnz = A_diag_nnz + A_offd_nnz; + HYPRE_Int N_nnz = A_diag_nnz + A_offd_nnz + A_ext_nnz; + P_i = hypre_CTAlloc(HYPRE_Int, P_nnz, memory_location); + P_n = hypre_CTAlloc(HYPRE_Int, (num_variables + 1), memory_location); + N_i = hypre_CTAlloc(HYPRE_BigInt, N_nnz, memory_location); + N_n = hypre_CTAlloc(HYPRE_Int, (num_variables + num_cols_A_offd + 1), memory_location); + + /* 1. Loop over matrix rows to extract positive and negative neighbors */ + P_n_max = 0; + N_n_max = 0; + kp = 0; + kn = 0; + for (i = 0; i < num_variables; i++) + { + cnt_p = kp; + cnt_n = kn; + + diag = A_diag_data[A_diag_i[i]]; + + // diagonal part + // Check for positive off-diagonal entry + for (j = A_diag_i[i]; j < A_diag_i[i + 1]; j++) + { + if ( (diag * A_diag_data[j] > 0) && (A_diag_j[j] != i)) + { + // position of positive neighbor. Store position for quick and easy access + P_i[kp++] = j; //A_diag_j[j]; + } + else if (A_diag_j[j] != i) + { + N_i[kn++] = A_diag_j[j] + A_col_starts[0]; + } + } + // Off-diagonal part + for (j = A_offd_i[i]; j < A_offd_i[i + 1]; j++) + { + if ( diag * A_offd_data[j] > 0 ) + { + // position of positive neighbor. Offset to distinguish from diag part columns + P_i[kp++] = j + A_diag_nnz; + } + else + { + N_i[kn++] = A_offd_colmap[A_offd_j[j]]; + } + } + cnt_p = kp - cnt_p; + cnt_n = kn - cnt_n; + P_n_max = P_n_max > cnt_p ? P_n_max : cnt_p; + N_n_max = N_n_max > cnt_n ? N_n_max : cnt_n; + // sort neighbor arrays + hypre_BigQsort0(N_i, N_n[i], (kn - 1)); + // update index pointer arrays + P_n[i + 1] = kp; + N_n[i + 1] = kn; + } + // add negative neighbors list of external rows + kn = N_n[num_variables]; + for (i = 0; i < num_cols_A_offd; i++) + { + diag = A_diag_data[A_diag_i[i]]; + + for (j = A_ext_i[i]; j < A_ext_i[i + 1]; j++) + { + if ( diag * A_ext_data[j] < 0 ) + { + // get negative neighbor list + N_i[kn++] = A_ext_j[j]; + } + } + jj = num_variables + i; + // sort neighbor arrays + hypre_BigQsort0(N_i, N_n[jj], (kn - 1)); + // update index pointer arrays + N_n[jj + 1] = kn; + } + hypre_CSRMatrixDestroy(A_ext); + /* 2. Loop to compute intersections of strong (negative) neighbors of common weak (positive) entry */ + B_i = hypre_CTAlloc(HYPRE_BigInt, (N_n_max * P_n[num_variables]), memory_location); + B_n = hypre_CTAlloc(HYPRE_Int, P_n[num_variables] + 1, memory_location); + if (method % 10 == 2) + { + B_a = hypre_CTAlloc(HYPRE_Complex, (N_n_max * P_n[num_variables]), memory_location); + B_a_sum = hypre_CTAlloc(HYPRE_Complex, P_n[num_variables], memory_location); + } + kn = 0; + kp = 0; + for (i = 0; i < num_variables; i++) + { + // loop over positive neighbor list + for (j = P_n[i]; j < P_n[i + 1]; j++) + { + jA = P_i[j] - A_diag_nnz; + // Local row (assumes square matrix) + if ( jA < 0 ) + { + k = A_diag_j[P_i[j]] ; + } + else + { + // external row. Access external variable range of N. + k = A_offd_j[jA] + num_variables; + } + // get intersection of N_i [i] and N_i[k] + m = N_n[i + 1] - N_n[i]; + n = N_n[k + 1] - N_n[k]; + hypre_IntersectTwoBigIntegerArrays(&N_i[N_n[i]], m, &N_i[N_n[k]], n, &B_i[kn], &kp); + // update position of intersection of strong neighbors + kn += kp; + B_n[j + 1] = kn; + + // WM: get sum of entries in A for the intersection + // WM: todo - switching around variables to use as my indices, which is kind of ugly... + // want the connection from positive connection to each column marked in B_i + if (method % 10 == 2) + { + m = k; + // if positive connection, m, is local + if ( m < num_variables ) + { + for (k = B_n[j]; k < B_n[j + 1]; k++) + { + HYPRE_BigInt big_col = B_i[k]; + if ( big_col >= A_col_starts[0] && big_col < A_col_starts[1]) + { + // neighbor of positive entry is in diag part + jA = A_diag_i[m]; + HYPRE_Int col = (HYPRE_Int) (big_col - A_col_starts[0]); + while (A_diag_j[jA] != col) { jA++; } + B_a[k] = A_diag_data[jA]; + B_a_sum[j] += A_diag_data[jA]; + } + else + { + // neighbor of positive entry is in offd part + jA = A_offd_i[m]; + HYPRE_Int col = hypre_BigBinarySearch( A_offd_colmap, big_col, num_cols_A_offd); + while (A_offd_j[jA] != col) { jA++; } + B_a[k] = A_offd_data[jA]; + B_a_sum[j] += A_offd_data[jA]; + } + } + } + else + { + m -= num_variables; + /* m is now a row index for A_ext? */ + /* it's a col index for A_offd, so these correspond to rows of A_ext, right? */ + for (k = B_n[j]; k < B_n[j + 1]; k++) + { + HYPRE_BigInt big_col = B_i[k]; + jA = A_ext_i[m]; + while (A_ext_j[jA] != big_col) { jA++; } + B_a[k] = A_ext_data[jA]; + B_a_sum[j] += A_ext_data[jA]; + } + } + } + } + } + hypre_TFree(N_i, memory_location); + hypre_TFree(N_n, memory_location); + /* 3. Fill in data for S_aux */ + kp = 0; + kn = 0; + for (i = 0; i < num_variables; i++) + { + diag = A_diag_data[A_diag_i[i]]; + + // diagonal part + for (j = A_diag_i[i]; j < A_diag_i[i + 1]; j++) + { + // Skip positive off-diagonal entry + if ( (diag * A_diag_data[j] > 0) && (A_diag_j[j] != i)) { continue; } + + S_aux_diag_data[kp] = A_diag_data[j]; + S_aux_diag_j[kp++] = A_diag_j[j]; + } + // Off-diagonal part + for (j = A_offd_i[i]; j < A_offd_i[i + 1]; j++) + { + // Skip positive off-diagonal entry + if ( diag * A_offd_data[j] > 0 ) { continue; } + + S_aux_offd_data[kn] = A_offd_data[j]; + S_aux_offd_j[kn++] = A_offd_j[j]; + } + + // update index pointer arrays + S_aux_diag_i[i + 1] = kp; + S_aux_offd_i[i + 1] = kn; + } + /* Distribute positive contributions to strong neighbors */ + for (i = 0; i < num_variables; i++) + { + // loop over positive neighbors and compute contributions to (strong) negative neighbors + for (j = P_n[i]; j < P_n[i + 1]; j++) + { + jA = P_i[j] - A_diag_nnz; + // positive entry is in diag part + if ( jA < 0 ) + { + pval = A_diag_data[P_i[j]] ; + } + else + { + // positive entry is in offd part + pval = A_offd_data[jA]; + } + + // compute distribution to negative neighbors + bval = -2.0 * pval / (HYPRE_Complex) (B_n[ j + 1 ] - B_n[ j ]); + + // Loop over negative connections to positive neighbor and distribute + for (k = B_n[ j ]; k < B_n[ j + 1 ]; k++) + { + HYPRE_BigInt big_col = B_i[k]; + + // WM: compute bval weigthed by connections in A + if (method % 10 == 2) + { + bval = -2.0 * pval * B_a[k] / B_a_sum[j]; + } + + if ( big_col >= A_col_starts[0] && big_col < A_col_starts[1]) + { + // neighbor is in diag part + jS = S_aux_diag_i[i] + 1; + jj = (HYPRE_Int) (big_col - A_col_starts[0]); + while (S_aux_diag_j[jS] != jj) { jS++; } + S_aux_diag_data[jS] -= bval; + } + else + { + // neighbor is in offd part + jj = hypre_BigBinarySearch( A_offd_colmap, big_col, num_cols_A_offd); + jS = S_aux_offd_i[i]; + while (S_aux_offd_j[jS] != jj) { jS++; } + S_aux_offd_data[jS] -= bval; + } + } + // update diagonal entry + S_aux_diag_data[S_aux_diag_i[i]] -= pval; + } + } + + hypre_CSRMatrixNumNonzeros(S_aux_diag) = S_aux_diag_i[num_variables]; + hypre_CSRMatrixNumNonzeros(S_aux_offd) = S_aux_offd_i[num_variables]; + hypre_CSRMatrixJ(S_aux_diag) = S_aux_diag_j; + hypre_CSRMatrixJ(S_aux_offd) = S_aux_offd_j; + + hypre_CSRMatrixMemoryLocation(S_aux_diag) = memory_location; + hypre_CSRMatrixMemoryLocation(S_aux_offd) = memory_location; + + // hypre_ParCSRMatrixCommPkg(S_aux) = NULL; + + *S_aux_ptr = S_aux; + + hypre_TFree(P_i, memory_location); + hypre_TFree(P_n, memory_location); + hypre_TFree(B_i, memory_location); + hypre_TFree(B_n, memory_location); + if (method % 10 == 2) + { + hypre_TFree(B_a, memory_location); + hypre_TFree(B_a_sum, memory_location); + } + + return (ierr); +} + +/*==========================================================================*/ +/*==========================================================================*/ +/** + Generates an auxiilliary matrix, S_aux, with M-matrix properties such that S_aux = A - B, + where B is obtained by distributing positive off-diagonal contributions to 'strong' neighbors. + Here, strong neighbors are defined by the pattern of a strength matrix S. + + NOTE:: Current implementation assumes A is a square matrix with a regular partitioning of + row and columns. + + {\bf Input files:} + _hypre_parcsr_ls.h + + @return Error code. + + @param A [IN] + coefficient matrix + @param S_ptr [OUT] + strength matrix + + @see */ +/*--------------------------------------------------------------------------*/ +HYPRE_Int +hypre_BoomerAMGCreateAuxMMatrixFromS(hypre_ParCSRMatrix *A, + hypre_ParCSRMatrix *S, + hypre_ParCSRMatrix **S_aux_ptr) +{ +#ifdef HYPRE_PROFILE + hypre_profile_times[HYPRE_TIMER_ID_CREATES] -= hypre_MPI_Wtime(); +#endif + + MPI_Comm comm = hypre_ParCSRMatrixComm(A); + HYPRE_MemoryLocation memory_location_P = hypre_ParCSRMatrixMemoryLocation(A); + hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A); + hypre_ParCSRCommHandle *comm_handle; + + /* diag part of A */ + hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A); + HYPRE_Complex *A_diag_data = hypre_CSRMatrixData(A_diag); + HYPRE_Int *A_diag_i = hypre_CSRMatrixI(A_diag); + HYPRE_Int *A_diag_j = hypre_CSRMatrixJ(A_diag); + /* off-diag part of A */ + hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A); + HYPRE_Complex *A_offd_data = hypre_CSRMatrixData(A_offd); + HYPRE_Int *A_offd_i = hypre_CSRMatrixI(A_offd); + HYPRE_Int *A_offd_j = hypre_CSRMatrixJ(A_offd); + + /* diag part of S */ + hypre_CSRMatrix *S_diag = hypre_ParCSRMatrixDiag(S); + HYPRE_Int *S_diag_i = hypre_CSRMatrixI(S_diag); + HYPRE_Int *S_diag_j = hypre_CSRMatrixJ(S_diag); + /* off-diag part of S */ + hypre_CSRMatrix *S_offd = hypre_ParCSRMatrixOffd(S); + HYPRE_Int *S_offd_i = hypre_CSRMatrixI(S_offd); + HYPRE_Int *S_offd_j = hypre_CSRMatrixJ(S_offd); + + HYPRE_BigInt *row_starts = hypre_ParCSRMatrixRowStarts(A); + HYPRE_Int num_variables = hypre_CSRMatrixNumRows(A_diag); + HYPRE_BigInt global_num_vars = hypre_ParCSRMatrixGlobalNumRows(A); + HYPRE_Int num_cols_A_diag = hypre_CSRMatrixNumCols(A_diag); + HYPRE_Int num_cols_A_offd = hypre_CSRMatrixNumCols(A_offd); + + HYPRE_Int num_cols_S_diag = hypre_CSRMatrixNumCols(S_diag); + HYPRE_Int num_cols_S_offd = hypre_CSRMatrixNumCols(S_offd); + + HYPRE_BigInt *A_col_starts = hypre_ParCSRMatrixColStarts(A); + HYPRE_BigInt *A_offd_colmap = hypre_ParCSRMatrixColMapOffd(A); + + HYPRE_BigInt *S_col_starts = hypre_ParCSRMatrixColStarts(S); + HYPRE_BigInt *S_offd_colmap = hypre_ParCSRMatrixColMapOffd(S); + + HYPRE_Int A_diag_nnz = A_diag_i[num_variables]; + HYPRE_Int A_offd_nnz = A_offd_i[num_variables]; + + HYPRE_Int S_diag_nnz = S_diag_i[num_variables]; + HYPRE_Int S_offd_nnz = S_offd_i[num_variables]; + + HYPRE_Int num_nonzeros_diag; + HYPRE_Int num_nonzeros_offd = 0; + HYPRE_Int num_cols_offd = 0; + + hypre_ParCSRMatrix *S_aux; + hypre_CSRMatrix *S_aux_diag; + HYPRE_Int *S_aux_diag_i = NULL; + HYPRE_Int *S_aux_diag_j = NULL; + HYPRE_Complex *S_aux_diag_data = NULL; + /* HYPRE_Real *S_diag_data; */ + hypre_CSRMatrix *S_aux_offd; + HYPRE_Int *S_aux_offd_i = NULL; + HYPRE_Int *S_aux_offd_j = NULL; + HYPRE_Complex *S_aux_offd_data = NULL; + HYPRE_BigInt *S_aux_offd_colmap = NULL; + /* HYPRE_Real *S_offd_data; */ + /* off processor portions of S */ + hypre_CSRMatrix *S_ext = NULL; + HYPRE_Int *S_ext_i = NULL; + HYPRE_Real *S_ext_data = NULL; + HYPRE_BigInt *S_ext_j = NULL; + HYPRE_Int S_ext_nnz = 0; + + HYPRE_Real diag, row_scale, row_sum; + HYPRE_Int i, j, k, jj, jA, jS, kp, kn, m, n, P_n_max, N_n_max, cnt_n, cnt_p; + + HYPRE_Int ierr = 0; + + HYPRE_Int *dof_func_offd; + HYPRE_Int num_sends; + HYPRE_Int *int_buf_data; + HYPRE_Int index, start; + + HYPRE_Int *P_i = NULL; + HYPRE_Int *P_n = NULL; + HYPRE_BigInt *N_i = NULL; + HYPRE_Int *N_n = NULL; + HYPRE_BigInt *B_i = NULL; + HYPRE_Int *B_n = NULL; + + HYPRE_Complex pval, bval; + + HYPRE_Int *prefix_sum_workspace; + + HYPRE_MemoryLocation memory_location = hypre_ParCSRMatrixMemoryLocation(A); + + HYPRE_Int num_procs, my_id; + hypre_MPI_Comm_size(comm, &num_procs); + hypre_MPI_Comm_rank(comm, &my_id); + + // Begin: + /* Allocate memory for auxilliary strength matrix */ + S_aux = hypre_ParCSRMatrixCreate(comm, global_num_vars, global_num_vars, + row_starts, row_starts, + num_cols_A_offd, A_diag_nnz, A_offd_nnz); + + S_aux_diag = hypre_ParCSRMatrixDiag(S_aux); + hypre_CSRMatrixI(S_aux_diag) = hypre_CTAlloc(HYPRE_Int, num_variables + 1, memory_location); + hypre_CSRMatrixJ(S_aux_diag) = hypre_CTAlloc(HYPRE_Int, A_diag_nnz, memory_location); + hypre_CSRMatrixData(S_aux_diag) = hypre_CTAlloc(HYPRE_Complex, A_diag_nnz, memory_location); + + S_aux_diag_i = hypre_CSRMatrixI(S_aux_diag); + S_aux_diag_j = hypre_CSRMatrixJ(S_aux_diag); + S_aux_diag_data = hypre_CSRMatrixData(S_aux_diag); + + S_aux_offd = hypre_ParCSRMatrixOffd(S_aux); + hypre_CSRMatrixI(S_aux_offd) = hypre_CTAlloc(HYPRE_Int, num_variables + 1, memory_location); + + S_aux_offd_i = hypre_CSRMatrixI(S_aux_offd); + + + if (num_cols_A_offd) + { + hypre_CSRMatrixJ(S_aux_offd) = hypre_CTAlloc(HYPRE_Int, A_offd_nnz, memory_location); + hypre_CSRMatrixData(S_aux_offd) = hypre_CTAlloc(HYPRE_Complex, A_offd_nnz, memory_location); + + S_aux_offd_j = hypre_CSRMatrixJ(S_aux_offd); + S_aux_offd_data = hypre_CSRMatrixData(S_aux_offd); + + S_aux_offd_colmap = hypre_TAlloc(HYPRE_BigInt, num_cols_A_offd, memory_location); + hypre_ParCSRMatrixColMapOffd(S_aux) = S_aux_offd_colmap; + +#ifdef HYPRE_USING_OPENMP + #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE +#endif + for (i = 0; i < num_cols_A_offd; i++) + { + S_aux_offd_colmap[i] = A_offd_colmap[i]; + } + } + // extract off-processor portion of A + if (num_procs > 1) + { + S_ext = hypre_ParCSRMatrixExtractBExt(S, A, 0); + S_ext_i = hypre_CSRMatrixI(S_ext); + S_ext_j = hypre_CSRMatrixBigJ(S_ext); + S_ext_nnz = S_ext_i[num_cols_A_offd]; + + } + /* allocate work arrays */ + HYPRE_Int P_nnz = S_diag_nnz + S_offd_nnz + num_variables; + HYPRE_Int N_nnz = S_diag_nnz + S_offd_nnz + S_ext_nnz + num_cols_S_offd + num_variables; + P_i = hypre_CTAlloc(HYPRE_Int, P_nnz, memory_location); + P_n = hypre_CTAlloc(HYPRE_Int, (num_variables + 1), memory_location); + N_i = hypre_CTAlloc(HYPRE_BigInt, N_nnz, memory_location); + N_n = hypre_CTAlloc(HYPRE_Int, (num_variables + num_cols_A_offd + 1), memory_location); + + /* 1. Loop over matrix rows to extract positive and negative neighbors */ + P_n_max = 0; + N_n_max = 0; + kp = 0; + kn = 0; + for (i = 0; i < num_variables; i++) + { + cnt_p = kp; + cnt_n = kn; + + diag = A_diag_data[A_diag_i[i]]; + + // diagonal part + // Check for positive off-diagonal entry + for (j = A_diag_i[i]; j < A_diag_i[i + 1]; j++) + { + if ( (diag * A_diag_data[j] > 0) && (A_diag_j[j] != i)) + { + // position of positive neighbor. Store position for quick and easy access + P_i[kp++] = j; //A_diag_j[j]; + } + } + // Off-diagonal part + for (j = A_offd_i[i]; j < A_offd_i[i + 1]; j++) + { + if ( diag * A_offd_data[j] > 0 ) + { + // position of positive neighbor. Offset to distinguish from diag part columns + P_i[kp++] = j + A_diag_nnz; + } + } + // Populate negative neighbor list from S. We need this so we can sort it. + // diagonal part + for (j = S_diag_i[i]; j < S_diag_i[i + 1]; j++) + { + N_i[kn++] = S_diag_j[j] + S_col_starts[0]; + } + // Off-diagonal part + for (j = S_offd_i[i]; j < S_offd_i[i + 1]; j++) + { + N_i[kn++] = S_offd_colmap[S_offd_j[j]]; + } + + cnt_p = kp - cnt_p; + cnt_n = kn - cnt_n; + P_n_max = P_n_max > cnt_p ? P_n_max : cnt_p; + N_n_max = N_n_max > cnt_n ? N_n_max : cnt_n; + // sort neighbor arrays + hypre_BigQsort0(N_i, N_n[i], (kn - 1)); + // update index pointer arrays + P_n[i + 1] = kp; + N_n[i + 1] = kn; + } + // add negative neighbors list of external rows + kn = N_n[num_variables]; + for (i = 0; i < num_cols_A_offd; i++) + { + for (j = S_ext_i[i]; j < S_ext_i[i + 1]; j++) + { + N_i[kn++] = S_ext_j[j]; + } + jj = num_variables + i; + // sort neighbor arrays + hypre_BigQsort0(N_i, N_n[jj], (kn - 1)); + // update index pointer arrays + N_n[jj + 1] = kn; + } + hypre_CSRMatrixDestroy(S_ext); + /* 2. Loop to compute intersections of strong (negative) neighbors of common weak (positive) entry */ + B_i = hypre_CTAlloc(HYPRE_BigInt, (N_n_max * P_n[num_variables]), memory_location); + B_n = hypre_CTAlloc(HYPRE_Int, (P_n_max * num_variables + 1), memory_location); + kn = 0; + kp = 0; + for (i = 0; i < num_variables; i++) + { + // loop over positive neighbor list + for (j = P_n[i]; j < P_n[i + 1]; j++) + { + jA = P_i[j] - A_diag_nnz; + // Local row (assumes square matrix) + if ( jA < 0 ) + { + k = A_diag_j[P_i[j]] ; + } + else + { + // external row. Access external variable range of N. + k = A_offd_j[jA] + num_variables; + } + // get intersection of N_i [i] and N_i[k] + m = N_n[i + 1] - N_n[i]; + n = N_n[k + 1] - N_n[k]; + hypre_IntersectTwoBigIntegerArrays(&N_i[N_n[i]], m, &N_i[N_n[k]], n, &B_i[kn], &kp); + // update position of intersection of strong neighbors + kn += kp; + B_n[j + 1] = kn; + } + } + hypre_TFree(N_i, memory_location); + hypre_TFree(N_n, memory_location); + /* 3. Fill in data for S_aux */ + kp = 0; + kn = 0; + for (i = 0; i < num_variables; i++) + { + // first insert diagonal entry + S_aux_diag_data[kp] = A_diag_data[A_diag_i[i] ]; + S_aux_diag_j[kp++] = A_diag_j[A_diag_i[i] ]; + // diagonal part + S_diag = hypre_ParCSRMatrixDiag(S); + S_diag_i = hypre_CSRMatrixI(S_diag); + S_diag_j = hypre_CSRMatrixJ(S_diag); + + for (j = S_diag_i[i]; j < S_diag_i[i + 1]; j++) + { + jA = A_diag_i[i] + 1; + jS = S_diag_j[j]; + while (A_diag_j[jA] != jS) { jA++; } + S_aux_diag_data[kp] = A_diag_data[jA]; + S_aux_diag_j[kp++] = A_diag_j[jA]; + } + // Off-diagonal part + S_offd = hypre_ParCSRMatrixOffd(S) ; + S_offd_i = hypre_CSRMatrixI(S_offd); + S_offd_j = hypre_CSRMatrixJ(S_offd); + + for (j = S_offd_i[i]; j < S_offd_i[i + 1]; j++) + { + jA = A_offd_i[i]; + jS = S_offd_j[j]; + while (jS != A_offd_j[jA]) { jA++; } + S_aux_offd_data[kn] = A_offd_data[jA]; + S_aux_offd_j[kn++] = A_offd_j[jA]; + } + // update index pointer arrays + S_aux_diag_i[i + 1] = kp; + S_aux_offd_i[i + 1] = kn; + } + /* Distribute positive contributions to strong neighbors */ + for (i = 0; i < num_variables; i++) + { + // loop over positive neighbors and compute contributions to (strong) negative neighbors + for (j = P_n[i]; j < P_n[i + 1]; j++) + { + jA = P_i[j] - A_diag_nnz; + // positive entry is in diag part + if ( jA < 0 ) + { + pval = A_diag_data[P_i[j]] ; + } + else + { + // positive entry is in offd part + pval = A_offd_data[jA]; + } + + // compute distribution to negative neighbors + bval = -2.0 * pval / (HYPRE_Complex) (B_n[ j + 1 ] - B_n[ j ]); + + // Loop over negative connections to positive neighbor and distribute + for (k = B_n[ j ]; k < B_n[ j + 1 ]; k++) + { + HYPRE_BigInt big_col = B_i[k]; + if ( big_col >= A_col_starts[0] && big_col < A_col_starts[1]) + { + // neighbor is in diag part + jS = S_aux_diag_i[i] + 1; + jj = (HYPRE_Int) (big_col - A_col_starts[0]); + while (S_aux_diag_j[jS] != jj) { jS++; } + S_aux_diag_data[jS] -= bval; + } + else + { + // neighbor is in offd part + jj = hypre_BigBinarySearch( A_offd_colmap, big_col, num_cols_A_offd); + jS = S_aux_offd_i[i]; + while (S_aux_offd_j[jS] != jj) { jS++; } + S_aux_offd_data[jS] -= bval; + } + } + // update diagonal entry + S_aux_diag_data[S_aux_diag_i[i]] -= pval; + } + } + + hypre_CSRMatrixNumNonzeros(S_aux_diag) = S_aux_diag_i[num_variables]; + hypre_CSRMatrixNumNonzeros(S_aux_offd) = S_aux_offd_i[num_variables]; + hypre_CSRMatrixJ(S_aux_diag) = S_aux_diag_j; + hypre_CSRMatrixJ(S_aux_offd) = S_aux_offd_j; + + hypre_CSRMatrixMemoryLocation(S_aux_diag) = memory_location; + hypre_CSRMatrixMemoryLocation(S_aux_offd) = memory_location; + + // hypre_ParCSRMatrixCommPkg(S_aux) = NULL; + + *S_aux_ptr = S_aux; + + hypre_TFree(P_i, memory_location); + hypre_TFree(P_n, memory_location); + hypre_TFree(B_i, memory_location); + hypre_TFree(B_n, memory_location); + + return (ierr); +} + +HYPRE_Int hypre_BoomerAMGCreateAuxS(hypre_ParCSRMatrix *A, hypre_ParCSRMatrix *S, + hypre_ParCSRMatrix **S_aux_ptr, HYPRE_Int method) +{ + int ierr; + if (method == 0) + { + ierr = hypre_BoomerAMGCreateAuxMMatrix(A, S_aux_ptr, method); + } + else if (method == 1) + { + ierr = hypre_BoomerAMGCreateAuxMMatrixFromS(A, S, S_aux_ptr); + } + else // default + { + ierr = hypre_BoomerAMGCreateAuxMMatrix(A, S_aux_ptr, method); + } + + return (ierr); +} + + +/* Compute the intersection of x and y, placing + * the intersection in z. Additionally, the array + * x_data is associated with x, i.e., the entries + * that we grab from x, we also grab from x_data. + * If x[k] is placed in z[m], then x_data[k] goes to + * output_x_data[m]. + * + * Assumptions: + * z is of length min(x_length, y_length) + * x and y are sorted + * x_length and y_length are similar in size, otherwise, + * looping over the smaller array and doing binary search + * in the longer array is faster. + * */ +HYPRE_Int +hypre_IntersectTwoIntegerArrays(HYPRE_Int *x, + HYPRE_Int x_length, + HYPRE_Int *y, + HYPRE_Int y_length, + HYPRE_Int *z, + HYPRE_Int *intersect_length) +{ + HYPRE_Int x_index = 0; + HYPRE_Int y_index = 0; + *intersect_length = 0; + + /* Compute Intersection, looping over each array */ + while ( (x_index < x_length) && (y_index < y_length) ) + { + if (x[x_index] > y[y_index]) + { + y_index = y_index + 1; + } + else if (x[x_index] < y[y_index]) + { + x_index = x_index + 1; + } + else + { + z[*intersect_length] = x[x_index]; + x_index = x_index + 1; + y_index = y_index + 1; + *intersect_length = *intersect_length + 1; + } + } + + return 1; +} + +HYPRE_Int +hypre_IntersectTwoBigIntegerArrays(HYPRE_BigInt *x, + HYPRE_Int x_length, + HYPRE_BigInt *y, + HYPRE_Int y_length, + HYPRE_BigInt *z, + HYPRE_Int *intersect_length) +{ + HYPRE_Int x_index = 0; + HYPRE_Int y_index = 0; + *intersect_length = 0; + + /* Compute Intersection, looping over each array */ + while ( (x_index < x_length) && (y_index < y_length) ) + { + if (x[x_index] > y[y_index]) + { + y_index = y_index + 1; + } + else if (x[x_index] < y[y_index]) + { + x_index = x_index + 1; + } + else + { + z[*intersect_length] = x[x_index]; + x_index = x_index + 1; + y_index = y_index + 1; + *intersect_length = *intersect_length + 1; + } + } + + return 1; +} + diff --git a/src/parcsr_ls/par_stats.c b/src/parcsr_ls/par_stats.c index 16ef79c628..988aa91a4d 100644 --- a/src/parcsr_ls/par_stats.c +++ b/src/parcsr_ls/par_stats.c @@ -1002,6 +1002,10 @@ hypre_BoomerAMGSetupStats( void *amg_vdata, hypre_printf(" operator = %f\n", operat_cmplxty); hypre_printf(" memory = %f\n", memory_cmplxty); hypre_printf("\n\n"); + + hypre_ParAMGDataGridComplexity(amg_data) = grid_cmplxty; + hypre_ParAMGDataOperatorComplexity(amg_data) = operat_cmplxty; + hypre_ParAMGDataMemoryComplexity(amg_data) = memory_cmplxty; } if (my_id == 0) diff --git a/src/parcsr_ls/protos.h b/src/parcsr_ls/protos.h index 354f2b4614..e7f202c39e 100644 --- a/src/parcsr_ls/protos.h +++ b/src/parcsr_ls/protos.h @@ -603,6 +603,8 @@ HYPRE_Int HYPRE_BoomerAMGSetIsolatedFPoints( HYPRE_Solver solver, HYPRE_Int num_ HYPRE_BigInt *isolated_fpt_index ); HYPRE_Int HYPRE_BoomerAMGSetFPoints( HYPRE_Solver solver, HYPRE_Int num_fpt, HYPRE_BigInt *fpt_index ); +HYPRE_Int HYPRE_BoomerAMGSetUseAuxStrengthMatrix(HYPRE_Solver solver, + HYPRE_Int use_aux_strength_mat); HYPRE_Int HYPRE_BoomerAMGSetCumNnzAP ( HYPRE_Solver solver, HYPRE_Real cum_nnz_AP ); HYPRE_Int HYPRE_BoomerAMGGetCumNnzAP ( HYPRE_Solver solver, HYPRE_Real *cum_nnz_AP ); @@ -1229,6 +1231,8 @@ HYPRE_Int hypre_BoomerAMGSetFPoints( void *data, HYPRE_Int isolated, HYPRE_Int n HYPRE_Int hypre_BoomerAMGSetCumNnzAP ( void *data, HYPRE_Real cum_nnz_AP ); HYPRE_Int hypre_BoomerAMGGetCumNnzAP ( void *data, HYPRE_Real *cum_nnz_AP ); +HYPRE_Int hypre_BoomerAMGSetUseAuxStrengthMatrix( void *data, HYPRE_Int use_aux_strength_mat); + /* par_amg_setup.c */ HYPRE_Int hypre_BoomerAMGSetup ( void *amg_vdata, hypre_ParCSRMatrix *A, hypre_ParVector *f, hypre_ParVector *u ); @@ -1643,6 +1647,9 @@ HYPRE_Int hypre_IntersectTwoArrays ( HYPRE_Int *x, HYPRE_Real *x_data, HYPRE_Int HYPRE_Int hypre_IntersectTwoBigArrays ( HYPRE_BigInt *x, HYPRE_Real *x_data, HYPRE_Int x_length, HYPRE_BigInt *y, HYPRE_Int y_length, HYPRE_BigInt *z, HYPRE_Real *output_x_data, HYPRE_Int *intersect_length ); +HYPRE_Int hypre_IntersectTwoBigIntegerArrays ( HYPRE_BigInt *x, HYPRE_Int x_length, + HYPRE_BigInt *y, HYPRE_Int y_length, HYPRE_BigInt *z, + HYPRE_Int *intersect_length ); HYPRE_Int hypre_SortedCopyParCSRData ( hypre_ParCSRMatrix *A, hypre_ParCSRMatrix *B ); HYPRE_Int hypre_BoomerAMG_MyCreateS ( hypre_ParCSRMatrix *A, HYPRE_Real strength_threshold, HYPRE_Real max_row_sum, HYPRE_Int num_functions, HYPRE_Int *dof_func, hypre_ParCSRMatrix **S_ptr ); @@ -1910,7 +1917,8 @@ HYPRE_Int hypre_BoomerAMGCreateSabsHost ( hypre_ParCSRMatrix *A, HYPRE_Real stre HYPRE_Int hypre_BoomerAMGCreate2ndSDevice( hypre_ParCSRMatrix *S, HYPRE_Int *CF_marker, HYPRE_Int num_paths, HYPRE_BigInt *coarse_row_starts, hypre_ParCSRMatrix **C_ptr); HYPRE_Int hypre_BoomerAMGMakeSocFromSDevice( hypre_ParCSRMatrix *A, hypre_ParCSRMatrix *S); - +HYPRE_Int hypre_BoomerAMGCreateAuxS(hypre_ParCSRMatrix *A, hypre_ParCSRMatrix *S, + hypre_ParCSRMatrix **S_aux_ptr, HYPRE_Int method); /* par_sv_interp.c */ HYPRE_Int hypre_BoomerAMGSmoothInterpVectors ( hypre_ParCSRMatrix *A, HYPRE_Int num_smooth_vecs, diff --git a/src/parcsr_mv/par_csr_communication.c b/src/parcsr_mv/par_csr_communication.c index bdaa8d8ffc..d28ea3229c 100644 --- a/src/parcsr_mv/par_csr_communication.c +++ b/src/parcsr_mv/par_csr_communication.c @@ -669,7 +669,7 @@ hypre_ParCSRCommHandleDestroy( hypre_ParCSRCommHandle *comm_handle ) if (!hypre_GetGpuAwareMPI()) { hypre_MemoryLocation act_send_memory_location = - hypre_GetActualMemLocation(hypre_ParCSRCommHandleSendMemoryLocation(comm_handle)); + hypre_GetActualMemLocation(hypre_ParCSRCommHandleSendMemoryLocation(comm_handle)); if ( act_send_memory_location == hypre_MEMORY_DEVICE || act_send_memory_location == hypre_MEMORY_UNIFIED ) @@ -679,7 +679,7 @@ hypre_ParCSRCommHandleDestroy( hypre_ParCSRCommHandle *comm_handle ) } hypre_MemoryLocation act_recv_memory_location = - hypre_GetActualMemLocation(hypre_ParCSRCommHandleRecvMemoryLocation(comm_handle)); + hypre_GetActualMemLocation(hypre_ParCSRCommHandleRecvMemoryLocation(comm_handle)); if ( act_recv_memory_location == hypre_MEMORY_DEVICE || act_recv_memory_location == hypre_MEMORY_UNIFIED ) diff --git a/src/test/Makefile b/src/test/Makefile index 292313ebf8..5bf88d33f5 100644 --- a/src/test/Makefile +++ b/src/test/Makefile @@ -71,6 +71,7 @@ HYPRE_DRIVERS =\ sstruct_fac.c\ ij_mm.c\ zboxloop.c +# cogent_proxy_2d_RIPS.c\ HYPRE_DRIVERS_CXX =\ cxx_ij.cxx\ @@ -129,6 +130,10 @@ distclean: clean ################################################################## # C +# +cogent_proxy_2d_RIPS: cogent_proxy_2d_RIPS.o + @echo "Building" $@ "... " + ${LINK_CC} -o $@ $< ${LFLAGS} ij: ij.o @echo "Building" $@ "... " diff --git a/src/test/cogent_proxy_2d_RIPS.c b/src/test/cogent_proxy_2d_RIPS.c new file mode 100644 index 0000000000..87f74f94f0 --- /dev/null +++ b/src/test/cogent_proxy_2d_RIPS.c @@ -0,0 +1,1405 @@ +/****************************************************************************** + * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other + * HYPRE Project Developers. See the top-level COPYRIGHT file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + ******************************************************************************/ + +/* + Example cogent_proxy_2d + + Verion 6 (2/6/2022) + + Compile with: make cogent_proxy_2d + + To see options: cogent_proxy_2d -help + + Description: This code solves a system corresponding to a discretization + of the system of equations + + (M + N_1)u - N_2 v = 1 + Mu + v = 0 + + on the unit square, where the matrices M, N_1 and N_2 are obtained from + the a second-order, centered-difference discretization of the operators + + M(w) = -a*w_xx + N_1(w) = -b*w_yy + N_2(w) = -c*w_yy + + where a, b and c are positive constants. The corresponding boundary + conditions are zero Dirichlet in x and periodic in y. + + The domain is split into an Nx x 1 processor grid. + Each processor's piece of the grid has nx x ny cells + We use cell-centered variables, and, therefore, the + nodes are not shared. Note that we have two variables, u and + v, and need only one part to describe the domain. + +*/ + +#include +#include "_hypre_utilities.h" +#include "HYPRE_sstruct_ls.h" +#include "HYPRE_krylov.h" +#include "_hypre_parcsr_mv.h" +#include "_hypre_parcsr_ls.h" + + + +int main (int argc, char *argv[]) +{ + int i, j; + + int myid, num_procs; + + int nx, ny, Nx, pi, pj; + HYPRE_Real hx, hy, hx2, hy2; + int ilower[2], iupper[2]; + + int solver_id; + + int object_type = HYPRE_PARCSR; + + HYPRE_SStructGrid grid; + HYPRE_SStructStencil stencil_v; + HYPRE_SStructStencil stencil_u; + HYPRE_SStructGraph graph; + HYPRE_SStructMatrix A; + HYPRE_SStructVector b; + HYPRE_SStructVector x; + + HYPRE_Solver par_solver; + + /* solver params */ + int max_iter = 10; + double tol = 1.e-6; + int print_level = 3; + int print_cgrid = 0; + + /* ILU params */ + int ILU_type = 11; + int ILU_fill_level = 0; + double ILU_drop_threshold = -1.; + double ILU_B_drop_threshold = -1.; + double ILU_EF_drop_threshold = -1.; + double ILU_Schur_drop_threshold = -1.; + int ILU_max_nnz_per_row = 1000; + int ILU_max_schur_iter = 5; + + /* MGR params */ + HYPRE_Int mgr_bsize = 2; + HYPRE_Int mgr_nlevels = 1; + HYPRE_Int mgr_non_c_to_f = 1; + HYPRE_Int mgr_frelax_method = 1; + HYPRE_Int *mgr_num_cindexes = NULL; + HYPRE_Int **mgr_cindexes = NULL; + HYPRE_Int mgr_relax_type = 0; + HYPRE_Int mgr_num_relax_sweeps = 1; + HYPRE_Int mgr_interp_type = 3; + HYPRE_Int mgr_num_interp_sweeps = 1; + HYPRE_Int mgr_gsmooth_type = 0; + HYPRE_Int mgr_num_gsmooth_sweeps = 0; + HYPRE_Int mgr_restrict_type = 0; + HYPRE_Int mgr_num_restrict_sweeps = 0; + HYPRE_Int mgr_cpoint = 0; + HYPRE_Real mgr_csolve_threshold = 0.25; + HYPRE_Real mgr_csolve_max_iter = 20; + HYPRE_Int mgr_use_non_galerkin = 0; + HYPRE_Int mgr_pmax = 0; + + HYPRE_Int strength_method = -1; + + HYPRE_Int semicoarsening_dir = 0; + HYPRE_Int semicoarsening_offset = 0; + + /* Initialize MPI */ + MPI_Init(&argc, &argv); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + + /* Set defaults */ + Nx = 1; + nx = 8; + ny = 8; + solver_id = 0; + HYPRE_Real a_fac = 1.e-5; + HYPRE_Real b_fac = 2.e2; + HYPRE_Real c_fac = 0.352; + // double R0 = 1.6; + // double axisym_fac = sqrt(2.*3.14159*R0); + double axisym_fac = 1.; + + /* Parse command line */ + { + int arg_index = 0; + int print_usage = 0; + + while (arg_index < argc) + { + if ( strcmp(argv[arg_index], "-Nx") == 0 ) + { + arg_index++; + Nx = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-nx") == 0 ) + { + arg_index++; + nx = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ny") == 0 ) + { + arg_index++; + ny = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-a") == 0 ) + { + arg_index++; + a_fac = atof(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-b") == 0 ) + { + arg_index++; + b_fac = atof(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-c") == 0 ) + { + arg_index++; + c_fac = atof(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-solver") == 0 ) + { + arg_index++; + solver_id = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-max_iter") == 0 ) + { + arg_index++; + max_iter = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-tol") == 0 ) + { + arg_index++; + tol = atof(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-print_level") == 0 ) + { + arg_index++; + print_level = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-print_cgrid") == 0 ) + { + arg_index++; + print_cgrid = atoi(argv[arg_index++]); + } + /* ILU options */ + else if ( strcmp(argv[arg_index], "-ilu_type") == 0 ) + { + arg_index++; + ILU_type = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ilu_fill_level") == 0 ) + { + arg_index++; + ILU_fill_level = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ilu_drop") == 0 ) + { + arg_index++; + ILU_drop_threshold = atof(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ilu_B_drop") == 0 ) + { + arg_index++; + ILU_B_drop_threshold = atof(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ilu_EF_drop") == 0 ) + { + arg_index++; + ILU_EF_drop_threshold = atof(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ilu_schur_drop") == 0 ) + { + arg_index++; + ILU_Schur_drop_threshold = atof(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ilu_max_nnz") == 0 ) + { + arg_index++; + ILU_max_nnz_per_row = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-ilu_max_schur_iter") == 0 ) + { + arg_index++; + ILU_max_schur_iter = atoi(argv[arg_index++]); + } + /* MGR options */ + else if ( strcmp(argv[arg_index], "-mgr_bsize") == 0 ) + { + /* mgr block size */ + arg_index++; + mgr_bsize = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_nlevels") == 0 ) + { + /* mgr number of coarsening levels */ + arg_index++; + mgr_nlevels = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_cpoint") == 0 ) + { + /* coarse point index in block system */ + arg_index++; + mgr_cpoint = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_non_c_to_f") == 0 ) + { + /* mgr intermediate coarse grid strategy */ + arg_index++; + mgr_non_c_to_f = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_frelax_method") == 0 ) + { + /* mgr F-relaxation strategy: single/ multi level */ + arg_index++; + mgr_frelax_method = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_relax_type") == 0 ) + { + /* relax type for "single level" F-relaxation */ + arg_index++; + mgr_relax_type = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_relax_sweeps") == 0 ) + { + /* number of relaxation sweeps */ + arg_index++; + mgr_num_relax_sweeps = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_interp_type") == 0 ) + { + /* interpolation type */ + arg_index++; + mgr_interp_type = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_interp_sweeps") == 0 ) + { + /* number of interpolation sweeps*/ + arg_index++; + mgr_num_interp_sweeps = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_gsmooth_type") == 0 ) + { + /* global smoother type */ + arg_index++; + mgr_gsmooth_type = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_gsmooth_sweeps") == 0 ) + { + /* number of global smooth sweeps*/ + arg_index++; + mgr_num_gsmooth_sweeps = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_restrict_type") == 0 ) + { + /* restriction type */ + arg_index++; + mgr_restrict_type = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_restrict_sweeps") == 0 ) + { + /* number of restriction sweeps*/ + arg_index++; + mgr_num_restrict_sweeps = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_csolve_threshold") == 0 ) + { + /* number of restriction sweeps*/ + arg_index++; + mgr_csolve_threshold = atof(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_csolve_max_iter") == 0 ) + { + /* number of restriction sweeps*/ + arg_index++; + mgr_csolve_max_iter = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_use_non_galerkin") == 0 ) + { + /* use non Galerkin coarse grid */ + arg_index++; + mgr_use_non_galerkin = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-mgr_pmax") == 0 ) + { + /* Pmax elements for non Galerkin case */ + arg_index++; + mgr_pmax = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-strength_method") == 0 ) + { + /* Strength of connection method */ + arg_index++; + strength_method = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-sdir") == 0 ) + { + /* Semicoarsening direction */ + arg_index++; + semicoarsening_dir = atoi(argv[arg_index++]); + } + else if ( strcmp(argv[arg_index], "-soffset") == 0 ) + { + /* Semicoarsening offset */ + arg_index++; + semicoarsening_offset = atoi(argv[arg_index++]); + } + /* end mgr options */ + else if ( strcmp(argv[arg_index], "-help") == 0 ) + { + print_usage = 1; + break; + } + else + { + arg_index++; + } + } + + if ((print_usage) && (myid == 0)) + { + printf("\n"); + printf("Usage: %s []\n", argv[0]); + printf("\n"); + printf(" -Nx : number of processors in x (default: 1)\n"); + printf(" -nx : grid size per processor in x (default: 8)\n"); + printf(" -ny : grid size per processor in y (default: 8)\n"); + printf(" -a : a factor (default: 1.e-5)\n"); + printf(" -b : b factor (default: 2.e2)\n"); + printf(" -solver : solver ID\n"); + printf(" -max_iter : max solver iterations (default 20) \n"); + printf(" -tol : solver convergence tolerance (default 1.e-7) \n"); + printf(" -print_level : print level (default 3) \n"); + printf(" 0 - hypre-ILU (default)\n"); + printf(" 1 - AMG \n"); + printf(" 2 - MGR \n"); + printf(" -ilu_type : ILU solver type (default 11) \n"); + printf(" -ilu_fill_level : ILU(k) fill level k (default 0) \n"); + printf(" -ilu_drop : ILU drop tolerance (default 1.e-2) \n"); + printf(" -ilu_B_drop : ILU B block drop tolerance (default 1.e-2) \n"); + printf(" -ilu_EF_drop : ILU E and F block drop tolerance (default 1.e-2) \n"); + printf(" -ilu_schur_drop : ILU Schur drop tolerance (default 1.e-2) \n"); + printf(" -ilu_max_nnz : ILU max number of nonzeros per row (default 1000) \n"); + printf(" -ilu_max_schur_iter : ILU max Schur iteratations (default 5) \n"); + /* MGR options */ + hypre_printf(" -mgr_bsize : set block size = val\n"); + hypre_printf(" -mgr_nlevels : set number of coarsening levels = val\n"); + hypre_printf(" to be kept till the coarsest grid = val\n"); + hypre_printf(" -mgr_non_c_to_f : set strategy for intermediate coarse grid \n"); + hypre_printf(" -mgr_non_c_to_f 0 : Allow some non Cpoints to be labeled \n"); + hypre_printf(" Cpoints on intermediate grid \n"); + hypre_printf(" -mgr_non_c_to_f 1 : set non Cpoints strictly to Fpoints \n"); + hypre_printf(" -mgr_frelax_method : set F-relaxation strategy \n"); + hypre_printf(" -mgr_frelax_method 0 : Use 'single-level smoother' strategy \n"); + hypre_printf(" for F-relaxation \n"); + hypre_printf(" -mgr_frelax_method 1 : Use a 'multi-level smoother' strategy \n"); + hypre_printf(" for F-relaxation \n"); + hypre_printf(" -mgr_csolve_threshold : AMG coarse solve strong threshold = val\n"); + hypre_printf(" -mgr_csolve_max_iter : AMG coarse solve max iterations = val\n"); + hypre_printf(" -strength_method : method for creating AMG strength matrix= val\n"); + /* end MGR options */ + printf("\n"); + } + + if (print_usage) + { + MPI_Finalize(); + return (0); + } + } + + /* Set up the Nx x 1 processor grid. pi, pj indicate + the position of the executing process in the processor grid. + The local problem size nx x ny. */ + + pj = 0; + pi = myid; + + /* Figure out the extents of each processor's piece of the grid. */ + ilower[0] = pi * nx; + ilower[1] = pj * ny; + + iupper[0] = ilower[0] + nx - 1; + iupper[1] = ilower[1] + ny - 1; + + /* Get the mesh size */ + + hx = 0.3 * axisym_fac / (Nx * nx); + hy = 3.77 * axisym_fac / ny; + hx2 = hx * hx; + hy2 = hy * hy; + + /* 1. Set up a grid - we have one part and two variables */ + { + int nparts = 1; + int part = 0; + int ndim = 2; + int periodic[2]; + + /* Create an empty 2D grid object */ + HYPRE_SStructGridCreate(MPI_COMM_WORLD, ndim, nparts, &grid); + + /* Add a new box to the grid */ + HYPRE_SStructGridSetExtents(grid, part, ilower, iupper); + + /* Set the variable type and number of variables on each part.*/ + { + int i; + int nvars = 2; + HYPRE_SStructVariable vartypes[2] = {HYPRE_SSTRUCT_VARIABLE_CELL, + HYPRE_SSTRUCT_VARIABLE_CELL + }; + + for (i = 0; i < nparts; i++) + { + HYPRE_SStructGridSetVariables(grid, i, nvars, vartypes); + } + } + + /* Make the second direction periodic */ + periodic[0] = 0; + periodic[1] = ny; + HYPRE_SStructGridSetPeriodic(grid, part, periodic); + + /* This is a collective call finalizing the grid assembly. + The grid is now ``ready to be used'' */ + HYPRE_SStructGridAssemble(grid); + } + + /* 2. Define the discretization stencils */ + { + int entry; + int stencil_size; + int var; + int ndim = 2; + + int M_offsets[5][2] = {{0, -1}, {-1, 0}, {0, 0}, {1, 0}, {0, 1}}; + int N_offsets[3][2] = {{0, -1}, {0, 0}, {0, 1}}; + int I_offsets[1][2] = {{0, 0}}; + + /* Stencil object for variable u (labeled as variable 0) */ + { + int stencil_size_uu = 5; + int stencil_size_uv = 3; + stencil_size = stencil_size_uu + stencil_size_uv; + + HYPRE_SStructStencilCreate(ndim, stencil_size, &stencil_u); + + /* The first stencil_size_uu entries are for the u-u connections */ + var = 0; /* connect to variable 0 */ + for (entry = 0; entry < stencil_size_uu; entry++) + { + HYPRE_SStructStencilSetEntry(stencil_u, entry, M_offsets[entry], var); + } + + /* The last stencil_size_uv entries are for the u-v connections */ + var = 1; /* connect to variable 1 */ + for (entry = stencil_size_uu, j = 0; entry < stencil_size; entry++, j++) + { + HYPRE_SStructStencilSetEntry(stencil_u, entry, N_offsets[j], var); + } + } + + /* Stencil object for variable v (variable 1) */ + { + int stencil_size_vu = 5; + int stencil_size_vv = 1; + stencil_size = stencil_size_vu + stencil_size_vv; + + HYPRE_SStructStencilCreate(ndim, stencil_size, &stencil_v); + + /* These are the v-u connections */ + var = 0; /* Connect to variable 0 */ + for (entry = 0; entry < stencil_size_vu; entry++) + { + HYPRE_SStructStencilSetEntry(stencil_v, entry, M_offsets[entry], var); + } + + /* These are the v-v connections */ + var = 1; /* Connect to variable 1 */ + for (entry = stencil_size_vu, j = 0; entry < stencil_size; entry++, j++) + { + HYPRE_SStructStencilSetEntry(stencil_v, entry, I_offsets[j], var); + } + } + } + + /* 3. Set up the Graph - this determines the non-zero structure + of the matrix and allows non-stencil relationships between the parts. */ + { + int var; + int part = 0; + + /* Create the graph object */ + HYPRE_SStructGraphCreate(MPI_COMM_WORLD, grid, &graph); + HYPRE_SStructGraphSetObjectType(graph, object_type); + + /* Assign the u-stencil we created to variable u (variable 0) */ + var = 0; + HYPRE_SStructGraphSetStencil(graph, part, var, stencil_u); + + /* Assign the v-stencil we created to variable v (variable 1) */ + var = 1; + HYPRE_SStructGraphSetStencil(graph, part, var, stencil_v); + + /* Assemble the graph */ + HYPRE_SStructGraphAssemble(graph); + } + + /* 4. Set up the SStruct Matrix */ + { + int nentries; + int nvalues; + int var; + int part = 0; + int num_local_cells = nx * ny; + double M_stencil[5]; + double N1_stencil[3]; + double N2_stencil[3]; + + M_stencil[0] = 0.; + M_stencil[1] = -a_fac * ( 1. / hx2 ); + M_stencil[2] = -a_fac * (-2. / hx2 ); + M_stencil[3] = -a_fac * ( 1. / hx2 ); + M_stencil[4] = 0.; + + N1_stencil[0] = -b_fac / hy2; + N1_stencil[1] = 2.*b_fac / hy2; + N1_stencil[2] = -b_fac / hy2; + + N2_stencil[0] = -c_fac / hy2; + N2_stencil[1] = 2.*c_fac / hy2; + N2_stencil[2] = -c_fac / hy2; + + /* Create an empty matrix object */ + HYPRE_SStructMatrixCreate(MPI_COMM_WORLD, graph, &A); + + /* Set the object type (by default HYPRE_SSTRUCT). This determines the + data structure used to store the matrix. If you want to use + unstructured solvers, e.g. BoomerAMG, the object type should be + HYPRE_PARCSR. If the problem is purely structured (with one part), you + may want to use HYPRE_STRUCT to access the structured solvers. */ + HYPRE_SStructMatrixSetObjectType(A, object_type); + + /* Indicate that the matrix coefficients are ready to be set */ + HYPRE_SStructMatrixInitialize(A); + + /* Each processor must set the stencil values for their boxes on each part. + In this example, we only set stencil entries and therefore use + HYPRE_SStructMatrixSetBoxValues. If we need to set non-stencil entries, + we have to use HYPRE_SStructMatrixSetValues. */ + + /* First set the u-stencil entries. Note that + HYPRE_SStructMatrixSetBoxValues can only set values corresponding + to stencil entries for the same variable. Therefore, we must set the + entries for each variable within a stencil with separate function calls. + For example, below the u-u connections and u-v connections are handled + in separate calls. */ + { + int i, j; + HYPRE_Real *u_values; + int u_u_indices[5] = {0, 1, 2, 3, 4}; + + var = 0; /* Set values for the u connections */ + + /* First the u-u connections */ + nentries = 5; + nvalues = nentries * num_local_cells; + u_values = (HYPRE_Real*) calloc(nvalues, sizeof(HYPRE_Real)); + + for (i = 0; i < nvalues; i += nentries) + { + for (j = 0; j < nentries; ++j) + { + u_values[i + j] = 0.; //M_stencil[j]; + } + + u_values[i] += 0.;//N1_stencil[0]; + u_values[i + 2] += 0.; //N1_stencil[1]; + u_values[i + 4] += 0.; //N1_stencil[2]; + } + + HYPRE_SStructMatrixAddToBoxValues(A, part, ilower, iupper, + var, nentries, + u_u_indices, u_values); + free(u_values); + + /* Next the u-v connections */ + int u_v_indices[3] = {5, 6, 7}; + nentries = 3; + nvalues = nentries * num_local_cells; + u_values = (HYPRE_Real*) calloc(nvalues, sizeof(HYPRE_Real)); + + for (i = 0; i < nvalues; i += nentries) + { + for (j = 0; j < nentries; ++j) + { + u_values[i + j] = -N2_stencil[j]; + } + } + + HYPRE_SStructMatrixAddToBoxValues(A, part, ilower, iupper, + var, nentries, + u_v_indices, u_values); + free(u_values); + } + /* Now set the v-stencil entries */ + { + int i, j; + HYPRE_Real *v_values; + int v_u_indices[5] = {0, 1, 2, 3, 4}; + + var = 1; /* the v connections */ + + /* the v-u connections */ + nentries = 5; + nvalues = nentries * num_local_cells; + v_values = (HYPRE_Real*) calloc(nvalues, sizeof(HYPRE_Real)); + + for (i = 0; i < nvalues; i += nentries) + { + for (j = 0; j < nentries; ++j) + { + v_values[i + j] = M_stencil[j]; + } + } + + HYPRE_SStructMatrixAddToBoxValues(A, part, ilower, iupper, + var, nentries, + v_u_indices, v_values); + + free(v_values); + + /* the v-v connections */ + int v_v_indices[1] = {5}; + nentries = 1; + nvalues = nentries * num_local_cells; + v_values = (HYPRE_Real*) calloc(nvalues, sizeof(HYPRE_Real)); + + for (i = 0; i < nvalues; i += nentries) + { + v_values[i] = 1.; + } + + HYPRE_SStructMatrixAddToBoxValues(A, part, ilower, iupper, + var, nentries, + v_v_indices, v_values); + + free(v_values); + } + } + + /* 5. Incorporate the zero Dirichlet boundary condition in the first coordinate: + at those boundaries, set the stencil entries that reach across it to zero.*/ + { + int bc_ilower[2]; + int bc_iupper[2]; + int nentries = 1; + int var; + int stencil_indices[1]; + + int part = 0; + + /* Recall: pi, pj describe position in the processor grid */ + if (pi == 0) + { + int nvalues = nentries * ny; /* number of stencil entries times the length + of one side of my grid box */ + + HYPRE_Real* values = (HYPRE_Real*) calloc(nvalues, sizeof(HYPRE_Real)); + for (j = 0; j < nvalues; j++) + { + values[j] = 0.0; + } + + /* Bottom plane of grid points */ + bc_ilower[0] = pi * nx; + bc_ilower[1] = pj * ny; + + bc_iupper[0] = bc_ilower[0]; + bc_iupper[1] = bc_ilower[1] + ny - 1; + + stencil_indices[0] = 1; + + var = 0; + HYPRE_SStructMatrixSetBoxValues(A, part, bc_ilower, bc_iupper, + var, nentries, + stencil_indices, values); + + var = 1; + HYPRE_SStructMatrixSetBoxValues(A, part, bc_ilower, bc_iupper, + var, nentries, + stencil_indices, values); + free(values); + } + + if (pi == Nx - 1) + { + int nvalues = nentries * ny; /* number of stencil entries times the length + of one side of my grid box */ + + HYPRE_Real* values = (HYPRE_Real*) calloc(nvalues, sizeof(HYPRE_Real)); + for (j = 0; j < nvalues; j++) + { + values[j] = 0.0; + } + + /* upper plane of grid points */ + bc_ilower[0] = pi * nx + nx - 1; + bc_ilower[1] = pj * ny; + + bc_iupper[0] = bc_ilower[0]; + bc_iupper[1] = bc_ilower[1] + ny - 1; + + stencil_indices[0] = 3; + + var = 0; + HYPRE_SStructMatrixSetBoxValues(A, part, bc_ilower, bc_iupper, + var, nentries, + stencil_indices, values); + + var = 1; + HYPRE_SStructMatrixSetBoxValues(A, part, bc_ilower, bc_iupper, + var, nentries, + stencil_indices, values); + free(values); + } + + /* Done with the boundary modifications, since the other direction + is periodic + */ + + } + + /* This is a collective call finalizing the matrix assembly. + The matrix is now ``ready to be used'' */ + HYPRE_SStructMatrixAssemble(A); + + + HYPRE_SStructMatrixPrint("A", A, 0); + + /* 5. Set up SStruct Vectors for b and x */ + { + int nvalues = nx * ny; + HYPRE_Real *values; + int part = 0; + int var; + + values = (HYPRE_Real*) calloc(nvalues, sizeof(HYPRE_Real)); + + /* Create an empty vector object */ + HYPRE_SStructVectorCreate(MPI_COMM_WORLD, grid, &b); + HYPRE_SStructVectorCreate(MPI_COMM_WORLD, grid, &x); + + /* Set the object type for the vectors + to be the same as was already set for the matrix */ + HYPRE_SStructVectorSetObjectType(b, object_type); + HYPRE_SStructVectorSetObjectType(x, object_type); + + /* Indicate that the vector coefficients are ready to be set */ + HYPRE_SStructVectorInitialize(b); + HYPRE_SStructVectorInitialize(x); + + /* Set the values for b */ + /* project out kernel for a consistent rhs (for 4th-order problem) */ + /* This is hardcoded to use semicoarsening in y */ + HYPRE_Real gamma = (double)1. / ny; + HYPRE_Real *pvalues = (HYPRE_Real*) calloc(nvalues, sizeof(HYPRE_Real)); + for (i = 0; i < nvalues; i ++) + { + values[i] = 1.0; //(double)rand() / (double)RAND_MAX; //1.0; + } + memcpy(pvalues, values, nvalues * sizeof(HYPRE_Real)); + for (HYPRE_Int j = 0; j < ny; j++) + { + for (i = 0; i < nx; i++) + { + HYPRE_Int idx = i + j * nx; + for (HYPRE_Int k = 0; k < ny; k++) + { + HYPRE_Int pos = idx + k * semicoarsening_offset; + HYPRE_Real val = pos >= nvalues ? pvalues[pos % nvalues] : pvalues[pos]; + values[idx] -= gamma * val; + } + } + } + free(pvalues); + for (i = 0; i < nvalues; i++) + { + printf("values[%d] = %f, gamma = %f \n", i, values[i], gamma); + } + + var = 0; + HYPRE_SStructVectorSetBoxValues(b, part, ilower, iupper, var, values); + + for (i = 0; i < nvalues; i ++) + { + values[i] = 0.0; + } + var = 1; + HYPRE_SStructVectorSetBoxValues(b, part, ilower, iupper, var, values); + + /* Set the values for the initial guess */ + // for (i = 0; i < nvalues; i ++) + // values[i] = (double)rand() / (double)RAND_MAX; + var = 0; + HYPRE_SStructVectorSetBoxValues(x, part, ilower, iupper, var, values); + + var = 1; + HYPRE_SStructVectorSetBoxValues(x, part, ilower, iupper, var, values); + + free(values); + + /* This is a collective call finalizing the vector assembly. + The vector is now ``ready to be used'' */ + HYPRE_SStructVectorAssemble(b); + HYPRE_SStructVectorAssemble(x); + } + + /* 6. Set up and use a solver + (Solver options can be found in the Reference Manual.) */ + { + HYPRE_Real final_res_norm; + int its; + + /* Get the object for the matrix and vectors. */ + HYPRE_ParCSRMatrix par_A; + HYPRE_ParVector par_b; + HYPRE_ParVector par_x; + HYPRE_SStructMatrixGetObject(A, (void **) &par_A); + HYPRE_SStructVectorGetObject(b, (void **) &par_b); + HYPRE_SStructVectorGetObject(x, (void **) &par_x); + + if (solver_id == 0) /* ILU */ + { + HYPRE_ILUCreate(&par_solver); + + /* + ILU_type options: + 0: Block-Jacobi with ILU(k) + 1: Block-Jacobi with ILUT + 10: GMRES with ILU(k) + 11: GMRES with ILUT + 20: Newton Schulz Hotelling with ILU(k) + 21: Newton Schulz Hotelling with ILUT + 30: Restricted Additive Schwarz with ILU(k) + 31: Restricted Additive Schwarz with ILUT + 40: DDPQ-GMRES with ILU(k) + 41: DDPQ-GMRES with ILUT + 50: GMRES with RAP-ILU(0) using MILU(0) for P + */ + + HYPRE_ILUSetType(par_solver, ILU_type); + + if (ILU_fill_level >= 0) + { + HYPRE_ILUSetLevelOfFill(par_solver, ILU_fill_level); + } + + if (ILU_drop_threshold >= 0.) + { + HYPRE_ILUSetDropThreshold(par_solver, ILU_drop_threshold); + } + + if (ILU_B_drop_threshold >= 0. && + ILU_EF_drop_threshold >= 0. && + ILU_Schur_drop_threshold >= 0. ) + { + double ILU_drop_threshold_array[3] = {ILU_B_drop_threshold, + ILU_EF_drop_threshold, + ILU_Schur_drop_threshold + }; + HYPRE_ILUSetDropThresholdArray(par_solver, ILU_drop_threshold_array); + } + + HYPRE_ILUSetMaxNnzPerRow(par_solver, ILU_max_nnz_per_row); + HYPRE_ILUSetMaxIter(par_solver, max_iter); + HYPRE_ILUSetSchurMaxIter(par_solver, ILU_max_schur_iter); + HYPRE_ILUSetTol(par_solver, tol); + HYPRE_ILUSetPrintLevel(par_solver, print_level); + + /* do the setup */ + HYPRE_ILUSetup(par_solver, par_A, par_b, par_x); + + /* do the solve */ + HYPRE_ILUSolve(par_solver, par_A, par_b, par_x); + + /* get some info */ + HYPRE_ILUGetNumIterations(par_solver, &its); + HYPRE_ILUGetFinalRelativeResidualNorm(par_solver, &final_res_norm); + + /* clean up */ + HYPRE_ILUDestroy(par_solver); + } + else if (solver_id == 1) /* AMG */ + { + HYPRE_BoomerAMGCreate(&par_solver); + HYPRE_BoomerAMGSetCoarsenType(par_solver, 6); + HYPRE_BoomerAMGSetOldDefault(par_solver); + HYPRE_BoomerAMGSetStrongThreshold(par_solver, 0.25); + HYPRE_BoomerAMGSetTol(par_solver, tol); + HYPRE_BoomerAMGSetPrintLevel(par_solver, 3); + HYPRE_BoomerAMGSetPrintFileName(par_solver, "polcor.out.log"); + HYPRE_BoomerAMGSetMaxIter(par_solver, max_iter); + // HYPRE_BoomerAMGSetUseAuxStrengthMatrix(par_solver, 1); + /* do the setup */ + HYPRE_BoomerAMGSetup(par_solver, par_A, par_b, par_x); + + /* do the solve */ + HYPRE_BoomerAMGSolve(par_solver, par_A, par_b, par_x); + + /* get some info */ + HYPRE_BoomerAMGGetNumIterations(par_solver, &its); + HYPRE_BoomerAMGGetFinalRelativeResidualNorm(par_solver, + &final_res_norm); + /* clean up */ + HYPRE_BoomerAMGDestroy(par_solver); + } + else if (solver_id == 2) /* MGR */ + { + HYPRE_MGRCreate(&par_solver); + + mgr_num_cindexes = hypre_CTAlloc(HYPRE_Int, mgr_nlevels, HYPRE_MEMORY_HOST); + for (i = 0; i < mgr_nlevels; i++) + { + /* assume 1 coarse index per level */ + mgr_num_cindexes[i] = 1; + } + mgr_cindexes = hypre_CTAlloc(HYPRE_Int*, mgr_nlevels, HYPRE_MEMORY_HOST); + for (i = 0; i < mgr_nlevels; i++) + { + mgr_cindexes[i] = hypre_CTAlloc(HYPRE_Int, mgr_num_cindexes[i], HYPRE_MEMORY_HOST); + } + for (i = 0; i < mgr_nlevels; i++) + { + /* assume coarse point is at index 0 */ + mgr_cindexes[i][0] = mgr_cpoint; + } + + /* set MGR data by block */ +#if 0 + { + HYPRE_BigInt rowstart = hypre_ParCSRMatrixFirstRowIndex(par_A); + HYPRE_BigInt rowend = hypre_ParCSRMatrixLastRowIndex(par_A); + HYPRE_Int fsize = (rowend - rowstart + 1) / 2 ; + HYPRE_BigInt next_block = rowstart + fsize; + // hypre_printf("%d: row_start = %d, next_block = %d, diff = %d \n",myid, rowstart, next_block,(rowend - rowstart + 1) ); + + HYPRE_BigInt idx_array[2] = {rowstart, next_block}; + HYPRE_MGRSetCpointsByContiguousBlock( par_solver, mgr_bsize, mgr_nlevels, idx_array, + mgr_num_cindexes, mgr_cindexes); + } +#else + { + int var, n; + HYPRE_Int* CF_index = (HYPRE_Int*) calloc(nx * ny * 2, sizeof(HYPRE_Int)); + HYPRE_Int* ptr = CF_index; + for (var = 0; var < 2; ++var) + { + for (n = 0; n < nx * ny; ++n) { *ptr++ = var; } + } + HYPRE_MGRSetCpointsByPointMarkerArray(par_solver, mgr_bsize, mgr_nlevels, mgr_num_cindexes, + mgr_cindexes, CF_index); + } +#endif + + // HYPRE_MGRSetBlockJacobiBlockSize(par_solver, 3); + + /* set intermediate coarse grid strategy */ + HYPRE_MGRSetNonCpointsToFpoints(par_solver, mgr_non_c_to_f); + /* set F relaxation strategy */ + HYPRE_MGRSetFRelaxMethod(par_solver, mgr_frelax_method); + /* set relax type for single level F-relaxation and post-relaxation */ + HYPRE_MGRSetRelaxType(par_solver, mgr_relax_type); + HYPRE_MGRSetNumRelaxSweeps(par_solver, mgr_num_relax_sweeps); + /* set interpolation type */ + HYPRE_MGRSetRestrictType(par_solver, mgr_restrict_type); + HYPRE_MGRSetNumRestrictSweeps(par_solver, mgr_num_restrict_sweeps); + HYPRE_MGRSetInterpType(par_solver, mgr_interp_type); + HYPRE_MGRSetNumInterpSweeps(par_solver, mgr_num_interp_sweeps); + /* set print level */ + HYPRE_MGRSetPrintLevel(par_solver, print_level); + /* set max iterations */ + HYPRE_MGRSetMaxIter(par_solver, max_iter); + HYPRE_MGRSetTol(par_solver, tol); + + HYPRE_MGRSetGlobalSmoothType(par_solver, mgr_gsmooth_type); + HYPRE_MGRSetMaxGlobalSmoothIters( par_solver, mgr_num_gsmooth_sweeps ); + + HYPRE_MGRSetCoarseGridMethod(par_solver, &mgr_use_non_galerkin); + HYPRE_MGRSetPMaxElmts(par_solver, mgr_pmax); // no truncation + //HYPRE_Int num_functions = 1; + //HYPRE_MGRSetLevelFRelaxNumFunctions(par_solver, &num_functions); + + /* Create the coarse grid solver */ + + HYPRE_Solver coarse_solver; + HYPRE_BoomerAMGCreate(&coarse_solver); + HYPRE_BoomerAMGSetCGCIts(coarse_solver, 1); + HYPRE_BoomerAMGSetInterpType(coarse_solver, 0); + HYPRE_BoomerAMGSetPostInterpType(coarse_solver, 0); + HYPRE_BoomerAMGSetCoarsenType(coarse_solver, 6); + HYPRE_BoomerAMGSetPMaxElmts(coarse_solver, 0); // no truncation + HYPRE_BoomerAMGSetCycleType(coarse_solver, 1); + HYPRE_BoomerAMGSetFCycle(coarse_solver, 0); + HYPRE_BoomerAMGSetNumSweeps(coarse_solver, 1); + HYPRE_BoomerAMGSetRelaxType(coarse_solver, 3); + HYPRE_BoomerAMGSetRelaxOrder(coarse_solver, 1); + HYPRE_BoomerAMGSetMaxLevels(coarse_solver, 25); + HYPRE_BoomerAMGSetStrongThreshold(coarse_solver, mgr_csolve_threshold); + HYPRE_BoomerAMGSetMaxIter(coarse_solver, mgr_csolve_max_iter); + HYPRE_BoomerAMGSetTol(coarse_solver, 1.e-12 ); + HYPRE_BoomerAMGSetPrintLevel(coarse_solver, 3); + // HYPRE_BoomerAMGSetUseAuxStrengthMatrix(coarse_solver, strength_method); + // HYPRE_BoomerAMGSetMinCoarseSize(coarse_solver, 100); + // HYPRE_BoomerAMGSetMaxCoarseSize(coarse_solver, 500); + /* set the MGR coarse solver. Comment out to use default Coarse Grid solver in MGR */ + HYPRE_MGRSetCoarseSolver( par_solver, HYPRE_BoomerAMGSolve, HYPRE_BoomerAMGSetup, coarse_solver); + + hypre_MGRSetFrelaxPrintLevel(par_solver, 3); + + /* setup MGR solver */ + HYPRE_MGRSetup(par_solver, par_A, par_b, par_x); + + if (print_cgrid) + { + HYPRE_Int* cgrid = (HYPRE_Int*) calloc(nx * ny, sizeof(HYPRE_Int)); + HYPRE_BoomerAMGGetGridHierarchy (coarse_solver, cgrid); + + FILE* fd = fopen ("cgrid", "w"); + + fprintf(fd, "%d %d\n", nx, ny); + for (i = 0; i < nx * ny; ++i) + { + fprintf(fd, "%d\n", cgrid[i]); + } + + fclose(fd); + free(cgrid); + } + /* MGR solve */ + HYPRE_MGRSolve(par_solver, par_A, par_b, par_x); + + HYPRE_MGRGetNumIterations(par_solver, &its); + HYPRE_MGRGetFinalRelativeResidualNorm(par_solver, &final_res_norm); + + /* free memory */ + if (mgr_num_cindexes) + { + hypre_TFree(mgr_num_cindexes, HYPRE_MEMORY_HOST); + } + mgr_num_cindexes = NULL; + + if (mgr_cindexes) + { + for ( i = 0; i < mgr_nlevels; i++) + { + if (mgr_cindexes[i]) + { + hypre_TFree(mgr_cindexes[i], HYPRE_MEMORY_HOST); + } + } + hypre_TFree(mgr_cindexes, HYPRE_MEMORY_HOST); + mgr_cindexes = NULL; + } + + HYPRE_BoomerAMGDestroy(coarse_solver); + HYPRE_MGRDestroy(par_solver); + } + else if (solver_id == 3) /* semicoarsening MGR (sMGR) */ + { + + HYPRE_MGRCreate(&par_solver); + + mgr_num_cindexes = hypre_CTAlloc(HYPRE_Int, mgr_nlevels, HYPRE_MEMORY_HOST); + for (i = 0; i < mgr_nlevels; i++) + { + /* assume 1 coarse index per level */ + mgr_num_cindexes[i] = 1; + } + mgr_cindexes = hypre_CTAlloc(HYPRE_Int*, mgr_nlevels, HYPRE_MEMORY_HOST); + for (i = 0; i < mgr_nlevels; i++) + { + mgr_cindexes[i] = hypre_CTAlloc(HYPRE_Int, mgr_num_cindexes[i], HYPRE_MEMORY_HOST); + } + for (i = 0; i < mgr_nlevels; i++) + { + /* assume coarse point is at index 0 */ + mgr_cindexes[i][0] = mgr_cpoint; + } + + /* set MGR data by marker array */ + int var, n; + HYPRE_Int* CF_index = (HYPRE_Int*) calloc(nx * ny * 2, sizeof(HYPRE_Int)); + HYPRE_Int* ptr = CF_index; + for (var = 0; var < 2; ++var) + { + for (n = 0; n < nx * ny; ++n) { *ptr++ = var; } + } + HYPRE_MGRSetCpointsByPointMarkerArray(par_solver, mgr_bsize, mgr_nlevels, mgr_num_cindexes, + mgr_cindexes, CF_index); + + /* set intermediate coarse grid strategy */ + HYPRE_MGRSetNonCpointsToFpoints(par_solver, mgr_non_c_to_f); + /* set F relaxation strategy */ + HYPRE_MGRSetFRelaxMethod(par_solver, 0); + /* set relax type for single level F-relaxation and post-relaxation */ + HYPRE_MGRSetRelaxType(par_solver, 0); + HYPRE_MGRSetNumRelaxSweeps(par_solver, 1); + /* set interpolation type */ + HYPRE_MGRSetRestrictType(par_solver, 0); + HYPRE_MGRSetInterpType(par_solver, 4); + HYPRE_MGRSetNumInterpSweeps(par_solver, mgr_num_interp_sweeps); + /* set print level */ + HYPRE_MGRSetPrintLevel(par_solver, print_level); + /* set max iterations */ + HYPRE_MGRSetMaxIter(par_solver, 1); + HYPRE_MGRSetTol(par_solver, 0.); + + hypre_MGRPrintCoarseSystem(par_solver, 1); + + + /* Create the inner MGR coarse grid solver */ + HYPRE_Solver mgr_coarse_solver; + HYPRE_MGRCreate(&mgr_coarse_solver); + HYPRE_MGRSetPrintLevel(mgr_coarse_solver, print_level); + hypre_MGRSetFrelaxPrintLevel(mgr_coarse_solver, 3); +#if 1 + /* set coarsening data for inner mgr solver */ + { + HYPRE_Int smgr_nblocks = semicoarsening_dir == 0 ? nx : ny; + HYPRE_Int smgr_nlevels = ceil(log2(smgr_nblocks)); + printf("smgr_nblocks = %d, smgr_nlevels = %d \n", smgr_nblocks, smgr_nlevels); + + /* Assume C-labels come first. So C F C F C F C ... */ + HYPRE_Int *lvl_num_cpts = (int*)calloc(smgr_nlevels, sizeof(int)); + HYPRE_Int **lvl_cpts = (int**)malloc(smgr_nlevels * sizeof(int*)); + HYPRE_Int num_c_pts = smgr_nblocks; + HYPRE_Int skip = 1; + for (i = 0; i < smgr_nlevels; i++) + { + num_c_pts = (num_c_pts + 1) / 2; // ensure we have ceiling of num_c_pts / 2 + lvl_cpts[i] = (int*)calloc(num_c_pts, sizeof(int)); + lvl_num_cpts[i] = num_c_pts; + printf("lvl_num_cpts[%d] = %d \n", i, num_c_pts); + skip *= 2; + int jc = 0; + int counter = 0; + /* populate coarse point labels for each level */ + while (counter < smgr_nblocks) + { + lvl_cpts[i][jc] = counter; + counter += skip; + jc++; + } + } + /* now generate marker labels for domain points */ + HYPRE_Int* smgr_CF_index = (HYPRE_Int*) calloc(nx * ny, sizeof(HYPRE_Int)); + ptr = smgr_CF_index; + HYPRE_Int ctr = 0; + for (j = 0; j < ny; j++) + { + for (i = 0; i < nx; i++) + { + *ptr++ = semicoarsening_dir == 0 ? i : j; + // printf("smgr_CF_index[%d] = %d \n", i+j*nx, smgr_CF_index[ctr++]); + } + } + + /* + for(i=0; i : set AMG cycles (1=V, 2=W, etc.)\n"); hypre_printf(" -cutf : set coarsening cut factor for dense rows\n"); hypre_printf(" -th : set AMG threshold Theta = val \n"); + /* RDF: Update help message for aux strength */ + hypre_printf(" -use_aux_str : use aux strength method (= 0,1,2,10,11,12) \n"); hypre_printf(" -tr : set AMG interpolation truncation factor = val \n"); hypre_printf(" -Pmx : set maximal no. of elmts per row for AMG interpolation (default: 4)\n"); hypre_printf(" -jtr : set truncation threshold for Jacobi interpolation = val \n"); @@ -4496,6 +4504,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetConvergeType(amg_solver, converge_type); HYPRE_BoomerAMGSetTol(amg_solver, tol); HYPRE_BoomerAMGSetStrongThreshold(amg_solver, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(amg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(amg_solver, seq_threshold); HYPRE_BoomerAMGSetRedundant(amg_solver, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(amg_solver, coarse_threshold); @@ -4825,6 +4834,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetMeasureType(amg_solver, measure_type); HYPRE_BoomerAMGSetTol(amg_solver, tol); HYPRE_BoomerAMGSetStrongThreshold(amg_solver, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(amg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(amg_solver, seq_threshold); HYPRE_BoomerAMGSetRedundant(amg_solver, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(amg_solver, coarse_threshold); @@ -5042,6 +5052,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetFPoints(pcg_precond, num_fpt, fpt_index); HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(pcg_precond, seq_threshold); HYPRE_BoomerAMGSetRedundant(pcg_precond, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(pcg_precond, coarse_threshold); @@ -5244,6 +5255,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetFPoints(pcg_precond, num_fpt, fpt_index); HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(pcg_precond, seq_threshold); HYPRE_BoomerAMGSetRedundant(pcg_precond, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(pcg_precond, coarse_threshold); @@ -5864,6 +5876,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetTruncFactor(pcg_precond, trunc_factor); HYPRE_BoomerAMGSetPMaxElmts(pcg_precond, P_max_elmts); HYPRE_BoomerAMGSetPostInterpType(pcg_precond, post_interp_type); @@ -5982,6 +5995,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetTruncFactor(pcg_precond, trunc_factor); HYPRE_BoomerAMGSetPMaxElmts(pcg_precond, P_max_elmts); HYPRE_BoomerAMGSetPostInterpType(pcg_precond, post_interp_type); @@ -6292,6 +6306,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetTruncFactor(pcg_precond, trunc_factor); HYPRE_BoomerAMGSetPMaxElmts(pcg_precond, P_max_elmts); HYPRE_BoomerAMGSetPostInterpType(pcg_precond, post_interp_type); @@ -6421,6 +6436,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetTruncFactor(pcg_precond, trunc_factor); HYPRE_BoomerAMGSetPMaxElmts(pcg_precond, P_max_elmts); HYPRE_BoomerAMGSetPostInterpType(pcg_precond, post_interp_type); @@ -6773,6 +6789,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(amg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(amg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(amg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(amg_precond, seq_threshold); HYPRE_BoomerAMGSetRedundant(amg_precond, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(amg_precond, coarse_threshold); @@ -6987,6 +7004,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(pcg_precond, seq_threshold); HYPRE_BoomerAMGSetRedundant(pcg_precond, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(pcg_precond, coarse_threshold); @@ -7368,6 +7386,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(pcg_precond, seq_threshold); HYPRE_BoomerAMGSetRedundant(pcg_precond, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(pcg_precond, coarse_threshold); @@ -7594,6 +7613,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(pcg_precond, seq_threshold); HYPRE_BoomerAMGSetRedundant(pcg_precond, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(pcg_precond, coarse_threshold); @@ -8029,6 +8049,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(pcg_precond, seq_threshold); HYPRE_BoomerAMGSetRedundant(pcg_precond, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(pcg_precond, coarse_threshold); @@ -8456,6 +8477,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(pcg_precond, seq_threshold); HYPRE_BoomerAMGSetRedundant(pcg_precond, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(pcg_precond, coarse_threshold); @@ -8832,6 +8854,7 @@ main( hypre_int argc, HYPRE_BoomerAMGSetIsolatedFPoints(pcg_precond, num_isolated_fpt, isolated_fpt_index); HYPRE_BoomerAMGSetMeasureType(pcg_precond, measure_type); HYPRE_BoomerAMGSetStrongThreshold(pcg_precond, strong_threshold); + HYPRE_BoomerAMGSetUseAuxStrengthMatrix(pcg_precond, use_aux_strength_matrix); HYPRE_BoomerAMGSetSeqThreshold(pcg_precond, seq_threshold); HYPRE_BoomerAMGSetRedundant(pcg_precond, redundant); HYPRE_BoomerAMGSetMaxCoarseSize(pcg_precond, coarse_threshold); diff --git a/src/utilities/_hypre_utilities.h b/src/utilities/_hypre_utilities.h index ccdce50c7f..f9ffc64c1f 100644 --- a/src/utilities/_hypre_utilities.h +++ b/src/utilities/_hypre_utilities.h @@ -3946,3 +3946,4 @@ HYPRE_Int hypre_mm_read_mtx_crd_size(FILE *f, HYPRE_Int *M, HYPRE_Int *N, HYPRE_ #endif #endif +