Add COCG method for complex symmetric linear systems #289

wsshin · 2021-01-25T05:21:21Z

This PR adds the conjugate orthogonal conjugate gradient (COCG) method for complex symmetric matrices A (i.e., transpose(A) == A rather than adjoint(A) == A). The PR implements cocg and cocg! functions (corresponding to cg and cg! for Hermitian matrices), supporting functions and types, and basic unit tests.

This PR is related with JuliaLang/julia#288.

… symmetric linear systems

wsshin · 2021-01-25T05:23:47Z

I haven't added any documentation pages for COCG, as IterativeSolvers.jl seems to use Documenter.jl, which I haven't used before. Please let me know if there are some basic guidelines as to how to create documentation pages using Documenter.jl.

coveralls · 2021-01-25T05:34:35Z

Pull Request Test Coverage Report for Build 508517998

0 of 63 (100.0%) changed or added relevant lines in 1 file are covered.
No unchanged relevant lines lost coverage.
Overall coverage decreased (-0.3%) to 97.198%

Totals
Change from base Build 444205317:	-0.3%
Covered Lines:	1873
Relevant Lines:	1927

💛 - Coveralls

coveralls · 2021-01-25T05:34:35Z

Pull Request Test Coverage Report for Build 652039845

Warning: This coverage report may be inaccurate.

This pull request's base commit is no longer the HEAD commit of its target branch. This means it includes changes from outside the original pull request, including, potentially, unrelated coverage changes.

For more information on this, see Tracking coverage changes with pull request builds.
To avoid this issue with future PRs, see these Recommended CI Configurations.
For a quick fix, rebase this PR at GitHub. Your next report should be accurate.

Details

17 of 17 (100.0%) changed or added relevant lines in 2 files are covered.
No unchanged relevant lines lost coverage.
Overall coverage increased (+0.3%) to 97.83%

Totals
Change from base Build 444205317:	0.3%
Covered Lines:	1939
Relevant Lines:	1982

💛 - Coveralls

haampie · 2021-01-25T09:05:31Z

Hi @wsshin, thanks for your contribution!

I think another way to do this is to dispatch to a different dot-function inside of CG based on a tag.

For instance:

struct InnerProduct end
struct UnconjugatedProduct end

_norm(::InnerProduct, xs) = norm(xs)
_dot(::InnerProduct, xs, ys) = dot(xs, ys)

_norm(::UnconjugatedProduct, xs) = sqrt(sum(x^2 for x in xs))
_dot(::UnconjugatedProduct, xs, ys) = sum(prod, zip(xs, ys))

and then add this tag InnerProduct / NoInnerProduct to the CG iterable type. That saves some duplication.

Is there a name for this non-conjugate innerproduct-like operation sum(prod, zip(xs, ys))?

haampie

See above, would be nice to reduce duplication a bit.

…product

wsshin · 2021-01-25T21:55:03Z

@haampie, thanks for the suggestion. Please let me know if the new commit fulfills your request.

wsshin

@haampie, I submitted new commits implementing your ideas. Please take a look and let me know what you think.

wsshin · 2021-01-26T14:03:25Z

Not sure why codecov generated a failure. Tried to remove it by adding more unit tests, but it didn't work. Please advise.

wsshin · 2021-01-28T21:41:12Z

@haampie, any thoughts on the new commit in this PR? Hope to finish this PR soon and move on to the implementation of COCR as mentioned in JuliaLang/julia#288.

haampie · 2021-01-28T21:42:32Z

Hi @wsshin , i'll review in the weekend if that's ok with you.

haampie · 2021-01-31T13:16:47Z

src/cg.jl

@@ -47,18 +57,19 @@ function iterate(it::CGIterable, iteration::Int=start(it))
    end

    # u := r + βu (almost an axpy)
-    β = it.residual^2 / it.prev_residual^2
+    ρ = _dot(it.r, it.conjugate_dot)


CG and preconditioned CG are split into separate methods because CG requires one fewer dot product. Can you check if you can drop this dot in your version somehow?

