# Custom dot product for GMRES

**URL:** <https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313>\
**Category:** Numerics\
**Tags:** question, package, linearalgebra\
**Created:** [April 22, 2025, 9:29pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313 "2025-04-22T21:29:08Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [April 22, 2025, 9:29pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/1 "2025-04-22T21:29:09Z")

</div>

Is it possible to define a custom inner product for gmres in Krylov.jl ?

---

<div class="post-metadata">

**Author:** ![vchuravy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vchuravy/32/8_2.png) [@vchuravy](https://discourse.julialang.org/u/vchuravy)\
**Post date:** [April 22, 2025, 10:10pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/2 "2025-04-22T22:10:51Z")

</div>

Yes you can implement your own `kdot`

> **[Custom workspaces · Krylov.jl](https://jso.dev/Krylov.jl/stable/custom_workspaces/#Methods-to-overload-for-compatibility-with-Krylov.jl)**
>
> Documentation for Krylov.jl.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 23, 2025, 1:02am UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/3 "2025-04-23T01:02:20Z")

</div>

The example implementations of `kdot` and `knorm` seem wrong? The element type is `FloatOrComplex`, but there are no complex conjugations?

That is, `kdot` should use `res += conj(_x[i,j]) * _y[i,j]` and `knorm` should use `res += abs2(_x[i,j])`, no? (Putting aside the risk of spurious overflow in an un-rescaled norm algorithm.)

Also, in `knorm`, the result should be initialized to `res = zero(real(T))` in order to give a real-type result for complex vectors.

(And why do `kdot` and `knorm` have an `n` argument that is not used? Also, aren’t the loops in the wrong order for locality?)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 23, 2025, 1:25am UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/4 "2025-04-23T01:25:14Z")

</div>

Another alternative might be rescale your problem to implicitly use the desired norm. It depends on what, precisely, you want to do, and why you want to use a different norm or inner product.

Suppose you want to work with a weighted norm \Vert x \Vert\_W = \sqrt{x^\* W x} = \sqrt{(Rx)^\* (Rx)} = \Vert Rx \Vert\_2, where W = R^\* R is a Hermitian positive-definite weight matrix and R is e.g. a Cholesky factor (or perhaps the matrix square root in cases where it is easy like a diagonal W).

Now, suppose you are solving Ax = b and you want GMRES to minimize the residual \Vert Ax-b\Vert\_W = \Vert RAx-Rb\Vert\_2 in your W norm, instead of in the usual Euclidean norm, over x in the Krylov space at each step. By inspection, however, that is ~~equivalent~~ _nearly_ equivalent to applying GMRES to the equations RAx = Rb (_except that the Krylov space changes_). That is, you pass Rb for your right-hand-side and RA as your matrix (perhaps doing the multiplication lazily as a composition, e.g. via LinearMaps.jl), much like a left preconditioner.

(If A is Hermitian/real-symmetric and you want to take advantage of Lanczos, then RA breaks the symmetry, but you can restore the symmetry by a further change of variables x = R^\* y, much like a right preconditioner.)

Note that this is not equivalent to changing _all_ your inner products: the basis for the Krylov space is still orthonormalized in the Euclidean inner product. But I’m not sure why you should care about that?

PS. If your norm is not L2-like, e.g. you want to use an L1 norm or some such radical change, then it is a very different algorithm than GMRES and you probably need to re-think and re-implement it from scratch.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [April 26, 2025, 8:21am UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/5 "2025-04-26T08:21:36Z")

</div>

@amontoison

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [April 26, 2025, 3:01pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/6 "2025-04-26T15:01:13Z")

</div>

@stevengj You are right, I forgot to handle the complex numbers with my example of custom workspace.  
It is not what I use internally by default in Krylov.jl.

I will fix that, the idea was just to illustrate what to define if you don’t have a classic `AbstractVector` where the operations norm, dot, etc… are not defined.

The argument `n` is here to directly call BLAS / rocBLAS / cuBLAS if you want.

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [April 26, 2025, 3:06pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/7 "2025-04-26T15:06:40Z")

</div>

@Veenty We need more details to understand what you want as inner product.  
If you just want to minimize the residual norm in a different metric than the L2 norm, using a left preconditioner (argument `M` in `Krylov.gmres`) could be what you want.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 26, 2025, 8:04pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/8 "2025-04-26T20:04:27Z")

