diag.v

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)
}