Add COCG method for complex symmetric linear systems #289
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change | ||||
|---|---|---|---|---|---|---|
| @@ -1,6 +1,6 @@ | ||||||
| import LinearAlgebra: ldiv!, \ | ||||||
|
|
||||||
| export Identity, ConjugatedDot, UnconjugatedDot | ||||||
|
|
||||||
| #### Type-handling | ||||||
| """ | ||||||
|
|
@@ -30,3 +30,16 @@ struct Identity end | |||||
| \(::Identity, x) = copy(x) | ||||||
| ldiv!(::Identity, x) = x | ||||||
| ldiv!(y, ::Identity, x) = copyto!(y, x) | ||||||
|
|
||||||
| """ | ||||||
| Conjugated and unconjugated dot products | ||||||
| """ | ||||||
| abstract type AbstractDot end | ||||||
| struct ConjugatedDot <: AbstractDot end | ||||||
| struct UnconjugatedDot <: AbstractDot end | ||||||
|
|
||||||
| _norm(x, ::ConjugatedDot) = norm(x) | ||||||
| _dot(x, y, ::ConjugatedDot) = dot(x, y) | ||||||
|
|
||||||
| _norm(x, ::UnconjugatedDot) = sqrt(sum(xₖ->xₖ^2, x)) | ||||||
| _dot(x, y, ::UnconjugatedDot) = transpose(@view(x[:])) * @view(y[:]) # allocating, but faster than sum(prod, zip(x,y)) | ||||||
|
||||||
| _dot(x, y, ::UnconjugatedDot) = transpose(@view(x[:])) * @view(y[:]) # allocating, but faster than sum(prod, zip(x,y)) | |
| _dot(x, y, ::UnconjugatedDot) = transpose(x) * y |
this shouldn't allocate, at least not on julia 1.5 and above; so let's drop that comment.
There was a problem hiding this comment.
I think it is allocating when x and y are not vectors; see my previous comment #289 (comment). The test in the comment was done with Julia 1.5.4. Could you confirm if it is non-allocating for x and y that are multi-dimensional arrays?
If it is allocating, how about changing the comment in the code from # allocating, but ... to # allocating for non-vectorial x and y, but ...?
There was a problem hiding this comment.
Sorry, but I'm not following why you're using view(x, :) in the first place. x and y are assumed to behave as vectors; we have no block methods.
There was a problem hiding this comment.
and the only thing that might allocate if julia is not smart enough is the array wrapper, but it's O(1) not O(n), so not worth documenting
There was a problem hiding this comment.
Sorry, but I'm not following why you're using
view(x, :)in the first place.xandyare assumed to behave as vectors; we have no block methods.
What if x and y are multi-dimensional arrays? It is quite common for iterative methods to handle such x and y with matrix-free operators A acting on them.
It is true that multi-dimensional arrays x and y are vector-like in a sense that vector operators such as dot(x, y) are defined on them. Other methods in IterativeMethods must have worked fine so far because they relied only on such vector operators. Unfortunately the unconjugated dot product is not one of such vector operators in Julia, so I had to turn those vector-like objects into actual vectors using the wrapper @view.
If the unconjugated dot product had been one of those operators required for making objects behave as vectors, I wouldn't have had to go through these complications...
and the only thing that might allocate if julia is not smart enough is the array wrapper, but it's O(1) not O(n), so not worth documenting
Fair enough. I will remove the comment in the next commit.
There was a problem hiding this comment.
The CG method only solves Ax=b right now where x and b are vectors, not AX=B for multiple right-hand sides. I think we have to reserve B::AbstractMatrix for the use case of block methods if they ever get implemented.
Orthogonal to that is how you internally represent your data. If you want X and B to be a tensor because that maps well to your PDE finite differences scheme, then the user should still wrap it in a vector type
struct MyTensorWithVectorSemantics{T} <: AbstractVector{T}
X::Array{T, 3}
end
and define the interface for it.
Or you don't and you call cg(my_operator, reshape(B, :)) where my_operator restores the shape to B's.
But generally we assume the interface of vectors in iterative methods, and later we can think about supporting block methods to solve AX=B for multiple right-hand sides, where we assume a matrix interface.
There was a problem hiding this comment.
Not sure if we should support multiple right-hand sides with the same cg() interface. The current cg() implementation interprets a non-vector b (including b::AbstractMatrix) as a single right-hand-side column vector, as long as * is defined between the given A and b. I think supporting such usage in the high-level interface in IterativeSolvers.jl was the conclusion people reached after the lengthy discussion in #2 about duck-typing.
That said, my current proposal
_norm(x, ::UnconjugatedDot) = sqrt(sum(xₖ->xₖ^2, x))
_dot(x, y, ::UnconjugatedDot) = transpose(@view(x[:])) * @view(y[:])is not the perfect solution that supports the most general type of x either: _norm assumes iterable x and _dot assumes reshapable x, and the Banach space element mentioned in #289 (comment) satisfies neither condition. The proposal just makes cocg() works for the most popular cases of x::AbstractArray without requiring the users to define _norm(x, ::UnconjugatedDot) and _dot(x, y, ::UnconjugatedDot) separately. For other types of x, the users must define _norm and _dot for the type.
(I think the best way to solve the present issue might be to define a lazy wrapper Conjugate like Transpose or Adjoint in the standard library and rewrite LinearAlgebra:dot to interpret dot(Conjugate(x), y) as the unconjugated dot product.)
Uh oh!
There was an error while loading. Please reload this page.