Good point, but I'm afraid we can't. It turns out that the saving of one dot product in the non-preconditioned CG comes from the fact that we calculate ρ = r⋅r and the norm of the residual vector can be calculated by sqrt(ρ) instead of performing the dot product again.

In COCG, we use the unconjugated dot product, so ρ = transpose(r)*r. Therefore, the norm of the residual vector cannot be calculated by sqrt(ρ), but needs to be calculated by sqrt(r⋅r).

In that case you should probably have

β = it.conjugate_dot isa Val{true} ? it.residual^2 / it.prev_residual^2 : _dot(it.r, it.conjugate_dot)

so that it only does the extra dot product in the unconjugated case.

And if you do that please change Val{true} into a named type to improve the readability a bit

haampie · 2021-01-31T13:19:24Z

src/cg.jl

 end

-mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number}
+mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number, boolT <: Union{Val{true},Val{false}}}


Can you turn it into two named structs instead of Val{true/false}, it would be easier to read.

I understand defining InnerProduct and UnconjugatedProduct was your original suggestion, but I thought passing a boolean flag in the top-level user interface like

cg!(x, A, b; conjugate_dot=true)

would be easier for the users than something like

cg!(x, A, b; innerprodtype=UnconjugatedProduct())

Introduction of a new type should be minimized, because it reduces portability. Also, I don't think using a boolean flag reduces readability in this case, because the flag is a keyword argument.

If we agree that passing a boolean flag conjugate_dot to the top interface is a slightly better choice, it would be more natural to do

CGIterable(A, x, r, ..., Val(conjugate_dot))

than

CGIterable(A, x, r, ..., conjugate_dot ? InnerProduct() : UnconjugatedProduct())

This is why I decided to include a boolean type in CGIterable than a completely new type as you suggested originally. Please let me know what you think

I'm thinking the high-level interface could just be cocg! instead of cg! with a flag. Would be nice if the methods would be type stable without relying on constant propagation of conjugate_dot.

Introduction of a new type should be minimized, because it reduces portability.

I'm not sure what you mean by "reducing portability," Wonseok?

A new type has the advantage of allowing more than two possibilities. For example, one could imagine a WeightedDot(B) that corresponds to the dot(x,B,y) dot product — this is useful in cases where A is Hermitian with respect to some nonstandard dot product (the alternative being to perform a change of variables via Cholesky factors of B, which may not always be convenient).

conjugate_dot=true has the disadvantage of potentially requiring a dynamic dispatch at the top level. (Unless the compiler does constant propagation, but I think it only does that if the function is short enough to be inlined, which I don't think cg! is?) That's not such a big deal here, however, because overhead of a single dynamic dispatch should be negligible for large problems where iterative solvers are used.

I'm not sure what you mean by "reducing portability," Wonseok?

Please ignore my comment on portability. I think it was irrelevant to the current situation. (For clarification and my future reference, I meant it is not desirable to define a package-specific collection type that aggregates various fields and to pass such a collection to a function in order to reduce the number of function arguments.)

A new type has the advantage of allowing more than two possibilities. For example, one could imagine a WeightedDot(B) that corresponds to the dot(x,B,y) dot product — this is useful in cases where A is Hermitian with respect to some nonstandard dot product (the alternative being to perform a change of variables via Cholesky factors of B, which may not always be convenient).

I really like this flexibility!

conjugate_dot=true has the disadvantage of potentially requiring a dynamic dispatch at the top level. (Unless the compiler does constant propagation, but I think it only does that if the function is short enough to be inlined, which I don't think cg! is?) That's not such a big deal here, however, because overhead of a single dynamic dispatch should be negligible for large problems where iterative solvers are used.

I didn't know that constant propagation can be a means to avoid dynamic dispatch. Glad to learn a new thing. Thanks for the information! (I find this useful on this topic.)

But before knowing about constant propagation, I thought the the overhead of a dynamic dispatch should be negligible when iterative solvers are useful, as you pointed out.

I have pushed a commit complying with this request.

haampie · 2021-01-31T13:22:45Z

src/cg.jl

@@ -2,20 +2,29 @@ import Base: iterate
 using Printf
 export cg, cg!, CGIterable, PCGIterable, cg_iterator!, CGStateVariables

-mutable struct CGIterable{matT, solT, vecT, numT <: Real}
+# Conjugated dot product
+_dot(x, ::Val{true}) = sum(abs2, x)  # for x::Complex, returns Real


Maybe better to just use _norm / norm; that's what was used before and norm is slightly more stable. Also I'm not sure if sum(abs2, x) works as efficiently on GPUs, norm is definitely optimized.

What do you mean by "norm is more stable"? If you are talking about the possibility of overflow, I think it doesn't matter because the calculated norm is squared again in the algorithm.

Also, for the unconjugated dot product, the norm(x) function would return sqrt(xᵀx), which is not a norm because xᵀx is complex. In COCG the quantity we use is xᵀx, not sqrt(xᵀx), so I don't feel that it is a good idea to store the square-rooted quantity and then to square it again when using it.

I verify that the GPU performance of sum does not degrade.

julia> using LinearAlgebra, CuArrays, BenchmarkTools julia> v = rand(ComplexF64, 10^9); julia> cv = cu(v); julia> @btime norm($v); 1.890 s (0 allocations: 0 bytes) julia> @btime sum(abs2, $v); 1.452 s (0 allocations: 0 bytes) # sum() is significantly faster than norm() on CPU julia> @btime norm($cv); 18.341 ms (1 allocation: 16 bytes) julia> @btime sum(abs2, $cv); 18.341 ms (361 allocations: 13.94 KiB) # sum is on par with norm() on GPU

The main advantage of falling back to norm and dot in the conjugated case is that they generalize better — if the user has defined some specialized Banach-space type (and corresponding self-adjoint linear operators A with overloaded *), they should have overloaded norm and dot, whereas sum(abs2, x) might no longer work.

(Even as simple a generalization as an array of arrays will fail with sum(abs2, x), whereas they work with dot and norm.)

The overhead of the additional sqrt should be negligible.

1.452 s (0 allocations: 0 bytes) # sum() is significantly faster than norm() on CPU

This is because it dispatches to BLAS which does a stable norm computation. But generally BLAS libs are free to choose how to implement norm.

The main advantage of falling back to norm and dot in the conjugated case is that they generalize better — if the user has defined some specialized Banach-space type (and corresponding self-adjoint linear operators A with overloaded *), they should have overloaded norm and dot, whereas sum(abs2, x) might no longer work.

(Even as simple a generalization as an array of arrays will fail with sum(abs2, x), whereas they work with dot and norm.)

This is a very convincing argument. Additionally, considering that a user-defined Banach-space type should overload norm and dot but not the unconjugated dot, it will be nice to be able to define the unconjugated dot using dot like

_dot(x, y, ::UnconjugatedDot) = dot(conj(x), y)

but of course this is inefficient as it conjugates x twice, not to mention extra allocations. I wish dotu was exported, as you suggested in https://github.com/JuliaLang/julia/issues/22227#issuecomment-306224429.

Users can overload _dot(x, y, ::UnconjugatedDot) for the Banach-space type they define, but what would be the best way to minimize such a requirement? I have pushed a commit implementing the unconjugated dot with sum, but if there is a better approach, please let me know.

You could implement the unconjugated dot as transpose(x) * y ... this should have no extra allocations since transpose is a "view" type by default?

Though transpose(x) * y wouldn't work if x and y are matrices (and A is some kind of multlinear operator); maybe it's best to stick with the sum(zip(x,y)) solution you have now.

Can we assume that getindex() is always defined for x and y? Then how about using transpose(@view(x[:])) * @view(y[:])? This uses allocations, but it seems much faster than sum(prod, zip(x,y)):

julia> x = rand(100,100); y = rand(100,100); julia> @btime sum(prod, zip($x,$y)); 12.564 μs (0 allocations: 0 bytes) julia> @btime transpose(@view($x[:])) * @view($y[:]); 1.966 μs (4 allocations: 160 bytes)

It is nearly on a par with dot:

julia> @btime dot($x,$y); 1.839 μs (0 allocations: 0 bytes)

Also, when x and y are AbstractVectors, the proposed method does not use allocations:

julia> x = rand(1000); y = rand(1000); julia> @btime transpose(@view($x[:])) * @view($y[:]); 108.091 ns (0 allocations: 0 bytes)

You could do dot(transpose(x'), x), maybe?

I didn't know that x' was nonallocating! This looked like a great nonallocating implementation, but it turns out that this is slower than the earlier proposed solution using @view:

julia> @btime transpose(@view($x[:])) * @view($y[:]); 267.794 ns (0 allocations: 0 bytes) julia> @btime dot(transpose($x'), $y); 1.547 μs (0 allocations: 0 bytes)

Also, neither x' nor transpose(x) works for x whose dimension is greater than 2.

I just pushed an implementation using @view for now.

wsshin · 2021-02-07T14:12:05Z

@haampie, please let me know if you don't agree with my opinions; I can accommodate yours.

The usage of the unconjugated dot product is not limited to implementing COCG from CG, so it is more appropriate to move its definition outside cg.jl. For example, the unconjugated dot product can be used in implementing COCR from CR.

wsshin · 2021-03-06T22:28:33Z

I think I have addressed all the concerns raised so far. If there is anything else, please let me know.

stevengj · 2021-03-08T20:30:03Z

src/common.jl

+"""
+Conjugated and unconjugated dot products
+"""
+abstract type DotType end


I would call this AbstractDot in keeping with the usual Julia naming conventions.

stevengj · 2021-03-08T20:31:39Z

src/cg.jl

    maxiter::Int
    mv_products::Int
+    dot_type::dotT


Maybe dotproduct — dot_type is a bit inapt since it is not a Type.

stevengj · 2021-03-08T20:36:33Z

src/common.jl

+_dot(x, y, ::ConjugatedDot) = dot(x, y)
+
+_norm(x, ::UnconjugatedDot) = sqrt(sum(xₖ^2 for xₖ in x))
+_dot(x, y, ::UnconjugatedDot) = sum(prod, zip(x, y))


sqrt(sum(xₖ -> xₖ^2, x)) for the norm, so that it can use pairwise summation (iterators are not indexable)?

Wow, this change makes a significant difference in performance:

julia> x = rand(1000); julia> @btime sum(xₖ^2 for xₖ in $x); 1.221 μs (0 allocations: 0 bytes) julia> @btime sum(xₖ->xₖ^2, $x); 85.675 ns (0 allocations: 0 bytes)

But it doesn't seem like indexibility is the only condition for the pairwise summation to kick in. I tried to achieve a similar enhancement by pairwise summation in _dot(x, y, ::UnconjugatedDot) = sum(prod, zip(x, y)) by replacing zip(x,y) with an indexible custom type, but it didn't work as expected. Looked into the code in detail, and it seems that mapreduce_impl will need to be implemented for the custom type replacing zip(x,y)...

stevengj · 2021-03-16T16:17:42Z

test/cg.jl

+        @test norm(A * x2 - b) / norm(b) ≤ reltol
+
+        # The following tests fails CI on Windows and Ubuntu due to a
+        # `SingularException(4)`


What's going on with this failure?

Actually the enclosing @testset of this @test is copied from L13-L43 of test/bicgstabl.jl. I didn't pay too much attention to the copied tests. Let me try to run the tests on a Windows box and get back to you.

I did pkg> test IterativeSolvers in a Windows box with the continue block (L35-L37 of test/bicgstabl.jl) commented out, and all the tests run fine. I will push the commit without the continue block in test/cg.jl. Let's see if the tests succeed...

No failure during CI (see below). I will probably open another PR to remove the continue block from test/bicgstabl.jl.

stevengj · 2021-03-16T16:19:58Z

test/cg.jl

+        # Do an exact LU decomp of a nearby matrix
+        F = lu(A + rand(T, n, n))
+        x3, his3 = cocg(A, b, Pl = F, maxiter = 100, reltol = reltol, log = true)
+        @test norm(A * x3 - b) / norm(b) ≤ reltol


For tests like this you can use ≈. It should be equivalent to:

@test A*x3 ≈ b rtol=reltol

(see isapprox).

Hmmm. I cannot seem to specify the keyword arguments of isapprox in its infix form:

julia> VERSION v"1.5.4-pre.0" julia> 0 ≈ 0 true julia> 0 ≈ 0, rtol=1e-8 ERROR: syntax: "0" is not a valid function argument name around REPL[23]:1

Is this a new capability introduced in version > 1.5.4?

Don't use the comma, and have it in a @test:

julia> @test 0 ≈ 0 rtol=1e-8 Test Passed

stevengj · 2021-03-19T00:46:48Z

src/cg.jl

+    cocg!(x, A, b; kwargs...) -> x, [history]
+
+Same as [`cg!`](@ref), but uses the unconjugated dot product instead of the usual,
+conjugated dot product.


You should say that this is for the case of complex-symmetric (not Hermitian) matrices A.

src/cg.jl

haampie

Can you also list this method in the docs? https://julialinearalgebra.github.io/IterativeSolvers.jl/dev/linear_systems/cg/ maybe as part of this page?

haampie · 2021-03-22T10:15:18Z

src/cg.jl

+
+Same as [`cg!`](@ref), but uses the unconjugated dot product (`xᵀy`) instead of the usual,
+conjugated dot product (`x'y`) in the algorithm.  It is for solving linear systems with
+matrices `A` that are complex-symmetric (`Aᵀ == A`) rahter than Hermitian (`A' == A`).


Suggested change

matrices `A` that are complex-symmetric (`Aᵀ == A`) rahter than Hermitian (`A' == A`).

matrices `A` that are complex-symmetric (`Aᵀ == A`) rather than Hermitian (`A' == A`).

Thanks for finding out the typo! Will correct this in the next commit.

haampie · 2021-03-22T10:18:21Z

src/common.jl

+_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))


Suggested change

_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.

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 ...?

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.

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

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.

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.

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.

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.)

@haampie and @stevengj, any thoughts on this?

wsshin · 2021-03-22T14:40:47Z

Can you also list this method in the docs? https://julialinearalgebra.github.io/IterativeSolvers.jl/dev/linear_systems/cg/ maybe as part of this page?

So this is the part I was curious about. I haven't created docs for Julia packages before. Is there a quick guide for creating docs?

haampie · 2021-03-22T15:18:16Z

You can do

$ cd path/to/IterativeSolvers/docs
$ julia
julia> using Revise

(@v1.5) pkg> activate .
 Activating environment at `~/.julia/dev/IterativeSolvers/docs/Project.toml`

(docs) pkg> instantiate

julia> include("make.jl") # repeatedly

and edit docs/* and open build/index.html in your browser

stevengj · 2021-03-22T18:38:50Z

@wsshin, see https://juliadocs.github.io/Documenter.jl/stable/ … you might add a note to https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/master/docs/src/linear_systems/cg.md and https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/master/docs/src/index.md

(The above commands are only only needed to build the doc pages locally; they will be built automatically when a PR is merged.)

wsshin · 2022-01-05T03:25:21Z

@haampie, it took a while for me to add the documentation for COCG. Hopefully the latest commits are good to be merged. Note that CI tests are failing for Julia 1, but the failures are from gmres, not from cg or cocg.

Add conjugate orthogonal conjugate gradient (COCG) method for complex…

c5b440e

… symmetric linear systems

haampie requested changes Jan 25, 2021

View reviewed changes

Remove cocg.jl; implement COCG as modified CG using unconjugated dot …

91208cd

…product

Add more unit tests for COCG

3337884

wsshin commented Jan 26, 2021

View reviewed changes

wsshin requested a review from haampie January 27, 2021 14:15

haampie requested changes Jan 31, 2021

View reviewed changes

wsshin requested a review from haampie February 13, 2021 16:10

wsshin mentioned this pull request Feb 22, 2021

COCG and COCR #288

Open

Define cocg() using DotType in cg()

d5b31f4

stevengj reviewed Mar 8, 2021

View reviewed changes

wsshin added 3 commits March 10, 2021 19:12

Rename DotType AbstractDot

97315be

Rename dot_type dotproduct

593e88c

Improve efficiency of unconjugated norm and dot

835a894

stevengj reviewed Mar 16, 2021

View reviewed changes

Utilize ≈ in test/cg.jl

9c47e7f

Remove unnecessary if ...; continue; end block from test/cg.jl

3cd7969

stevengj reviewed Mar 19, 2021

View reviewed changes

src/cg.jl Show resolved Hide resolved

Update docstring for cocg!()

d16fecb

haampie requested changes Mar 22, 2021

View reviewed changes

wsshin added 5 commits April 12, 2021 12:26

Fix typo in comment

04c3c16

Remove unnecessary comments

4800129

Merge branch 'master' into cocg

9a7fb26

Use Conjugate Gradient (no plural) as full name of CG

f3710c3

Add COCG in documentation

b457247

stevengj approved these changes Jan 5, 2022

View reviewed changes

wsshin requested a review from haampie January 16, 2022 16:12

wsshin requested a review from dkarrasch February 2, 2022 01:47

	matrices `A` that are complex-symmetric (`Aᵀ == A`) rahter than Hermitian (`A' == A`).
	matrices `A` that are complex-symmetric (`Aᵀ == A`) rather than Hermitian (`A' == A`).

	_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

Add COCG method for complex symmetric linear systems #289

Are you sure you want to change the base?

Add COCG method for complex symmetric linear systems #289

Conversation

wsshin commented Jan 25, 2021 • edited Loading

wsshin commented Jan 25, 2021

coveralls commented Jan 25, 2021

Pull Request Test Coverage Report for Build 508517998

💛 - Coveralls

coveralls commented Jan 25, 2021 • edited Loading

Pull Request Test Coverage Report for Build 652039845

Warning: This coverage report may be inaccurate.

Details

💛 - Coveralls

haampie commented Jan 25, 2021 • edited Loading

haampie left a comment

Choose a reason for hiding this comment

wsshin commented Jan 25, 2021

wsshin left a comment

Choose a reason for hiding this comment

wsshin commented Jan 26, 2021

wsshin commented Jan 28, 2021 • edited Loading

haampie commented Jan 28, 2021

haampie Jan 31, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

wsshin Jan 31, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

wsshin Jan 31, 2021 • edited Loading

Choose a reason for hiding this comment

stevengj Mar 5, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

wsshin Mar 6, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

wsshin commented Feb 7, 2021

wsshin commented Mar 6, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj Mar 16, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

haampie left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

haampie Mar 22, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

haampie Mar 23, 2021 • edited Loading

Choose a reason for hiding this comment

wsshin Mar 27, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

wsshin commented Mar 22, 2021

haampie commented Mar 22, 2021 • edited Loading

wsshin commented Jan 25, 2021 •

edited

Loading

coveralls commented Jan 25, 2021 •

edited

Loading

haampie commented Jan 25, 2021 •

edited

Loading

wsshin commented Jan 28, 2021 •

edited

Loading

haampie Jan 31, 2021 •

edited

Loading

wsshin Jan 31, 2021 •

edited

Loading

wsshin Jan 31, 2021 •

edited

Loading

stevengj Mar 5, 2021 •

edited

Loading

wsshin Mar 6, 2021 •

edited

Loading

stevengj Mar 16, 2021 •

edited

Loading

haampie Mar 22, 2021 •

edited

Loading

haampie Mar 23, 2021 •

edited

Loading

wsshin Mar 27, 2021 •

edited

Loading

haampie commented Mar 22, 2021 •

edited

Loading

stevengj commented Mar 22, 2021 •

edited

Loading

wsshin commented Jan 5, 2022 •

edited

Loading