diff --git a/CHANGELOG.md b/CHANGELOG.md
index ccd535926..0997af23b 100644
--- a/CHANGELOG.md
+++ b/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.
diff --git a/docs/src/man/performance.md b/docs/src/man/performance.md
index 159fbb9d5..3947c07e4 100644
--- a/docs/src/man/performance.md
+++ b/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
 
diff --git a/scripts/renames.jl b/scripts/renames.jl
index be91aed76..c32b37971 100644
--- a/scripts/renames.jl
+++ b/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,
diff --git a/src/PETSc.jl b/src/PETSc.jl
index 8844a9d85..7d281d556 100644
--- a/src/PETSc.jl
+++ b/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)
diff --git a/src/dmda.jl b/src/dmda.jl
index a87cdc884..dc65683e4 100644
--- a/src/dmda.jl
+++ b/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)
 
diff --git a/src/dmstag.jl b/src/dmstag.jl
index 6907962da..09fcbada0 100644
--- a/src/dmstag.jl
+++ b/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
diff --git a/src/ksp.jl b/src/ksp.jl
index 51a9f7657..d9168f089 100644
--- a/src/ksp.jl
+++ b/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
diff --git a/src/mat.jl b/src/mat.jl
index 962b5dcbb..6dc4ed3b4 100644
--- a/src/mat.jl
+++ b/src/mat.jl
@@ -51,7 +51,16 @@ MatPtr(::Type{PetscLib}, x...) where {PetscLib <: PetscLibType} =
 
 Base.size(m::AbstractPetscMat{PetscLib}) where {PetscLib} = LibPETSc.MatGetSize(PetscLib,m)
 Base.length(m::AbstractPetscMat{PetscLib}) where {PetscLib} = prod(size(m))
-Base.ndims(m::AbstractPetscMat{PetscLib}) where {PetscLib} = length(LibPETSc.MatGetSize(PetscLib,m))
+Base.ndims(::Type{<:AbstractPetscMat}) = 2
+Base.ndims(m::AbstractPetscMat) = ndims(typeof(m))
+
+# As for any matrix: the global size along dimensions 1 and 2, and 1 beyond them
+function Base.size(m::AbstractPetscMat, d::Integer)
+    d >= 1 || throw(ArgumentError("dimension must be ≥ 1, got $d"))
+    return d <= 2 ? Int(size(m)[d]) : 1
+end
+Base.axes(m::AbstractPetscMat) = map(Base.OneTo, size(m))
+Base.axes(m::AbstractPetscMat, d::Integer) = Base.OneTo(size(m, d))
 """
     type_name(A::AbstractPetscMat)
 
@@ -78,7 +87,6 @@ function set_type!(m::AbstractPetscMat{PetscLib}, type::Symbol) where {PetscLib}
     LibPETSc.MatSetType(getlib(PetscLib), m, String(type))
     return m
 end
-Base.axes(m::PetscMat{PetscLib}, i::Integer) where {PetscLib} = Base.OneTo(Base.size(m)[i])
 
 """
     M::PetscMat = PetscMat(petsclib, S::SparseMatrixCSC; with_arrays = false)
@@ -201,6 +209,10 @@ end
 
 PETSc dense array. This wraps a Julia `Matrix{PetscScalar}` object.
 
+PETSc uses the memory of `A` as its storage, so `A` must be a `Matrix` of exactly the
+library's scalar type; anything else throws an `ArgumentError`. To start from another
+array, convert it first: `PetscMat(petsclib, Matrix{PetscScalar}(B))`.
+
 Replaces v0.4's `MatSeqDense`: construction goes through the type
 (docs/src/man/naming.md §5.1).
 
