NVIDIA / NVIDIA/AMGX

[Issue] Device dotc uses unconjugated cuBLAS dotu and breaks complex GMRES orthogonalization

Open
#377 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

bug
Dominant language
Cuda
Stars
692
Forks
197
PR merge metrics
No merged PRs in 30d

Description

Describe the issue

AMGX v2.5.0 appears to use an unconjugated cuBLAS dot product in the device implementation of dotc, which causes incorrect behavior for complex-valued GMRES.

In src/blas.cu, the device implementation of dotc(...) calls:

Cublas::dot(size, a_raw, 1, b_raw, 1, &result);

For complex double precision, Cublas::dot(...) is eventually mapped to cublasZdotu, which computes an unconjugated dot product.

However, the expected semantics of dotc(x, y) for complex-valued vectors are:

$$
x^H y
$$

i.e. the first argument should be conjugated.

AMGX already provides Cublas::dotc(...), which maps to cublasZdotc.

There is also a related issue in the GMRES Arnoldi orthogonalization in src/solvers/gmres_solver.cu. The current code computes:

m_H(k, i) = dot(A,
                m_V_vectors[i + 1],
                m_V_vectors[k]);

After restoring the expected dotc(x,y) = x^H y semantics, this computes

$$
w^H v_k
$$

while Arnoldi / Modified Gram-Schmidt requires

$$
v_k^H w.
$$

Therefore the two arguments need to be swapped for the complex case.

With the original AMGX v2.5.0 implementation, the minimal 3x3 complex system below does not converge after 500 GMRES iterations.

With the following two changes:

- Cublas::dot(size, a_raw, 1, b_raw, 1, &result);
+ Cublas::dotc(size, a_raw, 1, b_raw, 1, &result);

in the device dotc implementations, and

- m_H(k, i) = dot(A, m_V_vectors[i + 1], m_V_vectors[k]);
+ m_H(k, i) = dot(A, m_V_vectors[k], m_V_vectors[i + 1]);

in GMRES, the same 3x3 system converges in 3 iterations to approximately machine precision.

Observed result
Original AMGX v2.5.0 Patched
Matrix size 3x3 complex 3x3 complex
Initial residual 1.270827e+01 1.270827e+01
Total iterations 500 3
Final residual 2.955383e-05 4.981581e-16
Status Did not converge Converged

The original solver terminates with:

Total Iterations: 500
Final Residual: 2.955383e-05
Total Reduction in Residual: 2.325560e-06
***Solver Did Not Converge

The patched solver gives:

Ini     1.270827e+01
0       2.975743e+00
1       1.027928e+00
2       4.981581e-16

Total Iterations: 3
Final Residual: 4.981581e-16
Total Reduction in Residual: 3.919954e-17

Because this is only a 3x3 system and the GMRES restart length is larger than the system dimension, the patched result is also consistent with the expected behavior of full GMRES for this small problem.


Environment information

  • OS: Ubuntu 24.04.4 LTS
  • CUDA runtime: CUDA 13.0
  • MPI version: N/A (CMAKE_NO_MPI=ON)
  • AMGX version: v2.5.0
  • AMGX commit hash: cc1cebdbb32b14d33762d4ddabcb2e23c1669f47
  • NVIDIA driver: 591.86
  • CUDA driver API reported by AMGX: 13.1
  • NVIDIA GPU: NVIDIA GeForce RTX 5070 Ti
  • Build type: Release
  • AMGX mode: dZZI
  • Relevant environment variables: none required for reproduction other than the library path used to load the corresponding AMGX build.

AMGX startup information:

AMGX version 2.5.0
Built on Aug 25 2026
Compiled with CUDA Runtime 13.0, using CUDA driver 13.1

AMGX solver configuration

The following configuration was used:

{
  "config_version": 2,
  "determinism_flag": 1,
  "exception_handling": 1,
  "solver": {
    "scope": "main",
    "solver": "GMRES",
    "preconditioner": {
      "scope": "noprec",
      "solver": "NOSOLVER"
    },
    "max_iters": 500,
    "gmres_n_restart": 10,
    "tolerance": 1e-8,
    "norm": "L2",
    "convergence": "RELATIVE_INI",
    "monitor_residual": 1,
    "store_res_history": 1,
    "print_solve_stats": 1,
    "obtain_timings": 1
  }
}

