Support of xAy and traceproduct

[jvdp1]

No description provided.

[Copilot]

🟡 Changes recommended

Critical bounds-safety and moderate parallelism and test-control-flow issues remain unresolved.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

Adds xAx, xAy, and block traceproduct operations for COO and CRS sparse matrices.

Changes:

• Exposes and implements the new sparse operations.
• Adds COO/CRS tests and usage examples.
• Documents behavior and complexity.

File summaries

File	Review
test/modtest_crs.f90	Adds CRS tests. Missing rectangular regressions; one failure can be masked by a later check.
test/modtest_coo.f90	Adds COO tests but lacks rectangular xAy coverage.
src/modsparse.f90	Exposes the new APIs; no issues found.
src/modsparse_gen.f90	Critical rectangular xAy out-of-bounds risk; OpenMP directives do not create parallel regions.
README.md	Adds usage examples; no issues found.
doc/documentation.md	Incorrectly describes implementation reuse and traceproduct complexity.

Review details

Suppressed comments (8)

doc/documentation.md:211

• This implementation does not reuse the format-specific matrix-vector product: xAx_gen delegates to xAy_gen, which traverses the stored entries directly. Update the implementation description so users are not given an incorrect execution model.

where *mat* is a `coosparse` or `crssparse` object and *x* has the same length as the matrix dimension. The result is computed as *dot(x, A*x)*, reusing the format-specific matrix-vector product, and therefore runs in O(nnz).

doc/documentation.md:219

• xAy_gen directly accumulates over COO/CRS storage rather than reusing the format-specific matrix-vector product. Please make this description match the implementation.

where *mat* is a `coosparse` or `crssparse` object, *x* has the length of the number of rows of *A*, and *y* has the length of the number of columns of *A* (both hold when *A* is square). `xAx` is the special case `xAx(x) = xAy(x, x)`. The result is computed as *dot(x, A*y)*, reusing the format-specific matrix-vector product, and therefore runs in O(nnz). Unlike `xAx`, `xAy` also works on non-square matrices.

doc/documentation.md:227

• The CRS implementation scans every row from min(r1,c1) through max(r2,c2), so a small off-diagonal symmetric block with distant row/column ranges can scan almost the entire matrix. Its cost is therefore not O(nnz) in the block as stated; either iterate only the union of the two relevant row ranges or document the actual span-based cost.

where *(r1, r2, c1, c2)* is a 1-based, inclusive tuple defining the square sub-block *A(r1:r2, c1:c2)* (i.e. `r2-r1+1 == c2-c1+1`), and *b* is a dense square matrix of that block size. The returned value is `trace(A(r1:r2, c1:c2) * b)`. On CRS the cost is O(nnz) in the block; on COO it is O(nnz) overall; an informative error is raised on unsupported formats. Setting *b* to the identity gives the plain sub-block trace.

src/modsparse_gen.f90:126

• This merge can evaluate x(sparse%ja(c))*y(i) on a non-square CRS matrix even though lsym is false. Column indices can exceed size(x) (or row indices size(y)), causing an out-of-bounds access; use an explicit conditional for the symmetric mirror.

      rv=rv+sparse%a(c)*(x(i)*y(sparse%ja(c))&
            +merge(x(sparse%ja(c))*y(i),0._wp,lsym.and.i.ne.sparse%ja(c)))

src/modsparse_gen.f90:210

• For a wide non-symmetric matrix, a valid off-diagonal square block may have c2 > nrow. Using the column range here then iterates past the CRS ia row array, although the mirror scan is only needed for symmetric matrices. Restrict non-symmetric scans to r1:r2.

    rmin=min(r1,c1)
    rmax=max(r2,c2)

test/modtest_coo.f90:1233

• All COO xAy cases use a square matrix, so the newly documented non-square behavior is not exercised. Add a rectangular case (for example, 2x3 with a nonzero in column 3) and compare against the dense bilinear form; this also catches out-of-range evaluation of the symmetric-only term.

  coo = coosparse(n, lupper = .true., unlog = sparse_unit)

test/modtest_crs.f90:553

• The CRS tests likewise only cover square xAy inputs despite the public non-square contract. Add a rectangular CRS case with unequal vector lengths and a stored column index greater than the row count.

  coo = coosparse(n, lupper = .true., unlog = sparse_unit)
  call coo%setsymmetric()
  call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, aspsd)
  call coo%add(2, 2, 202._wp)
  crs = coo

test/modtest_crs.f90:623

• These trace-product cases all use a square parent matrix. Add a wide rectangular CRS regression whose valid square block lies in columns beyond the row count (for example, rows 1:2 and columns 3:4 of a 2x4 matrix), since that path has distinct row-indexing behavior.

  result = crs%traceproduct(1, n, 1, n, b)

• Files reviewed: 6/6 changed files
• Comments generated: 3
• Review effort level: Balanced

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[Copilot]

Changed all four reduction loops in xAy_gen and traceproduct_gen to paired parallel do/end parallel do directives. Commit a042d74.