Giles method to differentiate solvers and minor edits to docs - #11
Merged
Conversation
Reference LAPACK's CMake build produces libblas.a/liblapack.a rather than librefblas.a, which the Make and Meson builds expect by default. Document the symlink/override workarounds. Tested (Make path only; Meson not installed locally): - Make + symlink (ln -s libblas.a librefblas.a): build succeeds - Make + BLAS_LIB override (BLAS_LIB="-L$LAPACKDIR -lblas"): build succeeds - Meson + -Dlibblas=blas -Dlibblas_path=$LAPACKDIR: build succeeds
Shreyas911
changed the title
sriharikrishna
self-requested a review
August 31, 2026 18:56
Correct and harden the black-box (Giles) derivatives of the eight
triangular-solve routines (d/s/c/z × trsm/trsv) across all four modes.
Source (BLAS/src):
- Wrap fixed-form lines exceeding column 72 that broke compilation with
gfortran defaults (strsm_b.f and the *trsm_bv.f IF/expression lines).
- Forward/vector-forward: pass DIAG='N' (not the caller's DIAG) to the
inner TRMM/TRMV. For a unit-triangular A the diagonal is constant, so
its derivative is zero; forwarding DIAG='U' added a spurious +X and
gave wrong tangents for every DIAG='U' case (adjoint already skipped
the diagonal, so this restores forward/reverse consistency).
- Vector modes (_dv/_bv): dimension the direction arrays by the runtime
nbdirs and drop the DIFFSIZES.inc/nbdirsmax cap, matching the
hand-modified dgemv_dv/dgemv_bv convention (nbdirs no longer limited
to 4).
- Remove the commented-out Tapenade-generated code from the hand-written
Giles sources.
Tests (BLAS/test):
- Rewrite the trs test drivers following the gemv/gemv_vector structure:
multi-size, sweep DIAG in {N,U}. Forward uses full-complex central
differences (holomorphic for TRANS='N'); reverse uses the real-part
Hermitian VJP identity. Verified against finite differences for all
eight routines in all four modes.
Build:
- BLAS/meson.build: fold the trsm/trsv sources into the per-mode file
lists in alphabetical order instead of a trailing catch-all block.
Co-Authored-By: Claude Code <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What
libblas.a, notlibrefblas.a, which the Make/Meson builds expect by default. Documents the naming fix (symlink,BLAS_LIBoverride, or-Dlibblas=blas) and adds a-jparallel-build note.DTRSM,STRSM,CTRSM,ZTRSM,DTRSV,STRSV,CTRSV,ZTRSV. Differentiates the defining relation (op(A)X=alpha*B/op(A)x=b) rather than the substitution algorithm, so each_d/_b/_dv/_bvis just calls to the routine's own undifferentiated BLAS.Fix
Real case: standard Giles adjoint,
_dvia perturb-and-resolve,_bvia transpose-solve on the seed + triangle-restricted outer product (DIAG='U'diagonal zeroed)._dv/_bvgather/scatter around Tapenade's direction-first array layout.Complex case (Wirtinger/CR-calculus convention):
A^H, notA^T.'N'/'C'map directly;'T'needsconj(A)Y=S ⟺ A·conj(Y)=conj(S).TRSM(2 side × 3 transa), 2 forTRSV.Bb = CONJG(ALPHA)*Yb(complexALPHAinTRSM)._d/_dvunaffected — passing the same trans flag toCTRMM/CTRMVhandles conjugation for free.Also fixes the pre-existing duplicate
seed_arraydeclaration in every scalar-mode test file (run_tapenade_blas.pyoutput, declared after executable statements — illegal Fortran). Vector-mode tests unaffected.DTRSMDTRSM,DTRMMSTRSMSTRSM,STRMMCTRSMCTRSM,CTRMMALPHA, 3-way transaZTRSMZTRSM,ZTRMMDTRSVDTRSV,DTRMVINCXSTRSVSTRSV,STRMVCTRSVCTRSV,CTRMVINCXZTRSVZTRSV,ZTRMVTesting
README/library naming:
makemake BLAS_LIB=...-Dlibblas=blas)Giles routines, FD/VJP (complex checks perturb Re+Im, use
conj(·)·(·)inner product):_d_b_dv_bvDTRSM/STRSM/CTRSM/ZTRSMDTRSV/STRSV/CTRSV/ZTRSVAll 32 combinations pass.
Notes for reviewers
run_tapenade_blas.pyrewrites top-level build files on every invocation regardless of--file/--mode— separate issue if useful.