@@ -597,8 +609,6 @@ function diagonal!(
 end
 
 
-LinearAlgebra.norm(M::PetscMat{PetscLib}, normtype::NormType = NORM_FROBENIUS) where {PetscLib} = LibPETSc.MatNorm(PetscLib, M, normtype)
-
 """
     mul!(y::PetscVec{PetscLib}, M::AbstractPetscMat{PetscLib}, x::PetscVec{PetscLib})
 
@@ -620,30 +630,52 @@ function Base.:*(
     return y
 end
 
+# `A'` and `transpose(A)` are lazy wrappers, so `A' * x` and `mul!(y, A', x)` reach
+# MatMultHermitianTranspose and MatMultTranspose without forming the transpose.
+LinearAlgebra.adjoint(A::AbstractPetscMat) = Adjoint(A)
+LinearAlgebra.transpose(A::AbstractPetscMat) = Transpose(A)
+
+const TransposedPetscMat{PetscLib} = Union{
+    Adjoint{<:Any, <:AbstractPetscMat{PetscLib}},
+    Transpose{<:Any, <:AbstractPetscMat{PetscLib}},
+}
+
+"""
+    mul!(y::PetscVec, A', x::PetscVec)
+    mul!(y::PetscVec, transpose(A), x::PetscVec)
+
+Compute `y = Aᴴx` or `y = Aᵀx` without forming the transpose, and return `y`.
+
+# External Links
+$(doc_external("Mat/MatMultHermitianTranspose"))
+$(doc_external("Mat/MatMultTranspose"))
+"""
 function LinearAlgebra.mul!(
     y::PetscVec{PetscLib},
-    M::Adjoint{AM},
+    M::TransposedPetscMat{PetscLib},
     x::PetscVec{PetscLib},
-) where {PetscLib, AM <: PetscMat{PetscLib}}
-    LibPETSc.MatMultHermitianTranspose(PetscLib, parent(M), x, y)
+) where {PetscLib}
+    multiply = M isa Adjoint ? LibPETSc.MatMultHermitianTranspose : LibPETSc.MatMultTranspose
+    capture_callback_errors(() -> multiply(PetscLib, parent(M), x, y))
     return y
 end
 
-function LinearAlgebra.mul!(
-    y::PetscVec{PetscLib},
-    M::Transpose{AM},
-    x::PetscVec{PetscLib},
-) where {PetscLib, AM <: PetscMat{PetscLib}}
-    LibPETSc.MatMultTranspose(PetscLib, parent(M), x, y)
+# The result has the layout of the columns of `A`, which MatCreateVecs returns first
+function Base.:*(
+    M::TransposedPetscMat{PetscLib},
+    x::AbstractPetscVec{PetscLib},
+) where {PetscLib}
+    y, _ = LibPETSc.MatCreateVecs(getlib(PetscLib), parent(M))
+    mul!(y, M, x)
     return y
 end
 
-function LinearAlgebra.issymmetric(A::PetscMat{PetscLib}; tol = 0.0) where {PetscLib} 
-    PetscReal = real(scalartype(PetscLib))        
+function LinearAlgebra.issymmetric(A::PetscMat{PetscLib}; tol = 0.0) where {PetscLib}
+    PetscReal = real(scalartype(PetscLib))
     return LibPETSc.MatIsSymmetric(PetscLib, A, PetscReal(tol))
 end
-function LinearAlgebra.ishermitian(A::PetscMat{PetscLib}; tol = 0.0) where {PetscLib} 
-    PetscReal = scalartype(PetscLib)
+function LinearAlgebra.ishermitian(A::PetscMat{PetscLib}; tol = 0.0) where {PetscLib}
+    PetscReal = real(scalartype(PetscLib))
     return LibPETSc.MatIsHermitian(PetscLib, A, PetscReal(tol))
 end
 
@@ -714,8 +746,10 @@ end
 """
     B = PetscMat(petsclib, rowptr, colval, nzval; comm = MPI.COMM_SELF, ncols = …)
 
-Create a PETSc SeqAIJ matrix directly on the CSR arrays `rowptr`, `colval` and
-`nzval`, which PETSc borrows rather than copies.
+Create a PETSc SeqAIJ matrix on the CSR arrays `rowptr`, `colval` and `nzval`.
+Arrays that are already `Vector`s of the library's integer and scalar types are
+borrowed rather than copied; 
+any other vector (another element type, a view, an offset array) is copied once.
 
 `rowptr` and `colval` are 0-based, PETSc's own base for bulk index arrays
 (docs/src/man/naming.md §12.1). The number of rows is `length(rowptr) - 1`;
@@ -732,9 +766,9 @@ $(doc_external("Mat/MatCreateSeqAIJWithArrays"))
 """
 function LibPETSc.PetscMat(
     petsclib::PetscLibType,
-    rowptr::Vector{<:Integer},
-    colval::Vector{<:Integer},
-    nzval::Vector;
+    rowptr::AbstractVector{<:Integer},
+    colval::AbstractVector{<:Integer},
+    nzval::AbstractVector;
     comm = MPI.COMM_SELF,
     ncols::Integer = isempty(colval) ? 0 : (maximum(colval) + 1),
 )
@@ -742,9 +776,9 @@ function LibPETSc.PetscMat(
     PetscInt = inttype(petsclib)
     PetscScalar = scalartype(petsclib)
 
-    row_ptr = convert(Vector{PetscInt}, rowptr)
-    col_idx = convert(Vector{PetscInt}, colval)
-    values = convert(Vector{PetscScalar}, nzval)
+    row_ptr = c_vector(PetscInt, rowptr)
+    col_idx = c_vector(PetscInt, colval)
+    values = c_vector(PetscScalar, nzval)
 
     mat = LibPETSc.MatCreateSeqAIJWithArrays(
         petsclib,
@@ -803,15 +837,17 @@ function destroy!(m::AbstractPetscMat{PetscLib}) where {PetscLib}
     return nothing
 end
 
-function LinearAlgebra.mul!(
-    y::PetscVec{PetscLib},
-    M::Transpose{PetscScalar, AM},
-    x::PetscVec{PetscLib},
-) where {PetscLib, PetscScalar, AM <: AbstractPetscMat{PetscLib}}
-    LibPETSc.MatMultTranspose(PetscLib, parent(M), x, y)
-    return y
-end
+"""
+    copyto!(M::AbstractPetscMat, S::SparseMatrixCSC)
 
+Write the nonzeros of `S` into `M` and return `M`. Call [`assemble!`](@ref) afterwards.
+
+`S` is indexed from the first row this rank owns: entry `S[i, j]` goes to global row
+`r + i` and global column `r + j`, where `r` is the number of rows before this rank's.
+In serial `r = 0`, so `S` is copied as it stands. In parallel `S` is the rank's diagonal
+block, and couplings to other ranks' unknowns are not written. Throws a
+`DimensionMismatch` when `S` reaches past this rank's rows or past the last column.
+"""
 function Base.copyto!(
     M::AbstractPetscMat{PetscLib},
     S::SparseMatrixCSC,
@@ -820,7 +856,13 @@ function Base.copyto!(
     PetscInt = PetscLib.PetscInt
     PetscScalar = PetscLib.PetscScalar
     row_start = row_rng[1]
-    _, n = size(S)
+    m, n = size(S)
+    m <= row_rng[2] - row_start && row_start + n <= size(M, 2) || throw(
+        DimensionMismatch(
+            "a $(m)x$(n) block starting at row $(row_start + 1) does not fit the " *
+            "$(row_rng[2] - row_start) rows this rank owns of a matrix with $(size(M, 2)) columns",
+        ),
+    )
     for j in 1:n
         for ii in S.colptr[j]:(S.colptr[j + 1] - 1)
             i = S.rowval[ii]
@@ -832,17 +874,19 @@ end
 
 """
     set_values!(
-        M::AbstractPetscMat{PetscLib},
-        rows_0b::Vector{MatStencil},
-        cols_0b::Vector{MatStencil},
-        rowvals::Array{PetscScalar},
+        M::AbstractPetscMat,
+        rows_0b::AbstractVector{MatStencil},
+        cols_0b::AbstractVector{MatStencil},
+        rowvals::AbstractArray,
         insertmode::InsertMode = INSERT_VALUES;
         num_rows = length(rows_0b),
         num_cols = length(cols_0b)
     )
 
 Set values of the matrix `M` with base-0 row and column indices `rows_0b` and
-`cols_0b`, inserting the values `rowvals`.
+`cols_0b`, inserting the values `rowvals`, read row by row in linear order.
+A `Vector{MatStencil}` and a `Vector` (or `Array`) of the library's scalar type 
+go to PETSc as they are; anything else is copied first.
 
 The `_0b` suffix says what the base is: a bulk index array handed to C keeps
 PETSc's base rather than being rebuilt on a hot path (docs/src/man/naming.md
@@ -856,38 +900,52 @@ $(doc_external("Mat/MatSetValuesStencil"))
 """
 function set_values!(
     M::AbstractPetscMat{PetscLib},
-    rows_0b::Vector{MatStencil},
-    cols_0b::Vector{MatStencil},
-    rowvals::Array{PetscScalar},
+    rows_0b::AbstractVector{MatStencil},
+    cols_0b::AbstractVector{MatStencil},
+    rowvals::AbstractArray,
     insertmode::InsertMode = INSERT_VALUES;
     num_rows = length(rows_0b),
     num_cols = length(cols_0b),
-) where {PetscLib, PetscScalar}
-    PetscScalar === PetscLib.PetscScalar || throw(
-        ArgumentError(
-            "values have element type $PetscScalar, but the library uses " *
-            "$(PetscLib.PetscScalar)",
-        ),
-    )
-    num_rows * num_cols <= length(rowvals) || throw(
-        DimensionMismatch(
-            "a $(num_rows)x$(num_cols) block needs $(num_rows * num_cols) values, " *
-            "got $(length(rowvals))",
-        ),
-    )
+) where {PetscLib}
+    PetscInt = PetscLib.PetscInt
+    check_block(num_rows, num_cols, length(rows_0b), length(cols_0b), length(rowvals))
     LibPETSc.MatSetValuesStencil(
         PetscLib,
         M,
-        num_rows,
-        rows_0b,
-        num_cols,
-        cols_0b,
-        rowvals,
+        PetscInt(num_rows),
+        c_vector(MatStencil, rows_0b),
+        PetscInt(num_cols),
+        c_vector(MatStencil, cols_0b),
+        c_vector(PetscLib.PetscScalar, rowvals),
         insertmode,
     )
     return M
 end
 
+# The array PETSc reads, of element type `T`: a `Vector{T}` as it is, any other
+# `Array{T}` reshaped without a copy, anything else copied in linear order, so
+# views, ranges and offset arrays all work.
+c_vector(::Type{T}, v::Vector{T}) where {T} = v
+c_vector(::Type{T}, v::Array{T}) where {T} = vec(v)
+c_vector(::Type{T}, v::AbstractArray) where {T} = copyto!(Vector{T}(undef, length(v)), v)
+
+# A num_rows x num_cols block needs that many indices and values
+function check_block(num_rows, num_cols, nrows, ncols, nvals)
+    num_rows <= nrows && num_cols <= ncols || throw(
+        DimensionMismatch(
+            "a $(num_rows)x$(num_cols) block needs $num_rows row and $num_cols column " *
+            "indices, got $nrows and $ncols",
+        ),
+    )
+    num_rows * num_cols <= nvals || throw(
+        DimensionMismatch(
+            "a $(num_rows)x$(num_cols) block needs $(num_rows * num_cols) values, " *
+            "got $nvals",
+        ),
+    )
+    return nothing
+end
+
 
 function Base.setindex!(
     M::AbstractPetscMat{PetscLib},
@@ -1110,7 +1168,7 @@ function Base.:*(M::MatShell{PetscLib}, x::AbstractVector) where {PetscLib}
     # PETSc works on this copy of `x` in place, so it must outlive the product
     x_copy = PetscScalar.(x)
     return GC.@preserve x_copy begin
-        petsc_x = LibPETSc.VecCreateSeqWithArray(petsclib, MPI.COMM_SELF, PetscInt(1), PetscInt(n), x_copy)
+        petsc_x = seq_vec_with_array(petsclib, MPI.COMM_SELF, PetscInt(1), PetscInt(n), x_copy)
         petsc_y = LibPETSc.VecCreateSeq(petsclib, MPI.COMM_SELF, PetscInt(m))
         mul!(petsc_y, M, petsc_x)
         result = petsc_y[:]
@@ -1148,17 +1206,20 @@ end
 
 """
     set_values!(
-        M::AbstractMat{PetscLib},
-        rows_0b::Vector{PetscInt},
-        cols_0b::Vector{PetscInt},
-        rowvals::Array{PetscScalar},
+        M::AbstractPetscMat,
+        rows_0b::AbstractVector{<:Integer},
+        cols_0b::AbstractVector{<:Integer},
+        rowvals::AbstractArray,
         insertmode::InsertMode = INSERT_VALUES;
         num_rows = length(rows_0b),
         num_cols = length(cols_0b)
     )
 
 Set values of the matrix `M` with base-0 row and column indices `rows_0b` and
-`cols_0b`, inserting the values `rowvals`.
+`cols_0b`, inserting the values `rowvals`, read row by row in linear order.
+A `Vector` of the library's integer type for the indices, and a `Vector` (or `Array`)
+of its scalar type for the values, go to PETSc as they are; anything else is
+copied first.
 
 The `_0b` suffix says what the base is: a bulk index array handed to C keeps
 PETSc's base rather than being rebuilt on a hot path (docs/src/man/naming.md
@@ -1172,52 +1233,76 @@ $(doc_external("Mat/MatSetValues"))
 """
 function set_values!(
     M::AbstractPetscMat{PetscLib},
-    rows_0b::Vector{PetscInt},
-    cols_0b::Vector{PetscInt},
-    rowvals::Array{PetscScalar},
+    rows_0b::AbstractVector{<:Integer},
+    cols_0b::AbstractVector{<:Integer},
+    rowvals::AbstractArray,
     insertmode::InsertMode = INSERT_VALUES;
     num_rows = length(rows_0b),
     num_cols = length(cols_0b),
-) where {PetscLib, PetscScalar, PetscInt}
-    PetscScalar === PetscLib.PetscScalar || throw(
-        ArgumentError(
-            "values have element type $PetscScalar, " *
-            "but the library uses $(PetscLib.PetscScalar)",
-        ),
-    )
-    PetscInt === PetscLib.PetscInt || throw(
-        ArgumentError(
-            "indices have element type $PetscInt, " *
-            "but the library uses $(PetscLib.PetscInt)",
-        ),
-    )
-    num_rows * num_cols <= length(rowvals) || throw(
-        DimensionMismatch(
-            "a $(num_rows)x$(num_cols) block needs $(num_rows * num_cols) values, " *
-            "got $(length(rowvals))",
-        ),
-    )
+) where {PetscLib}
+    PetscInt = PetscLib.PetscInt
+    check_block(num_rows, num_cols, length(rows_0b), length(cols_0b), length(rowvals))
     LibPETSc.MatSetValues(
         PetscLib,
         M,
         PetscInt(num_rows),
-        PetscInt.(rows_0b),
+        c_vector(PetscInt, rows_0b),
         PetscInt(num_cols),
-        PetscInt.(cols_0b),
-        rowvals,
+        c_vector(PetscInt, cols_0b),
+        c_vector(PetscLib.PetscScalar, rowvals),
         insertmode,
     )
     return M
 end
 
+"""
+    norm(A::AbstractPetscMat, p::Real = 2)
+    norm(A::AbstractPetscMat, normtype::NormType)
+
+The norm of the entries of `A` taken as one vector, as `norm` means for a Julia matrix.
+Only `p = 2`, the Frobenius norm, is available: PETSc computes no other entrywise norm.
+For the norms induced by a vector norm use [`opnorm`](@ref). The second form passes a
+PETSc `NormType` straight to `MatNorm`, where `NORM_1` and `NORM_INFINITY` are the
+induced norms.
+
+# External Links
+$(doc_external("Mat/MatNorm"))
+"""
 function LinearAlgebra.norm(
     M::AbstractPetscMat{PetscLib},
     normtype::NormType = NORM_FROBENIUS,
 ) where {PetscLib}
-    PetscReal = PetscLib.PetscReal
-    #r_val = Ref{PetscReal}()
-    r_val = LibPETSc.MatNorm(PetscLib, M, normtype)
-    return r_val
+    return LibPETSc.MatNorm(PetscLib, M, normtype)
+end
+
+function LinearAlgebra.norm(M::AbstractPetscMat, p::Real)
+    p == 2 || throw(
+        ArgumentError(
+            "PETSc computes only the p = 2 (Frobenius) entrywise norm of a matrix, " *
+            "got p = $p; use opnorm(A, 1) or opnorm(A, Inf) for the induced norms",
+        ),
+    )
+    return norm(M, NORM_FROBENIUS)
+end
+
+"""
+    opnorm(A::AbstractPetscMat, p::Real = 2)
+
+The operator norm of `A` induced by the vector `p`-norm: the largest column sum of
+`|A|` for `p = 1`, the largest row sum for `p = Inf`. PETSc does not compute the
+induced 2-norm, so `p = 2` throws an `ArgumentError`, as does any other `p`.
+
+# External Links
+$(doc_external("Mat/MatNorm"))
+"""
+function LinearAlgebra.opnorm(M::AbstractPetscMat, p::Real = 2)
+    p == 1 && return norm(M, NORM_1)
+    p == Inf && return norm(M, NORM_INFINITY)
+    throw(
+        ArgumentError(
+            "PETSc computes the induced matrix norm only for p = 1 and p = Inf, got p = $p",
+        ),
+    )
 end
 
 # ====
\ No newline at end of file
diff --git a/src/ts.jl b/src/ts.jl
index 61ee9bdca..9179224e9 100644
--- a/src/ts.jl
+++ b/src/ts.jl
@@ -374,9 +374,9 @@ for the argument of the matching setter and would be shadowed by it, and because
 comes back is the name of a PETSc implementation rather than a Julia type
 (docs/src/man/naming.md §3).
 
-Declared here, with methods on `TS`, `KSP`, `PetscVec`, `PetscMat` and `AbstractPetscDM`.
-On a `TS` the answer is a `Symbol` (or `nothing` before a type is set); the other
-methods still answer with a `String`.
+Declared here, with methods on `TS`, `KSP`, `SNES`, `PC`, `PetscVec`, `PetscMat` and
+`AbstractPetscDM`.
+Every method answers with a `Symbol`, or `nothing` while no type has been set.
 
 # External Links
 $(doc_external("TS/TSGetType"))
diff --git a/src/vec.jl b/src/vec.jl
index e3c058474..4d9667f77 100644
--- a/src/vec.jl
+++ b/src/vec.jl
@@ -111,7 +111,9 @@ A standard, sequentially-stored serial PETSc vector, wrapping the Julia vector
 `v`.
 
 This reuses the array `v` as storage, and so `v` should not be `resize!`-ed or
-otherwise have its length modified while the PETSc object exists.
+otherwise have its length modified while the PETSc object exists. For the same
+reason `v` must be a `Vector` of exactly the library's scalar type; anything else
+throws an `ArgumentError`. Convert other arrays first: `PetscVec(petsclib, Vector{PetscScalar}(w))`.
 The vector keeps `v` alive, so `PetscVec(petsclib, [1.0, 2.0])` is safe.
 
 This should only be need to be called for more advanced uses, for most simple
@@ -135,25 +137,31 @@ function LibPETSc.PetscVec(
         ),
     )
     PetscInt = petsclib.PetscInt
-    v = LibPETSc.VecCreateSeqWithArray(
-        petsclib,
-        comm,
-        PetscInt(blocksize),
-        PetscInt(length(array)),
-        array,
-    )
+    v = seq_vec_with_array(petsclib, comm, PetscInt(blocksize), PetscInt(length(array)), array)
     keep_alive!(v, array)
     finalizer(destroy!, v)
     return v
 end
 
 
+# A sequential Vec on the Julia array `array`, without a copy. PETSc starts it at
+# object state 0, which is also the state the norm cache slots get when PETSc first
+# allocates them, so after one norm every other norm type would read as cached and
+# zero. Assembling once moves the Vec past that state; on one process it costs nothing.
+function seq_vec_with_array(petsclib, comm, bs, n, array)
+    v = LibPETSc.VecCreateSeqWithArray(petsclib, comm, bs, n, array)
+    LibPETSc.VecAssemblyBegin(petsclib, v)
+    LibPETSc.VecAssemblyEnd(petsclib, v)
+    return v
+end
+
 # =============================================================================
 # Multiple dispatch to make PetscVec behave like Julia Vector
 # =============================================================================
 
 # Treat PETSc vectors as 1-D array-like objects for broadcasting/indexing
 Base.ndims(::Type{<:AbstractPetscVec}) = 1
+Base.ndims(v::AbstractPetscVec) = ndims(typeof(v))
 Base.IndexStyle(::Type{<:AbstractPetscVec}) = IndexLinear()
 Base.axes(v::AbstractPetscVec) = (Base.OneTo(length(v)),)
 # A Vec takes part in a broadcast as itself; the default would `collect` it
@@ -181,6 +189,23 @@ petsclib_of(::AbstractPetscVec{PetscLib}) where {PetscLib} = getlib(PetscLib)
 Base.size(v::AbstractPetscVec{PetscLib}) where {PetscLib} = LibPETSc.VecGetSize(PetscLib,v)
 Base.length(v::AbstractPetscVec{PetscLib}) where {PetscLib} = prod(size(v))
 Base.lastindex(v::AbstractPetscVec{PetscLib}) where {PetscLib} = length(v)
+
+# As for any vector: the global length along dimension 1, and 1 beyond it
+function Base.size(v::AbstractPetscVec, d::Integer)
+    d >= 1 || throw(ArgumentError("dimension must be ≥ 1, got $d"))
+    return d == 1 ? Int(length(v)) : 1
+end
+
+"""
+    similar(v::AbstractPetscVec)
+
+A new vector with the same layout as `v` (size, parallel distribution, type), its
+entries not set. There is no `similar(v, T)` or `similar(v, dims)`: every vector of a
+library holds that library's scalar type, and a different size has no layout to copy.
+
+# External Links
+$(doc_external("Vec/VecDuplicate"))
+"""
 Base.similar(v::AbstractPetscVec{PetscLib}) where {PetscLib} =  LibPETSc.VecDuplicate(getlib(PetscLib), v)
 """
     type_name(v::AbstractPetscVec)
@@ -227,7 +252,15 @@ function Base.getindex(v::AbstractPetscVec{PetscLib}, ::Colon) where {PetscLib}
     return getindex(v, 1:n)
 end
 
-Base.isapprox(v::AbstractPetscVec{PetscLib}, w::AbstractPetscVec{PetscLib}; kwargs...) where {PetscLib} = all(v[:] .≈ w[:]; kwargs...)
+"""
+    isapprox(v::AbstractPetscVec, w::AbstractPetscVec; atol, rtol, nans)
+
+Whether every entry of `v` is approximately equal to the matching entry of `w`, with
+the keywords of `isapprox` for numbers applied entry by entry. Unlike `isapprox` on two
+Julia arrays, this does not compare the norm of the difference.
+"""
+Base.isapprox(v::AbstractPetscVec{PetscLib}, w::AbstractPetscVec{PetscLib}; kwargs...) where {PetscLib} =
+    all(isapprox.(v[:], w[:]; kwargs...))
 
 function Base.setindex!(v::AbstractPetscVec{PetscLib}, val, i::Integer) where {PetscLib}
     PetscInt = inttype(PetscLib)
@@ -762,20 +795,25 @@ function ghost_update!(
 end
 
 """
-    v = PetscVec(petsclib, comm, array)
+    v = PetscVec(petsclib, comm, array::Vector)
 
-Creates a sequential PETSc vector of length `n` given a julia array `array`,
-on the communicator `comm`.
-The vector uses `array` as its storage and keeps it alive.
+Creates a sequential PETSc vector on the communicator `comm` that uses `array` as its
+storage and keeps it alive. The element type of `array` must be the library's scalar
+type, which throws an `ArgumentError` otherwise: PETSc reads the memory as it is.
 
 # External Links
 $(doc_external("Vec/VecCreateSeqWithArray"))
 """
 function LibPETSc.PetscVec(petsclib::PetscLib, comm, x::Vector) where {PetscLib <: PetscLibType}
     check_initialized(petsclib)
+    eltype(x) === petsclib.PetscScalar || throw(
+        ArgumentError(
+            "array has element type $(eltype(x)), but the library uses $(petsclib.PetscScalar)",
+        ),
+    )
     PetscInt = petsclib.PetscInt
 
-    v = LibPETSc.VecCreateSeqWithArray(petsclib, comm, PetscInt(1), PetscInt(length(x)), x)
+    v = seq_vec_with_array(petsclib, comm, PetscInt(1), PetscInt(length(x)), x)
     keep_alive!(v, x)
     finalizer(destroy!, v)
 
@@ -811,14 +849,28 @@ function ownership_range(vec::AbstractPetscVec{PetscLib}) where {PetscLib}
     return (r_lo + PetscInt(1)):r_hi
 end
 
-# Overload norm function
+"""
+    norm(v::AbstractPetscVec, p::Real = 2)
+    norm(v::AbstractPetscVec, normtype::NormType)
+
+The `p`-norm of `v`, for `p` = 1, 2 or `Inf`; other values throw an `ArgumentError`.
+The second form passes a PETSc `NormType` straight to `VecNorm`.
+
+# External Links
+$(doc_external("Vec/VecNorm"))
+"""
 function LinearAlgebra.norm(
     v::AbstractPetscVec{PetscLib},
     normtype::LibPETSc.NormType = LibPETSc.NORM_2,
 ) where {PetscLib}
-    PetscReal = PetscLib.PetscReal
-    r_val = LibPETSc.VecNorm(PetscLib, v, normtype)
-    return r_val
+    return LibPETSc.VecNorm(PetscLib, v, normtype)
+end
+
+function LinearAlgebra.norm(v::AbstractPetscVec, p::Real)
+    p == 1 && return norm(v, LibPETSc.NORM_1)
+    p == 2 && return norm(v, LibPETSc.NORM_2)
+    p == Inf && return norm(v, LibPETSc.NORM_INFINITY)
+    throw(ArgumentError("PETSc computes the vector norm only for p = 1, 2 and Inf, got p = $p"))
 end
 
 # ── GPU-aware array access helpers ────────────────────────────────────────────
diff --git a/test/core/errors.jl b/test/core/errors.jl
index baa9ed55d..e438d08e8 100644
--- a/test/core/errors.jl
+++ b/test/core/errors.jl
@@ -66,6 +66,18 @@ MPI.Initialized() || MPI.Init()
             wrong = PetscScalar === Float64 ? Float32 : Float64
             @test_throws ArgumentError PETSc.PetscVec(petsclib, wrong[1, 2, 3])
             @test_throws ArgumentError PETSc.PetscMat(petsclib, wrong[1 2; 3 4])
+            @test_throws "array has element type" PETSc.PetscVec(petsclib, MPI.COMM_SELF, wrong[1, 2, 3])
+            @test_throws ArgumentError PETSc.PetscVec(petsclib, MPI.COMM_SELF, wrong[1, 2, 3])
+        end
+
+        @testset "DMs have 1 to 3 dimensions ($PetscScalar)" begin
+            none = ntuple(_ -> PETSc.DM_BOUNDARY_NONE, 4)
+            @test_throws "a DMDA has 1, 2 or 3 dimensions, got 4" PETSc.DMDA(
+                petsclib, MPI.COMM_SELF, none, (3, 3, 3, 3), 1, 1,
+            )
+            @test_throws "a DMStag has 1, 2 or 3 dimensions, got 4" PETSc.DMStag(
+                petsclib, MPI.COMM_SELF, none, (3, 3, 3, 3), (1, 0, 0, 0, 0), 1,
+            )
         end
 
         @testset "size mismatch ($PetscScalar)" begin
diff --git a/test/vecmat/mat_vec_methods.jl b/test/vecmat/mat_vec_methods.jl
index f7f2e27b5..22ec0e282 100644
--- a/test/vecmat/mat_vec_methods.jl
+++ b/test/vecmat/mat_vec_methods.jl
@@ -1,11 +1,15 @@
 # test/vecmat/mat_vec_methods.jl
 # Base and Mat methods: copyto! between vectors, fill! on a matrix, row zeroing
-# for Dirichlet conditions, matrix options, the diagonal and isassembled.
+# for Dirichlet conditions, matrix options, the diagonal and isassembled; the
+# LinearAlgebra and array-interface contracts (transposed products, norms,
+# isapprox, size and axes); arguments of any vector type, offset arrays included.
 
 using Test
 using PETSc
 using MPI
+using LinearAlgebra: norm, opnorm, mul!, issymmetric, ishermitian
 using SparseArrays: sparse
+using OffsetArrays: OffsetVector
 
 MPI.Initialized() || MPI.Init()
 
@@ -90,5 +94,154 @@ MPI.Initialized() || MPI.Init()
         PETSc.destroy!(da)
     end
 
+    # a nonsymmetric matrix, so a product with A and one with its transpose differ
+    nonsymmetric() = PETSc.PetscMat(petsclib, sparse(PetscScalar[1 2; 3 4]))
+
+    @testset "transposed products" begin
+        A = nonsymmetric()
+        x = PETSc.PetscVec(petsclib, PetscScalar[1, 1])
+        y = PETSc.PetscVec(petsclib, PetscScalar[0, 0])
+        @test mul!(y, A', x) === y
+        @test y[:] == PetscScalar[4, 6]
+        @test mul!(y, transpose(A), x) === y
+        @test y[:] == PetscScalar[4, 6]
+        @test (A' * x)[:] == PetscScalar[4, 6]
+        @test (transpose(A) * x)[:] == PetscScalar[4, 6]
+        @test (A * x)[:] == PetscScalar[3, 7]
+        foreach(PETSc.destroy!, (y, x, A))
+    end
+
+    @testset "issymmetric, ishermitian" begin
+        L, A = laplacian(), nonsymmetric()
+        @test issymmetric(L) && ishermitian(L)
+        @test !issymmetric(A) && !ishermitian(A)
+        foreach(PETSc.destroy!, (A, L))
+    end
+
+    @testset "norms" begin
+        v = PETSc.PetscVec(petsclib, PetscScalar[3, -4])
+        @test norm(v) == 5
+        @test norm(v, 1) == 7
+        @test norm(v, 2) == 5
+        @test norm(v, Inf) == 4
+        @test norm(v, LibPETSc.NORM_1) == 7
+        @test_throws "only for p = 1, 2 and Inf, got p = 3" norm(v, 3)
+
+        A = nonsymmetric()
+        @test norm(A) ≈ sqrt(30)
+        @test norm(A, 2) ≈ sqrt(30)
+        @test opnorm(A, 1) == 6       # largest column sum
+        @test opnorm(A, Inf) == 7     # largest row sum
+        @test_throws "use opnorm(A, 1) or opnorm(A, Inf)" norm(A, 1)
+        @test_throws "only for p = 1 and p = Inf, got p = 2" opnorm(A)
+        foreach(PETSc.destroy!, (A, v))
+    end
+
+    @testset "isapprox keywords" begin
+        v = PETSc.PetscVec(petsclib, PetscScalar[1, 2])
+        w = PETSc.PetscVec(petsclib, PetscScalar[1, 2 + 1e-6])
+        @test isapprox(v, v)
+        @test !isapprox(v, w)
+        @test isapprox(v, w; rtol = 1e-5)
+        @test isapprox(v, w; atol = 1e-5)
+        foreach(PETSc.destroy!, (w, v))
+    end
+
+    @testset "size, axes, ndims" begin
+        v = PETSc.PetscVec(petsclib, PetscScalar[1, 2, 3])
+        @test size(v, 1) === 3
+        @test size(v, 2) === 1
+        @test_throws ArgumentError size(v, 0)
+        @test ndims(v) == 1
+
+        A = PETSc.PetscMat(petsclib, sparse(PetscScalar[1 2 3; 4 5 6]))
+        @test size(A, 1) === 2
+        @test size(A, 2) === 3
+        @test size(A, 3) === 1
+        @test_throws ArgumentError size(A, 0)
+        @test axes(A) == (Base.OneTo(2), Base.OneTo(3))
+        @test axes(A, 2) == Base.OneTo(3)
+        @test ndims(A) == 2
+        @test ndims(typeof(A)) == 2
+        foreach(PETSc.destroy!, (A, v))
+    end
+
+    @testset "copyto! from a SparseMatrixCSC" begin
+        S = sparse(PetscScalar[1 0 2; 0 3 0])
+        M = PETSc.PetscMat(petsclib, 2, 3, 3)
+        @test copyto!(M, S) === M
+        PETSc.assemble!(M)
+        @test M[:, :] == Matrix(S)
+        @test_throws "does not fit" copyto!(M, sparse(PetscScalar[1 2 3 4; 5 6 7 8]))
+        PETSc.destroy!(M)
+    end
+
+    @testset "set_values! takes any vector" begin
+        PetscInt = petsclib.PetscInt
+        reference = PetscScalar[1 2; 3 4]
+
+        A = PETSc.PetscMat(petsclib, 2, 2, 2)
+        PETSc.set_values!(A, PetscInt[0, 1], PetscInt[0, 1], PetscScalar[1, 2, 3, 4])
+        PETSc.assemble!(A)
+        @test A[:, :] == reference
+
+        # other integer and scalar types, a view and offset vectors give the same matrix
+        B = PETSc.PetscMat(petsclib, 2, 2, 2)
+        rows = OffsetVector(Int32[0, 1], -1:0)
+        cols = view([0, 1, 7], 1:2)
+        vals = OffsetVector(Float32[1, 2, 3, 4], 0:3)
+        @test PETSc.set_values!(B, rows, cols, vals) === B
+        PETSc.assemble!(B)
+        @test B[:, :] == reference
+
+        @test_throws "needs 3 row and 2 column indices" PETSc.set_values!(
+            B, PetscInt[0, 1], PetscInt[0, 1], PetscScalar[1, 2, 3, 4, 5, 6]; num_rows = 3,
+        )
+        @test_throws "needs 4 values" PETSc.set_values!(B, PetscInt[0, 1], PetscInt[0, 1], PetscScalar[1, 2, 3])
+
+        # the MatStencil form, on a 1D DMDA: the x index is the third field
+        da = PETSc.DMDA(petsclib, comm, (PETSc.DM_BOUNDARY_NONE,), (3,), 1, 1)
+        S = PETSc.PetscMat(da)
+        stencils = OffsetVector([LibPETSc.MatStencil(0, 0, i, 0) for i in 0:1], 0:1)
+        PETSc.set_values!(S, stencils, stencils, OffsetVector(PetscScalar[1, 2, 3, 4], 0:3))
+        PETSc.assemble!(S)
+        @test S[1:2, 1:2] == reference
+        foreach(PETSc.destroy!, (S, da, B, A))
+    end
+
+    @testset "PetscMat from CSR arrays of any vector type" begin
+        expected = PetscScalar[5 0; 0 6]
+        A = PETSc.PetscMat(petsclib, [0, 1, 2], [0, 1], PetscScalar[5, 6])
+        B = PETSc.PetscMat(
+            petsclib,
+            OffsetVector(Int32[0, 1, 2], 0:2),
+            view([0, 1], :),
+            Float32[5, 6],
+        )
+        @test A[:, :] == expected
+        @test B[:, :] == expected
+        foreach(PETSc.destroy!, (B, A))
+    end
+
+    @testset "ksp \\ b for any vector" begin
+        L = laplacian()
+        ksp = PETSc.KSP(L; ksp_type = "preonly", pc_type = "lu")
+        x = ksp \ PetscScalar[1, 2, 3]
+        b = OffsetVector(PetscScalar[1, 2, 3], 0:2)
+        y = ksp \ b
+        @test axes(y) == axes(b)
+        @test parent(y) ≈ x
+        @test ksp \ Float32[1, 2, 3] ≈ x
+        PETSc.destroy!(ksp)
+        PETSc.destroy!(L)
+    end
+
+    @testset "dof_slot takes any integer" begin
+        dm = PETSc.DMStag(petsclib, comm, (PETSc.DM_BOUNDARY_NONE,), (3,), (1, 1), 1)
+        loc = PETSc.element_location(dm)
+        @test PETSc.dof_slot(dm, loc, Int32(0)) == PETSc.dof_slot(dm, loc, 0)
+        PETSc.destroy!(dm)
+    end
+
     PETSc.finalize(petsclib)
 end