No preconditioner is used so that the reproducer isolates the behavior of GMRES itself.


Matrix Data

The issue can be reproduced using the following 3x3 complex-valued Matrix Market system:

%%MatrixMarket matrix coordinate complex general
%%NVAMG 1 1 rhs solution
3 3 7
1 1 4.0 0.0
1 2 1.0 1.0
2 1 1.0 -1.0
2 2 3.0 0.0
2 3 1.0 2.0
3 2 1.0 -2.0
3 3 5.0 0.0
3
7.0 5.0
6.0 -4.5
-5.0 -2.5
3
0.0 0.0
0.0 0.0
0.0 0.0

The matrix is

$$
A=
\begin{bmatrix}
4 & 1+i & 0\
1-i & 3 & 1+2i\
0 & 1-2i & 5
\end{bmatrix}.
$$

The right-hand side is

$$
b=
\begin{bmatrix}
7+5i\
6-4.5i\
-5-2.5i
\end{bmatrix}.
$$

A complete reproducer containing the matrix, solver configuration, original log, patched log, environment information, and patch is attached to this issue.


Reproduction steps

AMGX was built without MPI.

Example configuration:

mkdir build
cd build
cmake \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_NO_MPI=ON \
    -DCMAKE_CUDA_ARCHITECTURES=120 \
    ..
make -j8

Run the reproducer using the device complex-double mode:

LD_LIBRARY_PATH="$PWD/build" \
./build/examples/amgx_capi \
    -mode dZZI \
    -m A.mtx \
    -c gmres_complex.json
Original AMGX v2.5.0

The original implementation runs all 500 iterations and does not satisfy the requested tolerance:

Total Iterations: 500
Final Residual: 2.955383e-05
***Solver Did Not Converge
Patched version

Changing the device dotc calls from:

Cublas::dot(...)

to:

Cublas::dotc(...)

and changing the GMRES projection from:

dot(A, m_V_vectors[i + 1], m_V_vectors[k])

to:

dot(A, m_V_vectors[k], m_V_vectors[i + 1])

produces:

Total Iterations: 3
Final Residual: 4.981581e-16

No other solver changes were made for this comparison.


Additional context

For complex vectors, the distinction between the unconjugated and conjugated dot products is fundamental.

For example, for

$$
x=[1,i]^T,
$$

the unconjugated product gives

$$
x^T x = 1+i^2=0,
$$

while the Hermitian inner product gives

$$
x^H x=|1|^2+|i|^2=2.
$$

Therefore an unconjugated dot product cannot be used as the complex inner product required by Arnoldi orthogonalization.

For real-valued inputs this problem is hidden because

$$
x^T y=y^T x,
$$

whereas for complex-valued inputs

$$
x^H y
$$

and

$$
y^H x
$$

are generally different.

The minimal reproducer shows a significant behavioral difference:

Original v2.5.0:
    500 iterations
    final residual = 2.955383e-05
    did not converge

Patched:
    3 iterations
    final residual = 4.981581e-16
    converged

I can also prepare a pull request with the two changes above and a regression test if this diagnosis is confirmed.

amgx-v2.5.0-complex-gmres-reproducer.zip

Contributor guide

No contributing guide indexed for this repository

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Inspect the device dotc implementation in src/blas.cu and the Arnoldi projection in src/solvers/gmres_solver.cu, then build the provided 3x3 complex reproducer with the documented CMake and amgx_capi commands. Verify that complex GMRES converges in about three iterations with the expected residual, and add regression coverage if the project’s test structure supports it.

Written by the indexing model from the issue text.

Assessment

Tech stack
cpp
Domain
hpc
Issue type
Bug
Difficulty
3/5
Estimated time
1-2 days
Activity status
Active
Clarity
Clearly specified
Newbie friendliness
72/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.