Vec/Mat: fix transposed products and norms, take any vector

CHANGELOG.md

@@ -1,5 +1,23 @@
 # Changelog
 
+## Unreleased
+
+### Fixed
+
+- An MPI run started against a stale precompile cache (after an update, or a change to a `dev`ed PETSc.jl) could hang for good at 0% CPU. The rank that precompiled PETSc called `MPI.Init` inside the precompile process, which then waited for peers that never came. PETSc no longer initializes MPI while it loads; `initialize` does it, as before.
+- `A' * x`, `transpose(A) * x` and `mul!(y, A', x)` on a PETSc matrix threw a `MethodError`: the transposed products existed but could not be reached. They now run `MatMultHermitianTranspose` and `MatMultTranspose`, without forming the transpose.
+- On a `PetscVec` made on a Julia array (`PetscVec(petsclib, x)`), every norm type after the first came back as 0: `norm(v)` followed by `norm(v, LibPETSc.NORM_1)` gave 5 and then 0. PETSc saw the other norms as already cached, because such a vector starts in the same object state as PETSc's empty norm cache. The vector is now moved past that state when it is created.
+- `isapprox(v, w; rtol = 1e-8)` on two `PetscVec`s threw a `MethodError`; `atol`, `rtol` and `nans` now apply to each entry.
+- `DMDA` and `DMStag` with more than three dimensions failed with an `UndefVarError`; they throw an `ArgumentError` saying a DM has 1 to 3 dimensions.
+- `PetscVec(petsclib, comm, x)` with an `x` of the wrong element type throws an `ArgumentError` naming both types, as `PetscVec(petsclib, x)` does, instead of a `MethodError` from the generated layer.
+- `copyto!(A, S::SparseMatrixCSC)` has a docstring saying how `S` is placed (from the first row the rank owns; in parallel, the rank's diagonal block), and throws a `DimensionMismatch` when `S` does not fit.
+
+### Added
+
+- `set_values!` (both 0-based forms), `ksp \ b`, the CSR `PetscMat(petsclib, rowptr, colval, nzval)` constructor and `dof_slot` accept any vector or integer type: other element types, views, ranges and offset arrays. Arrays already of the library's types still go to PETSc without a copy; anything else is copied once. `ksp \ b` returns `x` with the axes of `b`.
+- `norm(v, p)` for `p` = 1, 2 and `Inf`, and `norm(A, 2)` (Frobenius), beside the `NormType` methods; `opnorm(A, 1)` and `opnorm(A, Inf)` for the induced matrix norms.
+- `size(v, d)`, `size(A, d)`, `axes(A)` and `axes(A, d)` for any PETSc vector or matrix, with `size(_, d) == 1` beyond the last dimension as for a Julia array; `ndims` on a matrix no longer calls into PETSc.
+
 ## v0.5.2
 
 0.5.2 is what a staggered-grid solver needs from PETSc.jl: DMStag gets locations keyed by axis, stencils built from 1-based indices, assembly and halo exchange, and field views that hand a callback concretely typed arrays. TS gets what a long run needs, and Vec and Mat get the Base verbs that were missing. The rest is leaks closed in the generated layer and two performance traps, one of which made a 4-rank DMStag solve 25 times slower. Working code is most likely to notice the three behaviour changes below; the rest is additive.

docs/src/man/performance.md

@@ -61,7 +61,7 @@ julia --project=. -e 'using PETSc'
 mpiexec -n 4 julia --project=. run.jl
 ```
 
-Ranks that start against a stale cache all try to precompile at once and wait on each other's lock. The run then sits at 0% CPU with no error.
+Ranks that start against a stale cache take turns: one precompiles while the others wait on its lock, so every rank pays for the compile before the run starts.
 
 ## Profiling a solve
 

scripts/renames.jl

@@ -462,6 +462,10 @@ const INTERNAL = Set{Symbol}([
     :aliases_dest,
     :check_undistributed,
     :PetscVecStyle,
+    :TransposedPetscMat,
+    :seq_vec_with_array,
+    :c_vector,
+    :check_block,
     :checkout_arrays,
     :checked_array,
     :flags,

src/PETSc.jl

@@ -4,8 +4,6 @@ module PETSc
 
 using MPI, LinearAlgebra, SparseArrays, OffsetArrays, Preferences
 
-MPI.Initialized() || MPI.Init()
-
 function petsc_link(fname)
 """
 [`$fname`](https://petsc.org/release/docs/manualpages/$fname.html)

src/dmda.jl

@@ -75,6 +75,7 @@ function DMDA(
     options...,
 ) where {PetscLib, N}
     check_initialized(getlib(PetscLib))
+    1 <= N <= 3 || throw(ArgumentError("a DMDA has 1, 2 or 3 dimensions, got $N"))
     PetscInt = inttype(PetscLib)
     stencil_type = stencil_type_enum(DMDAStencilType, stencil_type)
 

src/dmstag.jl

@@ -63,6 +63,7 @@ function DMStag(
     options...,
 ) where {PetscLib, N, N1}
     check_initialized(getlib(PetscLib))
+    1 <= N <= 3 || throw(ArgumentError("a DMStag has 1, 2 or 3 dimensions, got $N"))
     N1 == N + 1 || throw(
         DimensionMismatch(
             "dof_per_node has length $N1, " *
@@ -409,13 +410,13 @@ function global_indices(dm::DMStag{PetscLib, N}) where {PetscLib, N}
 end
 
 """
-    slot::Int = dof_slot(dm::DMStag, loc::LibPETSc.DMStagStencilLocation, dof::Int) 
+    slot::Int = dof_slot(dm::DMStag, loc::LibPETSc.DMStagStencilLocation, dof::Integer)
 
 Returns the location `slot` for a degree of freedom `dof` at a given stencil location `loc` in the DMStag `dm`.
 Note that the returned `slot` is 1-based for Julia compatibility.
 `dof` is PETSc's 0-based component number at that location, as in [`stencil`](@ref).
 """
-function dof_slot(dm::DMStag{PetscLib}, loc::LibPETSc.DMStagStencilLocation, dof::Int) where {PetscLib} 
+function dof_slot(dm::DMStag{PetscLib}, loc::LibPETSc.DMStagStencilLocation, dof::Integer) where {PetscLib}
     slot = LibPETSc.DMStagGetLocationSlot(getlib(PetscLib), dm, loc, PetscLib.PetscInt(dof))
     return slot+1
 end

src/ksp.jl

@@ -166,8 +166,10 @@ end
 """
     KSP(petsclib, comm::MPI.Comm, A::SparseMatrixCSC; options...)
 
-Create a [`KSP`](@ref) with the sparse matrix `A` using the `petsclib`. If
-`petsclib` is not given, the default library will be used`.
+Create a [`KSP`](@ref) for the sparse matrix `A`, on `comm`, using the library
+`petsclib`. `A` is copied into a new PETSc matrix that the solver owns, so changing
+`A` afterwards does not change the solver. The keyword arguments are PETSc options,
+as for `KSP(A::AbstractPetscMat; options...)`.
 """
 KSP(petsclib, comm, S::SparseMatrixCSC; kwargs...) 
 
@@ -219,16 +221,9 @@ function Base.:\(ksp::KSP, b::PetscVec{PetscLib}) where {PetscLib}
     return x
 end
 
-function Base.:\(
-    ksp::KSP{PetscLib},
-    b::Vector{PetscScalar},
-) where {PetscLib, PetscScalar}
-    PetscScalar === PetscLib.PetscScalar || throw(
-        ArgumentError(
-            "right-hand side has element type $PetscScalar, " *
-            "but the library uses $(PetscLib.PetscScalar)",
-        ),
-    )
+# Solves into a copy of `b` in the library's scalar type, so any vector works, 
+# and returns `x` with the axes of `b`
+function Base.:\(ksp::KSP{PetscLib}, b::AbstractVector) where {PetscLib}
     c = comm(ksp)
     MPI.Comm_size(c) == 1 || throw(
         ArgumentError(
@@ -237,13 +232,14 @@ function Base.:\(
         ),
     )
     PetscInt = PetscLib.PetscInt
+    PetscScalar = PetscLib.PetscScalar
 
     # PETSc works on this copy of `b` in place, so it must outlive the solve
-    b_copy = PetscScalar.(b)
+    b_copy = copyto!(Vector{PetscScalar}(undef, length(b)), b)
     x = GC.@preserve b_copy begin
-        petsc_b = LibPETSc.VecCreateSeqWithArray(getlib(PetscLib), c, PetscInt(1), PetscInt(length(b)), b_copy)
+        petsc_b = seq_vec_with_array(getlib(PetscLib), c, PetscInt(1), PetscInt(length(b)), b_copy)
         petsc_x = ksp \ petsc_b
-        x = petsc_x[:]
+        x = copyto!(similar(b, PetscScalar), petsc_x[:])
         destroy!(petsc_b)
         destroy!(petsc_x)
         x