</div>

> [@amontoison](#):
>
> The argument `n` is here to directly call BLAS / rocBLAS / cuBLAS if you want.

Can’t you get that from the size of the arrays? Why would you need to pass it as a separate argument?

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [April 26, 2025, 8:25pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/9 "2025-04-26T20:25:28Z")

</div>

It depends on your custom workspace.  
You can also imagine that it is a `Ptr{Cvoid}` or `CuPtr{Cvoid}`, and you would need to use `unsafe_wrap` on it with the length.  
Also, why would I want to recompute the length of the vector if it is already available?  
You can just ignore this argument if you don’t need it.

In any case, these functions `kdot`, `knorm`, `kaxpy`, etc., are intended for advanced users and specific applications.  
You should define your own versions only if the default dispatch can’t handle your type:

- [Krylov.jl/src/krylov\_utils.jl at main · JuliaSmoothOptimizers/Krylov.jl · GitHub](https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/main/src/krylov_utils.jl#L309-L345)

Almost nobody asked for this feature before, but I needed to make this custom API available for `Oceananigans.jl`, where they use “Fields” (3D offset arrays with specific halo regions) as vectors for Krylov workspaces.

This custom API could also be useful for doing distributed Krylov solvers with a custom “MPI” vector type but I don’t want to diverge from the initial topic.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 26, 2025, 11:04pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/10 "2025-04-26T23:04:55Z")

</div>

> [@amontoison](#):
>
> Also, why would I want to recompute the length of the vector if it is already available?

Because (a) passing/storing redundant data is a notorious magnet for bugs and (b) it’s not “computed”, it’s usually stored in any high-level container data structure (e.g. in a `Vector`).

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [April 26, 2025, 11:57pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/11 "2025-04-26T23:57:55Z")

</div>

@stevengj I agree with your points.  
For (a), I have some consistency checks at the beginning for the lengths of the vectors. I just compute properties of the vectors outside the call to these methods.  
We could discuss the design, but most of these things were developed back when Julia 0.6 (or earlier) was current, and they haven’t evolved much since then.  
These functions were initially created because Julia didn’t dispatch to BLAS automatically, so we handled it manually with these “k” routines.

For (b), I remember that the length of a `CuVector` wasn’t stored in the structure and was recomputed at every `length` call. I just checked, and it seems that was changed at the end of 2022. However, the GPU support in Krylov.jl was added way earlier, back when CuArrays.jl was still around.  
So, quite a long time ago.

@Veenty Please explain what you are trying to achieve. 🙂

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [April 27, 2025, 12:24am UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/12 "2025-04-27T00:24:15Z")

</div>

This is good. But why is this topic put under `optimization` category?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 27, 2025, 11:46am UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/13 "2025-04-27T11:46:29Z")

</div>

I edited the category, thanks.

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [April 27, 2025, 8:59pm UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/14 "2025-04-27T20:59:27Z")

</div>

I’m solving an integral equation problem and wanted to see if constructing the Krylov spaces in the natural dot product of the space would give better results. It all boils down into doing the naive dot product with a diagonal matrix in the middle.  
I did try the suggestion of @stevengj of using that matrix instead as a left preconditioner, but it was giving me pretty bad results.

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [April 28, 2025, 3:39am UTC](https://discourse.julialang.org/t/custom-dot-product-for-gmres/128313/15 "2025-04-28T03:39:58Z")

</div>

Thanks for your explanation @Veenty.  
Did you pass the diagonal matrix D as left preconditioner or sqrt.(D) ?  
You need `sqrt.(D)` for recovering the norm that you want.

If you already provide the correct matrix, I suggest to try to plug a custom `kdot` and `knorm` to quickly test if it helps.