From 3cac332e990a5403011d7eb69f91b6e704e59479 Mon Sep 17 00:00:00 2001 From: Weston Ortiz Date: Fri, 5 Apr 2024 09:48:10 -0600 Subject: [PATCH] Add Amesos2, Tpetra matrix format, update legacy build script (#456) - Adds Tpetra Matrix format - Refactors Epetra to use generic form - Adds Amesos2 solver - Adds Tpetra based Stratimikos solver --- .github/workflows/build-unit.yml | 5 + CMakeLists.txt | 69 +- .../matrix_storage_format.rst | 3 + .../solution_algorithm.rst | 16 + include/dp_comm.h | 10 + include/linalg/sparse_matrix.h | 117 ++ include/linalg/sparse_matrix_epetra.h | 65 + include/linalg/sparse_matrix_tpetra.h | 64 + include/linalg/vector.h | 0 include/mm_mp_const.h | 2 +- include/rf_solver.h | 2 + include/rf_solver_const.h | 2 +- include/sl_amesos2_interface.h | 21 + include/sl_epetra_interface.h | 49 - include/sl_epetra_util.h | 27 - include/sl_stratimikos_interface.h | 7 + include/sl_util_structs.h | 34 +- scripts/build-goma-dependencies.sh | 1645 +---------------- src/ac_conti.c | 18 +- src/ac_hunt.c | 26 +- src/adapt/resetup_problem.c | 20 +- src/bc/rotate.c | 30 +- src/dp_comm.c | 112 ++ src/dp_ghost.cpp | 1 - src/dp_vif.c | 2 + src/globals.c | 4 + .../sparse_matrix.cpp} | 275 ++- src/linalg/sparse_matrix_epetra.cpp | 180 ++ src/linalg/sparse_matrix_tpetra.cpp | 213 +++ src/metis_decomp.c | 4 +- src/mm_fill.c | 7 +- src/mm_fill_ls.c | 8 +- src/mm_fill_species.c | 6 +- src/mm_fill_stress.c | 4 +- src/mm_input.c | 45 +- src/mm_input_bc.c | 2 +- src/mm_sol_nonlinear.c | 47 +- src/rf_solve.c | 36 +- src/rf_solve_segregated.c | 26 +- src/sl_amesos2_interface.cpp | 132 ++ src/sl_amesos_interface.cpp | 15 +- src/sl_aztecoo_interface.cpp | 6 +- src/sl_epetra_interface.cpp | 121 -- src/sl_matrix_util.c | 21 +- src/sl_stratimikos_interface.cpp | 270 ++- src/sl_util.c | 8 + 46 files changed, 1703 insertions(+), 2074 deletions(-) create mode 100644 include/linalg/sparse_matrix.h create mode 100644 include/linalg/sparse_matrix_epetra.h create mode 100644 include/linalg/sparse_matrix_tpetra.h create mode 100644 include/linalg/vector.h create mode 100644 include/sl_amesos2_interface.h delete mode 100644 include/sl_epetra_interface.h delete mode 100644 include/sl_epetra_util.h rename src/{sl_epetra_util.cpp => linalg/sparse_matrix.cpp} (70%) create mode 100644 src/linalg/sparse_matrix_epetra.cpp create mode 100644 src/linalg/sparse_matrix_tpetra.cpp create mode 100644 src/sl_amesos2_interface.cpp delete mode 100644 src/sl_epetra_interface.cpp diff --git a/.github/workflows/build-unit.yml b/.github/workflows/build-unit.yml index 530f4e615..99f4b195f 100644 --- a/.github/workflows/build-unit.yml +++ b/.github/workflows/build-unit.yml @@ -12,10 +12,15 @@ jobs: runs-on: ubuntu-latest container: westonortiz/goma-libs + defaults: + run: + shell: bash -l -e {0} + steps: - uses: actions/checkout@v3 - name: Check versions run: | + echo $PETSC_DIR mpirun --version cmake --version - name: CMake Configure Goma diff --git a/CMakeLists.txt b/CMakeLists.txt index a185893aa..77c069389 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -64,6 +64,11 @@ else() endif() endif() +option(GOMA_MATRIX_GO_LONG_LONG OFF) +if(GOMA_MATRIX_GO_LONG_LONG) + list(APPEND GOMA_COMPILE_DEFINITIONS GOMA_MATRIX_GO_LONG_LONG) +endif() + option(ENABLE_CCACHE OFF) if(ENABLE_CCACHE) message(STATUS "Searching for ccache.") @@ -108,6 +113,47 @@ if(${stratimikos_package_index} GREATER_EQUAL 0) endif() endif() +list(FIND Trilinos_PACKAGE_LIST Tpetra tpetra_package_index) +if(${tpetra_package_index} GREATER_EQUAL 0) + option(ENABLE_TPETRA "ENABLE_TPETRA" ON) + if((NOT ${HAVE_TPETRA_INST_INT_INT}) AND (NOT ${GOMA_MATRIX_GO_LONG_LONG})) + message( + WARNING + "Trilinos not built with Tpetra with int int, enable GOMA_MATRIX_GO_LONG_LONG (experimental), or rebuild Trilinos with Tpetra int int" + ) + set(ENABLE_TPETRA OFF) + endif() + if(ENABLE_TPETRA) + option(GOMA_MATRIX_GO_LONG_LONG "GOMA_MATRIX_GO_LONG_LONG" OFF) + if(GOMA_MATRIX_GO_LONG_LONG) + message( + WARNING + "GOMA_MATRIX_GO_LONG_LONG is enabled, using long long for Tpetra this is experimental" + ) + list(APPEND GOMA_COMPILE_DEFINITIONS GOMA_MATRIX_GO_LONG_LONG) + endif() + message(STATUS "TRILINOS: Tpetra found, enabling in Goma") + list(APPEND GOMA_COMPILE_DEFINITIONS GOMA_ENABLE_TPETRA) + list(FIND Trilinos_PACKAGE_LIST Amesos2 amesos2_package_index) + if(${amesos2_package_index} GREATER_EQUAL 0) + option(ENABLE_AMESOS2 "ENABLE_AMESOS2" ON) + if(ENABLE_AMESOS2) + message(STATUS "TRILINOS: Amesos2 found, enabling in Goma") + list(APPEND GOMA_COMPILE_DEFINITIONS GOMA_ENABLE_AMESOS2) + endif() + endif() + endif() +endif() + +list(FIND Trilinos_PACKAGE_LIST Epetra epetra_package_index) +if(${epetra_package_index} GREATER_EQUAL 0) + option(ENABLE_EPETRA "ENABLE_EPETRA" ON) + if(ENABLE_EPETRA) + message(STATUS "TRILINOS: Epetra found, enabling in Goma") + list(APPEND GOMA_COMPILE_DEFINITIONS GOMA_ENABLE_EPETRA) + endif() +endif() + # Exodus find_package(SEACASExodus REQUIRED HINTS ${Trilinos_DIR}/../SEACASExodus) message(STATUS "SEACASExodus_DIR = ${SEACASExodus_DIR}") @@ -439,6 +485,9 @@ set(GOMA_INCLUDES include/el_quality.h include/exo_conn.h include/exo_struct.h + include/linalg/sparse_matrix.h + include/linalg/sparse_matrix_tpetra.h + include/linalg/sparse_matrix_epetra.h include/load_field_variables.h include/loca_const.h include/loca_util_const.h @@ -536,13 +585,12 @@ set(GOMA_INCLUDES include/shell_tfmp_struct.h include/shell_tfmp_util.h include/sl_amesos_interface.h + include/sl_amesos2_interface.h include/sl_aux.h include/sl_auxutil.h include/sl_aztecoo_interface.h include/sl_eggroll_def.h include/sl_eggroll.h - include/sl_epetra_interface.h - include/sl_epetra_util.h include/sl_lu.h include/sl_petsc.h include/sl_petsc_complex.h @@ -591,6 +639,9 @@ set(GOMA_SOURCES src/el_quality.c src/exo_conn.c src/globals.c + src/linalg/sparse_matrix.cpp + src/linalg/sparse_matrix_tpetra.cpp + src/linalg/sparse_matrix_epetra.cpp src/load_field_variables.c src/loca_bord.c src/loca_eigen_c2f.F @@ -684,6 +735,7 @@ set(GOMA_SOURCES src/bc/rotate_coordinates.c src/shell_tfmp_util.c src/sl_amesos_interface.cpp + src/sl_amesos2_interface.cpp src/sl_aux.c src/sl_auxutil.c src/sl_aztecoo_interface.cpp @@ -694,8 +746,6 @@ set(GOMA_SOURCES src/sl_eggroll05.c src/sl_eggrollutil.c src/sl_eggrollwrap.c - src/sl_epetra_interface.cpp - src/sl_epetra_util.cpp src/sl_leastsquares.cpp src/sl_lu.c src/sl_lu_fill.c @@ -948,14 +998,11 @@ if(ENABLE_PETSC AND PETSC_FOUND) ${PETSC_STATIC_LIBRARY_DIRS}) endif() -if ($Trilinos_VERSION VERSION_LESS 14.0.0) - set(GOMA_TRILINOS_LIBRARIES - ${Trilinos_LIBRARIES} - ${Trilinos_TPL_LIBRARIES} - ${Trilinos_EXTRA_LD_FLAGS}) +if($Trilinos_VERSION VERSION_LESS 14.0.0) + set(GOMA_TRILINOS_LIBRARIES ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES} + ${Trilinos_EXTRA_LD_FLAGS}) else() - set(GOMA_TRILINOS_LIBRARIES - Trilinos::all_libs) + set(GOMA_TRILINOS_LIBRARIES Trilinos::all_libs) endif() set(GOMA_TPL_LIBRARIES diff --git a/docs/problem_description_file/solver_specifications/matrix_storage_format.rst b/docs/problem_description_file/solver_specifications/matrix_storage_format.rst index 91a0cd3b0..e9c65d1cd 100644 --- a/docs/problem_description_file/solver_specifications/matrix_storage_format.rst +++ b/docs/problem_description_file/solver_specifications/matrix_storage_format.rst @@ -21,6 +21,9 @@ vbr selected when an Aztec iterative solver is chosen. epetra Compressed Sparse Row format using the Epetra library from Trilinos +tpetra + FECRSMatrix Compressed Sparse Row format using the Tpetra library from Trilinos + currently only works using Stratimikos and Amesos2 ------------ Examples diff --git a/docs/problem_description_file/solver_specifications/solution_algorithm.rst b/docs/problem_description_file/solver_specifications/solution_algorithm.rst index f70967f13..96a6c9aad 100644 --- a/docs/problem_description_file/solver_specifications/solution_algorithm.rst +++ b/docs/problem_description_file/solver_specifications/solution_algorithm.rst @@ -98,6 +98,21 @@ amesos Of these four options, we currently recommend **mumps**. All options can be run in parallel. +amesos2 + Allows access to direct solver options implemented in parallel. Please see + the user-notes below for Goma build options that must be exercised. With + this option, you must add an additional input card to specify the parallel + direct solvers: + + :: + + Amesos2 Solver Package = {SuperLUDist | Mumps | KLU2} + # Optional Amesos2 File to set parameters for the solver + # xml or yaml supported through Trilinos + # Amesos2 File = params_file.(xml | yaml) + + Of these four options, we currently recommend **Mumps**. + All options can be run in parallel. stratimikos Interface to Trilino's Stratimikos package requires: @@ -105,6 +120,7 @@ stratimikos :: Matrix storage format = epetra + # Or tpetra Allows block solvers, see also ref:`Stratimikos File` petsc diff --git a/include/dp_comm.h b/include/dp_comm.h index fb7de5f47..201dacc94 100644 --- a/include/dp_comm.h +++ b/include/dp_comm.h @@ -35,6 +35,16 @@ EXTERN void exchange_dof(Comm_Ex *, /* cx - ptr to communications exchange info double *, /* x - local processor dof-based vector */ int); +EXTERN void exchange_dof_int(Comm_Ex *, /* cx - ptr to communications exchange info */ + Dpi *, /* dpi - distributed processing info */ + int *, /* x - local processor dof-based vector */ + int); + +EXTERN void exchange_dof_long_long(Comm_Ex *, /* cx - ptr to communications exchange info */ + Dpi *, /* dpi - distributed processing info */ + long long *, /* x - local processor dof-based vector */ + int); + EXTERN void exchange_node(Comm_Ex *cx, /* cx - ptr to communications exchange info */ Dpi *d, /* dpi - distributed processing info */ double *a); /* x - local processor node-based vector */ diff --git a/include/linalg/sparse_matrix.h b/include/linalg/sparse_matrix.h new file mode 100644 index 000000000..aa66558d8 --- /dev/null +++ b/include/linalg/sparse_matrix.h @@ -0,0 +1,117 @@ +#ifndef GOMA_SPARSE_MATRIX +#define GOMA_SPARSE_MATRIX + +#ifdef __cplusplus +extern "C" { +#endif +#ifdef GOMA_MATRIX_GO_LONG_LONG +typedef long long GomaGlobalOrdinal; +#define MPI_GOMA_ORDINAL MPI_LONG_LONG +#else +typedef int GomaGlobalOrdinal; +#define MPI_GOMA_ORDINAL MPI_INT +#endif + +#ifdef DISABLE_CPP +#define CPP_ALREADY_DISABLED +#endif +#define DISABLE_CPP +#include "dp_types.h" +#include "exo_struct.h" +#include "mm_eh.h" +#include "rf_node_const.h" +#include "sl_util_structs.h" +#ifndef CPP_ALREADY_DISABLED +#undef DISABLE_CPP +#endif + +enum GomaSparseMatrixType { + GOMA_SPARSE_MATRIX_TYPE_EPETRA, + GOMA_SPARSE_MATRIX_TYPE_AZTEC_MRS, + GOMA_SPARSE_MATRIX_TYPE_TPETRA, + GOMA_SPARSE_MATRIX_TYPE_PETSC +}; + +struct g_GomaSparseMatrix { + enum GomaSparseMatrixType type; + void *data; + GomaGlobalOrdinal *global_ids; + GomaGlobalOrdinal n_rows; + GomaGlobalOrdinal n_cols; + GomaGlobalOrdinal nnz; + // Create matrix with given rows and columns + // coo_rows and coo_cols are the COO representation of the matrix ordered by row + // local_nnz is the number of non-zero entries in the local partition + // max_nz_per_row is the maximum number of non-zero entries in any row + goma_error (*create_graph)(struct g_GomaSparseMatrix *matrix, + GomaGlobalOrdinal n_rows, + GomaGlobalOrdinal *row_list, + GomaGlobalOrdinal n_cols, + GomaGlobalOrdinal *col_list, + GomaGlobalOrdinal local_nnz, + GomaGlobalOrdinal max_nz_per_row, + GomaGlobalOrdinal *coo_rows, + GomaGlobalOrdinal *coo_cols); + // optional, run after create graph to finalize structure + goma_error (*complete_graph)(struct g_GomaSparseMatrix *matrix); + // Insert values into matrix row replacing existing values + goma_error (*insert_row_values)(struct g_GomaSparseMatrix *matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices); + // Add values into a matrix row, sums existing values + goma_error (*sum_into_row_values)(struct g_GomaSparseMatrix *matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices); + // set matrix non-zeros to specified scalar value (commonly re-zero for next assembly) + goma_error (*put_scalar)(struct g_GomaSparseMatrix *matrix, double scalar); + // row sum scaling, compute row sum scale, scale matrix and b, and return scaling vector + goma_error (*row_sum_scaling)(struct g_GomaSparseMatrix *matrix, double *b, double *scale); + // Zeros a global row + goma_error (*zero_global_row)(struct g_GomaSparseMatrix *matrix, GomaGlobalOrdinal global_row); + // Zeros a global row and sets diagonal to 1.0 + goma_error (*zero_global_row_set_diag)(struct g_GomaSparseMatrix *matrix, + GomaGlobalOrdinal global_row); + // delete the allocated matrix; + goma_error (*destroy)(struct g_GomaSparseMatrix *matrix); +}; +typedef struct g_GomaSparseMatrix *GomaSparseMatrix; + +goma_error GomaSparseMatrix_CreateFromFormat(GomaSparseMatrix *matrix, char *matrix_format); +goma_error GomaSparseMatrix_Create(GomaSparseMatrix *matrix, enum GomaSparseMatrixType type); +// populates: +// - global_ids +// - n_rows +// - n_cols +goma_error GomaSparseMatrix_SetProblemGraph( + GomaSparseMatrix matrix, + int num_internal_dofs, + int num_boundary_dofs, + int num_external_dofs, + int local_nodes, + NODE_INFO_STRUCT **Nodes, + int MaxVarPerNode, + int *Matilda, + int Inter_Mask[MAX_NUM_MATRICES][MAX_VARIABLE_TYPES][MAX_VARIABLE_TYPES], + Exo_DB *exo, + Dpi *dpi, + Comm_Ex *cx, + int imtrx, + int Debug_Flag, + struct GomaLinearSolverData *ams); + +goma_error GomaSparseMatrix_LoadLec(GomaSparseMatrix matrix, + int ielem, + struct Local_Element_Contributions *lec, + double resid_vector[]); + +goma_error GomaSparseMatrix_Destroy(GomaSparseMatrix *matrix); + +#ifdef __cplusplus +} +#endif + +#endif // GOMA_SPARSE_MATRIX \ No newline at end of file diff --git a/include/linalg/sparse_matrix_epetra.h b/include/linalg/sparse_matrix_epetra.h new file mode 100644 index 000000000..6c68fa93e --- /dev/null +++ b/include/linalg/sparse_matrix_epetra.h @@ -0,0 +1,65 @@ +#ifndef GOMA_SPARSE_MATRIX_EPETRA +#define GOMA_SPARSE_MATRIX_EPETRA +#ifdef GOMA_ENABLE_EPETRA +#ifdef __cplusplus +#include "Teuchos_RCP.hpp" +#include +#include +#include + +#include "linalg/sparse_matrix.h" + +using LO = int; +using GO = GomaGlobalOrdinal; + +struct EpetraSparseMatrix { + Teuchos::RCP matrix; + Teuchos::RCP row_map; + Teuchos::RCP crs_graph; + EpetraSparseMatrix() = default; +}; + +extern "C" { +#endif + +goma_error GomaSparseMatrix_Epetra_Create(GomaSparseMatrix *matrix); + +goma_error g_epetra_create_graph(GomaSparseMatrix matrix, + GomaGlobalOrdinal n_rows, + GomaGlobalOrdinal *row_list, + GomaGlobalOrdinal n_cols, + GomaGlobalOrdinal *col_list, + GomaGlobalOrdinal local_nnz, + GomaGlobalOrdinal max_per_row, + GomaGlobalOrdinal *coo_rows, + GomaGlobalOrdinal *coo_cols); + +goma_error g_epetra_complete_graph(GomaSparseMatrix matrix); + +goma_error g_epetra_insert_row_values(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices); + +goma_error g_epetra_sum_into_row_values(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices); + +goma_error g_epetra_put_scalar(GomaSparseMatrix matrix, double scalar); + +goma_error g_epetra_row_sum_scaling(GomaSparseMatrix matrix, double *b, double *scale); + +goma_error g_epetra_zero_row(GomaSparseMatrix matrix, GomaGlobalOrdinal global_row); + +goma_error g_epetra_zero_row_set_diag(GomaSparseMatrix matrix, GomaGlobalOrdinal global_row); + +goma_error g_epetra_destroy(GomaSparseMatrix matrix); + +#ifdef __cplusplus +} +#endif +#endif +#endif // GOMA_SPARSE_MATRIX_EPETRA \ No newline at end of file diff --git a/include/linalg/sparse_matrix_tpetra.h b/include/linalg/sparse_matrix_tpetra.h new file mode 100644 index 000000000..d5f4319a4 --- /dev/null +++ b/include/linalg/sparse_matrix_tpetra.h @@ -0,0 +1,64 @@ +#ifndef GOMA_SPARSE_MATRIX_TPETRA +#define GOMA_SPARSE_MATRIX_TPETRA +#ifdef GOMA_ENABLE_TPETRA +#ifdef __cplusplus +#include "Teuchos_RCP.hpp" +#include + +#include "linalg/sparse_matrix.h" + +using LO = int; +using GO = GomaGlobalOrdinal; + +struct TpetraSparseMatrix { + Teuchos::RCP> matrix; + Teuchos::RCP> row_map; + Teuchos::RCP> col_map; + Teuchos::RCP> crs_graph; + TpetraSparseMatrix() = default; +}; + +extern "C" { +#endif + +goma_error GomaSparseMatrix_Tpetra_Create(GomaSparseMatrix *matrix); + +goma_error g_tpetra_create_graph(GomaSparseMatrix matrix, + GomaGlobalOrdinal n_rows, + GomaGlobalOrdinal *row_list, + GomaGlobalOrdinal n_cols, + GomaGlobalOrdinal *col_list, + GomaGlobalOrdinal local_nnz, + GomaGlobalOrdinal max_per_row, + GomaGlobalOrdinal *coo_rows, + GomaGlobalOrdinal *coo_cols); + +goma_error g_tpetra_complete_graph(GomaSparseMatrix matrix); + +goma_error g_tpetra_insert_row_values(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices); + +goma_error g_tpetra_sum_into_row_values(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices); + +goma_error g_tpetra_put_scalar(GomaSparseMatrix matrix, double scalar); + +goma_error g_tpetra_row_sum_scaling(GomaSparseMatrix matrix, double *b, double *scale); + +goma_error g_tpetra_zero_row(GomaSparseMatrix matrix, GomaGlobalOrdinal global_row); + +goma_error g_tpetra_zero_row_set_diag(GomaSparseMatrix matrix, GomaGlobalOrdinal global_row); + +goma_error g_tpetra_destroy(GomaSparseMatrix matrix); + +#ifdef __cplusplus +} +#endif +#endif +#endif // GOMA_SPARSE_MATRIX_TPETRA \ No newline at end of file diff --git a/include/linalg/vector.h b/include/linalg/vector.h new file mode 100644 index 000000000..e69de29bb diff --git a/include/mm_mp_const.h b/include/mm_mp_const.h index efb05089b..3daf2486b 100644 --- a/include/mm_mp_const.h +++ b/include/mm_mp_const.h @@ -564,7 +564,7 @@ extern int Num_Var_Init_Mat[MAX_NUMBER_MATLS]; /* number of variables to overwri #define SOLID_DIFFUSION_ELECTRONEUTRALITY 911 #define SOLID_DIFFUSION 912 #define GAS_DIFFUSION 913 -#define FULL 914 +#define METAL_CORROSION_FULL 914 #define ANNIHILATION_ELECTRONEUTRALITY 915 #define ANNIHILATION 916 #define NET_CHARGE 917 /* refer to F multiplied by sum of ci zi */ diff --git a/include/rf_solver.h b/include/rf_solver.h index 068c3833b..66ecc7728 100644 --- a/include/rf_solver.h +++ b/include/rf_solver.h @@ -83,10 +83,12 @@ extern String_line Matrix_Relative_Threshold; /* Trilinos 2 */ extern String_line Matrix_Absolute_Threshold; /* Trilinos 2 */ extern String_line Amesos_Package; +extern String_line Amesos2_Package; extern String_line AztecOO_Solver; extern String_line Stratimikos_File[MAX_NUM_MATRICES]; +extern String_line Amesos2_File[MAX_NUM_MATRICES]; /* extern * A new Aztec 2.0 option. There are more and difft options and our diff --git a/include/rf_solver_const.h b/include/rf_solver_const.h index 9e7e48a09..6d3d3aa33 100644 --- a/include/rf_solver_const.h +++ b/include/rf_solver_const.h @@ -61,7 +61,7 @@ #define STRATIMIKOS 10 #define PETSC_SOLVER 11 #define PETSC_COMPLEX_SOLVER 12 - +#define AMESOS2 13 /* * FORTRAN BLAS functions. Inside C, use "DCOPY" and the preprocessor to * make it look like the FORTRAN name for this routine. diff --git a/include/sl_amesos2_interface.h b/include/sl_amesos2_interface.h new file mode 100644 index 000000000..cddf7cb3c --- /dev/null +++ b/include/sl_amesos2_interface.h @@ -0,0 +1,21 @@ +#ifndef GOMA_SL_AMESOS2_INTERFACE_H_ +#define GOMA_SL_AMESOS2_INTERFACE_H_ + +#ifdef __cplusplus +extern "C" { +#endif + +#include "rf_fem_const.h" +#include "rf_io_const.h" + +int amesos2_solve(struct GomaLinearSolverData *ams, + double *x_, + double *b_, + char *amesos2_solver, + char *amesos2_file); + +#ifdef __cplusplus +} // end of extern "C" +#endif + +#endif /* GOMA_SL_AMESOS2_INTERFACE_H_ */ diff --git a/include/sl_epetra_interface.h b/include/sl_epetra_interface.h deleted file mode 100644 index a0b9e107c..000000000 --- a/include/sl_epetra_interface.h +++ /dev/null @@ -1,49 +0,0 @@ -/* - * sl_epetra_interface.h - * - * Created on: Oct 8, 2014 - * Author: wortiz - */ - -#ifndef INCLUDE_SL_EPETRA_INTERFACE_H_ -#define INCLUDE_SL_EPETRA_INTERFACE_H_ - -#if defined(__cplusplus) && !defined(DISABLE_CPP) -#include "Epetra_RowMatrix.h" -typedef Epetra_RowMatrix C_Epetra_RowMatrix_t; -#else -typedef struct Epetra_RowMatrix C_Epetra_RowMatrix_t; -#endif - -#if defined(__cplusplus) && !defined(DISABLE_CPP) -extern "C" { -#endif - -C_Epetra_RowMatrix_t *EpetraCreateRowMatrix(int NumberProcRows); - -void EpetraInsertGlobalRowMatrix(C_Epetra_RowMatrix_t *AMatrix, - int GlobalRow, - int NumEntries, - const double *Values, - const int *Indices); - -void EpetraSumIntoGlobalRowMatrix(C_Epetra_RowMatrix_t *AMatrix, - int GlobalRow, - int NumEntries, - const double *Values, - const int *Indices); - -void EpetraPutScalarRowMatrix(C_Epetra_RowMatrix_t *AMatrix, double scalar); - -void EpetraFillCompleteRowMatrix(C_Epetra_RowMatrix_t *AMatrix); - -int EpetraExtractRowValuesRowMatrix( - C_Epetra_RowMatrix_t *AMatrix, int GlobalRow, int size, double *values, int *indices); - -void EpetraDeleteRowMatrix(C_Epetra_RowMatrix_t *AMatrix); - -#if defined(__cplusplus) && !defined(DISABLE_CPP) -} // end of extern "C" -#endif - -#endif /* INCLUDE_SL_EPETRA_INTERFACE_H_ */ diff --git a/include/sl_epetra_util.h b/include/sl_epetra_util.h deleted file mode 100644 index af43cd827..000000000 --- a/include/sl_epetra_util.h +++ /dev/null @@ -1,27 +0,0 @@ -/* - * sl_epetra_util.h - * - */ - -#ifndef INCLUDE_SL_EPETRA_UTIL_H_ -#define INCLUDE_SL_EPETRA_UTIL_H_ - -#include "dpi.h" -#include "exo_struct.h" -#ifdef __cplusplus -extern "C" { -#endif - -void EpetraCreateGomaProblemGraph(struct GomaLinearSolverData *ams, Exo_DB *exo, Dpi *dpi); - -void EpetraLoadLec(int ielem, struct GomaLinearSolverData *ams, double resid_vector[]); - -void EpetraRowSumScale(struct GomaLinearSolverData *ams, double *b, double *scale); - -void EpetraSetDiagonalOnly(struct GomaLinearSolverData *ams, int GlobalRow); - -#ifdef __cplusplus -} // end of extern "C" -#endif - -#endif /* INCLUDE_SL_EPETRA_UTIL_H_ */ diff --git a/include/sl_stratimikos_interface.h b/include/sl_stratimikos_interface.h index d39b877f3..f6e7be484 100644 --- a/include/sl_stratimikos_interface.h +++ b/include/sl_stratimikos_interface.h @@ -9,6 +9,13 @@ extern "C" { #include "rf_fem_const.h" #include "rf_io_const.h" +int stratimikos_solve_tpetra(struct GomaLinearSolverData *ams, + double *x_, + double *b_, + int *iterations, + char stratimikos_file[MAX_NUM_MATRICES][MAX_CHAR_IN_INPUT], + int imtrx); + int stratimikos_solve(struct GomaLinearSolverData *ams, double *x_, double *b_, diff --git a/include/sl_util_structs.h b/include/sl_util_structs.h index 38218901b..e5b8ed4fd 100644 --- a/include/sl_util_structs.h +++ b/include/sl_util_structs.h @@ -41,7 +41,6 @@ #endif #include "az_aztec.h" -#include "sl_epetra_interface.h" /* * NUM_ALSS - Number of Aztec Linear Solver Systems. Here we have one. @@ -51,20 +50,6 @@ #define JAC 0 #define NUM_ALSS 1 -struct Matrix_Data { - struct GomaLinearSolverData *ams; - double *x; /* Solution vector */ - double *x_old; /* Solution vector , previous last time step */ - double *x_older; /* Solution vector , previous prev time step */ - double *x_oldest; - double *xdot; /* xdot of current solution */ - double *xdot_old; /* xdot_old of current solution */ - double *xdot_older; - double *x_update; /* last update vector */ - double *resid_vector; /* Residual vector */ - double *scale; -}; - struct GomaLinearSolverData { int proc_config[AZ_PROC_SIZE]; @@ -124,10 +109,23 @@ struct GomaLinearSolverData { int solveSetup; - C_Epetra_RowMatrix_t *RowMatrix; /* This is a Epetra_RowMatrix object */ - int *GlobalIDs; /* Pointer to global ids of DOFs (only available with epetra) */ - void *PetscMatrixData; + void *GomaMatrixData; + void *SolverData; + void (*DestroySolverData)(struct GomaLinearSolverData *ams); }; +struct Matrix_Data { + struct GomaLinearSolverData *ams; + double *x; /* Solution vector */ + double *x_old; /* Solution vector , previous last time step */ + double *x_older; /* Solution vector , previous prev time step */ + double *x_oldest; + double *xdot; /* xdot of current solution */ + double *xdot_old; /* xdot_old of current solution */ + double *xdot_older; + double *x_update; /* last update vector */ + double *resid_vector; /* Residual vector */ + double *scale; +}; #endif diff --git a/scripts/build-goma-dependencies.sh b/scripts/build-goma-dependencies.sh index 7404b82e0..3b2a99034 100755 --- a/scripts/build-goma-dependencies.sh +++ b/scripts/build-goma-dependencies.sh @@ -239,30 +239,36 @@ NETCDF_VERSION="c-4.9.0" NETCDF_VERSION_ONLY="4.9.0" NETCDF_MD5="26cfd1b2e32d511bd82477803976365d" -TRILINOS_VERSION="14.4.0" -TRILINOS_VERSION_DASH="14-4-0" -TRILINOS_MD5="334f9c3700c72f6ed5658eaa783ffccd" +PNETCDF_VERSION="1.13.0" +PNETCDF_MD5="31b94d39462b1f1f2293f735c9819bf2" +PNETCDF_URL="https://parallel-netcdf.github.io/Release/pnetcdf-$PNETCDF_VERSION.tar.gz" -MUMPS_VERSION="5.5.1" +TRILINOS_VERSION="15.1.0" +TRILINOS_VERSION_DASH="15-1-0" +TRILINOS_MD5="79237697af4fc42eaaf70f23104a8e12" + +MUMPS_VERSION="5.6.2" MUMPS_URL="https://graal.ens-lyon.fr/MUMPS/MUMPS_$MUMPS_VERSION.tar.gz" -MUMPS_MD5="da26c4b43d53a9a6096775245cee847f" +MUMPS_MD5="adc38b6cca34dfcc8f7036b7d39cf260" -OPENMPI_VERSION="4.1.4" -OPENMPI_MD5="f057e12aabaf7dd5a6a658180fca404e" +OPENMPI_VERSION="4.1.6" +OPENMPI_MD5="c9b1c974cfc23c77c0fbdb965cd58a1c" OPENMPI_ARCHIVE_URL="https://download.open-mpi.org/release/open-mpi/v4.1/openmpi-$OPENMPI_VERSION.tar.bz2" OPENMPI_EXTRA_CONFIGURE_FLAGS="" CMAKE_VERSION="3.24.2" CMAKE_MD5="f6f79ec66c91bbc075757a70205128ca" +#SUITESPARSE_VERSION="7.7.0" +#SUITESPARSE_MD5="e659373ed5e9b961d2fcb6d67d250783" SUITESPARSE_VERSION="5.13.0" SUITESPARSE_MD5="e9e7bc594b77ae4b58d943cdc286d679" MATIO_VERSION="1.5.21" MATIO_MD5="afeb5d21b234699fd5b9dc4564afe1ca" -SCALAPACK_VERSION="2.2.0" -SCALAPACK_MD5="2397d36790d1445383bc3cdb1e18ca5f" +SCALAPACK_VERSION="0234af94c6578c53ac4c19f2925eb6e5c4ad6f0f" +SCALAPACK_MD5="221dab094055b73fa7bed645c426548c" LAPACK_VERSION="3.8.0" LAPACK_MD5="96591affdbf58c450d45c1daa540dbd2" @@ -294,12 +300,12 @@ PARMETIS_EXTRACT_NAME="petsc-pkg-parmetis-5777d7ec2084" ARCHIVE_NAMES=("arpack-ng-$ARPACK_NG_VERSION.tar.gz" \ "hdf5-${HDF5_VERSION}.tar.bz2" \ +"pnetcdf-$PNETCDF_VERSION.tar.gz" \ "netcdf-${NETCDF_VERSION}.tar.gz" \ "metis-$METIS_VERSION.tar.gz" \ "parmetis-$PARMETIS_VERSION.tar.gz" \ "sparse.tar.gz" \ "superlu_dist-$SUPERLU_DIST_VERSION.tar.gz" \ -"y12m.tar.gz" \ "Trilinos-trilinos-release-$TRILINOS_VERSION_DASH.tar.gz" \ "MUMPS_$MUMPS_VERSION.tar.gz" \ "SuiteSparse-$SUITESPARSE_VERSION.tar.gz" \ @@ -312,12 +318,12 @@ ARCHIVE_NAMES=("arpack-ng-$ARPACK_NG_VERSION.tar.gz" \ #meaning each y12m tar has a unique MD5SUM. ARCHIVE_MD5SUMS=("$ARPACK_NG_MD5" \ "${HDF5_MD5}" \ +"$PNETCDF_MD5" \ "${NETCDF_MD5}" \ "$METIS_MD5" \ "$PARMETIS_MD5" \ "1566d914d1035ac17b73fe9bc0eed02a" \ "$SUPERLU_DIST_MD5" \ -"SKIP" \ $TRILINOS_MD5 \ $MUMPS_MD5 \ "$SUITESPARSE_MD5" \ @@ -327,12 +333,12 @@ $MUMPS_MD5 \ ARCHIVE_URLS=("$ARPACK_NG_URL" \ "https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.12/hdf5-${HDF5_VERSION}/src/hdf5-${HDF5_VERSION}.tar.bz2" \ +"$PNETCDF_URL" \ "https://downloads.unidata.ucar.edu/netcdf-c/${NETCDF_VERSION_ONLY}/netcdf-${NETCDF_VERSION}.tar.gz" \ "https://bitbucket.org/petsc/pkg-metis/get/v$METIS_VERSION.tar.gz" \ "https://bitbucket.org/petsc/pkg-parmetis/get/v$PARMETIS_VERSION.tar.gz" \ "http://downloads.sourceforge.net/project/sparse/sparse/sparse1.4b/sparse1.4b.tar.gz" \ "http://codeload.github.com/xiaoyeli/superlu_dist/tar.gz/v$SUPERLU_DIST_VERSION" \ -"http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz\\&filename=y12m%2Fy12m.f" \ "https://github.com/trilinos/Trilinos/archive/trilinos-release-$TRILINOS_VERSION_DASH.tar.gz" \ "$MUMPS_URL" \ "https://github.com/DrTimothyAldenDavis/SuiteSparse/archive/refs/tags/v$SUITESPARSE_VERSION.tar.gz" \ @@ -344,12 +350,12 @@ ARCHIVE_URLS=("$ARPACK_NG_URL" \ # When in reality it isn't ARCHIVE_DIR_NAMES=("arpack-ng-${ARPACK_NG_VERSION}" \ "hdf5-${HDF5_VERSION}" \ +"pnetcdf-$PNETCDF_VERSION" \ "netcdf-${NETCDF_VERSION}" \ "metis-$METIS_VERSION" \ "parmetis-$PARMETIS_VERSION" \ "sparse" \ "superlu_dist-$SUPERLU_DIST_VERSION" \ -"y12m" \ "Trilinos-trilinos-release-$TRILINOS_VERSION_DASH" \ "MUMPS_$MUMPS_VERSION" \ "SuiteSparse-$SUITESPARSE_VERSION" \ @@ -359,12 +365,12 @@ ARCHIVE_DIR_NAMES=("arpack-ng-${ARPACK_NG_VERSION}" \ ARCHIVE_HOMEPAGES=("https://github.com/opencollab/arpack-ng/" \ "https://www.hdfgroup.org/" \ +"https://parallel-netcdf.github.io/" \ "https://www.unidata.ucar.edu/software/netcdf/" \ "http://glaros.dtc.umn.edu/gkhome/metis/metis/overview" \ "http://glaros.dtc.umn.edu/gkhome/metis/metis/overview" \ "http://sparse.sourceforge.net/" \ "https://github.com/xiaoyeli/superlu_dist" \ -"http://doi.org/10.1007/3-540-10874-2" \ "https://trilinos.org/" \ "http://mumps.enseeiht.fr/" \ "http://faculty.cse.tamu.edu/davis/suitesparse.html" \ @@ -374,12 +380,12 @@ ARCHIVE_HOMEPAGES=("https://github.com/opencollab/arpack-ng/" \ ARCHIVE_REAL_NAMES=("arpack-ng" \ "HDF5" \ +"PNetCDF" \ "NetCDF" \ "METIS" \ "ParMETIS" \ "Sparse" \ "SuperLU_DIST" \ -"y12m" \ "Trilinos" \ "MUMPS" \ "SuiteSparse-$SUITESPARSE_VERSION" \ @@ -572,9 +578,9 @@ function setMathVars { export SCALAPACK_LIBRARY_NAME="scalapack" export SCALAPACK_LIBRARY_NAME_ARG="-L${SCALAPACK_LIBRARY_DIR} -lscalapack" SCALAPACK_INCLUDE_DIR="${GOMA_LIB}/scalapack-$SCALAPACK_VERSION/include" - ARCHIVE_NAMES+=("scalapack-$SCALAPACK_VERSION.tgz") + ARCHIVE_NAMES+=("$SCALAPACK_VERSION.tar.gz") ARCHIVE_MD5SUMS+=("$SCALAPACK_MD5") - ARCHIVE_URLS+=("http://www.netlib.org/scalapack/scalapack-$SCALAPACK_VERSION.tgz") + ARCHIVE_URLS+=("https://github.com/Reference-ScaLAPACK/scalapack/archive/$SCALAPACK_VERSION.tar.gz") ARCHIVE_DIR_NAMES+=("scalapack-$SCALAPACK_VERSION") ARCHIVE_HOMEPAGES+=("http://www.netlib.org/scalapack/") ARCHIVE_REAL_NAMES+=("ScaLAPACK") @@ -593,11 +599,11 @@ function setMathVars { USING_MKLS="OFF" #Add packages that otherwise come preinstalled in intel compiler. ARCHIVE_NAMES+=("lapack-$LAPACK_VERSION.tar.gz" \ - "scalapack-$SCALAPACK_VERSION.tgz") + "$SCALAPACK_VERSION.tar.gz") ARCHIVE_MD5SUMS+=("$LAPACK_MD5" \ "$SCALAPACK_MD5" ) ARCHIVE_URLS+=("https://github.com/Reference-LAPACK/lapack/archive/refs/tags/v$LAPACK_VERSION.tar.gz" \ - "http://www.netlib.org/scalapack/scalapack-$SCALAPACK_VERSION.tgz" ) + "https://github.com/Reference-ScaLAPACK/scalapack/archive/$SCALAPACK_VERSION.tar.gz" ) ARCHIVE_DIR_NAMES+=("lapack-$LAPACK_VERSION" \ "scalapack-$SCALAPACK_VERSION" ) ARCHIVE_HOMEPAGES+=("http://www.netlib.org/lapack/" \ @@ -626,11 +632,11 @@ function setMathVars { USING_MKLS="OFF" #Add packages that otherwise come preinstalled in intel compiler. ARCHIVE_NAMES+=("OpenBLAS-$OPENBLAS_VERSION.tar.gz" \ - "scalapack-$SCALAPACK_VERSION.tgz") + "$SCALAPACK_VERSION.tar.gz") ARCHIVE_MD5SUMS+=("$OPENBLAS_MD5" \ "$SCALAPACK_MD5" ) ARCHIVE_URLS+=("$OPENBLAS_URL" \ - "http://www.netlib.org/scalapack/scalapack-$SCALAPACK_VERSION.tgz" ) + "https://github.com/Reference-ScaLAPACK/scalapack/archive/$SCALAPACK_VERSION.tar.gz" ) ARCHIVE_DIR_NAMES+=("OpenBLAS-$OPENBLAS_VERSION" \ "scalapack-$SCALAPACK_VERSION" ) ARCHIVE_HOMEPAGES+=("https://github.com/xianyi/OpenBLAS/" \ @@ -685,9 +691,9 @@ function setMathVars { SCALAPACK_INCLUDE_DIR="${GOMA_LIB}/scalapack-$SCALAPACK_VERSION/include" export NON_INTEL_BLAS_LIBRARY="${BLAS_LIBRARY_DIR}/lib${BLAS_LIBRARY_NAME}.a" export NON_INTEL_BLAS_LINK="-L${BLAS_LIBRARY_DIR} -l${BLAS_LIBRARY_NAME} ${FORTRAN_LIBS}" - ARCHIVE_NAMES+=("scalapack-$SCALAPACK_VERSION.tgz") + ARCHIVE_NAMES+=("$SCALAPACK_VERSION.tar.gz") ARCHIVE_MD5SUMS+=("$SCALAPACK_MD5" ) - ARCHIVE_URLS+=("http://www.netlib.org/scalapack/scalapack-$SCALAPACK_VERSION.tgz" ) + ARCHIVE_URLS+=("https://github.com/Reference-ScaLAPACK/scalapack/archive/$SCALAPACK_VERSION.tar.gz" ) ARCHIVE_DIR_NAMES+=("scalapack-$SCALAPACK_VERSION" ) ARCHIVE_HOMEPAGES+=("http://www.netlib.org/scalapack/") ARCHIVE_REAL_NAMES+=("ScaLAPACK") @@ -992,9 +998,36 @@ else exit 1 fi fi + +# PNetCDF cd $GOMA_LIB +if [ -e pnetcdf-${PNETCDF_VERSION}/lib/libpnetcdf.a ] +then + log_echo "PNetCDF already built." +else + if ! [ -e pnetcdf-${NETCDF_VERSION}/.goma-extracted ] + then + mv pnetcdf-${PNETCDF_VERSION} tmpdir + mkdir pnetcdf-${PNETCDF_VERSION} + mv tmpdir pnetcdf-${PNETCDF_VERSION}/src + touch pnetcdf-${PNETCDF_VERSION}/.goma-extracted + fi + cd pnetcdf-${PNETCDF_VERSION}/src + CC=${MPI_C_COMPILER} ./configure --disable-fortran --enable-shared=off --disable-cxx --prefix=$GOMA_LIB/pnetcdf-${PNETCDF_VERSION} 2>&1 | tee -a $COMPILE_LOG + make -j$MAKE_JOBS 2>&1 | tee -a $COMPILE_LOG + make install 2>&1 | tee -a $COMPILE_LOG + if [ -e $GOMA_LIB/pnetcdf-${PNETCDF_VERSION}/lib/libpnetcdf.a ] + then + log_echo "Built PNetCDF $PNETCDF_VERSION" + else + log_echo "Failed to build PNetCDF $PNETCDF_VERSION" + exit 1 + fi +fi + -#netcdf +# NetCDF +cd $GOMA_LIB if [ -e netcdf-${NETCDF_VERSION}/lib/libnetcdf.a ] then log_echo "netcdf already built" @@ -1007,17 +1040,19 @@ else touch netcdf-${NETCDF_VERSION}/.goma-extracted fi cd $GOMA_LIB/netcdf-${NETCDF_VERSION}/src - export CPPFLAGS=-I$GOMA_LIB/hdf5-${HDF5_VERSION}/include + export CPPFLAGS="-I$GOMA_LIB/hdf5-${HDF5_VERSION}/include -I$GOMA_LIB/pnetcdf-${PNETCDF_VERSION}/include" + export LDFLAGS="-L${GOMA_LIB}/hdf5-${HDF5_VERSION}/lib -L${GOMA_LIB}/pnetcdf-${PNETCDF_VERSION}/lib" log_echo $CPPFLAGS log_echo $LDFLAGS - CC=${MPI_C_COMPILER} CFLAGS="-I${GOMA_LIB}/hdf5-${HDF5_VERSION}/include" \ + CC=${MPI_C_COMPILER} CFLAGS="-I${GOMA_LIB}/hdf5-${HDF5_VERSION}/include -I${GOMA_LIB}/pnetcdf-${PNETCDF_VERSION}/include" \ CPP="${MPI_C_COMPILER} -E" \ - CPPFLAGS="-I${GOMA_LIB}/hdf5-${HDF5_VERSION}/include" \ - LDFLAGS="-L${GOMA_LIB}/hdf5-${HDF5_VERSION}/lib" \ + CPPFLAGS="-I${GOMA_LIB}/hdf5-${HDF5_VERSION}/include -I${GOMA_LIB}/pnetcdf-${PNETCDF_VERSION}/include" \ + LDFLAGS="-L${GOMA_LIB}/hdf5-${HDF5_VERSION}/lib -L${GOMA_LIB}/pnetcdf-${PNETCDF_VERSION}/lib" \ ./configure \ --prefix=$GOMA_LIB/netcdf-${NETCDF_VERSION} \ --enable-shared=off \ + --enable-pnetcdf \ --disable-dap 2>&1 | tee -a $COMPILE_LOG make -j$MAKE_JOBS 2>&1 | tee -a $COMPILE_LOG @@ -1029,10 +1064,10 @@ else log_echo "Failed to build NetCDF $NETCDF_VERSION" exit 1 fi - cd ../.. fi #parmetis patch +cd $GOMA_LIB read -d '' PARMETIS_PATCH << "EOF" 55c55,58 < CONFIG_FLAGS += -DCMAKE_C_COMPILER=$(cc) @@ -1456,48 +1491,6 @@ fi # Otherwise Goma dynamically links to UMFPACK when intel is sourced and disabled export LD_LIBRARY_PATH="${GOMA_LIB}/SuiteSparse-$SUITESPARSE_VERSION/UMFPACK/Lib:$LD_LIBRARY_PATH" -#make y12m -cd $GOMA_LIB/y12m -if [ -e liby12m.a ] -then - log_echo "y12m already built" -else -# The y12m we use does not include a makefile, so this is a simple one for compiling a static library. - -cat > makefile << EOF -LIB=y12m -FFLAGS=-O -FC=$MPI_F77_COMPILER -OBJ = \ - y12m.o \ - y12mae.o \ - y12maf.o \ - y12mbe.o \ - y12mbf.o \ - y12mce.o \ - y12mcf.o \ - y12mde.o \ - y12mdf.o \ - y12mfe.o \ - y12mge.o \ - y12mhe.o \ - y12cck.o -lib\$(LIB).a: \$(OBJ) - ar ru lib\$(LIB).a \$? - ranlib lib\$(LIB).a -EOF - - make 2>&1 | tee -a $COMPILE_LOG - if [ -e $GOMA_LIB/y12m/liby12m.a ] - then - log_echo "Built y12m" - else - log_echo "Failed to build y12m" - exit 1 - fi -fi - - if [[ "$MATH_LIBRARIES" == "intel" ]] && [[ ! "$SCALAPACK_LIBRARY_NAME" = "scalapack" ]]; then log_echo "Not building scalapack because intel MKL used" else @@ -1511,1379 +1504,14 @@ else mv src-scalapack $GOMA_LIB/scalapack-$SCALAPACK_VERSION/src fi cd $GOMA_LIB/scalapack-$SCALAPACK_VERSION/src - # Use PETSc's parallel make changes -cat > Makefile.objs < Makefile.parallel < SLmake.inc <&1 | tee -a $COMPILE_LOG + cmake --build build -j $MAKE_JOBS 2>&1 | tee -a $COMPILE_LOG + cmake --install build 2>&1 | tee -a $COMPILE_LOG cd $GOMA_LIB if [ -f $GOMA_LIB/scalapack-$SCALAPACK_VERSION/lib/libscalapack.a ]; then log_echo "Build scalapack $SCALAPACK_VERSION" @@ -3018,123 +1646,6 @@ TRILINOS_INSTALL=$GOMA_LIB/trilinos-$TRILINOS_VERSION if [ -e $TRILINOS_INSTALL/bin/aprepro ]; then log_echo "Trilinos is already built!" else -cd $GOMA_LIB/Trilinos-trilinos-release-$TRILINOS_VERSION_DASH/packages/aztecoo/src -cat << "EOF" > az_aztec_h.patch -731c731 -< extern char *AZ_allocate(unsigned int iii); ---- -> extern char *AZ_allocate(size_t iii); -1024c1024 -< extern double *AZ_manage_memory(unsigned int size, int action, int type, ---- -> extern double *AZ_manage_memory(size_t size, int action, int type, -1198c1198 -< extern char *AZ_realloc(void *ptr, unsigned int size); ---- -> extern char *AZ_realloc(void *ptr, size_t size); -EOF -patch -f az_aztec.h < az_aztec_h.patch -rm az_aztec_h.patch - -cat << "EOF" > az_util_c.patch -843c843 -< (void) AZ_manage_memory((unsigned int) 0, AZ_CLEAR, label, (char *) NULL, ---- -> (void) AZ_manage_memory((size_t) 0, AZ_CLEAR, label, (char *) NULL, -851c851 -< double *AZ_manage_memory(unsigned int input_size, int action, int type, ---- -> double *AZ_manage_memory(size_t input_size, int action, int type, -936c936 -< int size; ---- -> size_t size; -941c941 -< long int size; ---- -> size_t size, aligned_size; -945,946c945,946 -< long int aligned_str_mem, aligned_j, aligned_size; -< double *dtmp; ---- -> unsigned int aligned_str_mem, aligned_j; -> double *dtmp; -949c949 -< size = (long int) input_size; ---- -> size = input_size; -997,998c997 -< dtmp = (double *) AZ_allocate((unsigned int) (aligned_str_mem+aligned_j+ -< aligned_size) ); ---- -> dtmp = (double *) AZ_allocate(aligned_str_mem+aligned_j+aligned_size); -1183,1184c1182 -< dtmp = (double *) AZ_realloc((char *) dtmp,(unsigned int) -< aligned_str_mem+aligned_j+aligned_size); ---- -> dtmp = (double *) AZ_realloc((char *) dtmp, aligned_str_mem+aligned_j+aligned_size); -1872c1870 -< int size; ---- -> size_t size; -1902c1900 -< char *AZ_allocate(unsigned int isize) { ---- -> char *AZ_allocate(size_t isize) { -1917,1918c1915,1917 -< int *size_ptr, i; -< unsigned int size; ---- -> int i; -> size_t size; -> size_t *size_ptr; -1952c1951 -< size_ptr = (int *) ptr; ---- -> size_ptr = (size_t *) ptr; -1992c1991 -< int *iptr, size, i; ---- -> int i; -1993a1993,1994 -> size_t size; -> size_t* iptr; -2023c2024 -< iptr = (int *) ptr; ---- -> iptr = (size_t *) ptr; -2073c2074 -< char *AZ_realloc(void *vptr, unsigned int new_size) { ---- -> char *AZ_realloc(void *vptr, size_t new_size) { -2076c2077 -< int i, *iptr, size, *new_size_ptr; ---- -> int i; -2079d2079 -< int newmsize, smaller; -2080a2081,2082 -> size_t size, newmsize, smaller; -> size_t *iptr, *new_size_ptr; -2108c2110 -< iptr = (int *) ptr; ---- -> iptr = (size_t *) ptr; -2128c2130 -< new_size_ptr = (int *) new_ptr; ---- -> new_size_ptr = (size_t *) new_ptr; -2175c2177 -< char *AZ_allocate(unsigned int size) { ---- -> char *AZ_allocate(size_t size) { -2185c2187 -< char *AZ_realloc(void *ptr, unsigned int size) { ---- -> char *AZ_realloc(void *ptr, size_t size) { -EOF -patch -f az_util.c < az_util_c.patch -rm az_util_c.patch cd $GOMA_LIB/trilinos-$TRILINOS_VERSION-Temp cmake \ -D CMAKE_AR=/usr/bin/ar \ @@ -3148,7 +1659,6 @@ cd $GOMA_LIB/trilinos-$TRILINOS_VERSION-Temp -D TPL_ENABLE_Boost:BOOL=OFF \ -D Trilinos_ENABLE_Triutils:BOOL=ON \ -D Trilinos_ENABLE_SEACAS:BOOL=ON \ --D Trilinos_GOMA_ENABLE_AMESOS:BOOL=ON \ -D Trilinos_ENABLE_Epetra:BOOL=ON \ -D Trilinos_ENABLE_Xpetra:BOOL=ON \ -D Trilinos_ENABLE_Ifpack:BOOL=ON \ @@ -3159,8 +1669,9 @@ cd $GOMA_LIB/trilinos-$TRILINOS_VERSION-Temp -D Trilinos_ENABLE_AztecOO:BOOL=ON \ -D Trilinos_ENABLE_Stratimikos:BOOL=ON \ -D Trilinos_ENABLE_Teko:BOOL=ON \ --D Trilinos_GOMA_ENABLE_AMESOS2:BOOL=ON \ -D Trilinos_ENABLE_Belos:BOOL=ON \ +-D Trilinos_ENABLE_Amesos2:BOOL=ON \ +-D Trilinos_ENABLE_Sacado:BOOL=ON \ -D Trilinos_ENABLE_EpetraExt:BOOL=ON \ -D Trilinos_ENABLE_Thyra:BOOL=ON \ -D Trilinos_ENABLE_ThyraTpetraAdapters:BOOL=ON \ @@ -3173,9 +1684,10 @@ cd $GOMA_LIB/trilinos-$TRILINOS_VERSION-Temp -D TPL_HDF5_INCLUDE_DIRS:PATH="$GOMA_LIB/hdf5-${HDF5_VERSION}/include" \ -D HDF5_LIBRARY_DIRS:PATH="$GOMA_LIB/hdf5-${HDF5_VERSION}/lib" \ -D HDF5_LIBRARY_NAMES:STRING="hdf5_hl;hdf5;z;dl" \ - -D Netcdf_LIBRARY_DIRS:PATH="$GOMA_LIB/netcdf-${NETCDF_VERSION}/lib" \ + -D Netcdf_LIBRARY_DIRS:PATH="$GOMA_LIB/netcdf-${NETCDF_VERSION}/lib;$GOMA_LIB/pnetcdf-${PNETCDF_VERSION}/lib" \ -D TPL_ENABLE_Netcdf:BOOL=ON \ - -D TPL_Netcdf_INCLUDE_DIRS:PATH="$GOMA_LIB/netcdf-${NETCDF_VERSION}/include" \ + -D TPL_Netcdf_INCLUDE_DIRS:PATH="$GOMA_LIB/netcdf-${NETCDF_VERSION}/include;$GOMA_LIB/pnetcdf-${PNETCDF_VERSION}/include" \ + -D Netcdf_LIBRARY_NAMES:STRING="netcdf;pnetcdf" \ -D Matio_LIBRARY_DIRS:PATH=$GOMA_LIB/matio-$MATIO_VERSION/lib \ -D Matio_INCLUDE_DIRS:PATH=$GOMA_LIB/matio-$MATIO_VERSION/include \ -D TPL_ENABLE_MPI:BOOL=ON \ @@ -3204,11 +1716,7 @@ cd $GOMA_LIB/trilinos-$TRILINOS_VERSION-Temp -D SuperLUDist_INCLUDE_DIRS:PATH=$GOMA_LIB/superlu_dist-$SUPERLU_DIST_VERSION/include \ -D TPL_ENABLE_ParMETIS:BOOL=ON \ -D ParMETIS_LIBRARY_DIRS:PATH="$GOMA_LIB/parmetis-$PARMETIS_VERSION/lib;$GOMA_LIB/metis-$METIS_VERSION/lib" \ - -D ParMETIS_INCLUDE_DIRS:PATH="$GOMA_LIB/parmetis-$PARMETIS_VERSION/include;$GOMA_LIB/metis-$METIS_VERSION/include" \ -D TPL_ParMETIS_INCLUDE_DIRS:PATH="$GOMA_LIB/parmetis-$PARMETIS_VERSION/include;$GOMA_LIB/metis-$METIS_VERSION/include" \ - -D TPL_ENABLE_y12m:BOOL=ON \ - -D y12m_LIBRARY_NAMES:STRING="y12m" \ - -D y12m_LIBRARY_DIRS:PATH=$GOMA_LIB/y12m \ -D TPL_ENABLE_MUMPS:BOOL=ON \ -D MUMPS_LIBRARY_NAMES:STRING="dmumps;mumps_common;pord;$BLACS_LIBRARY_NAME" \ -D MUMPS_LIBRARY_DIRS:PATH="$GOMA_LIB/MUMPS_$MUMPS_VERSION/lib;$GOMA_LIB/MUMPS_$MUMPS_VERSION/PORD/lib;$SCALAPACK_LIBRARY_DIR" \ @@ -3229,8 +1737,11 @@ cd $GOMA_LIB/trilinos-$TRILINOS_VERSION-Temp -D Amesos_ENABLE_KLU:BOOL=ON \ -D Amesos_ENABLE_UMFPACK:BOOL=ON \ -D Amesos_ENABLE_MUMPS:BOOL=ON \ +-D Tpetra_INST_INT_INT:BOOL=ON \ $EXTRA_ARGS \ $GOMA_LIB/Trilinos-trilinos-release-$TRILINOS_VERSION_DASH 2>&1 | tee -a $COMPILE_LOG +# can set LONG LONG but only one is valid +# -D Tpetra_INST_INT_LONG_LONG:BOOL=ON \ make -j$MAKE_JOBS 2>&1 | tee -a $COMPILE_LOG make install 2>&1 | tee -a $COMPILE_LOG diff --git a/src/ac_conti.c b/src/ac_conti.c index 353dba9e7..3d8f8a41a 100644 --- a/src/ac_conti.c +++ b/src/ac_conti.c @@ -34,6 +34,7 @@ #include "el_geom.h" #include "el_quality.h" #include "exo_struct.h" +#include "linalg/sparse_matrix.h" #include "mm_as.h" #include "mm_as_structs.h" #include "mm_augc_util.h" @@ -56,6 +57,7 @@ #include "rf_io.h" #include "rf_io_const.h" #include "rf_io_structs.h" +#include "rf_masks.h" #include "rf_mp.h" #include "rf_node_const.h" #include "rf_solve.h" @@ -547,7 +549,21 @@ void continue_problem(Comm_Ex *cx, /* array of communications structures */ pg->matrices[pg->imtrx].resid_vector = resid_vector; /* Allocate sparse matrix */ - if (strcmp(Matrix_Format, "msr") == 0) { + if ((strcmp(Matrix_Format, "tpetra") == 0) || (strcmp(Matrix_Format, "epetra") == 0)) { + err = check_compatible_solver(); + GOMA_EH(err, "Incompatible matrix solver for tpetra, tpetra supports stratimikos"); + check_parallel_error("Matrix format / Solver incompatibility"); + GomaSparseMatrix goma_matrix; + goma_error err = GomaSparseMatrix_CreateFromFormat(&goma_matrix, Matrix_Format); + GOMA_EH(err, "GomaSparseMatrix_CreateFromFormat"); + int local_nodes = Num_Internal_Nodes + Num_Border_Nodes + Num_External_Nodes; + err = GomaSparseMatrix_SetProblemGraph( + goma_matrix, num_internal_dofs[pg->imtrx], num_boundary_dofs[pg->imtrx], + num_external_dofs[pg->imtrx], local_nodes, Nodes, MaxVarPerNode, Matilda, Inter_Mask, exo, + dpi, cx, pg->imtrx, Debug_Flag, ams[JAC]); + GOMA_EH(err, "GomaSparseMatrix_SetProblemGraph"); + ams[JAC]->GomaMatrixData = goma_matrix; + } else if (strcmp(Matrix_Format, "msr") == 0) { log_msg("alloc_MSR_sparse_arrays..."); alloc_MSR_sparse_arrays(&ija, &a, &a_old, 0, node_to_fill, exo, dpi); /* diff --git a/src/ac_hunt.c b/src/ac_hunt.c index e82963f78..c8c97e456 100644 --- a/src/ac_hunt.c +++ b/src/ac_hunt.c @@ -34,6 +34,7 @@ #include "el_geom.h" #include "el_quality.h" #include "exo_struct.h" +#include "linalg/sparse_matrix.h" #include "mm_as.h" #include "mm_as_structs.h" #include "mm_augc_util.h" @@ -56,13 +57,13 @@ #include "rf_io.h" #include "rf_io_const.h" #include "rf_io_structs.h" +#include "rf_masks.h" #include "rf_mp.h" #include "rf_node_const.h" #include "rf_solve.h" #include "rf_solver.h" #include "rf_util.h" #include "sl_auxutil.h" -#include "sl_epetra_util.h" #include "sl_matrix_util.h" #ifdef GOMA_ENABLE_PETSC #include "sl_petsc.h" @@ -539,14 +540,25 @@ void hunt_problem(Comm_Ex *cx, /* array of communications structures */ /* Allocate sparse matrix */ - if (strcmp(Matrix_Format, "epetra") == 0) { + if ((strcmp(Matrix_Format, "tpetra") == 0) || (strcmp(Matrix_Format, "epetra") == 0)) { err = check_compatible_solver(); - GOMA_EH(err, - "Incompatible matrix solver for epetra, epetra supports amesos and aztecoo solvers."); + GOMA_EH(err, "Incompatible matrix solver for matrix format %s", Matrix_Format); check_parallel_error("Matrix format / Solver incompatibility"); - ams[JAC]->RowMatrix = - EpetraCreateRowMatrix(num_internal_dofs[pg->imtrx] + num_boundary_dofs[pg->imtrx]); - EpetraCreateGomaProblemGraph(ams[JAC], exo, dpi); + GomaSparseMatrix goma_matrix; + goma_error err = GomaSparseMatrix_CreateFromFormat(&goma_matrix, Matrix_Format); + GOMA_EH(err, "GomaSparseMatrix_CreateFromFormat"); + + int local_nodes = Num_Internal_Nodes + Num_Border_Nodes + Num_External_Nodes; + err = GomaSparseMatrix_SetProblemGraph( + goma_matrix, num_internal_dofs[pg->imtrx], num_boundary_dofs[pg->imtrx], + num_external_dofs[pg->imtrx], local_nodes, Nodes, MaxVarPerNode, Matilda, Inter_Mask, exo, + dpi, cx, pg->imtrx, Debug_Flag, ams[JAC]); + GOMA_EH(err, "GomaSparseMatrix_SetProblemGraph"); + + ams[JAC]->GomaMatrixData = goma_matrix; + ams[JAC]->npu = num_internal_dofs[pg->imtrx] + num_boundary_dofs[pg->imtrx]; + ams[JAC]->npu_plus = num_universe_dofs[pg->imtrx]; + #ifdef GOMA_ENABLE_PETSC #if PETSC_USE_COMPLEX } else if (strcmp(Matrix_Format, "petsc_complex") == 0) { diff --git a/src/adapt/resetup_problem.c b/src/adapt/resetup_problem.c index 48dffc106..6d364e414 100644 --- a/src/adapt/resetup_problem.c +++ b/src/adapt/resetup_problem.c @@ -10,19 +10,20 @@ #include "dp_map_comm_vec.h" #include "dp_types.h" #include "dpi.h" +#include "linalg/sparse_matrix.h" #include "mm_as.h" #include "mm_as_structs.h" #include "mm_eh.h" #include "mm_unknown_map.h" #include "rf_fem.h" #include "rf_fem_const.h" +#include "rf_io.h" +#include "rf_masks.h" #include "rf_mp.h" #include "rf_node_const.h" #include "rf_solver.h" #include "rf_util.h" #include "rf_vars_const.h" -#include "sl_epetra_interface.h" -#include "sl_epetra_util.h" #include "sl_util_structs.h" int resetup_problem(Exo_DB *exo, /* ptr to the finite element mesh database */ @@ -188,12 +189,17 @@ int resetup_problem(Exo_DB *exo, /* ptr to the finite element mesh database */ } int resetup_matrix(struct GomaLinearSolverData **ams, Exo_DB *exo, Dpi *dpi) { - if (strcmp(Matrix_Format, "epetra") == 0) { + if ((strcmp(Matrix_Format, "tpetra") == 0) || (strcmp(Matrix_Format, "epetra") == 0)) { for (pg->imtrx = 0; pg->imtrx < upd->Total_Num_Matrices; pg->imtrx++) { - EpetraDeleteRowMatrix(ams[pg->imtrx]->RowMatrix); - ams[pg->imtrx]->RowMatrix = - EpetraCreateRowMatrix(num_internal_dofs[pg->imtrx] + num_boundary_dofs[pg->imtrx]); - EpetraCreateGomaProblemGraph(ams[pg->imtrx], exo, dpi); + GomaSparseMatrix goma_matrix = ams[pg->imtrx]->GomaMatrixData; + GomaSparseMatrix_Destroy(&goma_matrix); + GomaSparseMatrix_CreateFromFormat(&goma_matrix, Matrix_Format); + ams[pg->imtrx]->GomaMatrixData = goma_matrix; + int local_nodes = exo->num_nodes; + GomaSparseMatrix_SetProblemGraph(goma_matrix, num_internal_dofs[pg->imtrx], + num_boundary_dofs[pg->imtrx], num_external_dofs[pg->imtrx], + local_nodes, Nodes, MaxVarPerNode, Matilda, Inter_Mask, exo, + dpi, cx[pg->imtrx], pg->imtrx, Debug_Flag, ams[JAC]); } pg->imtrx = 0; ams[pg->imtrx]->solveSetup = 0; diff --git a/src/bc/rotate.c b/src/bc/rotate.c index 197a4e41c..bb1c18058 100644 --- a/src/bc/rotate.c +++ b/src/bc/rotate.c @@ -27,6 +27,7 @@ #include "el_geom.h" #include "exo_struct.h" #include "gds/gds_vector.h" +#include "linalg/sparse_matrix.h" #include "load_field_variables.h" #include "mm_as.h" #include "mm_as_alloc.h" @@ -69,7 +70,6 @@ int dup_blks_list[MAX_MAT_PER_SS + 1]; int rotation_allocated = FALSE; #define GOMA_BC_ROTATE_C -#include "sl_epetra_interface.h" /*********** R O U T I N E S I N T H I S F I L E ************************* * @@ -782,7 +782,8 @@ void rotate_mesh_eqn(int id, /* Elemental stiffness matrix row index * } } } - } else if (strcmp(Matrix_Format, "epetra") == 0) { + } else if (ams->GomaMatrixData != NULL) { + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; if (I < (DPI_ptr->num_internal_nodes + DPI_ptr->num_boundary_nodes)) { /* * Find the global equation number @@ -796,8 +797,8 @@ void rotate_mesh_eqn(int id, /* Elemental stiffness matrix row index * */ for (j = 0; j < rot->d_vector_n; j++) { double sum_val; - int global_row; - int global_col; + GomaGlobalOrdinal global_row; + GomaGlobalOrdinal global_col; J = rot->d_vector_J[j]; if (Dolphin[pg->imtrx][I][MESH_DISPLACEMENT1] > 0 && @@ -813,10 +814,9 @@ void rotate_mesh_eqn(int id, /* Elemental stiffness matrix row index * sum_val += rot->d_vector_dx[ldir][b][j] * lec->R[LEC_R_INDEX(peqn_mesh[ldir], id)]; } - global_row = ams->GlobalIDs[index_eqn]; - global_col = ams->GlobalIDs[index_var]; - EpetraSumIntoGlobalRowMatrix(ams->RowMatrix, global_row, 1, &sum_val, - &global_col); + global_row = matrix->global_ids[index_eqn]; + global_col = matrix->global_ids[index_var]; + matrix->sum_into_row_values(matrix, global_row, 1, &sum_val, &global_col); } } } @@ -1133,13 +1133,14 @@ void rotate_momentum_eqn(int id, /* Elemental stiffness matrix row ind } } - } else if (strcmp(Matrix_Format, "epetra") == 0) { + } else if (ams->GomaMatrixData != NULL) { + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; if (I < (DPI_ptr->num_internal_nodes + DPI_ptr->num_boundary_nodes)) { // Direct translation from MSR for (j = 0; j < rotation[I][eq][kdir]->d_vector_n; j++) { double sum_val; - int global_row; - int global_col; + GomaGlobalOrdinal global_row; + GomaGlobalOrdinal global_col; int ktype, ndof, index_eqn, index_var; J = rotation[I][eq][kdir]->d_vector_J[j]; if (Dolphin[pg->imtrx][I][R_MOMENTUM1] > 0 && @@ -1164,10 +1165,9 @@ void rotate_momentum_eqn(int id, /* Elemental stiffness matrix row ind sum_val += rotation[I][eq][kdir]->d_vector_dx[ldir][b][j] * lec->R[LEC_R_INDEX(upd->ep[pg->imtrx][R_MESH1 + ldir], id)]; } - global_row = ams->GlobalIDs[index_eqn]; - global_col = ams->GlobalIDs[index_var]; - EpetraSumIntoGlobalRowMatrix(ams->RowMatrix, global_row, 1, &sum_val, - &global_col); + global_row = matrix->global_ids[index_eqn]; + global_col = matrix->global_ids[index_var]; + matrix->sum_into_row_values(matrix, global_row, 1, &sum_val, &global_col); } /* end of Baby_dolphin */ } } diff --git a/src/dp_comm.c b/src/dp_comm.c index 5f715e588..5e9199e48 100644 --- a/src/dp_comm.c +++ b/src/dp_comm.c @@ -100,6 +100,62 @@ void exchange_dof(Comm_Ex *cx, Dpi *dpi, double *x, int imtrx) safer_free((void **)&ptr_send_list); #endif /* PARALLEL */ } + +void exchange_dof_int(Comm_Ex *cx, Dpi *dpi, int *x, int imtrx) + +/************************************************************ + * + * exchange_dof(): + * + * send/recv appropriate pieces of a dof-based double array + ************************************************************/ +{ + COMM_NP_STRUCT *np_base, *np_ptr; + int *ptr_send_list, *ptr_recv_list; + register int *ptrd; + register int *ptr_int, i; + int p; + int num_neighbors = dpi->num_neighbors; + int total_num_send_unknowns; + + if (num_neighbors == 0) + return; + +#ifdef PARALLEL + total_num_send_unknowns = ptr_dof_send[imtrx][dpi->num_neighbors]; + np_base = alloc_struct_1(COMM_NP_STRUCT, dpi->num_neighbors); + ptrd = malloc(total_num_send_unknowns * sizeof(int)); + ptr_send_list = ptrd; + + /* + * gather up the list of send unknowns + */ + ptr_int = list_dof_send[imtrx]; + for (i = total_num_send_unknowns; i > 0; i--) { + *ptrd++ = x[*ptr_int++]; + } + + /* + * store base address for the start of the external degrees of freedom + * in this vector + */ + ptr_recv_list = x + num_internal_dofs[imtrx] + num_boundary_dofs[imtrx]; + + np_ptr = np_base; + for (p = 0; p < dpi->num_neighbors; p++) { + np_ptr->neighbor_ProcID = cx[p].neighbor_name; + np_ptr->send_message_buf = (void *)(ptr_send_list + ptr_dof_send[imtrx][p]); + np_ptr->send_message_length = sizeof(int) * cx[p].num_dofs_send; + np_ptr->recv_message_buf = (void *)ptr_recv_list; + np_ptr->recv_message_length = sizeof(int) * cx[p].num_dofs_recv; + ptr_recv_list += cx[p].num_dofs_recv; + np_ptr++; + } + exchange_neighbor_proc_info(dpi->num_neighbors, np_base); + safer_free((void **)&np_base); + safer_free((void **)&ptr_send_list); +#endif /* PARALLEL */ +} /********************************************************************/ /********************************************************************/ /********************************************************************/ @@ -196,6 +252,62 @@ void exchange_node(Comm_Ex *cx, Dpi *dpi, double *x) safer_free((void **)&ptr_send_list); #endif } + +void exchange_dof_long_long(Comm_Ex *cx, Dpi *dpi, long long *x, int imtrx) + +/************************************************************ + * + * exchange_dof(): + * + * send/recv appropriate pieces of a dof-based double array + ************************************************************/ +{ + COMM_NP_STRUCT *np_base, *np_ptr; + long long *ptr_send_list, *ptr_recv_list; + register long long *ptrd; + register int *ptr_int, i; + int p; + int num_neighbors = dpi->num_neighbors; + int total_num_send_unknowns; + + if (num_neighbors == 0) + return; + +#ifdef PARALLEL + total_num_send_unknowns = ptr_dof_send[imtrx][dpi->num_neighbors]; + np_base = alloc_struct_1(COMM_NP_STRUCT, dpi->num_neighbors); + ptrd = malloc(total_num_send_unknowns * sizeof(long long)); + ptr_send_list = ptrd; + + /* + * gather up the list of send unknowns + */ + ptr_int = list_dof_send[imtrx]; + for (i = total_num_send_unknowns; i > 0; i--) { + *ptrd++ = x[*ptr_int++]; + } + + /* + * store base address for the start of the external degrees of freedom + * in this vector + */ + ptr_recv_list = x + num_internal_dofs[imtrx] + num_boundary_dofs[imtrx]; + + np_ptr = np_base; + for (p = 0; p < dpi->num_neighbors; p++) { + np_ptr->neighbor_ProcID = cx[p].neighbor_name; + np_ptr->send_message_buf = (void *)(ptr_send_list + ptr_dof_send[imtrx][p]); + np_ptr->send_message_length = sizeof(long long) * cx[p].num_dofs_send; + np_ptr->recv_message_buf = (void *)ptr_recv_list; + np_ptr->recv_message_length = sizeof(long long) * cx[p].num_dofs_recv; + ptr_recv_list += cx[p].num_dofs_recv; + np_ptr++; + } + exchange_neighbor_proc_info(dpi->num_neighbors, np_base); + safer_free((void **)&np_base); + safer_free((void **)&ptr_send_list); +#endif /* PARALLEL */ +} /********************************************************************/ /********************************************************************/ /********************************************************************/ diff --git a/src/dp_ghost.cpp b/src/dp_ghost.cpp index c939da475..502bf5899 100644 --- a/src/dp_ghost.cpp +++ b/src/dp_ghost.cpp @@ -26,7 +26,6 @@ #include // not needed except to avoid including as a C file -#include "sl_epetra_interface.h" #include "dp_ghost.h" diff --git a/src/dp_vif.c b/src/dp_vif.c index f0b5cc894..c87b82266 100644 --- a/src/dp_vif.c +++ b/src/dp_vif.c @@ -691,7 +691,9 @@ void noahs_ark(void) { ddd_add_member(n, Matrix_Relative_Threshold, MAX_CHAR_IN_INPUT, MPI_CHAR); ddd_add_member(n, Matrix_Absolute_Threshold, MAX_CHAR_IN_INPUT, MPI_CHAR); ddd_add_member(n, Amesos_Package, MAX_CHAR_IN_INPUT, MPI_CHAR); + ddd_add_member(n, Amesos2_Package, MAX_CHAR_IN_INPUT, MPI_CHAR); ddd_add_member(n, Stratimikos_File, MAX_CHAR_IN_INPUT * MAX_NUM_MATRICES, MPI_CHAR); + ddd_add_member(n, Amesos2_File, MAX_CHAR_IN_INPUT * MAX_NUM_MATRICES, MPI_CHAR); ddd_add_member(n, &Linear_Solver, 1, MPI_INT); diff --git a/src/globals.c b/src/globals.c index 57738fc09..b6b5f9631 100644 --- a/src/globals.c +++ b/src/globals.c @@ -203,10 +203,14 @@ String_line Matrix_Absolute_Threshold; /* Trilinos 2 */ String_line Amesos_Package; +String_line Amesos2_Package; + String_line AztecOO_Solver; String_line Stratimikos_File[MAX_NUM_MATRICES]; +String_line Amesos2_File[MAX_NUM_MATRICES]; + /* * A new Aztec 2.0 option. There are more and difft options and our * previous options probably ought to be revised to reflect the newer diff --git a/src/sl_epetra_util.cpp b/src/linalg/sparse_matrix.cpp similarity index 70% rename from src/sl_epetra_util.cpp rename to src/linalg/sparse_matrix.cpp index 266609766..675a2bcfc 100644 --- a/src/sl_epetra_util.cpp +++ b/src/linalg/sparse_matrix.cpp @@ -1,63 +1,78 @@ -#ifndef __cplusplus -#define __cplusplus -#endif - -#if defined(PARALLEL) && !defined(EPETRA_MPI) -#define EPETRA_MPI -#endif - -#include -#include -#include +#include #include -#include "Epetra_ConfigDefs.h" -#include "Epetra_CrsMatrix.h" -#include "Epetra_Map.h" -#include "Epetra_Vector.h" -#include "dpi.h" -#include "exo_struct.h" - -#ifdef EPETRA_MPI -#else -#include "Epetra_SerialComm.h" +#include "linalg/sparse_matrix.h" +#ifdef GOMA_ENABLE_TPETRA +#include "linalg/sparse_matrix_tpetra.h" +#endif +#ifdef GOMA_ENABLE_EPETRA +#include "linalg/sparse_matrix_epetra.h" #endif extern "C" { +#define DISABLE_CPP #include "dp_comm.h" #include "dp_types.h" -#include "el_geom.h" #include "mm_as.h" #include "mm_as_const.h" -#include "mm_as_structs.h" #include "mm_eh.h" #include "mm_fill_util.h" #include "mm_mp.h" -#include "mm_mp_structs.h" #include "mm_unknown_map.h" -#include "rf_fem.h" -#include "rf_fem_const.h" -#include "rf_io.h" #include "rf_masks.h" #include "rf_node_const.h" -#include "rf_solve.h" -#include "rf_vars_const.h" #include "sl_util_structs.h" -#include "std.h" +#undef DISABLE_CPP } -#include "sl_epetra_util.h" - -#include "sl_epetra_interface.h" +extern "C" goma_error GomaSparseMatrix_CreateFromFormat(GomaSparseMatrix *matrix, + char *matrix_format) { + if (strcmp(matrix_format, "tpetra") == 0) { + return GomaSparseMatrix_Create(matrix, GOMA_SPARSE_MATRIX_TYPE_TPETRA); + } else if (strcmp(matrix_format, "epetra") == 0) { + return GomaSparseMatrix_Create(matrix, GOMA_SPARSE_MATRIX_TYPE_EPETRA); + } + return GOMA_ERROR; +} -extern "C" { +extern "C" goma_error GomaSparseMatrix_Create(GomaSparseMatrix *matrix, + enum GomaSparseMatrixType type) { + *matrix = (GomaSparseMatrix)malloc(sizeof(struct g_GomaSparseMatrix)); + switch (type) { +#ifdef GOMA_ENABLE_TPETRA + case GOMA_SPARSE_MATRIX_TYPE_TPETRA: + return GomaSparseMatrix_Tpetra_Create(matrix); + break; +#endif +#ifdef GOMA_ENABLE_EPETRA + case GOMA_SPARSE_MATRIX_TYPE_EPETRA: + return GomaSparseMatrix_Epetra_Create(matrix); + break; +#endif + default: + GOMA_EH(GOMA_ERROR, "Unknown matrix type, GomaSparseMatrix_Create"); + return GOMA_ERROR; + break; + } + return GOMA_SUCCESS; +} -/** - * Create the goma problem graph in the epetra matrix ams->RowMatrix - * @param ams ams structure containing appropriate RowMatrix - * @param exo exodus file for this processor - */ -void EpetraCreateGomaProblemGraph(struct GomaLinearSolverData *ams, Exo_DB *exo, Dpi *dpi) { +extern "C" goma_error GomaSparseMatrix_SetProblemGraph( + GomaSparseMatrix matrix, + int num_internal_dofs, + int num_boundary_dofs, + int num_external_dofs, + int local_nodes, + NODE_INFO_STRUCT **Nodes, + int MaxVarPerNode, + int *Matilda, + int Inter_Mask[MAX_NUM_MATRICES][MAX_VARIABLE_TYPES][MAX_VARIABLE_TYPES], + Exo_DB *exo, + Dpi *dpi, + Comm_Ex *cx, + int imtrx, + int Debug_Flag, + struct GomaLinearSolverData *ams) { int j, inode, i1, i2, eb1; int iunknown, inter_unknown, inter_node, row_num_unknowns, col_num_unknowns; int irow_index = 0; @@ -68,49 +83,46 @@ void EpetraCreateGomaProblemGraph(struct GomaLinearSolverData *ams, Exo_DB *exo, int inode_varType[MaxVarPerNode], inode_matID[MaxVarPerNode]; int inter_node_varType[MaxVarPerNode], inter_node_matID[MaxVarPerNode]; int nnz = 0; - int total_nodes = Num_Internal_Nodes + Num_Border_Nodes + Num_External_Nodes; - std::vector Indices; - std::vector Values; - int NumMyRows = num_internal_dofs[pg->imtrx] + num_boundary_dofs[pg->imtrx]; - int NumExternal = num_external_dofs[pg->imtrx]; - int NumMyCols = NumMyRows + NumExternal; + GomaGlobalOrdinal NumMyRows = num_internal_dofs + num_boundary_dofs; + GomaGlobalOrdinal NumExternal = num_external_dofs; + GomaGlobalOrdinal NumMyCols = NumMyRows + NumExternal; + matrix->n_rows = NumMyRows; + matrix->n_cols = NumMyCols; - /* - * This is kind of hacky the way it is being done, - * modified from the amesos interface for gomamsr to epetra - * - * Creates an array to be communicated of ints cast to double to use the exchange_dof - * communicator, then casts doubles back to int and sets those as global ids - * for the epetra array - * - * TODO: replace with non-double conversion routine - */ + GomaGlobalOrdinal RowOffset; + MPI_Scan(&NumMyRows, &RowOffset, 1, MPI_GOMA_ORDINAL, MPI_SUM, MPI_COMM_WORLD); + RowOffset -= NumMyRows; - // get the row map from the row matrix - Epetra_Map RowMap = ams->RowMatrix->RowMatrixRowMap(); + std::vector GlobalIDs(NumMyCols); - // get the global elements for this processor - int *MyGlobalElements = RowMap.MyGlobalElements(); + for (int i = 0; i < NumMyRows; i++) { + GlobalIDs[i] = i + RowOffset; + } + matrix->global_ids = (GomaGlobalOrdinal *)malloc(sizeof(GomaGlobalOrdinal) * NumMyCols); - double *dblColGIDs = new double[NumMyCols]; - ams->GlobalIDs = (int *)malloc(sizeof(int) * NumMyCols); +#ifdef GOMA_MATRIX_GO_LONG_LONG + exchange_dof_long_long(cx, dpi, GlobalIDs.data(), imtrx); +#else + exchange_dof_int(cx, dpi, GlobalIDs.data(), imtrx); +#endif - // copy global id's and convert to double for boundary exchange - for (int i = 0; i < NumMyRows; i++) - dblColGIDs[i] = (double)MyGlobalElements[i]; + for (size_t i = 0; i < GlobalIDs.size(); i++) { + matrix->global_ids[i] = GlobalIDs[i]; + } - exchange_dof(cx[pg->imtrx], dpi, dblColGIDs, pg->imtrx); + std::vector rows(GlobalIDs.begin(), GlobalIDs.begin() + NumMyRows); + std::vector cols(GlobalIDs.begin(), GlobalIDs.end()); + std::vector coo_rows; + std::vector coo_cols; - // convert back to int with known global id's from all processors - for (int j = 0; j < NumMyCols; j++) { - ams->GlobalIDs[j] = (int)dblColGIDs[j]; - } + int max_nz_per_row = 0; + int row_nz; /* * loop over all of the nodes on this processor */ - for (inode = 0; inode < total_nodes; inode++) { + for (inode = 0; inode < local_nodes; inode++) { nv = Nodes[inode]->Nodal_Vars_Info[pg->imtrx]; /* * Fill the vector list which points to the unknowns defined at this @@ -129,14 +141,12 @@ void EpetraCreateGomaProblemGraph(struct GomaLinearSolverData *ams, Exo_DB *exo, * Loop over the unknowns defined at this row node */ for (iunknown = 0; iunknown < row_num_unknowns; iunknown++) { + row_nz = 0; /* * Retrieve the var type of the current unknown */ rowVarType = inode_varType[iunknown]; - Indices.clear(); - Values.clear(); - /* * Loop over the nodes which are determined to have an interaction * with the current row node @@ -209,20 +219,26 @@ void EpetraCreateGomaProblemGraph(struct GomaLinearSolverData *ams, Exo_DB *exo, * Determine the equation number of the current unknown */ icol_index = nodeCol->First_Unknown[pg->imtrx] + inter_unknown; - Indices.push_back(ams->GlobalIDs[icol_index]); - Values.push_back(0); + coo_rows.push_back(GlobalIDs[irow_index]); + coo_cols.push_back(GlobalIDs[icol_index]); + nnz++; + row_nz++; } } } - EpetraInsertGlobalRowMatrix(ams->RowMatrix, ams->GlobalIDs[irow_index], Indices.size(), - &Values[0], &Indices[0]); - nnz += Indices.size(); + if (row_nz > max_nz_per_row) { + max_nz_per_row = row_nz; + } irow_index++; } } - EpetraFillCompleteRowMatrix(ams->RowMatrix); - EpetraPutScalarRowMatrix(ams->RowMatrix, 0); + matrix->create_graph(matrix, NumMyRows, rows.data(), NumMyCols, cols.data(), nnz, max_nz_per_row, + coo_rows.data(), coo_cols.data()); + + if (matrix->complete_graph != NULL) { + // matrix->complete_graph(matrix); + } /* * Add ams values that are needed elsewhere @@ -234,13 +250,11 @@ void EpetraCreateGomaProblemGraph(struct GomaLinearSolverData *ams, Exo_DB *exo, ; ams->npn_plus = dpi->num_universe_nodes; - ams->npu = num_internal_dofs[pg->imtrx] + num_boundary_dofs[pg->imtrx]; - ams->npu_plus = num_universe_dofs[pg->imtrx]; - - delete[] dblColGIDs; + ams->npu = num_internal_dofs + num_boundary_dofs; + ams->npu_plus = num_internal_dofs + num_boundary_dofs + num_external_dofs; int64_t num_unknowns; - int64_t my_unknowns = num_universe_dofs[pg->imtrx]; + int64_t my_unknowns = ams->npu_plus; int64_t num_nzz_global; int64_t my_nnz = nnz; @@ -250,27 +264,20 @@ void EpetraCreateGomaProblemGraph(struct GomaLinearSolverData *ams, Exo_DB *exo, DPRINTF(stdout, "\n%-30s= %ld\n", "Number of unknowns", num_unknowns); DPRINTF(stdout, "\n%-30s= %ld\n", "Number of matrix nonzeroes", num_nzz_global); + return GOMA_SUCCESS; } -/** - * Load local element contributions into the local matrix on this processor - * - * Modified from MSR version in mm_fill load_lec - * - * @param exo ptr to EXODUS II finite element mesh db - * @param ielem Element number we are working on - * @param ams Matrix contianer - * @param x Solution vector - * @param resid_vector residual vector - */ -void EpetraLoadLec(int ielem, struct GomaLinearSolverData *ams, double resid_vector[]) { +extern "C" goma_error GomaSparseMatrix_LoadLec(GomaSparseMatrix matrix, + int ielem, + struct Local_Element_Contributions *lec, + double resid_vector[]) { int e, v, i, j, pe, pv; int dofs; int gnn, row_index, ke, kv, nvdof; int col_index, ledof; int je_new; struct Element_Indices *ei_ptr; - std::vector Indices; + std::vector Indices; std::vector Values; for (e = V_FIRST; e < V_LAST; e++) { @@ -320,7 +327,7 @@ void EpetraLoadLec(int ielem, struct GomaLinearSolverData *ams, double resid_vec ei_ptr->Baby_Dolphin[v][j], ei_ptr->matID_ledof[ledof], pg->imtrx); GOMA_EH(col_index, "Bad var index."); - Indices.push_back(ams->GlobalIDs[col_index]); + Indices.push_back(matrix->global_ids[col_index]); Values.push_back(lec->J[LEC_J_INDEX(pe, pv, i, j)]); } } @@ -339,14 +346,14 @@ void EpetraLoadLec(int ielem, struct GomaLinearSolverData *ams, double resid_vec GOMA_EH(GOMA_ERROR, "LEC Indexing error"); } GOMA_EH(col_index, "Bad var index."); - Indices.push_back(ams->GlobalIDs[col_index]); + Indices.push_back(matrix->global_ids[col_index]); Values.push_back(lec->J[LEC_J_INDEX(pe, pv, i, j)]); } } } } - EpetraSumIntoGlobalRowMatrix(ams->RowMatrix, ams->GlobalIDs[row_index], - Indices.size(), &Values[0], &Indices[0]); + matrix->sum_into_row_values(matrix, matrix->global_ids[row_index], Indices.size(), + &Values[0], &Indices[0]); Indices.clear(); Values.clear(); } @@ -403,7 +410,7 @@ void EpetraLoadLec(int ielem, struct GomaLinearSolverData *ams, double resid_vec } } GOMA_EH(col_index, "Bad var index."); - Indices.push_back(ams->GlobalIDs[col_index]); + Indices.push_back(matrix->global_ids[col_index]); Values.push_back(lec->J[LEC_J_INDEX(pe, pv, i, j)]); } } @@ -432,14 +439,14 @@ void EpetraLoadLec(int ielem, struct GomaLinearSolverData *ams, double resid_vec GOMA_EH(GOMA_ERROR, "LEC Indexing error"); } GOMA_EH(col_index, "Bad var index."); - Indices.push_back(ams->GlobalIDs[col_index]); + Indices.push_back(matrix->global_ids[col_index]); Values.push_back(lec->J[LEC_J_INDEX(pe, pv, i, j)]); } } } } - EpetraSumIntoGlobalRowMatrix(ams->RowMatrix, ams->GlobalIDs[row_index], - Indices.size(), &Values[0], &Indices[0]); + matrix->sum_into_row_values(matrix, matrix->global_ids[row_index], Indices.size(), + &Values[0], &Indices[0]); Indices.clear(); Values.clear(); } @@ -448,50 +455,18 @@ void EpetraLoadLec(int ielem, struct GomaLinearSolverData *ams, double resid_vec } } } + return GOMA_SUCCESS; } -/** - * Inverse row sum scale, used by goma for scaling matrix and residual - * @param ams Aztec structure containing RowMatrix - * @param b b from Ax = b - * @param scale array for scale values to be placed for further usage (b[i] will equal old b[i] / - * scale[i]) - */ -void EpetraRowSumScale(struct GomaLinearSolverData *ams, double *b, double *scale) { - Epetra_Vector vector_scale(ams->RowMatrix->RowMatrixRowMap(), false); - ams->RowMatrix->InvRowSums(vector_scale); - ams->RowMatrix->LeftScale(vector_scale); - int local; - for (int i = 0; i < ams->RowMatrix->NumMyRows(); i++) { - local = ams->RowMatrix->RowMatrixRowMap().LID(ams->GlobalIDs[i]); - b[i] *= vector_scale[local]; - scale[i] = 1 / vector_scale[local]; +extern "C" goma_error GomaSparseMatrix_Destroy(GomaSparseMatrix *matrix) { + if (*matrix == NULL) { + return GOMA_SUCCESS; } -} - -/** - * Set the GlobalRow to zero and set the diagonal column to 1 in that row - * - * @param ams Aztec_Linear_Solver_System matrix struct containing epetra matrix - * @param GlobalRow global row to set to diagonal only - */ -void EpetraSetDiagonalOnly(struct GomaLinearSolverData *ams, int GlobalRow) { - Epetra_CrsMatrix *CrsMatrix = dynamic_cast(ams->RowMatrix); - int size = CrsMatrix->NumGlobalEntries(GlobalRow); - double *values = new double[size]; - int *indices = new int[size]; - int NumEntries; - CrsMatrix->ExtractGlobalRowCopy(GlobalRow, size, NumEntries, values, indices); - for (int i = 0; i < NumEntries; i++) { - if (indices[i] == GlobalRow) { - values[i] = 1; - } else { - values[i] = 0; - } + if ((*matrix)->destroy != NULL) { + (*matrix)->destroy(*matrix); } - CrsMatrix->ReplaceGlobalValues(GlobalRow, NumEntries, values, indices); - delete[] indices; - delete[] values; -} -} -/* End extern "C" */ + free((*matrix)->global_ids); + free(*matrix); + *matrix = NULL; + return GOMA_SUCCESS; +} \ No newline at end of file diff --git a/src/linalg/sparse_matrix_epetra.cpp b/src/linalg/sparse_matrix_epetra.cpp new file mode 100644 index 000000000..4f6d6050b --- /dev/null +++ b/src/linalg/sparse_matrix_epetra.cpp @@ -0,0 +1,180 @@ +#ifdef GOMA_ENABLE_EPETRA +#include "std.h" +#include +#include +#include +#include +#include +#include +#include +#ifdef EPETRA_MPI +#include "Epetra_MpiComm.h" +#else +#include "Epetra_SerialComm.h" +#endif +#include +#include +#include +#include +#include + +extern "C" { +#include "mm_eh.h" +#include "rf_io.h" +} +#include "linalg/sparse_matrix.h" +#include "linalg/sparse_matrix_epetra.h" + +using Teuchos::RCP; + +extern "C" goma_error GomaSparseMatrix_Epetra_Create(GomaSparseMatrix *matrix) { + EpetraSparseMatrix *tmp = new EpetraSparseMatrix(); + (*matrix)->type = GOMA_SPARSE_MATRIX_TYPE_EPETRA; + (*matrix)->data = reinterpret_cast(tmp); + (*matrix)->create_graph = g_epetra_create_graph; + (*matrix)->complete_graph = g_epetra_complete_graph; + (*matrix)->insert_row_values = g_epetra_insert_row_values; + (*matrix)->sum_into_row_values = g_epetra_sum_into_row_values; + (*matrix)->put_scalar = g_epetra_put_scalar; + (*matrix)->row_sum_scaling = g_epetra_row_sum_scaling; + (*matrix)->zero_global_row = g_epetra_zero_row; + (*matrix)->zero_global_row_set_diag = g_epetra_zero_row_set_diag; + (*matrix)->destroy = g_epetra_destroy; + return GOMA_SUCCESS; +} + +extern "C" goma_error g_epetra_create_graph(GomaSparseMatrix matrix, + GomaGlobalOrdinal n_rows, + GomaGlobalOrdinal *row_list, + GomaGlobalOrdinal n_cols, + GomaGlobalOrdinal *col_list, + GomaGlobalOrdinal local_nnz, + GomaGlobalOrdinal max_per_row, + GomaGlobalOrdinal *coo_rows, + GomaGlobalOrdinal *coo_cols) { + auto *tmp = static_cast(matrix->data); +#ifdef EPETRA_MPI + Epetra_MpiComm comm(MPI_COMM_WORLD); +#else + Epetra_SerialComm comm; +#endif + GomaGlobalOrdinal global_n_rows; + GomaGlobalOrdinal global_n_cols; + MPI_Allreduce(&n_rows, &global_n_rows, 1, MPI_GOMA_ORDINAL, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&n_cols, &global_n_cols, 1, MPI_GOMA_ORDINAL, MPI_SUM, MPI_COMM_WORLD); + + tmp->row_map = Teuchos::rcp(new Epetra_Map(global_n_rows, n_rows, row_list, 0, comm)); + tmp->crs_graph = Teuchos::rcp(new Epetra_CrsGraph(Copy, *(tmp->row_map), max_per_row)); + + GomaGlobalOrdinal current_row = row_list[0]; + std::vector indices(max_per_row); + std::vector values(max_per_row); + GomaGlobalOrdinal row_count = 0; + for (GomaGlobalOrdinal i = 0; i < local_nnz; i++) { + if (coo_rows[i] != current_row) { + tmp->crs_graph->InsertGlobalIndices(current_row, row_count, indices.data()); + current_row = coo_rows[i]; + row_count = 0; + } + indices[row_count] = coo_cols[i]; + row_count++; + } + tmp->crs_graph->InsertGlobalIndices(current_row, row_count, indices.data()); + + tmp->crs_graph->FillComplete(); + + tmp->matrix = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *(tmp->crs_graph))); + tmp->matrix->FillComplete(); + tmp->matrix->PutScalar(0.0); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_epetra_complete_graph(GomaSparseMatrix matrix) { + auto *tmp = static_cast(matrix->data); + tmp->matrix->FillComplete(); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_epetra_insert_row_values(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices) { + auto *tmp = static_cast(matrix->data); + GO global_row_t = static_cast(global_row); + tmp->matrix->InsertGlobalValues(global_row_t, num_entries, values, indices); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_epetra_sum_into_row_values(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices) { + auto *tmp = static_cast(matrix->data); + GO global_row_t = static_cast(global_row); + tmp->matrix->SumIntoGlobalValues(global_row_t, num_entries, values, indices); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_epetra_put_scalar(GomaSparseMatrix matrix, double scalar) { + auto *tmp = static_cast(matrix->data); + tmp->matrix->PutScalar(scalar); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_epetra_row_sum_scaling(GomaSparseMatrix matrix, double *b, double *scale) { + auto *tmp = static_cast(matrix->data); + tmp->matrix->RowMatrixRowMap(); + Epetra_Vector vector_scale(tmp->matrix->RowMatrixRowMap(), false); + tmp->matrix->InvRowSums(vector_scale); + tmp->matrix->LeftScale(vector_scale); + int local; + for (int i = 0; i < tmp->matrix->NumMyRows(); i++) { + local = tmp->matrix->RowMatrixRowMap().LID(matrix->global_ids[i]); + b[i] *= vector_scale[local]; + scale[i] = 1 / vector_scale[local]; + } + return GOMA_SUCCESS; +} + +extern "C" goma_error g_epetra_zero_row(GomaSparseMatrix matrix, GomaGlobalOrdinal global_row) { + auto *tmp = static_cast(matrix->data); + + int size = tmp->matrix->NumGlobalEntries(global_row); + std::vector values(size); + std::vector indices(size); + int NumEntries; + tmp->matrix->ExtractGlobalRowCopy(global_row, size, NumEntries, values.data(), indices.data()); + for (int i = 0; i < NumEntries; i++) { + values[i] = 0; + } + tmp->matrix->ReplaceGlobalValues(global_row, NumEntries, values.data(), indices.data()); + + return GOMA_SUCCESS; +} + +extern "C" goma_error g_epetra_zero_row_set_diag(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row) { + auto *tmp = static_cast(matrix->data); + int size = tmp->matrix->NumGlobalEntries(global_row); + std::vector values(size); + std::vector indices(size); + int NumEntries; + tmp->matrix->ExtractGlobalRowCopy(global_row, size, NumEntries, values.data(), indices.data()); + for (int i = 0; i < NumEntries; i++) { + if (indices[i] == global_row) { + values[i] = 1; + } else { + values[i] = 0; + } + } + tmp->matrix->ReplaceGlobalValues(global_row, NumEntries, values.data(), indices.data()); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_epetra_destroy(GomaSparseMatrix matrix) { + delete static_cast(matrix->data); + return GOMA_SUCCESS; +} +#endif // GOMA_ENABLE_EPETRA \ No newline at end of file diff --git a/src/linalg/sparse_matrix_tpetra.cpp b/src/linalg/sparse_matrix_tpetra.cpp new file mode 100644 index 000000000..36d49515c --- /dev/null +++ b/src/linalg/sparse_matrix_tpetra.cpp @@ -0,0 +1,213 @@ +#ifdef GOMA_ENABLE_TPETRA +#include "Tpetra_computeRowAndColumnOneNorms.hpp" +#include "std.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +extern "C" { +#include "mm_eh.h" +#include "rf_io.h" +} +#include "linalg/sparse_matrix.h" +#include "linalg/sparse_matrix_tpetra.h" + +using Teuchos::RCP; + +extern "C" goma_error GomaSparseMatrix_Tpetra_Create(GomaSparseMatrix *matrix) { + TpetraSparseMatrix *tmp = new TpetraSparseMatrix(); + (*matrix)->type = GOMA_SPARSE_MATRIX_TYPE_TPETRA; + (*matrix)->data = reinterpret_cast(tmp); + (*matrix)->create_graph = g_tpetra_create_graph; + (*matrix)->complete_graph = g_tpetra_complete_graph; + (*matrix)->insert_row_values = g_tpetra_insert_row_values; + (*matrix)->sum_into_row_values = g_tpetra_sum_into_row_values; + (*matrix)->put_scalar = g_tpetra_put_scalar; + (*matrix)->row_sum_scaling = g_tpetra_row_sum_scaling; + (*matrix)->zero_global_row = g_tpetra_zero_row; + (*matrix)->zero_global_row_set_diag = g_tpetra_zero_row_set_diag; + (*matrix)->destroy = g_tpetra_destroy; + return GOMA_SUCCESS; +} + +extern "C" goma_error g_tpetra_create_graph(GomaSparseMatrix matrix, + GomaGlobalOrdinal n_rows, + GomaGlobalOrdinal *row_list, + GomaGlobalOrdinal n_cols, + GomaGlobalOrdinal *col_list, + GomaGlobalOrdinal local_nnz, + GomaGlobalOrdinal max_per_row, + GomaGlobalOrdinal *coo_rows, + GomaGlobalOrdinal *coo_cols) { + auto *tmp = static_cast(matrix->data); + RCP> comm(new Teuchos::MpiComm(MPI_COMM_WORLD)); + GomaGlobalOrdinal global_n_rows; + GomaGlobalOrdinal global_n_cols; + MPI_Allreduce(&n_rows, &global_n_rows, 1, MPI_GOMA_ORDINAL, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&n_cols, &global_n_cols, 1, MPI_GOMA_ORDINAL, MPI_SUM, MPI_COMM_WORLD); + + Teuchos::ArrayView row_list_view(row_list, n_rows); + Teuchos::ArrayView col_list_view(col_list, n_cols); + + tmp->row_map = Teuchos::rcp(new Tpetra::Map(global_n_rows, row_list_view, 0, comm)); + tmp->col_map = Teuchos::rcp(new Tpetra::Map(global_n_rows, row_list_view, 0, comm)); + tmp->crs_graph = + Teuchos::rcp(new Tpetra::FECrsGraph(tmp->row_map, tmp->col_map, max_per_row)); + + GomaGlobalOrdinal current_row = row_list[0]; + std::vector indices(max_per_row); + std::vector values(max_per_row); + GomaGlobalOrdinal row_count = 0; + for (GomaGlobalOrdinal i = 0; i < local_nnz; i++) { + if (coo_rows[i] != current_row) { + tmp->crs_graph->insertGlobalIndices( + current_row, Teuchos::ArrayView(indices.data(), row_count)); + current_row = coo_rows[i]; + row_count = 0; + } + indices[row_count] = coo_cols[i]; + values[row_count] = coo_cols[i]; + row_count++; + } + tmp->crs_graph->insertGlobalIndices( + current_row, Teuchos::ArrayView(indices.data(), row_count)); + + if (Debug_Flag > 2) { + RCP fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); + tmp->crs_graph->describe(*fos, Teuchos::VERB_EXTREME); + } + tmp->crs_graph->fillComplete(); + + tmp->matrix = Teuchos::rcp(new Tpetra::FECrsMatrix((tmp->crs_graph))); + tmp->matrix->beginAssembly(); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_tpetra_complete_graph(GomaSparseMatrix matrix) { + auto *tmp = static_cast(matrix->data); + tmp->matrix->fillComplete(); + tmp->matrix->resumeFill(); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_tpetra_insert_row_values(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices) { + auto *tmp = static_cast(matrix->data); + GO global_row_t = static_cast(global_row); + tmp->matrix->insertGlobalValues(global_row_t, + Teuchos::ArrayView(indices, num_entries), + Teuchos::ArrayView(values, num_entries)); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_tpetra_sum_into_row_values(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row, + GomaGlobalOrdinal num_entries, + double *values, + GomaGlobalOrdinal *indices) { + auto *tmp = static_cast(matrix->data); + GO global_row_t = static_cast(global_row); + tmp->matrix->sumIntoGlobalValues(global_row_t, num_entries, values, indices); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_tpetra_put_scalar(GomaSparseMatrix matrix, double scalar) { + auto *tmp = static_cast(matrix->data); + tmp->matrix->setAllToScalar(scalar); + return GOMA_SUCCESS; +} + +extern "C" goma_error g_tpetra_row_sum_scaling(GomaSparseMatrix matrix, double *b, double *scale) { + auto *tmp = static_cast(matrix->data); + bool ended_assembly = false; + if (!tmp->matrix->isFillComplete()) { + tmp->matrix->endAssembly(); + ended_assembly = true; + } + RCP> b_vec = + Teuchos::rcp(new Tpetra::Vector(tmp->row_map)); + auto n_rows = tmp->matrix->getLocalNumRows(); + + using crs_t = Tpetra::CrsMatrix; + typename crs_t::local_inds_host_view_type indices; + typename crs_t::values_host_view_type values; + for (size_t i = 0; i < n_rows; i++) { + tmp->matrix->getLocalRowView(i, indices, values); + double row_sum = 0; + for (size_t j = 0; j < values.size(); j++) { + row_sum += std::abs(values[j]); + } + if (row_sum == 0) { + GOMA_WH_MANY(GOMA_ERROR, "Row sum is zero setting to 1.0, g_tpetra_row_sum_scaling"); + row_sum = 1.0; + } + b_vec->replaceLocalValue(i, 1.0 / row_sum); + } + tmp->matrix->leftScale(*b_vec); + if (ended_assembly) { + tmp->matrix->beginAssembly(); + } + // tmp->matrix->beginAssembly(); + auto data = b_vec->getData(); + for (size_t i = 0; i < n_rows; i++) { + auto local = tmp->matrix->getRowMap()->getLocalElement(matrix->global_ids[i]); + scale[local] = 1 / data[local]; + b[local] *= data[local]; + } + return GOMA_SUCCESS; +} + +extern "C" goma_error g_tpetra_zero_row(GomaSparseMatrix matrix, GomaGlobalOrdinal global_row) { + auto *tmp = static_cast(matrix->data); + using crs_t = Tpetra::CrsMatrix; + typename crs_t::nonconst_global_inds_host_view_type Indices; + typename crs_t::nonconst_values_host_view_type Values; + size_t NumEntries; + tmp->matrix->getGlobalRowCopy(global_row, Indices, Values, NumEntries); + if (NumEntries == Teuchos::OrdinalTraits::invalid()) { + GOMA_EH(GOMA_ERROR, "Global row does not exist on this processor, g_tptra_zero_row"); + } + for (size_t i = 0; i < NumEntries; i++) { + Values[i] = 0; + } + + return GOMA_SUCCESS; +} + +extern "C" goma_error g_tpetra_zero_row_set_diag(GomaSparseMatrix matrix, + GomaGlobalOrdinal global_row) { + auto *tmp = static_cast(matrix->data); + using crs_t = Tpetra::CrsMatrix; + typename crs_t::nonconst_global_inds_host_view_type Indices; + typename crs_t::nonconst_values_host_view_type Values; + size_t NumEntries; + tmp->matrix->getGlobalRowCopy(global_row, Indices, Values, NumEntries); + if (NumEntries == Teuchos::OrdinalTraits::invalid()) { + GOMA_EH(GOMA_ERROR, "Global row does not exist on this processor, g_tptra_zero_row"); + } + for (size_t i = 0; i < NumEntries; i++) { + if (Indices[i] == global_row) { + Values[i] = 1; + } else { + Values[i] = 0; + } + } + return GOMA_SUCCESS; +} + +extern "C" goma_error g_tpetra_destroy(GomaSparseMatrix matrix) { + delete static_cast(matrix->data); + return GOMA_SUCCESS; +} +#endif // GOMA_ENABLE_TPETRA \ No newline at end of file diff --git a/src/metis_decomp.c b/src/metis_decomp.c index 79ec3ed88..fef0004a1 100644 --- a/src/metis_decomp.c +++ b/src/metis_decomp.c @@ -455,12 +455,12 @@ goma_error goma_metis_decomposition(char **filenames, int n_files) { int *partitions = malloc(sizeof(int) * monolith->num_elems); idx_t n_parts = Num_Proc; - if (Decompose_Type == 1 || (Decompose_Type == 0 && n_parts < 8)) { + if (Decompose_Type == 1 || Decompose_Type == 0) { DPRINTF(stdout, "\nInternal METIS decomposition using Recursive Bisection.\n\n"); METIS_PartGraphRecursive(&monolith->num_elems, &n_con, elem_adj_pntr, elem_adj_list, vwgt, NULL, NULL, &n_parts, NULL, NULL, options, &edgecut, partitions); - } else if (Decompose_Type == 2 || (Decompose_Type == 0 && n_parts >= 8)) { + } else if (Decompose_Type == 2) { DPRINTF(stdout, "\nInternal METIS decomposition using KWAY.\n\n"); METIS_PartGraphKway(&monolith->num_elems, &n_con, elem_adj_pntr, elem_adj_list, vwgt, NULL, NULL, &n_parts, NULL, NULL, options, &edgecut, partitions); diff --git a/src/mm_fill.c b/src/mm_fill.c index fe5c863f9..f19853d1d 100644 --- a/src/mm_fill.c +++ b/src/mm_fill.c @@ -40,6 +40,7 @@ #include "el_elm_info.h" #include "el_geom.h" #include "exo_struct.h" +#include "linalg/sparse_matrix.h" #include "load_field_variables.h" #include "md_timer.h" #include "mm_as.h" @@ -93,7 +94,6 @@ #include "rf_node_const.h" #include "rf_solver.h" #include "rf_solver_const.h" -#include "sl_epetra_util.h" #include "sl_util.h" #include "sl_util_structs.h" #include "std.h" @@ -5000,8 +5000,9 @@ static void load_lec(Exo_DB *exo, /* ptr to EXODUS II finite element mesh db */ fprintf(rrrr, "\nGlobal_NN Proc_NN Equation idof Proc_SolnNum ResidValue\n"); #endif - if (strcmp(Matrix_Format, "epetra") == 0) { - EpetraLoadLec(ielem, ams, resid_vector); + if (ams->GomaMatrixData != NULL) { + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; + GomaSparseMatrix_LoadLec(matrix, ielem, lec, resid_vector); } #ifdef GOMA_ENABLE_PETSC #if PETSC_USE_COMPLEX diff --git a/src/mm_fill_ls.c b/src/mm_fill_ls.c index 9c25969b0..a9e866b91 100644 --- a/src/mm_fill_ls.c +++ b/src/mm_fill_ls.c @@ -36,6 +36,7 @@ #include "exo_conn.h" #include "exo_struct.h" #include "exodusII.h" +#include "linalg/sparse_matrix.h" #include "load_field_variables.h" #include "mm_as_alloc.h" #include "mm_fill_aux.h" @@ -89,7 +90,6 @@ #include "sl_util.h" #define GOMA_MM_FILL_LS_C -#include "sl_epetra_util.h" struct Extended_Shape_Fcn_Basics *xfem; /* This is a global structure for the basic pieces needed for XFEM */ @@ -10590,7 +10590,7 @@ void check_xfem_contribution( } } } - } else if (strcmp(Matrix_Format, "epetra") == 0) { + } else if (ams->GomaMatrixData != NULL) { for (irow = 0; irow < N; irow++) { eqn = idv[pg->imtrx][irow][0]; if (eqn == R_MASS || eqn == R_ENERGY) { @@ -10599,8 +10599,8 @@ void check_xfem_contribution( eps = eps_standard; } if (fabs(xfem->active_vol[irow]) < eps * xfem->tot_vol[irow]) { - - EpetraSetDiagonalOnly(ams, ams->GlobalIDs[irow]); + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; + matrix->zero_global_row_set_diag(matrix, matrix->global_ids[irow]); resid[irow] = x[irow] - x_old_static[irow]; if (FALSE && xfem->active_vol[irow] != 0.) /* debugging */ diff --git a/src/mm_fill_species.c b/src/mm_fill_species.c index c775c9e6f..cf8812198 100644 --- a/src/mm_fill_species.c +++ b/src/mm_fill_species.c @@ -2471,7 +2471,7 @@ void mass_flux_surf_SULFIDATION(dbl mass_flux[MAX_CONC], c_H2S = fv->c[0]; /* H2S is the 1st diffusing species */ c_O2 = fv->c[1]; /* O2 is the 2nd diffusion species */ mass_flux[wspec] = nu * k1 * exp(-E1 / R / T) * c_H2S * sqrt(c_O2); - } else if (mode == FULL) { + } else if (mode == METAL_CORROSION_FULL) { fprintf(stderr, "The full model has not yet implemented - awaits future efforts\n"); exit(1); } else if (mode == ANNIHILATION_ELECTRONEUTRALITY) { @@ -2542,7 +2542,7 @@ void mass_flux_surf_SULFIDATION(dbl mass_flux[MAX_CONC], } d_mass_flux[wspec][TEMPERATURE] = nu * k1 * (E1 / R / T / T) * exp(-E1 / R / T) * c_H2S * c_O2; - } else if (mode == FULL) { + } else if (mode == METAL_CORROSION_FULL) { fprintf(stderr, "The full model has not yet implemented - awaits future efforts\n"); exit(1); } else if (mode == ANNIHILATION_ELECTRONEUTRALITY) { @@ -7312,7 +7312,7 @@ void compute_leak_velocity(double *vnorm, } else if (mode == GAS_DIFFUSION) { StoiCoef[wspec] = -1.0; /* 1 mole of Cu2S produced */ /* per mole of H2S consumped */ - } else if (mode == FULL) { + } else if (mode == METAL_CORROSION_FULL) { fprintf(stderr, "The full model has not yet implemented - awaits future efforts\n"); exit(1); } diff --git a/src/mm_fill_stress.c b/src/mm_fill_stress.c index 777b228b1..9cb1de717 100644 --- a/src/mm_fill_stress.c +++ b/src/mm_fill_stress.c @@ -6006,8 +6006,8 @@ void load_neighbor_pointers(Exo_DB *exo, } } } - } else if (strcmp(Matrix_Format, "epetra") == 0) { - GOMA_EH(GOMA_ERROR, "load_neighbor_pointers unsupported by epetra"); + } else { + GOMA_EH(GOMA_ERROR, "load_neighbor_pointers unsupported by %s", Matrix_Format); } } /***************************************************************************/ diff --git a/src/mm_input.c b/src/mm_input.c index e8cfa9c51..2c8db9eca 100644 --- a/src/mm_input.c +++ b/src/mm_input.c @@ -5832,6 +5832,7 @@ void rd_solver_specs(FILE *ifp, char *input) { strcpy(Matrix_Absolute_Threshold, "0"); strcpy(Matrix_Reorder, "none"); strcpy(Amesos_Package, "KLU"); + strcpy(Amesos2_Package, "KLU2"); /* Read in Solver specifications */ @@ -5956,6 +5957,9 @@ void rd_solver_specs(FILE *ifp, char *input) { } else if (strcmp(Matrix_Solver, "amesos") == 0) { Linear_Solver = AMESOS; is_Solver_Serial = FALSE; + } else if (strcmp(Matrix_Solver, "amesos2") == 0) { + Linear_Solver = AMESOS2; + is_Solver_Serial = FALSE; } else if (strcmp(Matrix_Solver, "aztecoo") == 0) { Linear_Solver = AZTECOO; is_Solver_Serial = FALSE; @@ -5996,6 +6000,10 @@ void rd_solver_specs(FILE *ifp, char *input) { strcpy(Matrix_Format, "petsc_complex"); /* save string for aztec use */ strcpy(search_string, "Matrix storage format"); snprintf(echo_string, MAX_CHAR_ECHO_INPUT, eoformat, search_string, Matrix_Format); + } else if (strcmp(Matrix_Solver, "amesos2") == 0) { + strcpy(Matrix_Format, "tpetra"); /* save string for aztec use */ + strcpy(search_string, "Matrix storage format"); + snprintf(echo_string, MAX_CHAR_ECHO_INPUT, eoformat, search_string, Matrix_Format); } else if (strcmp(Matrix_Solver, "front") != 0) { /* Read in Matrix Format only if we are going to assemble one */ @@ -6089,6 +6097,27 @@ void rd_solver_specs(FILE *ifp, char *input) { ECHO(echo_string, echo_file); } + for (int i = 1; i < MAX_NUM_MATRICES; i++) { + strcpy(Stratimikos_File[i], Stratimikos_File[0]); + } + + strcpy(search_string, "Amesos2 File"); + iread = look_for_optional(ifp, search_string, input, '='); + if (iread == 1) { + read_string(ifp, input, '\n'); + strip(input); + strcpy(Amesos2_File[0], input); + snprintf(echo_string, MAX_CHAR_ECHO_INPUT, eoformat, search_string, input); + ECHO(echo_string, echo_file); + } else { + // Set stratimikos.xml as defualt stratimikos file + strcpy(Amesos2_File[0], ""); + } + + for (int i = 1; i < MAX_NUM_MATRICES; i++) { + strcpy(Amesos2_File[i], Amesos2_File[0]); + } + strcpy(search_string, "Preconditioner"); iread = look_for_optional(ifp, search_string, input, '='); @@ -6435,6 +6464,19 @@ void rd_solver_specs(FILE *ifp, char *input) { ECHO(echo_string, echo_file); } + strcpy(search_string, "Amesos2 Solver Package"); + iread = look_for_optional(ifp, search_string, input, '='); + if (iread == 1) { + read_string(ifp, input, '\n'); + strip(input); + strcpy(Amesos2_Package, input); + snprintf(echo_string, MAX_CHAR_ECHO_INPUT, eoformat, search_string, input); + ECHO(echo_string, echo_file); + } else { + snprintf(echo_string, MAX_CHAR_ECHO_INPUT, def_form, search_string, "KLU", default_string); + ECHO(echo_string, echo_file); + } + /* first initialize modified newton parameter to false */ modified_newton = FALSE; @@ -8099,9 +8141,6 @@ void rd_eq_specs(FILE *ifp, char *input, const int mn) { snprintf(echo_string, MAX_CHAR_ECHO_INPUT, "Stratimikos file = %s for matrix %d", input, mtrx_index1); ECHO(echo_string, echo_file); - } else { - // Set stratimikos.xml as defualt stratimikos file - strcpy(Stratimikos_File[imtrx], "stratimikos.xml"); } iread = look_forward_optional_until(ifp, "Normalized Residual Tolerance", "MATRIX", input, '='); diff --git a/src/mm_input_bc.c b/src/mm_input_bc.c index fe9c28e3d..7855bdfb7 100644 --- a/src/mm_input_bc.c +++ b/src/mm_input_bc.c @@ -1270,7 +1270,7 @@ void rd_bc_specs(FILE *ifp, char *input) { } else if (!strcmp(input, "GAS_DIFFUSION")) { BC_Types[ibc].BC_Data_Int[2] = GAS_DIFFUSION; } else if (!strcmp(input, "FULL")) { - BC_Types[ibc].BC_Data_Int[2] = FULL; + BC_Types[ibc].BC_Data_Int[2] = METAL_CORROSION_FULL; } else if (!strcmp(input, "ANNIHILATION_ELECTRONEUTRALITY")) { BC_Types[ibc].BC_Data_Int[2] = ANNIHILATION_ELECTRONEUTRALITY; } else if (!strcmp(input, "ANNIHILATION")) { diff --git a/src/mm_sol_nonlinear.c b/src/mm_sol_nonlinear.c index 3dd870cf8..ec2fb4203 100644 --- a/src/mm_sol_nonlinear.c +++ b/src/mm_sol_nonlinear.c @@ -17,9 +17,9 @@ */ #include "bc_contact.h" +#include "linalg/sparse_matrix.h" #include "mm_eh.h" #include "mm_mp_const.h" -#include "sl_epetra_interface.h" #include "sl_util_structs.h" #define GOMA_MM_SOL_NONLINEAR_C @@ -71,6 +71,7 @@ #include "rf_solver.h" #include "rf_solver_const.h" #include "rf_util.h" +#include "sl_amesos2_interface.h" #include "sl_amesos_interface.h" #include "sl_auxutil.h" #include "sl_aztecoo_interface.h" @@ -767,8 +768,9 @@ int solve_nonlinear_problem(struct GomaLinearSolverData *ams, init_vec_value(resid_vector, 0.0, numProcUnknowns); init_vec_value(delta_x, 0.0, numProcUnknowns); /* Zero matrix values */ - if (strcmp(Matrix_Format, "epetra") == 0) { - EpetraPutScalarRowMatrix(ams->RowMatrix, 0.0); + if (ams->GomaMatrixData != NULL) { + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; + matrix->put_scalar(matrix, 0.0); } else { init_vec_value(a, 0.0, ams->nnz); } @@ -1415,6 +1417,18 @@ int solve_nonlinear_problem(struct GomaLinearSolverData *ams, amesos_solve(Amesos_Package, ams, delta_x, resid_vector, 1, pg->imtrx); strcpy(stringer, " 1 "); break; + case AMESOS2: + + if (ams->GomaMatrixData != NULL) { + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; + if (matrix->type != GOMA_SPARSE_MATRIX_TYPE_TPETRA) { + GOMA_EH(GOMA_ERROR, " Sorry, only Tpetra matrix formats are currently supported with " + "the Amesos2 solver suite\n"); + } + } + amesos2_solve(ams, delta_x, resid_vector, Amesos2_Package, Amesos2_File[pg->imtrx]); + strcpy(stringer, " 1 "); + break; case AZTECOO: if (strcmp(Matrix_Format, "epetra") == 0) { @@ -1474,6 +1488,15 @@ int solve_nonlinear_problem(struct GomaLinearSolverData *ams, check_parallel_error("Error in solve - stratimikos"); } aztec_stringer(AZ_normal, iterations, &stringer[0]); + } else if (strcmp(Matrix_Format, "tpetra") == 0) { + int iterations; + int err = stratimikos_solve_tpetra(ams, delta_x, resid_vector, &iterations, + Stratimikos_File, pg->imtrx); + if (err) { + GOMA_EH(err, "Error in stratimikos solve"); + check_parallel_error("Error in solve - stratimikos"); + } + aztec_stringer(AZ_normal, iterations, &stringer[0]); } else { GOMA_EH(GOMA_ERROR, "Sorry, only Epetra matrix formats are currently supported with the Stratimikos " @@ -1626,6 +1649,15 @@ int solve_nonlinear_problem(struct GomaLinearSolverData *ams, } else { aztec_stringer(AZ_normal, iterations, &stringer[0]); } + } else if (strcmp(Matrix_Format, "tpetra") == 0) { + int iterations; + int err = stratimikos_solve_tpetra(ams, &wAC[iAC][0], &bAC[iAC][0], &iterations, + Stratimikos_File, pg->imtrx); + if (err) { + GOMA_EH(err, "Error in stratimikos solve"); + check_parallel_error("Error in solve - stratimikos"); + } + aztec_stringer(AZ_normal, iterations, &stringer[0]); } else { GOMA_EH( GOMA_ERROR, @@ -3449,6 +3481,15 @@ static int soln_sens(double lambda, /* parameter */ } else { aztec_stringer(AZ_normal, iterations, &stringer[0]); } + } else if (strcmp(Matrix_Format, "tpetra") == 0) { + int iterations; + int err = stratimikos_solve_tpetra(ams, x_sens, resid_vector_sens, &iterations, + Stratimikos_File, pg->imtrx); + if (err) { + GOMA_EH(err, "Error in stratimikos solve"); + check_parallel_error("Error in solve - stratimikos"); + } + aztec_stringer(AZ_normal, iterations, &stringer[0]); } else { GOMA_EH(GOMA_ERROR, "Sorry, only Epetra matrix formats are currently supported with the Stratimikos " diff --git a/src/rf_solve.c b/src/rf_solve.c index a08c9f85a..d558f6ac1 100644 --- a/src/rf_solve.c +++ b/src/rf_solve.c @@ -40,6 +40,7 @@ #include "el_elm_info.h" #include "el_geom.h" #include "exo_struct.h" +#include "linalg/sparse_matrix.h" #include "mm_as.h" #include "mm_as_structs.h" #include "mm_augc_util.h" @@ -71,14 +72,13 @@ #include "rf_io.h" #include "rf_io_const.h" #include "rf_io_structs.h" +#include "rf_masks.h" #include "rf_mp.h" #include "rf_node_const.h" #include "rf_solve_segregated.h" #include "rf_solver.h" #include "rf_util.h" #include "sl_auxutil.h" -#include "sl_epetra_interface.h" -#include "sl_epetra_util.h" #include "sl_matrix_util.h" #include "sl_petsc.h" #include "sl_petsc_complex.h" @@ -718,14 +718,21 @@ void solve_problem(Exo_DB *exo, /* ptr to the finite element mesh database */ /* Allocate sparse matrix */ - if (strcmp(Matrix_Format, "epetra") == 0) { + ams[JAC]->GomaMatrixData = NULL; + if ((strcmp(Matrix_Format, "tpetra") == 0) || (strcmp(Matrix_Format, "epetra") == 0)) { err = check_compatible_solver(); - GOMA_EH(err, - "Incompatible matrix solver for epetra, epetra supports amesos and aztecoo solvers."); + GOMA_EH(err, "Incompatible matrix solver for tpetra, tpetra supports stratimikos"); check_parallel_error("Matrix format / Solver incompatibility"); - ams[JAC]->RowMatrix = - EpetraCreateRowMatrix(num_internal_dofs[pg->imtrx] + num_boundary_dofs[pg->imtrx]); - EpetraCreateGomaProblemGraph(ams[JAC], exo, dpi); + GomaSparseMatrix goma_matrix; + goma_error err = GomaSparseMatrix_CreateFromFormat(&goma_matrix, Matrix_Format); + GOMA_EH(err, "GomaSparseMatrix_CreateFromFormat"); + int local_nodes = Num_Internal_Nodes + Num_Border_Nodes + Num_External_Nodes; + err = GomaSparseMatrix_SetProblemGraph( + goma_matrix, num_internal_dofs[pg->imtrx], num_boundary_dofs[pg->imtrx], + num_external_dofs[pg->imtrx], local_nodes, Nodes, MaxVarPerNode, Matilda, Inter_Mask, exo, + dpi, cx[pg->imtrx], pg->imtrx, Debug_Flag, ams[JAC]); + GOMA_EH(err, "GomaSparseMatrix_SetProblemGraph"); + ams[JAC]->GomaMatrixData = goma_matrix; #ifdef GOMA_ENABLE_PETSC #if PETSC_USE_COMPLEX } else if (strcmp(Matrix_Format, "petsc_complex") == 0) { @@ -788,8 +795,6 @@ void solve_problem(Exo_DB *exo, /* ptr to the finite element mesh database */ ams[JAC]->nnz = ija[num_internal_dofs[pg->imtrx] + num_boundary_dofs[pg->imtrx]] - 1; ams[JAC]->nnz_plus = ija[num_universe_dofs[pg->imtrx]]; - ams[JAC]->RowMatrix = NULL; - } else if (strcmp(Matrix_Format, "vbr") == 0) { log_msg("alloc_VBR_sparse_arrays..."); alloc_VBR_sparse_arrays(ams[JAC], exo, dpi); @@ -2650,6 +2655,7 @@ void solve_problem(Exo_DB *exo, /* ptr to the finite element mesh database */ for (i = 0; i < MAX_NUMBER_MATLS; i++) { for (n = 0; n < MAX_MODES; n++) { safer_free((void **)&(ve_glob[i][n]->gn)); + safer_free((void **)&(ve_glob[i][n]->time_const_st)); safer_free((void **)&(ve_glob[i][n])); } safer_free((void **)&(vn_glob[i])); @@ -2664,14 +2670,10 @@ void solve_problem(Exo_DB *exo, /* ptr to the finite element mesh database */ } if (last_call) { - if (strcmp(Matrix_Format, "epetra") == 0) { - EpetraDeleteRowMatrix(ams[JAC]->RowMatrix); - if (ams[JAC]->GlobalIDs != NULL) { - free(ams[JAC]->GlobalIDs); - } - } + GomaSparseMatrix goma_matrix = ams[JAC]->GomaMatrixData; + GomaSparseMatrix_Destroy(&goma_matrix); #ifdef GOMA_ENABLE_PETSC - else if (strcmp(Matrix_Format, "petsc") == 0) { + if (strcmp(Matrix_Format, "petsc") == 0) { err = goma_petsc_free_matrix(ams[JAC]); GOMA_EH(err, "free petsc matrix"); } diff --git a/src/rf_solve_segregated.c b/src/rf_solve_segregated.c index aec0757df..61e9f19e9 100644 --- a/src/rf_solve_segregated.c +++ b/src/rf_solve_segregated.c @@ -25,6 +25,7 @@ #include "dpi.h" #include "el_geom.h" #include "exo_struct.h" +#include "linalg/sparse_matrix.h" #include "mm_as.h" #include "mm_as_structs.h" #include "mm_augc_util.h" @@ -50,13 +51,12 @@ #include "rf_io.h" #include "rf_io_const.h" #include "rf_io_structs.h" +#include "rf_masks.h" #include "rf_mp.h" #include "rf_node_const.h" #include "rf_solve.h" #include "rf_solver.h" #include "rf_util.h" -#include "sl_epetra_interface.h" -#include "sl_epetra_util.h" #include "sl_matrix_util.h" #include "sl_petsc.h" #include "sl_petsc_complex.h" @@ -552,16 +552,26 @@ void solve_problem_segregated(Exo_DB *exo, /* ptr to the finite element mesh dat a = malloc(upd->Total_Num_Matrices * sizeof(double *)); a_old = malloc(upd->Total_Num_Matrices * sizeof(double *)); - if (strcmp(Matrix_Format, "epetra") == 0) { + if ((strcmp(Matrix_Format, "tpetra") == 0) || (strcmp(Matrix_Format, "epetra") == 0)) { err = check_compatible_solver(); - GOMA_EH(err, "Incompatible matrix solver for epetra, epetra supports amesos and " - "aztecoo solvers."); + GOMA_EH(err, "Incompatible matrix solver for tpetra, tpetra supports stratimikos"); check_parallel_error("Matrix format / Solver incompatibility"); + for (pg->imtrx = 0; pg->imtrx < upd->Total_Num_Matrices; pg->imtrx++) { - ams[pg->imtrx]->RowMatrix = - EpetraCreateRowMatrix(num_internal_dofs[pg->imtrx] + num_boundary_dofs[pg->imtrx]); - EpetraCreateGomaProblemGraph(ams[pg->imtrx], exo, dpi); + GomaSparseMatrix goma_matrix; + goma_error err = GomaSparseMatrix_CreateFromFormat(&goma_matrix, Matrix_Format); + GOMA_EH(err, "GomaSparseMatrix_CreateFromFormat"); + int local_nodes = Num_Internal_Nodes + Num_Border_Nodes + Num_External_Nodes; + err = GomaSparseMatrix_SetProblemGraph( + goma_matrix, num_internal_dofs[pg->imtrx], num_boundary_dofs[pg->imtrx], + num_external_dofs[pg->imtrx], local_nodes, Nodes, MaxVarPerNode, Matilda, Inter_Mask, exo, + dpi, cx[pg->imtrx], pg->imtrx, Debug_Flag, ams[JAC]); + GOMA_EH(err, "GomaSparseMatrix_SetProblemGraph"); + ams[pg->imtrx]->GomaMatrixData = goma_matrix; + ams[pg->imtrx]->npu = num_internal_dofs[pg->imtrx] + num_boundary_dofs[pg->imtrx]; + ams[pg->imtrx]->npu_plus = num_universe_dofs[pg->imtrx]; } + #ifdef GOMA_ENABLE_PETSC #if !(PETSC_USE_COMPLEX) } else if (strcmp(Matrix_Format, "petsc") == 0) { diff --git a/src/sl_amesos2_interface.cpp b/src/sl_amesos2_interface.cpp new file mode 100644 index 000000000..4af582f42 --- /dev/null +++ b/src/sl_amesos2_interface.cpp @@ -0,0 +1,132 @@ +#ifdef GOMA_ENABLE_AMESOS2 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "linalg/sparse_matrix.h" +#include "linalg/sparse_matrix_tpetra.h" +#include "sl_amesos2_interface.h" + +using MAT = Tpetra::CrsMatrix; +using VEC = Tpetra::Vector; +using MV = Tpetra::MultiVector; + +struct Amesos2_Solver_Data { + Teuchos::RCP> solver; + + Amesos2_Solver_Data() { solver = Teuchos::null; } +}; + +extern "C" void amesos2_solver_destroy(struct GomaLinearSolverData *ams) { + auto solver_data = static_cast(ams->SolverData); + delete solver_data; +} + +extern "C" { +int amesos2_solve(struct GomaLinearSolverData *ams, + double *x_, + double *b_, + char *amesos2_solver, + char *amesos2_file) { + using Teuchos::RCP; + auto matrix = static_cast(ams->GomaMatrixData); + auto *tpetra_data = static_cast(matrix->data); + bool success = true; + bool verbose = true; + + if (ams->SolverData == NULL) { + ams->SolverData = new Amesos2_Solver_Data(); + ams->DestroySolverData = amesos2_solver_destroy; + } + auto solver_data = static_cast(ams->SolverData); + + try { + if (!tpetra_data->matrix->isFillComplete()) { + tpetra_data->matrix->endAssembly(); + } + + RCP tpetra_x = rcp(new VEC(tpetra_data->matrix->getDomainMap())); + RCP tpetra_b = rcp(new VEC(tpetra_data->matrix->getRangeMap())); + + for (size_t i = 0; i < tpetra_x->getLocalLength(); i++) { + tpetra_x->replaceGlobalValue(matrix->global_ids[i], x_[i]); + } + for (size_t i = 0; i < tpetra_b->getLocalLength(); i++) { + tpetra_b->replaceGlobalValue(matrix->global_ids[i], b_[i]); + } + + if (solver_data->solver.is_null()) { + Teuchos::RCP crs_matrix = Teuchos::rcp_dynamic_cast(tpetra_data->matrix); + solver_data->solver = Amesos2::create("Mumps", crs_matrix); + if (amesos2_file != NULL && strlen(amesos2_file) > 0) { + std::filesystem::path path(amesos2_file); + Teuchos::RCP amesos2_params; + if (path.extension() == ".yaml") { + amesos2_params = Teuchos::getParametersFromYamlFile(amesos2_file); + } else { + amesos2_params = Teuchos::getParametersFromXmlFile(amesos2_file); + } + solver_data->solver->setParameters(amesos2_params); + } + + solver_data->solver->symbolicFactorization(); + } + solver_data->solver->numericFactorization(); + + solver_data->solver->solve(tpetra_x.ptr(), tpetra_b.ptr()); + + Teuchos::RCP outstream = Teuchos::VerboseObjectBase::getDefaultOStream(); + + /* Convert solution vector */ + int NumMyRows = matrix->n_rows; + + auto x_data = tpetra_x->getData(); + for (int i = 0; i < NumMyRows; i++) { + x_[i] = x_data[i]; + } + tpetra_data->matrix->beginAssembly(); + } + TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success) + if (success) { + return 0; + } + + return 1; +} + +} /* End extern "C" */ + +#else /* GOMA_ENABLE_AMESOS2 */ + +#include "mpi.h" + +extern "C" { +#include "mm_eh.h" +#include "std.h" + +int amesos2_solve(struct GomaLinearSolverData *ams, + double *x_, + double *b_, + int *iterations, + char *stratimikos_file) { + GOMA_EH(GOMA_ERROR, "Not built with Amesos2 support!"); + return -1; +} +} +#endif /* GOMA_ENABLE_AMESOS2 */ diff --git a/src/sl_amesos_interface.cpp b/src/sl_amesos_interface.cpp index 437db2ddf..3f2fcda29 100644 --- a/src/sl_amesos_interface.cpp +++ b/src/sl_amesos_interface.cpp @@ -48,7 +48,7 @@ #include "Epetra_SerialComm.h" #endif -#if 0 +#if 1 #include "EpetraExt_RowMatrixOut.h" #include "EpetraExt_VectorOut.h" #endif @@ -60,6 +60,8 @@ #include "Epetra_Map.h" #include "Epetra_Vector.h" #include "Trilinos_Util.h" +#include "linalg/sparse_matrix.h" +#include "linalg/sparse_matrix_epetra.h" #include "sl_util_structs.h" static void GomaMsr2EpetraCsr(struct GomaLinearSolverData *ams, Epetra_CrsMatrix *A, int newmatrix); @@ -84,7 +86,7 @@ void amesos_solve(char *choice, Amesos A_Factory; /* Convert to Epetra format */ - if (ams->RowMatrix == 0) { + if (ams->GomaMatrixData == NULL) { if (!ams->solveSetup) { if (A[imtrx] != nullptr) delete A[imtrx]; @@ -92,15 +94,18 @@ void amesos_solve(char *choice, } GomaMsr2EpetraCsr(ams, A[imtrx], !ams->solveSetup); } else { - A[imtrx] = dynamic_cast(ams->RowMatrix); + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; + EpetraSparseMatrix *epetra_matrix = static_cast(matrix->data); + A[imtrx] = dynamic_cast(epetra_matrix->matrix.get()); } const Epetra_Map &map = A[imtrx]->RowMatrixRowMap(); Epetra_Vector x(Copy, map, x_); Epetra_Vector b(Copy, map, b_); #if 0 - EpetraExt::RowMatrixToMatrixMarketFile("Jnext.mm", *A); - EpetraExt::VectorToMatrixMarketFile("bnext.mm", b); + EpetraExt::RowMatrixToMatrixMarketFile("Jep.mm", *A[imtrx]); + EpetraExt::VectorToMatrixMarketFile("xep.mm", x); + EpetraExt::VectorToMatrixMarketFile("bep.mm", b); MPI_Finalize(); exit(0); #endif diff --git a/src/sl_aztecoo_interface.cpp b/src/sl_aztecoo_interface.cpp index 51b31c2fd..f12f00351 100644 --- a/src/sl_aztecoo_interface.cpp +++ b/src/sl_aztecoo_interface.cpp @@ -11,6 +11,8 @@ * This software is distributed under the GNU General Public License. * * See LICENSE file. * \************************************************************************/ +#include "linalg/sparse_matrix.h" +#include "linalg/sparse_matrix_epetra.h" #if defined(PARALLEL) && !defined(EPETRA_MPI) #define EPETRA_MPI #endif @@ -49,6 +51,8 @@ extern "C" { * @param b_ residual vector */ void aztecoo_solve_epetra(struct GomaLinearSolverData *ams, double *x_, double *b_) { + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; + EpetraSparseMatrix *epetra_matrix = static_cast(matrix->data); /* Initialize MPI communications */ #ifdef EPETRA_MPI @@ -58,7 +62,7 @@ void aztecoo_solve_epetra(struct GomaLinearSolverData *ams, double *x_, double * #endif /* Define internal variables */ - Epetra_RowMatrix *A = ams->RowMatrix; + Epetra_RowMatrix *A = epetra_matrix->matrix.get(); Epetra_LinearProblem Problem; Epetra_Map map = A->RowMatrixRowMap(); diff --git a/src/sl_epetra_interface.cpp b/src/sl_epetra_interface.cpp deleted file mode 100644 index a557011d8..000000000 --- a/src/sl_epetra_interface.cpp +++ /dev/null @@ -1,121 +0,0 @@ -#ifndef __cplusplus -#define __cplusplus -#endif - -#if defined(PARALLEL) && !defined(EPETRA_MPI) -#define EPETRA_MPI -#endif - -#include - -#include "Epetra_CrsMatrix.h" -#include "Epetra_DataAccess.h" -#include "Epetra_Map.h" -#include "mpi.h" - -class Epetra_RowMatrix; - -#ifdef EPETRA_MPI -#include "Epetra_MpiComm.h" -#else -#include "Epetra_SerialComm.h" -#endif - -#include "sl_epetra_interface.h" - -/** - * Create a row matrix with a row map of size NumberProcRows - * @param NumberProcRows the number of rows on this processor - * @return C compatible pointer to row matrix - */ -C_Epetra_RowMatrix_t *EpetraCreateRowMatrix(int NumberProcRows) { -#ifdef EPETRA_MPI - Epetra_MpiComm comm(MPI_COMM_WORLD); -#else - Epetra_SerialComm comm; -#endif - Epetra_RowMatrix *AMatrix; - Epetra_Map RowMap(-1, NumberProcRows, 0, comm); - - AMatrix = new Epetra_CrsMatrix(Copy, RowMap, 0); - - return AMatrix; -} - -/** - * Insert values into an epetra row matrix (C interface) - * - * Insert if not filled - * - * Replace if filled - * - * @param AMatrix Matrix to insert values - * @param GlobalRow Global Row - * @param NumEntries Number of Entries - * @param Values Values to insert - * @param Indices Indices of values - */ -void EpetraInsertGlobalRowMatrix(C_Epetra_RowMatrix_t *AMatrix, - int GlobalRow, - int NumEntries, - const double *Values, - const int *Indices) { - Epetra_CrsMatrix *CrsMatrix = dynamic_cast(AMatrix); - if (CrsMatrix->Filled()) { - CrsMatrix->ReplaceGlobalValues(GlobalRow, NumEntries, Values, Indices); - } else { - CrsMatrix->InsertGlobalValues(GlobalRow, NumEntries, Values, Indices); - } -} - -/** - * Sum into global values for epetra row matrix (C interface) - * @param AMatrix Matrix to sum into values - * @param GlobalRow Global Row - * @param NumEntries Number of Entries - * @param Values Values to sum - * @param Indices Indices for values - */ -void EpetraSumIntoGlobalRowMatrix(C_Epetra_RowMatrix_t *AMatrix, - int GlobalRow, - int NumEntries, - const double *Values, - const int *Indices) { - Epetra_CrsMatrix *CrsMatrix = dynamic_cast(AMatrix); - int ierr = CrsMatrix->SumIntoGlobalValues(GlobalRow, NumEntries, Values, Indices); - if (ierr) - printf("Error in sum into epetra %d\n", ierr); -} - -/** - * Set an Epetra Row Matrix to a specified scalar value for all NZ entries (C interface) - * @param AMatrix Row matrix to set - * @param scalar scalar value to set - */ -void EpetraPutScalarRowMatrix(C_Epetra_RowMatrix_t *AMatrix, double scalar) { - Epetra_CrsMatrix *CrsMatrix = dynamic_cast(AMatrix); - CrsMatrix->PutScalar(scalar); -} - -/** - * Call fill complete on an Epetra Row Matrix (C interface) - * @param AMatrix Matrix to fill completely - */ -void EpetraFillCompleteRowMatrix(C_Epetra_RowMatrix_t *AMatrix) { - Epetra_CrsMatrix *CrsMatrix = dynamic_cast(AMatrix); - CrsMatrix->FillComplete(); -} - -int EpetraExtractRowValuesRowMatrix( - C_Epetra_RowMatrix_t *AMatrix, int GlobalRow, int size, double *values, int *indices) { - Epetra_CrsMatrix *CrsMatrix = dynamic_cast(AMatrix); - int NumEntries; - CrsMatrix->ExtractGlobalRowCopy(GlobalRow, size, NumEntries, values, indices); - return NumEntries; -} - -/** - * Call deconstructor for row matrix (C interface) - * @param AMatrix Row Matrix to delete - */ -void EpetraDeleteRowMatrix(C_Epetra_RowMatrix_t *AMatrix) { delete AMatrix; } diff --git a/src/sl_matrix_util.c b/src/sl_matrix_util.c index 6b4965fdf..9a9107581 100644 --- a/src/sl_matrix_util.c +++ b/src/sl_matrix_util.c @@ -23,6 +23,7 @@ #include "dpi.h" #include "el_geom.h" #include "exo_struct.h" +#include "linalg/sparse_matrix.h" #include "mm_as.h" #include "mm_as_structs.h" #include "mm_eh.h" @@ -35,7 +36,6 @@ #include "rf_solver.h" #include "rf_solver_const.h" #include "rf_vars_const.h" -#include "sl_epetra_util.h" #include "sl_matrix_util.h" #include "sl_util.h" #include "sl_util_structs.h" @@ -444,8 +444,9 @@ void row_sum_scaling_scale(struct GomaLinearSolverData *ams, double b[], double } else if (strcmp(Matrix_Format, "vbr") == 0) { row_sum_scale_VBR(ams->npn, ams->val, ams->bpntr, ams->bindx, ams->indx, ams->rpntr, ams->cpntr, b, scale); - } else if (strcmp(Matrix_Format, "epetra") == 0) { - row_sum_scale_epetra(ams, b, scale); + } else if (ams->GomaMatrixData != NULL) { + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; + matrix->row_sum_scaling(matrix, b, scale); #ifdef GOMA_ENABLE_PETSC #if PETSC_USE_COMPLEX } else if (strcmp(Matrix_Format, "petsc_complex") == 0) { @@ -661,10 +662,6 @@ void row_sum_scale_VBR(int N, } } /* END of routine row_sum_scale_VBR */ -void row_sum_scale_epetra(struct GomaLinearSolverData *ams, double *b, double *scale) { - EpetraRowSumScale(ams, b, scale); -} - /******************************************************************************/ /******************************************************************************/ /******************************************************************************/ /* @@ -820,7 +817,15 @@ void vector_scaling(const int N, double b[], double scale[]) { * @return 0 if compatible, -1 if not */ int check_compatible_solver(void) { - if (strcmp(Matrix_Format, "epetra") == 0) { + if (strcmp(Matrix_Format, "tpetra") == 0) { + switch (Linear_Solver) { + case STRATIMIKOS: + case AMESOS2: + return GOMA_SUCCESS; + default: + return GOMA_ERROR; + } + } else if (strcmp(Matrix_Format, "epetra") == 0) { switch (Linear_Solver) { case AZTECOO: case AMESOS: diff --git a/src/sl_stratimikos_interface.cpp b/src/sl_stratimikos_interface.cpp index 451f73fb2..e3a96ee3d 100644 --- a/src/sl_stratimikos_interface.cpp +++ b/src/sl_stratimikos_interface.cpp @@ -1,5 +1,7 @@ +#include #ifdef GOMA_ENABLE_STRATIMIKOS +#include #include #include #include @@ -17,6 +19,7 @@ #include "Teuchos_VerboseObject.hpp" #include "Teuchos_VerbosityLevel.hpp" #include "Teuchos_XMLParameterListCoreHelpers.hpp" +#include "Teuchos_YamlParameterListCoreHelpers.hpp" #include "Teuchos_config.h" #include "Thyra_EpetraLinearOp.hpp" #include "Thyra_EpetraThyraWrappers.hpp" @@ -27,7 +30,20 @@ #include "Thyra_OperatorVectorTypes.hpp" #include "Thyra_SolveSupportTypes.hpp" #include "Thyra_VectorBase.hpp" -#include "sl_epetra_interface.h" + +#include "Thyra_LinearOpTester.hpp" +#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp" +#include "Thyra_LinearOpWithSolveTester.hpp" + +#ifdef GOMA_ENABLE_TPETRA +#include "Thyra_TpetraLinearOp.hpp" +#include "Thyra_TpetraThyraWrappers.hpp" +#include "Thyra_TpetraVector.hpp" +#include "linalg/sparse_matrix_tpetra.h" +#endif + +#include "EpetraExt_RowMatrixOut.h" +#include "EpetraExt_VectorOut.h" #ifdef GOMA_ENABLE_TEKO // Teko-Package includes @@ -42,45 +58,42 @@ #include "Epetra_Map.h" #include "Epetra_RowMatrix.h" #include "Epetra_Vector.h" +#include "linalg/sparse_matrix.h" +#include "linalg/sparse_matrix_epetra.h" #include "sl_stratimikos_interface.h" #include "sl_util_structs.h" -extern "C" { - -int stratimikos_solve(struct GomaLinearSolverData *ams, - double *x_, - double *b_, - int *iterations, - char stratimikos_file[MAX_NUM_MATRICES][MAX_CHAR_IN_INPUT], - int imtrx) { - using Teuchos::RCP; - bool success = true; - bool verbose = true; - static RCP solverParams_static[MAX_NUM_MATRICES]; - static bool param_set[MAX_NUM_MATRICES] = {false}; - static bool param_echo[MAX_NUM_MATRICES] = {false}; +struct Stratimikos_Solver_Data { + Teuchos::RCP> solver; + Teuchos::RCP solverParams; + Teuchos::RCP> solverFactory; + Teuchos::RCP> A; - try { - Epetra_Map map = ams->RowMatrix->RowMatrixRowMap(); + Stratimikos_Solver_Data() { + solver = Teuchos::null; + solverParams = Teuchos::null; + solverFactory = Teuchos::null; + } +}; - // Assign A with false so it doesn't get garbage collected. - RCP epetra_A = Teuchos::rcp(ams->RowMatrix, false); - RCP epetra_x = Teuchos::rcp(new Epetra_Vector(Copy, map, x_)); - RCP epetra_b = Teuchos::rcp(new Epetra_Vector(Copy, map, b_)); +extern "C" void stratimikos_solver_destroy(struct GomaLinearSolverData *ams) { + auto solver_data = static_cast(ams->SolverData); + solver_data->solver = Teuchos::null; + solver_data->solverParams = Teuchos::null; + delete solver_data; +} - RCP> A = Thyra::epetraLinearOp(epetra_A); - RCP> x = Thyra::create_Vector(epetra_x, A->domain()); - RCP> b = Thyra::create_Vector(epetra_b, A->range()); +using Teuchos::RCP; - Teuchos::RCP outstream = Teuchos::VerboseObjectBase::getDefaultOStream(); +static void stratimikos_solve_setup(RCP> A, + Stratimikos_Solver_Data *solver_data, + std::string stratimikos_file, + bool echo_params) { - // Get parameters from file - if (!param_set[imtrx]) { - param_set[imtrx] = true; - solverParams_static[imtrx] = Teuchos::getParametersFromXmlFile(stratimikos_file[imtrx]); - } + if (solver_data->solver.is_null()) { + // Set up solver (only once per matrix) - RCP solverParams = solverParams_static[imtrx]; + RCP solverParams = solver_data->solverParams; // Set up base builder Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder; @@ -94,23 +107,193 @@ int stratimikos_solve(struct GomaLinearSolverData *ams, // set up solver factory using base/params RCP> solverFactory = linearSolverBuilder.createLinearSolveStrategy(""); + solver_data->solverFactory = solverFactory; - if (!param_echo[imtrx]) { - std::string echo_file(stratimikos_file[imtrx]); + if (echo_params) { + std::string echo_file(stratimikos_file); echo_file = "echo_" + echo_file; linearSolverBuilder.writeParamsFile(*solverFactory, echo_file); - param_echo[imtrx] = true; +#ifdef GOMA_STRATIMIKOS_WRITE_VALID_PARAMS + if (imtrx == 0) { + auto valid_params = linearSolverBuilder.getValidParameters(); + Teuchos::writeParameterListToYamlFile(*valid_params, "stratimikos_valid_params.yaml"); + } +#endif } + Teuchos::RCP outstream = Teuchos::VerboseObjectBase::getDefaultOStream(); // set output stream solverFactory->setOStream(outstream); // set solver verbosity solverFactory->setDefaultVerbLevel(Teuchos::VERB_NONE); - RCP> solver = Thyra::linearOpWithSolve(*solverFactory, A); + solver_data->solver = solverFactory->createOp(); + Thyra::initializeOp(*(solver_data->solverFactory), A, solver_data->solver.ptr()); + } else { + Thyra::initializeAndReuseOp(*(solver_data->solverFactory), A, solver_data->solver.ptr()); + } +} - Thyra::SolveStatus status = Thyra::solve(*solver, Thyra::NOTRANS, *b, x.ptr()); +extern "C" { +#ifdef GOMA_ENABLE_TPETRA +int stratimikos_solve_tpetra(struct GomaLinearSolverData *ams, + double *x_, + double *b_, + int *iterations, + char stratimikos_file[MAX_NUM_MATRICES][MAX_CHAR_IN_INPUT], + int imtrx) { + using Teuchos::RCP; + auto matrix = static_cast(ams->GomaMatrixData); + auto *tpetra_data = static_cast(matrix->data); + bool success = true; + bool verbose = true; + static bool param_echo[MAX_NUM_MATRICES] = {false}; + + if (ams->SolverData == NULL) { + ams->SolverData = new Stratimikos_Solver_Data(); + ams->DestroySolverData = stratimikos_solver_destroy; + } + auto solver_data = static_cast(ams->SolverData); + + try { + RCP> tpetra_A = tpetra_data->matrix; + if (!tpetra_data->matrix->isFillComplete()) { + tpetra_data->matrix->endAssembly(); + } + + RCP> tpetra_x = + rcp(new Tpetra::Vector(tpetra_A->getDomainMap())); + RCP> tpetra_b = + rcp(new Tpetra::Vector(tpetra_A->getRangeMap())); + + for (size_t i = 0; i < tpetra_x->getLocalLength(); i++) { + tpetra_x->replaceGlobalValue(matrix->global_ids[i], x_[i]); + } + for (size_t i = 0; i < tpetra_b->getLocalLength(); i++) { + tpetra_b->replaceGlobalValue(matrix->global_ids[i], b_[i]); + } +#if 0 + Tpetra::MatrixMarket::Writer>::writeSparseFile("A.mm", + tpetra_A); + Tpetra::MatrixMarket::Writer>::writeDenseFile("x.mm", tpetra_x); + Tpetra::MatrixMarket::Writer>::writeDenseFile("b.mm", tpetra_b); +#endif + + solver_data->A = Thyra::createConstLinearOp( + Teuchos::rcp_dynamic_cast>(tpetra_A)); + + RCP> x = Thyra::createVector(tpetra_x); + RCP> b = Thyra::createVector(tpetra_b); + + Teuchos::RCP outstream = Teuchos::VerboseObjectBase::getDefaultOStream(); + + // Get parameters from file + if (solver_data->solverParams.is_null()) { + std::filesystem::path path(stratimikos_file[imtrx]); + if (path.extension() == ".yaml") { + solver_data->solverParams = Teuchos::getParametersFromYamlFile(stratimikos_file[imtrx]); + } else { + solver_data->solverParams = Teuchos::getParametersFromXmlFile(stratimikos_file[imtrx]); + } + } + + stratimikos_solve_setup(solver_data->A, solver_data, stratimikos_file[imtrx], + param_echo[imtrx]); + param_echo[imtrx] = true; + + Thyra::SolveStatus status = + Thyra::solve(*(solver_data->solver), Thyra::NOTRANS, *b, x.ptr()); + + *iterations = 1; + if (!status.extraParameters.is_null()) { + try { + *iterations = status.extraParameters.get()->get("Iteration Count"); + } catch (const Teuchos::Exceptions::InvalidParameter &excpt) { + } + } + + /* Convert solution vector */ + int NumMyRows = matrix->n_rows; + + auto x_data = tpetra_x->getData(); + for (int i = 0; i < NumMyRows; i++) { + x_[i] = x_data[i]; + } + x = Teuchos::null; + Thyra::uninitializeOp(*(solver_data->solverFactory), solver_data->solver.ptr()); + } + TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success) + tpetra_data->matrix->beginAssembly(); + + if (success) { + return 0; + } else { + return -1; + } +} +#else /* GOMA_ENABLE_TPETRA */ +int stratimikos_solve_tpetra(struct GomaLinearSolverData *ams, + double *x_, + double *b_, + int *iterations, + char stratimikos_file[MAX_NUM_MATRICES][MAX_CHAR_IN_INPUT], + int imtrx) { + GOMA_EH(GOMA_ERROR, "Not built with Tpetra Stratimikos support!"); + return -1; +} +#endif /* GOMA_ENABLE_TPETRA */ + +int stratimikos_solve(struct GomaLinearSolverData *ams, + double *x_, + double *b_, + int *iterations, + char stratimikos_file[MAX_NUM_MATRICES][MAX_CHAR_IN_INPUT], + int imtrx) { + using Teuchos::RCP; + bool success = true; + bool verbose = true; + static bool param_echo[MAX_NUM_MATRICES] = {false}; + + GomaSparseMatrix matrix = (GomaSparseMatrix)ams->GomaMatrixData; + EpetraSparseMatrix *epetra_matrix = static_cast(matrix->data); + + if (ams->SolverData == NULL) { + ams->SolverData = new Stratimikos_Solver_Data(); + ams->DestroySolverData = stratimikos_solver_destroy; + } + auto solver_data = static_cast(ams->SolverData); + try { + Epetra_Map map = epetra_matrix->matrix->RowMatrixRowMap(); + + // Assign A with false so it doesn't get garbage collected. + RCP epetra_A = epetra_matrix->matrix; + RCP epetra_x = Teuchos::rcp(new Epetra_Vector(Copy, map, x_)); + RCP epetra_b = Teuchos::rcp(new Epetra_Vector(Copy, map, b_)); + + solver_data->A = Thyra::epetraLinearOp(epetra_A); + RCP> x = Thyra::create_Vector(epetra_x, solver_data->A->domain()); + RCP> b = + Thyra::create_Vector(epetra_b, solver_data->A->range()); + + Teuchos::RCP outstream = Teuchos::VerboseObjectBase::getDefaultOStream(); + + // Get parameters from file + if (solver_data->solverParams.is_null()) { + std::filesystem::path path(stratimikos_file[imtrx]); + if (path.extension() == ".yaml") { + solver_data->solverParams = Teuchos::getParametersFromYamlFile(stratimikos_file[imtrx]); + } else { + solver_data->solverParams = Teuchos::getParametersFromXmlFile(stratimikos_file[imtrx]); + } + } + + stratimikos_solve_setup(solver_data->A, solver_data, stratimikos_file[imtrx], + param_echo[imtrx]); + param_echo[imtrx] = true; + + Thyra::SolveStatus status = + Thyra::solve(*(solver_data->solver), Thyra::NOTRANS, *b, x.ptr()); x = Teuchos::null; @@ -145,15 +328,26 @@ int stratimikos_solve(struct GomaLinearSolverData *ams, #include "mpi.h" +#include "sl_stratimikos_interface.h" extern "C" { #include "mm_eh.h" #include "std.h" +int stratimikos_solve_tpetra(struct GomaLinearSolverData *ams, + double *x_, + double *b_, + int *iterations, + char stratimikos_file[MAX_NUM_MATRICES][MAX_CHAR_IN_INPUT], + int imtrx) { + GOMA_EH(GOMA_ERROR, "Not built with stratimikos support!"); + return -1; +} -int stratimikos_solve(struct Aztec_Linear_Solver_System *ams, +int stratimikos_solve(struct GomaLinearSolverData *ams, double *x_, double *b_, int *iterations, - char *stratimikos_file) { + char stratimikos_file[MAX_NUM_MATRICES][MAX_CHAR_IN_INPUT], + int imtrx) { GOMA_EH(GOMA_ERROR, "Not built with stratimikos support!"); return -1; } diff --git a/src/sl_util.c b/src/sl_util.c index b34080638..b42e1b3f2 100644 --- a/src/sl_util.c +++ b/src/sl_util.c @@ -337,6 +337,11 @@ void sl_free(unsigned int option_mask, struct GomaLinearSolverData *ams[]) { */ void free_ams(struct GomaLinearSolverData *a) { + if (a->DestroySolverData) { + a->DestroySolverData(a); + a->DestroySolverData = NULL; + a->SolverData = NULL; + } safer_free((void **)&(a->data_org)); safer_free((void **)&(a->val)); safer_free((void **)&(a->indx)); @@ -423,6 +428,9 @@ void set_aztec_options_params(int options[], double params[]) { } else if (strcmp(Matrix_Solver, "amesos") == 0) { Linear_Solver = AMESOS; options[AZ_solver] = -1; + } else if (strcmp(Matrix_Solver, "amesos2") == 0) { + Linear_Solver = AMESOS2; + options[AZ_solver] = -1; } else if (strcmp(Matrix_Solver, "stratimikos") == 0) { Linear_Solver = STRATIMIKOS; options[AZ_solver] = -1;