Skip to content

Add BLAS triangular solve derivatives #2825

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jlperla wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Add BLAS triangular solve derivatives #2825

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b11bdaf
feat: add BLAS triangular solve derivatives
jlperla May 14, 2026
8f3163c
fix: simplify triangular BLAS helper
jlperla May 14, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
137 changes: 124 additions & 13 deletions enzyme/Enzyme/BlasDerivatives.td
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ def Rows : MagicInst; // given a transpose, normal rows, normal cols get the tru
def Concat : MagicInst;

def ShadowNoInc : MagicInst;
def Dep : MagicInst;

class Binop<string _s, list<string> _tys> {
string s = _s;
Expand Down Expand Up @@ -89,8 +90,13 @@ class IntMatchers<string _inty, string _outty, list<string> _before, list<string

def side_to_trans : IntMatchers<
"charType", "charType",
["'l'", "'L'", "'R'", "'r'"],
["'n'", "'N'", "'T'", "'t'"]
["'l'", "'L'", "141", "0", "'R'", "'r'", "142", "1"],
["'n'", "'N'", "'N'", "'N'", "'T'", "'t'", "'T'", "'T'"]
>;
def side_to_rtrans : IntMatchers<
"charType", "charType",
["'l'", "'L'", "141", "0", "'R'", "'r'", "142", "1"],
["'t'", "'T'", "'T'", "'T'", "'N'", "'n'", "'N'", "'N'"]
>;

def is_diag_int : IntMatchers<
Expand All @@ -107,6 +113,7 @@ def is_nonunit : MagicInst;
def First : MagicInst;
def Lookup : MagicInst;
def LoadLookup : MagicInst;
def mat_ld : MagicInst;

class FAdd<string _tmp=""> {
string unused = _tmp;
Expand Down Expand Up @@ -163,6 +170,10 @@ class MemcpyMatAdd<bit _shouldZero> {
bit shouldZero = _shouldZero;
}

class MatAdd<bit _shouldZero> {
bit shouldZero = _shouldZero;
}

class For<string idx_, bit offset_> {
string idx = idx_;
bit offset = offset_;
Expand Down Expand Up @@ -672,6 +683,45 @@ def trtrs : CallBlasPattern<(Op $layout, $uplo, $trans, $diag, $n, $nrhs, $A, $l
]
>;

def trsv : CallBlasPattern<(Op $layout, $uplo, $trans, $diag, $n, $A, $lda, $x, $incx),
["x"],
[cblas_layout, uplo, trans, diag, len, mld<["uplo", "n", "n"]>, vinc<["n"]>],
[
/* A */
(Seq<["tri", "triangular", "n"], [], 1>
(BlasCall<"lacpy"> $layout, $uplo, $n, $n, (Shadow $A), use<"tri">, $n),
(BlasCall<"ger">
$layout,
$n,
$n,
Constant<"-1">,
(Rows $trans,
(Concat (Shadow $x)),
(Concat $x)),
(Rows $trans,
(Concat $x),
(Concat (Shadow $x))),
use<"tri">, $n),

(BlasCall<"copy">
(ISelect (is_nonunit $diag), ConstantInt<0>, $n),
(First (Shadow $A)), (Add $lda, ConstantInt<1>),
use<"tri">, (Add $n, ConstantInt<1>)),

(BlasCall<"lacpy"> $layout, $uplo, $n, $n, use<"tri">, $n, (Shadow $A))
),
/* x */ (BlasCall<"trsv"> $layout, $uplo, transpose<"trans">, $diag, $n, $A, (ld $A, Char<"N">, $lda, $n, $n), (Shadow $x)),
]
,
(Seq<["tmp", "vector", "n"], [], 0>
(BlasCall<"copy"> $n, (Dep (Shadow $A), (Dep (Shadow $x), $x)), use<"tmp">, ConstantInt<1>),
(BlasCall<"trmv"> $layout, $uplo, $trans, $diag, $n, (Shadow $A), (Dep (Shadow $x), use<"tmp">), ConstantInt<1>),
(BlasCall<"axpy"> $n, Constant<"-1">, (Dep (Shadow $A), use<"tmp">), ConstantInt<1>, (Shadow $x)),
(BlasCall<"axpy"> (ISelect (is_nonunit $diag), ConstantInt<0>, $n), Constant<"1">, (Dep (Shadow $A), $x), (Shadow $x)),
(BlasCall<"trsv"> $layout, $uplo, $trans, $diag, $n, $A, (ld $A, Char<"N">, $lda, $n, $n), (Shadow $x))
)
>;

def spr2 : CallBlasPattern<(Op $layout, $uplo, $n, $alpha, $x, $incx, $y, $incy, $ap),
["ap"],
[cblas_layout, uplo, len, fp, vinc<["n"]>, vinc<["n"]>, ap<["n"]>],
Expand Down Expand Up @@ -718,33 +768,78 @@ def trsm: CallBlasPattern<(Op $layout, $side, $uplo, $transa, $diag, $m, $n, $al
[cblas_layout, side, uplo, trans, diag, len, len, fp, mld<["side", "m", "n"]>, mld<["m","n"]>],
[
/* alpha */ (AssertingInactiveArg),
/* A */ (AssertingInactiveArg),
/* B */ (AssertingInactiveArg),
/* A */ (Seq<["tmp", "product", "m", "n"], [], 1>
(BlasCall<"lacpy"> $layout, Char<"G">, $m, $n, (Shadow $B), use<"tmp">, (mat_ld $layout, $m, $n)),
(BlasCall<"trsm"> $layout, $side, $uplo, transpose<"transa">, $diag, $m, $n, Constant<"1.0">, $A, (ld $A, (side_to_trans $side), $lda, $n, $m), use<"tmp">, (mat_ld $layout, $m, $n)),
(Seq<["tri", "side_square", "side", "m", "n"], [], 1>
(BlasCall<"lacpy"> $layout, $uplo, (ISelect (is_left $side), $m, $n), (ISelect (is_left $side), $m, $n), (Shadow $A), use<"tri">, (ISelect (is_left $side), $m, $n)),
(BlasCall<"gemm">
$layout,
(side_to_trans $side),
(side_to_rtrans $side),
(ISelect (is_left $side), $m, $n),
(ISelect (is_left $side), $m, $n),
(ISelect (is_left $side), $n, $m),
Constant<"-1">,
(Rows (side_to_trans $side),
(Rows $transa,
(Concat use<"tmp">, (mat_ld $layout, $m, $n)),
(Concat $B, $ldb)),
(Rows $transa,
(Concat $B, $ldb),
(Concat use<"tmp">, (mat_ld $layout, $m, $n)))),
(Rows (side_to_trans $side),
(Rows $transa,
(Concat $B, $ldb),
(Concat use<"tmp">, (mat_ld $layout, $m, $n))),
(Rows $transa,
(Concat use<"tmp">, (mat_ld $layout, $m, $n)),
(Concat $B, $ldb))),
Constant<"1">,
use<"tri">, (ISelect (is_left $side), $m, $n)),
(BlasCall<"copy">
(ISelect (is_nonunit $diag), ConstantInt<0>, (ISelect (is_left $side), $m, $n)),
(First (Shadow $A)), (Add $lda, ConstantInt<1>),
use<"tri">, (Add (ISelect (is_left $side), $m, $n), ConstantInt<1>)),
(BlasCall<"lacpy"> $layout, $uplo, (ISelect (is_left $side), $m, $n), (ISelect (is_left $side), $m, $n), use<"tri">, (ISelect (is_left $side), $m, $n), (Shadow $A))
)
),
/* B */ (BlasCall<"trsm"> $layout, $side, $uplo, transpose<"transa">, $diag, $m, $n, $alpha, $A, (ld $A, (side_to_trans $side), $lda, $n, $m), (Shadow $B)),
]
,
(Seq<["tmp", "product", "m", "n"], [], 0>
(Dep (Shadow $alpha), (AssertingInactiveArg)),
(BlasCall<"lascl"> $layout, Char<"G">, ConstantInt<0>, ConstantInt<0>, Constant<"1.0">, $alpha, $m, $n, (Shadow $B), Alloca<1>),
(BlasCall<"lacpy"> $layout, Char<"G">, $m, $n, (Dep (Shadow $A), (Dep (Shadow $B), $B)), $ldb, use<"tmp">, (mat_ld $layout, $m, $n)),
(BlasCall<"trmm"> $layout, $side, $uplo, $transa, $diag, $m, $n, Constant<"1.0">, (Shadow $A), (Dep (Shadow $B), use<"tmp">), (mat_ld $layout, $m, $n)),
(BlasCall<"lascl"> $layout, Char<"G">, ConstantInt<0>, ConstantInt<0>, Constant<"1.0">, Constant<"-1.0">, $m, $n, (Dep (Shadow $A), use<"tmp">), (mat_ld $layout, $m, $n), Alloca<1>),
(MatAdd<true> $layout, $m, $n, (Shadow $B), (Concat (Dep (Shadow $A), use<"tmp">), (mat_ld $layout, $m, $n))),
(MatAdd<false> $layout, (ISelect (is_nonunit $diag), ConstantInt<0>, $m), $n, (Shadow $B), (Concat (Dep (Shadow $A), $B), $ldb)),
(BlasCall<"trsm"> $layout, $side, $uplo, $transa, $diag, $m, $n, Constant<"1.0">, $A, (ld $A, (side_to_trans $side), $lda, $n, $m), (Shadow $B))
)
>;

def uplo_to_normal : IntMatchers<
"charType", "charType",
["'l'", "'L'", "'U'", "'u'"],
["'n'", "'N'", "'T'", "'t'"]
["'l'", "'L'", "122", "0", "'U'", "'u'", "121", "1"],
Comment thread
jlperla marked this conversation as resolved.
Outdated
["'n'", "'N'", "'N'", "'N'", "'T'", "'t'", "'T'", "'T'"]
>;
def uplo_to_trans : IntMatchers<
"charType", "charType",
["'l'", "'L'", "'U'", "'u'"],
["'t'", "'T'", "'N'", "'n'"]
["'l'", "'L'", "122", "0", "'U'", "'u'", "121", "1"],
["'t'", "'T'", "'T'", "'T'", "'N'", "'n'", "'N'", "'N'"]
>;
def uplo_to_side : IntMatchers<
"charType", "charType",
["'l'", "'L'", "'U'", "'u'"],
["'L'", "'L'", "'R'", "'R'"]
["'l'", "'L'", "122", "0", "'U'", "'u'", "121", "1"],
["'L'", "'L'", "'L'", "'L'", "'R'", "'R'", "'R'", "'R'"]
>;
def uplo_to_rside : IntMatchers<
"charType", "charType",
["'l'", "'L'", "'U'", "'u'"],
["'R'", "'R'", "'L'", "'L'"]
["'l'", "'L'", "122", "0", "'U'", "'u'", "121", "1"],
["'R'", "'R'", "'R'", "'R'", "'L'", "'L'", "'L'", "'L'"]
>;


def potrf: CallBlasPattern<(Op $layout, $uplo, $n, $A, $lda, $info),
["A"],
[cblas_layout, uplo, len, mld<["uplo", "n", "n"]>, info],
Expand Down Expand Up @@ -888,4 +983,20 @@ def potrs: CallBlasPattern<(Op $layout, $uplo, $n, $nrhs, $A, $lda, $B, $ldb, $i
),
(BlasCall<"potrs"> $layout, $uplo, $n, $nrhs, $A, (ld $A, Char<"N">, $lda, $n, $n), (Shadow $B), Alloca<1>)
]
,
(Seq<["tmp", "product", "n", "nrhs"], [], 0>
(BlasCall<"lacpy"> $layout, Char<"G">, $n, $nrhs, (Dep (Shadow $A), (Dep (Shadow $B), $B)), $ldb, use<"tmp">, (mat_ld $layout, $n, $nrhs)),
(BlasCall<"trmm"> $layout, Char<"L">, $uplo, (uplo_to_trans $uplo), Char<"N">, $n, $nrhs, Constant<"1.0">, (Dep (Shadow $A), $A), (ld $A, Char<"N">, $lda, $n, $n), (Dep (Shadow $B), use<"tmp">), (mat_ld $layout, $n, $nrhs)),
(BlasCall<"trmm"> $layout, Char<"L">, $uplo, (uplo_to_normal $uplo), Char<"N">, $n, $nrhs, Constant<"1.0">, (Shadow $A), (Dep (Shadow $B), use<"tmp">), (mat_ld $layout, $n, $nrhs)),
(BlasCall<"lascl"> $layout, Char<"G">, ConstantInt<0>, ConstantInt<0>, Constant<"1.0">, Constant<"-1.0">, $n, $nrhs, (Dep (Shadow $A), (Dep (Shadow $B), use<"tmp">)), (mat_ld $layout, $n, $nrhs), Alloca<1>),
(MatAdd<true> $layout, $n, $nrhs, (Shadow $B), (Concat (Dep (Shadow $A), use<"tmp">), (mat_ld $layout, $n, $nrhs))),

(BlasCall<"lacpy"> $layout, Char<"G">, $n, $nrhs, (Dep (Shadow $A), (Dep (Shadow $B), $B)), $ldb, use<"tmp">, (mat_ld $layout, $n, $nrhs)),
(BlasCall<"trmm"> $layout, Char<"L">, $uplo, (uplo_to_trans $uplo), Char<"N">, $n, $nrhs, Constant<"1.0">, (Shadow $A), (Dep (Shadow $B), use<"tmp">), (mat_ld $layout, $n, $nrhs)),
(BlasCall<"trmm"> $layout, Char<"L">, $uplo, (uplo_to_normal $uplo), Char<"N">, $n, $nrhs, Constant<"1.0">, (Dep (Shadow $A), $A), (ld $A, Char<"N">, $lda, $n, $n), (Dep (Shadow $B), use<"tmp">), (mat_ld $layout, $n, $nrhs)),
(BlasCall<"lascl"> $layout, Char<"G">, ConstantInt<0>, ConstantInt<0>, Constant<"1.0">, Constant<"-1.0">, $n, $nrhs, (Dep (Shadow $A), (Dep (Shadow $B), use<"tmp">)), (mat_ld $layout, $n, $nrhs), Alloca<1>),
(MatAdd<true> $layout, $n, $nrhs, (Shadow $B), (Concat (Dep (Shadow $A), use<"tmp">), (mat_ld $layout, $n, $nrhs))),

(BlasCall<"potrs"> $layout, $uplo, $n, $nrhs, $A, (ld $A, Char<"N">, $lda, $n, $n), (Shadow $B), Alloca<1>)
)
>;
130 changes: 129 additions & 1 deletion enzyme/Enzyme/Utils.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2068,6 +2068,129 @@ Function *getOrInsertDifferentialFloatMemcpyMat(
return F;
}

Function *getOrInsertDifferentialFloatMemcpyMatLayout(
Comment thread
jlperla marked this conversation as resolved.
Outdated
Module &Mod, Type *elementType, PointerType *PT, IntegerType *IT,
IntegerType *CT, unsigned dstalign, unsigned srcalign, bool zeroSrc) {
assert(elementType->isFPOrFPVectorTy());
#if LLVM_VERSION_MAJOR < 17
#if LLVM_VERSION_MAJOR >= 15
if (Mod.getContext().supportsTypedPointers()) {
#endif
#if LLVM_VERSION_MAJOR >= 13
if (!PT->isOpaquePointerTy())
#endif
assert(PT->getPointerElementType() == elementType);
#if LLVM_VERSION_MAJOR >= 15
}
#endif
#endif
std::string name = "__enzyme_dmemcpy_" + tofltstr(elementType) +
"_mat_layout_" +
std::to_string(cast<IntegerType>(CT)->getBitWidth()) +
"_" +
std::to_string(cast<IntegerType>(IT)->getBitWidth()) +
(zeroSrc ? "_zero" : "");
FunctionType *FT = FunctionType::get(
Type::getVoidTy(Mod.getContext()), {CT, IT, IT, PT, IT, PT, IT}, false);

Function *F = cast<Function>(Mod.getOrInsertFunction(name, FT).getCallee());

if (!F->empty())
return F;

F->setLinkage(Function::LinkageTypes::InternalLinkage);
#if LLVM_VERSION_MAJOR >= 16
F->setOnlyAccessesArgMemory();
#else
F->addFnAttr(Attribute::ArgMemOnly);
#endif
F->addFnAttr(Attribute::NoUnwind);
F->addFnAttr(Attribute::AlwaysInline);
F->addParamAttr(3, Attribute::NoAlias);
F->addParamAttr(5, Attribute::NoAlias);

BasicBlock *entry = BasicBlock::Create(F->getContext(), "entry", F);
BasicBlock *init = BasicBlock::Create(F->getContext(), "init.idx", F);
BasicBlock *body = BasicBlock::Create(F->getContext(), "for.body", F);
BasicBlock *initend = BasicBlock::Create(F->getContext(), "init.end", F);
BasicBlock *end = BasicBlock::Create(F->getContext(), "for.end", F);

auto layout = F->arg_begin();
layout->setName("layout");
auto M = layout + 1;
M->setName("M");
auto N = M + 1;
N->setName("N");
auto dst = N + 1;
dst->setName("dst");
auto ldst = dst + 1;
ldst->setName("ldst");
auto src = ldst + 1;
src->setName("src");
auto lsrc = src + 1;
lsrc->setName("lsrc");

{
IRBuilder<> B(entry);
Value *l0 = B.CreateICmpEQ(M, ConstantInt::get(IT, 0));
Value *l1 = B.CreateICmpEQ(N, ConstantInt::get(IT, 0));
B.CreateCondBr(B.CreateOr(l0, l1), end, init);
}

PHINode *j;
{
IRBuilder<> B(init);
j = B.CreatePHI(IT, 2, "j");
j->addIncoming(ConstantInt::get(IT, 0), entry);
B.CreateBr(body);
}

{
IRBuilder<> B(body);
PHINode *i = B.CreatePHI(IT, 2, "i");
i->addIncoming(ConstantInt::get(IT, 0), init);

Value *srci = lookup_with_layout(B, elementType, layout, src, lsrc, i, j);
Value *dsti = lookup_with_layout(B, elementType, layout, dst, ldst, i, j);
LoadInst *srcl = B.CreateLoad(elementType, srci, "src.i.l");
LoadInst *dstl = B.CreateLoad(elementType, dsti, "dst.i.l");
auto res = B.CreateFAdd(srcl, dstl);
StoreInst *dsts = B.CreateStore(res, dsti);
StoreInst *srcs = nullptr;
if (zeroSrc)
srcs = B.CreateStore(Constant::getNullValue(res->getType()), srci);
if (dstalign) {
dsts->setAlignment(Align(dstalign));
dstl->setAlignment(Align(dstalign));
}
if (srcalign) {
if (zeroSrc)
srcs->setAlignment(Align(srcalign));
srcl->setAlignment(Align(srcalign));
}

Value *nexti =
B.CreateAdd(i, ConstantInt::get(IT, 1), "i.next", true, true);
i->addIncoming(nexti, body);
B.CreateCondBr(B.CreateICmpEQ(nexti, M), initend, body);
}

{
IRBuilder<> B(initend);
Value *nextj =
B.CreateAdd(j, ConstantInt::get(IT, 1), "j.next", true, true);
j->addIncoming(nextj, initend);
B.CreateCondBr(B.CreateICmpEQ(nextj, N), end, init);
}

{
IRBuilder<> B(end);
B.CreateRetVoid();
}

return F;
}

// TODO implement differential memmove
Function *
getOrInsertDifferentialFloatMemmove(Module &M, Type *T, unsigned dstalign,
Expand Down Expand Up @@ -3567,6 +3690,7 @@ llvm::Optional<BlasInfo> extractBLAS(llvm::StringRef in)
"dot", "scal", "axpy", "gemv", "gemm", "spmv", "syrk", "nrm2",
"trmm", "trmv", "symm", "potrf", "potrs", "copy", "spmv", "syr2k",
"potrs", "getrf", "getrs", "trtrs", "getri", "symv", "lacpy", "trsv",
"trsm",
};
const char *floatType[] = {"s", "d", "c", "z"};
const char *prefixes[] = {"" /*Fortran*/, "cblas_"};
Expand Down Expand Up @@ -4100,7 +4224,11 @@ SmallVector<llvm::Value *, 1> get_blas_row(llvm::IRBuilder<> &B,
if (!cublas) {

if (!byRef) {
cond = B.CreateICmpEQ(trans, ConstantInt::get(trans->getType(), 111));
auto isCblasN =
B.CreateICmpEQ(trans, ConstantInt::get(trans->getType(), 111));
auto isN = B.CreateICmpEQ(trans, ConstantInt::get(trans->getType(), 'N'));
auto isn = B.CreateICmpEQ(trans, ConstantInt::get(trans->getType(), 'n'));
cond = B.CreateOr(isCblasN, B.CreateOr(isN, isn));
} else {
auto isn = B.CreateICmpEQ(trans, ConstantInt::get(trans->getType(), 'n'));
auto isN = B.CreateICmpEQ(trans, ConstantInt::get(trans->getType(), 'N'));
Expand Down
5 changes: 5 additions & 0 deletions enzyme/Enzyme/Utils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -811,6 +811,11 @@ llvm::Function *getOrInsertDifferentialFloatMemcpyMat(
llvm::IntegerType *IT, llvm::IntegerType *CT, unsigned dstalign,
unsigned srcalign, bool zeroSrc);

llvm::Function *getOrInsertDifferentialFloatMemcpyMatLayout(
llvm::Module &M, llvm::Type *elementType, llvm::PointerType *PT,
llvm::IntegerType *IT, llvm::IntegerType *CT, unsigned dstalign,
unsigned srcalign, bool zeroSrc);

/// Create function for type that performs the derivative memmove on floating
/// point memory
llvm::Function *getOrInsertDifferentialFloatMemmove(
Expand Down
28 changes: 28 additions & 0 deletions enzyme/test/Enzyme/ForwardMode/blas/cblas_trsm_f.ll
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
; RUN: %opt < %s %newLoadEnzyme -passes="enzyme" -enzyme-preopt=false -enzyme-detect-readthrow=0 -S | FileCheck %s

target datalayout = "e-m:e-i64:64-f80:128-n8:16:32:64-S128"
target triple = "x86_64-unknown-linux-gnu"

declare dso_local void @__enzyme_fwddiff(...)

declare void @cblas_dtrsm(i32, i32, i32, i32, i32, i32, i32, double, double*, i32, double*, i32)

define void @f(i32 %layout, double* %A, double* %B) {
entry:
call void @cblas_dtrsm(i32 %layout, i32 141, i32 121, i32 111, i32 131, i32 4, i32 3, double 2.000000e+00, double* %A, i32 4, double* %B, i32 4)
ret void
}

define void @active(i32 %layout, double* %A, double* %dA, double* %B, double* %dB) {
entry:
call void (...) @__enzyme_fwddiff(void (i32, double*, double*)* @f, metadata !"enzyme_const", i32 %layout, metadata !"enzyme_dup", double* %A, double* %dA, metadata !"enzyme_dup", double* %B, double* %dB)
ret void
}

; CHECK-LABEL: define internal void @fwddiffef(
; CHECK: call void @cblas_dtrsm
; CHECK: icmp eq i32 %layout, 101
; CHECK: __enzyme_dmemcpy_double_mat_layout_32_32_zero
; CHECK: call void @cblas_dtrsm
; CHECK-NOT: call void bitcast
; CHECK: define internal void @__enzyme_dmemcpy_double_mat_layout_32_32_zero(i32 %layout,
Loading