Skip to content

Instantly share code, notes, and snippets.

@hholst80
Last active December 24, 2023 00:59
Show Gist options
  • Save hholst80/2da72eb771b04ad4430f52cdba071bea to your computer and use it in GitHub Desktop.
Save hholst80/2da72eb771b04ad4430f52cdba071bea to your computer and use it in GitHub Desktop.
import strconv
import math { max, min }
// m-by-n matrix of float 64
// assume C storage (row-major)
pub fn diag(m i64, n i64, lda i64, mut arr []f64, k i64) {
// zero the matrix; potentially faster with vmemset if lda=n
assert lda >= n
zero := i64(0)
for i := zero; i < m; i++ {
for j := zero; j < n; j++ {
arr[i * lda + j] = f64(0) // i+j*lda for fortran column-major storage
}
}
// This early breakout is technically not needed since nelem will be zero.
if (m <= -k) || (k >= n) {
return
}
mut offset := max(zero, k) + lda*max(zero,-k)
nelem := max(zero, min(min(m, n - k), min(n, m + k)))
for i := zero; i < nelem; i++ {
arr[offset] = f64(1)
offset += lda + 1
}
}
fn disp(m i64, n i64, lda i64, arr []f64) {
assert lda >= n
for i := i64(0); i < m; i++ {
for j := i64(0); j < n; j++ {
v := arr[i * lda + j]
unsafe {
s := strconv.v_sprintf('%10f', v)
print(s)
}
}
print('\n')
}
}
fn main() {
mut arr := []f64{len: 100}
lda := i64(10)
m := i64(4)
n := i64(6)
k := i64(3)
diag(m, n, lda, mut arr, k)
disp(m, n, lda, arr)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment