CG and GMRES Incompatible with Systems of PDEs

[ChrisRackauckas]

CG and GMRES Incompatible with Systems of PDEs

[lopezm94]

Do you mean when b is a matrix? For me the program stops because of KrylovSubspace (not dot), which would require some change in the internals.

julia> x = rand(20,20); A = rand(20,20); b = rand(20,20);

julia> cg(A,b)
ERROR: DimensionMismatch("dimensions must match")
 in promote_shape at ./operators.jl:211

julia> cg!(x,A,b)
ERROR: MethodError: `init!` has no method matching init!(::IterativeSolvers.KrylovSubspace{Float64,Array{Float64,2}}, ::Array{Float64,2})
Closest candidates are:
  init!{T}(::IterativeSolvers.KrylovSubspace{T,OpT}, ::Array{T,1})

[ChrisRackauckas]

Yes, that's the error which stops it from working. Either cg needs to be changed in IterativeSolvers, or an alternative version needs to be implemented just for this.

[lopezm94]

I don't really know what the future plan for KrylovSubspace should be for handling matrices.

@jiahao

[ChrisRackauckas]

Now that I think about, if I'm not mistaken, this can just be done with a comprehension over the iterative solvers. For example, when solving AX=B with X a matrix, the values of X[i,:] only depend on B[i,:] (think about the way it multiplies). So actually this is just shorthand for n independent Ax=b problems, where the A happens to be the same each time. So a comprehension or loop over the columns with cg calls solves it.

(Would there be any efficiency gains by handling it somewhat similarly inside the cg loop? I can't see a reason why there would be)

Should there be type annotation on the b and x in cg then to clarify that it's for vectors? And would this wrapper for matrices be welcomed in IterativeSolvers.jl itself, or should I make it and keep it here?

[lopezm94]

Does the stopping criteria change for matrices? Yeah what you are saying could be done easily in a patch, but at the moment the pull request are merging slowly (might take some time).

[ChrisRackauckas]

I think the stopping criteria can be done either way. Theoretically there is none so it's just what's practical. In practice, by doing a loop over cg it will go until the tolerance is satisfied for each, whereas if it was implemented in the loop it could check the convergence of the whole B-AX in some suitable norm. By the Triangle Inequality the second form would be a weaker condition, though one could always adjust the tolerances to account for it.

However, in the case of PDEs it may be more practical to allow the different columns to have the errors/tolerances applied separately because they may differ vastly in their values (one case where this may come up is a PDE where the components are kept in the columns of X, but some substances have much higher values, and thus would dominate the norm allowing for a lot of error in the other values. This could be terrible in cases where relative tolerance dominates).

[ChrisRackauckas]

It works. However, I wasn't able to use views, so it has some extra allocations. Noted in #33