1# This script demonstrates solving a symmetric positive definite linear system using PETSc and HPDDM preconditioner 2# It must be run with exactly 4 MPI processes 3# Example run: 4# mpiexec -n 4 python3 hpddm.py -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_1_eps_threshold_absolute 0.1 -ksp_monitor 5# For more options, see ${PETSC_DIR}/src/ksp/ksp/tutorials/ex76.c 6 7import sys 8import petsc4py 9# Initialize PETSc with command-line arguments 10petsc4py.init(sys.argv) 11from petsc4py import PETSc 12 13# Get the rank of the current process 14rank = PETSc.COMM_WORLD.getRank() 15# Ensure that the script is run with exactly 4 processes 16if PETSc.COMM_WORLD.getSize() != 4: 17 if rank == 0: 18 print("This example requires 4 processes") 19 quit() 20 21# Load directory for input data 22load_dir = PETSc.Options().getString("load_dir", "${DATAFILESPATH}/matrices/hpddm/GENEO") 23 24# Load an index set (IS) from binary file 25sizes = PETSc.IS().load(PETSc.Viewer().createBinary(f"{load_dir}/sizes_{rank}_4.dat", "r", comm = PETSc.COMM_SELF)) 26# Get indices from the loaded IS 27idx = sizes.getIndices() 28 29# Create a PETSc matrix object 30A = PETSc.Mat().create() 31# Set the local and global sizes of the matrix 32A.setSizes([[idx[0], idx[2]], [idx[1], idx[3]]]) 33# Configure matrix using runtime options 34A.setFromOptions() 35# Load matrix A from binary file 36A = A.load(PETSc.Viewer().createBinary(f"{load_dir}/A.dat", "r", comm = PETSc.COMM_WORLD)) 37 38# Load an index set (IS) from binary file 39aux_IS = PETSc.IS().load(PETSc.Viewer().createBinary(f"{load_dir}/is_{rank}_4.dat", "r", comm = PETSc.COMM_SELF)) 40# Set the block size of the index set 41aux_IS.setBlockSize(A.getBlockSize()) 42# Load the Neumann matrix of the current process 43aux_Mat = PETSc.Mat().load(PETSc.Viewer().createBinary(f"{load_dir}/Neumann_{rank}_4.dat", "r", comm = PETSc.COMM_SELF)) 44 45# Create and configure the linear solver (KSP) and preconditioner (PC) 46ksp = PETSc.KSP(PETSc.COMM_WORLD).create() 47pc = ksp.getPC() 48# Use HPDDM as the preconditioner 49pc.setType(PETSc.PC.Type.HPDDM) 50# Set the index set (local-to-global numbering) and auxiliary matrix 51pc.setHPDDMAuxiliaryMat(aux_IS, aux_Mat) 52# Inform HPDDM that the auxiliary matrix is the local Neumann matrix 53pc.setHPDDMHasNeumannMat(True) 54# Apply any command-line options (which may override options from the source file) 55ksp.setFromOptions() 56# Set the system matrix (Amat = Pmat) 57ksp.setOperators(A) 58 59# Create RHS (b) and solution (x) vectors, set random values to b, and solve the system 60b, x = A.createVecs() 61b.setRandom() 62ksp.solve(b, x) 63 64# Output grid and operator complexities on rank 0 65gc, oc = pc.getHPDDMComplexities() 66if rank == 0: 67 print("grid complexity = ", gc, ", operator complexity = ", oc, sep = "") 68