Skip to content

Add UMFPACK sparse solver as PARDISO alternative for non-Intel platforms#35

Open
mahopan wants to merge 1 commit intonuc-astro:masterfrom
mahopan:feature/umfpack-sparse-solver-minimal
Open

Add UMFPACK sparse solver as PARDISO alternative for non-Intel platforms#35
mahopan wants to merge 1 commit intonuc-astro:masterfrom
mahopan:feature/umfpack-sparse-solver-minimal

Conversation

@mahopan
Copy link
Copy Markdown

@mahopan mahopan commented Feb 21, 2026

Motivation

WinNet requires Intel MKL PARDISO for sparse linear solves in netsolve(), which limits portability to Intel-compiler platforms. This PR adds SuiteSparse/UMFPACK as a drop-in alternative, enabling WinNet to build and run on Apple Silicon Macs (and other non-Intel platforms) using gfortran + Homebrew.

Changes

No existing WinNet source files are modified. Three new files are added:

src/external_tools/umfpack_interface.f90

Fortran 2003 iso_c_binding wrapper for the UMFPACK C API. Wraps umfpack_di_symbolic, umfpack_di_numeric, umfpack_di_solve, umfpack_di_defaults, and the memory-free routines.

src/external_tools/pardiso_shim.f90

Drop-in replacement for Intel MKL's pardiso() subroutine, matching the exact API signature that netsolve() calls. Provides two compile-time backends:

Backend Flag Dependencies
UMFPACK (default) -DUSE_UMFPACK SuiteSparse (-lumfpack -lamd -lsuitesparseconfig)
Dense LAPACK -DUSE_DENSE_SOLVER BLAS/LAPACK only (no external sparse solver)

Key technical detail: CSC/CSR transpose

WinNet's netsolve() constructs the Jacobian in CSC format (ia = column pointers, ja = row indices) but passes it to PARDISO using PARDISO's CSR calling convention (ia = row pointers, ja = column indices). PARDISO therefore interprets the CSC data as CSR and effectively solves A^T · x = b. This is not a bug — the Newton-Raphson iteration converges correctly with A^T — but any replacement solver must reproduce this transposed solve:

  • UMFPACK backend: uses umfpack_di_solve(UMFPACK_At, ...) (transpose solve)
  • Dense LAPACK backend: transposes during matrix assembly (A_dense(j, row))

Makefile.example_mac_arm

Build template for Apple Silicon Macs (aarch64):

  • gfortran via h5fc (Homebrew HDF5)
  • OpenBLAS (replacing Intel MKL BLAS/LAPACK)
  • SuiteSparse/UMFPACK (replacing Intel MKL PARDISO)
  • Auto-detects sparse solver backend; configurable via SPARSE_SOLVER={UMFPACK|PARDISO|DENSE}

Build Instructions (Apple Silicon Mac)

brew install gcc hdf5 openblas suite-sparse
cp Makefile.example_mac_arm Makefile
make

Test Results

Example Status Notes
neutron-decay 30 iterations
NSE (1910 species) 1 iteration
CCSN explosive Si-burning (2919 species) 2874 iterations, 52.6s

The CCSN example reproduces the expected results from Reichert et al. (2023) Fig. 14:

  • Iron-group peak at A ≈ 56: X(⁵⁶Fe) = 0.60, X(⁴He) = 0.24, X(⁵⁸Ni) = 0.083
  • Σ X = 1.000000 (mass conservation)
  • Complete trajectory: T₉ = 6.0 → 5.8 GK (NSE→network) → 0.01 GK (freeze-out)

Scope

This PR is intentionally minimal — it adds a new solver option without touching any existing code. The original Intel/PARDISO build path is completely unaffected.

Motivation:
  WinNet requires Intel MKL PARDISO for its sparse linear solves, which
  limits portability to Intel-compiler platforms. This PR adds
  SuiteSparse/UMFPACK as a drop-in alternative, enabling WinNet to build
  and run on Apple Silicon Macs (and other platforms) using gfortran +
  Homebrew's SuiteSparse.

New files:
  src/external_tools/umfpack_interface.f90
    Fortran 2003 iso_c_binding wrapper for the UMFPACK C API.
    Wraps umfpack_di_symbolic, umfpack_di_numeric, umfpack_di_solve,
    umfpack_di_defaults, and the free routines.

  src/external_tools/pardiso_shim.f90
    Drop-in replacement for Intel MKL's pardiso() subroutine.
    Provides two backends selectable at compile time:
    - UMFPACK (default on non-Intel): links against SuiteSparse
    - Dense LAPACK (fallback): uses dgesv, no external deps beyond BLAS

    Key implementation detail — CSC/CSR transpose handling:
    WinNet's netsolve() (pardiso_class.f90) constructs the Jacobian
    matrix in CSC format (ia = column pointers, ja = row indices) but
    passes it directly to PARDISO, which interprets ia/ja as CSR format
    (ia = row pointers, ja = column indices). The result is that PARDISO
    effectively solves A^T*x = b rather than A*x = b. This is not a bug
    in WinNet — the NR iteration converges correctly with A^T — but any
    replacement solver must reproduce this behavior. The UMFPACK backend
    uses umfpack_di_solve(UMFPACK_At, ...) for the transpose solve, and
    the dense LAPACK backend transposes during matrix assembly.

  Makefile.example_mac_arm
    Build template for Apple Silicon Macs (aarch64) using:
    - gfortran (via h5fc from Homebrew HDF5)
    - OpenBLAS (replacing Intel MKL BLAS/LAPACK)
    - SuiteSparse/UMFPACK (replacing Intel MKL PARDISO)
    Auto-detects sparse solver backend based on compiler and available
    libraries. New SPARSE_SOLVER variable: UMFPACK | PARDISO | DENSE.

Changes to existing files: NONE
  No modifications to any existing WinNet source files or Makefiles.
  The original Intel/PARDISO build path is completely unaffected.

Build instructions (Apple Silicon Mac):
  brew install gcc hdf5 openblas suite-sparse
  cp Makefile.example_mac_arm Makefile
  make

Tested examples:
  - neutron-decay: pass (30 iterations)
  - NSE 1910-species: pass
  - CCSN explosive Si-burning parametrized (2919 species):
    Runs to completion (2874 iter, 52.6s). Final mass fractions
    reproduce Reichert+2023 Fig. 14 iron-group peak at A~56
    with X(56Fe)=0.60, X(4He)=0.24, sum(X)=1.000000.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant