reproducer.chpl

use LinearAlgebra;

type eType = real;

config const N = 8;

proc asum1(a : [?asize] ?T) : T // can overflow
{
	return + reduce abs(a);
}

proc damian(inout a : [?asize] eType, w : [?wsize] eType, v : [?vsize] eType) : int
{
	proc fgh(c : [?asize] eType, scale : eType) : (eType, eType, eType)
	{
		c /= scale;
		{
			const i = asize.dim(1).low;
			const zero = 0.0:eType;
			// const s = + reduce(c[..] * c[..]);
			const s = dot(c, c);
			const t = sqrt(s);
			const f = c[i];
			const g = if f < zero then +t else -t;
			const z = f - g;
	
			c[i] = z;
	
			return (z * scale, g, f * g - s);
		}
	}

	const m = asize.dim(1).high;
	const n = asize.dim(2).high;
	var rv : [1..n] eType;
	const zero = 0.0:eType;
	var scale = zero;
	var g = zero;

	for i in 1..n do
	{
		const l = i + 1;
		const cols = l .. n;
		const rows = i .. m;

		rv[i] = scale * g;
		g = zero;

		if i <= m then // left hand reduction
		{
			var r = i..m;
			var aci = a[r, i];
			var s = asum1(aci);

			if s != zero then
			{
				var h : eType;

				( a[i, i], g, h ) = fgh(aci, s);

				forall j in l .. n do
				{
					var acj = a[r, j];
					const z = (+ reduce(aci * acj)) / h;

					acj += z * aci;
					a[r, j] = acj;
				}
			}
			scale = s;
		}
		w[i] = scale * g;

		scale = zero;
		g = zero;

		if i <= m && i != n then // right hand reduction
		{
			var r = l..n;
			var ari = a[i, r];
			var s = asum1(ari);

			if s != zero then
			{
				var h : eType;

				( a[l, l], g, h ) = fgh(ari, s);

				forall j in l .. m do
				{
					var arj = a[j, r];
					const z = (+ reduce(ari * arj)) / h;

					arj += z * ari;
					a[j, r] = arj;
				}
			}
			scale = s;
		}
	}
	return 0;
}

proc show(in x : [?asize] eType)
{
	const m = asize.dim(1).high;
	const n = asize.dim(2).high;

	for i in 1..m do
	{
		writeln(x[i, 1..n]);
	}
}

proc main()
{
	const n = 5;
	const m = 8;
	var y : [1..m, 1..n] int;
	var x : [1..m, 1..n] real;
	var v : [1..n, 1..n] real;
	var w : [1..n] real;

	// Example 1 Matrix from Golub and Reinsch 1970 

	y[1, 1..n] = [ 22, 10, 2, 3, 7];
	y[2, 1..n] = [ 14, 7, 10, 0, 8];
	y[3, 1..n] = [ -1, 12, -1, -11, 3];
	y[4, 1..n] = [ -3, -2, 13, -2, 4];
	y[5, 1..n] = [ 9, 8, 1, -2, 4];
	y[6, 1..n] = [ 9, 1, -7, 5, -1];
	y[7, 1..n] = [ 2, -6, 6, 5, 1];
	y[8, 1..n] = [ 4, 5, 0, -2, 2];
	for i in 1..m do
	{
		for j in 1..n do
		{
			x[i, j] = y[i, j] : real;
		}
	}
	show(x);
	damian(x, w, v);
	show(x);
}