JuliaSIMD / JuliaSIMD/Polyester.jl

Using @threads works, using @batch has no effect

Open
#138 4 comments 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
285
Forks
20
PR merge metrics
No merged PRs in 30d

Description

The following function
```
using LinearAlgebra, Polyester
function invMatU(X ::Matrix{tw}) ::Matrix{tw} where {tw}
m,n = size(X)
Y = zeros(tw,m,n)
_1 = one(tw)
@inbounds Threads.@threads for j = n : -1 : 1 # backward-substitution;
if j≤m
Y[j,j] = inv(X[j,j])
end
@inbounds for i = min(j-1,m) : -1 : 1
Yij = X[i,j] * (j≤m ? Y[j,j] : _1)
@inbounds for k = i+1 : min(j-1,m)
Yij += X[i,k]*Y[k,j]
end
if i≤n
Y[i,j] = -Yij/X[i,i]
end
end
end
return Y
end;
```
computes the inverse of an upper-triangular matrix. If X is not square, it is assumed that it represents a larger square matrix, by extending (`hcat`or `vcat`) by the identity matrix. We can see it works correctly, since
```
m,n=5,10; tw=Int;
X=rand(tw.(0:10),m,n);
triu!(X);
for i=1:min(m,n) X[i,i] = rand([-1,1]) end
@time Xi=invMatU(X);
X_,Xi_=(vcat(t,Matrix{tw}(I,n,n)[m+1:end,:]) for t∈(X,Xi));
display.((X,Xi,X_,Xi_,X_*Xi_==I));
```
returns
```
5×10 Matrix{Int64}:
-1 2 0 1 9 5 4 7 7 1
0 -1 0 6 6 2 2 0 10 1
0 0 1 10 5 3 2 10 3 3
0 0 0 1 1 3 9 0 8 7
0 0 0 0 -1 0 5 8 7 10
5×10 Matrix{Int64}:
-1 -2 0 13 -8 -30 -69 71 -21 -8
0 -1 0 6 0 -16 -52 0 -38 -41
0 0 1 -10 -5 27 113 30 112 117
0 0 0 1 1 -3 -14 -8 -15 -17
0 0 0 0 -1 0 5 8 7 10
10×10 Matrix{Int64}:
-1 2 0 1 9 5 4 7 7 1
0 -1 0 6 6 2 2 0 10 1
0 0 1 10 5 3 2 10 3 3
0 0 0 1 1 3 9 0 8 7
0 0 0 0 -1 0 5 8 7 10
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 1
10×10 Matrix{Int64}:
-1 -2 0 13 -8 -30 -69 71 -21 -8
0 -1 0 6 0 -16 -52 0 -38 -41
0 0 1 -10 -5 27 113 30 112 117
0 0 0 1 1 -3 -14 -8 -15 -17
0 0 0 0 -1 0 5 8 7 10
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 1
true
```
However, replacing `Threads.@threads` by `Polyester.@batch` results in matrix `Xi` being zero. So something isn't right here.

Btw, this was run with `Polyester v0.7.10` on `Julia 1.9.2` in `Linux Mint 20`.

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.