# Maple version of SciGMark
# specialized version
# author Laurentiu Dragan

Constants := module () # {{{
	export RESOLUTION_DEFAULT, RANDOM_SEED, FFT_SIZE, SOR_SIZE, 
		SPARSE_SIZE_M, SPARSE_SIZE_nz, LU_SIZE, POLY_SIZE;
	RESOLUTION_DEFAULT := 2.0;
	RANDOM_SEED := 101010;

	FFT_SIZE := 4096;
	SOR_SIZE := 400;
	SPARSE_SIZE_M := 4000;
	SPARSE_SIZE_nz := 20000;
	LU_SIZE := 400;
	POLY_SIZE := 160;
end module: # }}}

Stopwatch := proc() # {{{
	local m;
	m := module() 
		export init, strt, stp, rd;
		local reset, seconds, running, total, last_time;

		seconds := proc()
			return time();
		end proc;

		reset := proc()
			running := false; 
			last_time := 0.0; 
			total := 0.0; 
		end proc;

		init := proc()
			reset();
		end proc;

		strt := proc() 
			if not running then
				running := true;
				total := 0.0;
				last_time := seconds(); 
			end if;
		end proc;

		stp := proc()  
			if running then
				total := total + seconds() - last_time; 
				running := false;
			end if;
			return total; 
		end proc;

		rd := proc()   
			if running then
				total := total + seconds() - last_time;
				last_time := seconds();
			end if;
			return total;
		end proc;

		
	end module: 
	m:-init();
	return m;
end proc: # }}}

Random := proc(seed) # {{{
	local m;
	m := module()
		export init, nextDouble;  
		local dm1, haveRange, seed, initialize, m, i, j, m1, m2, mdig; 
		mdig := 32;
		m1 := (2^(mdig - 2)) + ((2^(mdig - 2)) - 1);
		m2 := 2^iquo(mdig, 2);
		dm1 := 1.0 / m1;
		haveRange := false;
		initialize := proc(s)
			local jseed,k0, k1, j0, j1, iloop;
			seed := s;
			m := array(0..16);
			jseed := min(abs(seed), m1);
			if irem(jseed, 2) = 0 then 
				jseed := jseed - 1;
			end if;
			k0 := irem(9069, m2);
			k1 := iquo(9069, m2);
			j0 := irem(jseed, m2);
			j1 := iquo(jseed, m2);
			for iloop from 0 to 16 do 
				jseed := j0 * k0;
				j1 := irem(iquo(jseed,m2) + j0*k1 + j1*k0, 
						iquo(m2,2));
				j0 := irem(jseed, m2);
				m[iloop] := j0 + m2 * j1;
			end do;
			i := 4;
			j := 16;
		end proc;

		nextDouble := proc() 
			local k, nextValue;

			k := m[i] - m[j];
			if k < 0 then k := k + m1; end if;
			m[j] := k;

			if i = 0 then
				i := 16;
			else
				i := i - 1;
			end if;

			if j = 0 then
				j := 16;
			else
				j := j - 1;
			end if;

			if haveRange then
				return left + dm1 * k * width;
			else
				return dm1 * k;
			end if;
		end proc;

		init := proc(seed)
			initialize(seed);
		end proc:
	end module: 
	m:-init(seed);
	return m;
end proc: #}}}

smFFT := proc() # {{{
local m;
m := module()
	export num_flops, transform, inverse, test;
	
	local makeRandom, log2, transform_internal, bitreverse;
	
	num_flops := (N) -> (5.0 * N - 2) * log2(N) + 2 * (N + 1);

	transform := proc(data::array) 
		transform_internal(data, -1);
	end proc;

	inverse := proc(data::array)
		local nd, n, norm, i;
		transform_internal(data, +1);	
		nd := op(2, eval(data));
		nd := op(2, nd) - op(1, nd) + 1;
		n := nd / 2;
		norm := 1 / n;
		for i from 0 to nd - 1 do 
			data[i] := data[i] * norm;
		end do;
		#for i from 0 to nd - 1 do print(data[i]); end do;
	end proc;

	test := proc(data::array)
		local nd, cpy, diff, d, i;
		nd := op(2, eval(data));
		nd := op(2, nd) - op(1, nd) + 1;
		# cpy := array(nd);
		cpy := copy(data);
		transform(data);
		inverse(data);
		diff := 0.0;
		for i from 0 to nd - 1 do 
			d := data[i] - cpy[i];
			diff := diff + d * d;
		end do;
		return sqrt(diff / nd); 
	end proc;

	makeRandom := proc(n::integer)
		local nd, data, i, r;
		nd := 2 * n;
		data = array(0..nd-1);
		r := Random(11);
		for i from 0 to nd - 1 do 
			data[i] := r:-nextDouble();
		end do;
		return data;
	end proc;

	log2 := proc(n::integer)
		local l;
		l := log[2](n);
		if l <> floor(l) then 
			print("FFT: Data length is not a power of 2!: ", n);
		end if;
		return l; 
	end proc;

	transform_internal := proc(data::array, direction::integer)
		local n, logn, dual, bit, w_real, w_imag, theta, s, t, s2, b, 
				i, j, wd_real, wd_imag, a, tmp_real, tmp_imag, 
				z1_real, z1_imag;
		n := op(2, eval(data));
		n := op(2, n) - op(1, n) + 1; # number of elements
		if n = 0 then return; end if;
		n := n/2;
		if n = 1 then return; end if;
		logn := log2(n);

		bitreverse(data);
		
		dual := 1;
		for bit from 0 to logn - 1 do 
			w_real := 1.0;
			w_imag := 0.0;

			theta := evalhf(2.0 * direction * Pi / (2.0 * dual));
			s := sin(theta);
			t := sin(theta / 2.0);
			s2 := 2.0 * t * t;

			for b from 0 by 2 * dual to n-1 do
				i := 2 * b ;
				j := 2 * (b + dual);

				wd_real := data[j] ;
				wd_imag := data[j+1] ;
					
				data[j]	:= data[i] - wd_real;
				data[j+1] := data[i+1] - wd_imag;
				data[i]	:= data[i] + wd_real;
				data[i+1] := data[i+1] + wd_imag;
			end do;
			
			for a from 1 to dual - 1 do 
				tmp_real := w_real - s * w_imag - s2 * w_real;
				tmp_imag := w_imag + s * w_real - s2 * w_imag;
				w_real := tmp_real;
				w_imag := tmp_imag;
				for b from 0 by 2*dual to n - 1 do
					i := 2*(b + a);
					j := 2*(b + a + dual);

					z1_real := data[j];
					z1_imag := data[j+1];

					wd_real := w_real * z1_real 
							- w_imag * z1_imag;
					wd_imag := w_real * z1_imag 
							+ w_imag * z1_real;

					data[j]	:= data[i] - wd_real;
					data[j+1] := data[i+1] - wd_imag;
					data[i]	:= data[i] + wd_real;
					data[i+1] := data[i+1] + wd_imag;
				end do;
			end do;
			dual := dual *2;
		end do;
	end proc;

	bitreverse := proc(data::array)
		local n, i, j, ii, jj, k, tmp_real, tmp_imag;
		n := op(2, eval(data));
		n := (op(2, n) - op(1, n) + 1)/2;
		i := 0;
		j := 0;
		for i from 0 to n-2 do
			ii := 2*i;
			jj := 2*j;
			k := n / 2 ;

			if i < j then
				tmp_real := data[ii];
				tmp_imag := data[ii+1];
				data[ii] := data[jj];
				data[ii+1] := data[jj+1];
				data[jj] := tmp_real;
				data[jj+1] := tmp_imag;
			end if;

			while k <= j do 
				j := j - k ;
				k := k / 2 ; 
			end do;
			j := j + k;
		end do;
	end proc;

end module:
return m;
end proc: # }}}

MonteCarlo := module () #{{{
	export num_flops, integrate;
	local SEED;
	SEED := 113;

	num_flops := proc(Num_samples)
		return (Num_samples)* 4.0;
	end proc;

	integrate := proc (Num_samples)
		local R, under_curve, count, x, y, nsm1;
		R := Random(SEED);

		under_curve := 0;
		nsm1 := Num_samples - 1;
		for count from 0 to nsm1 do
			x := R:-nextDouble();
			y := R:-nextDouble();

			if x*x + y*y <= 1.0 then
				 under_curve := under_curve + 1;
			end if;
		end do;

		return (under_curve / Num_samples) * 4.0;
	end proc; 
end module: # }}}

SOR := proc() # {{{
local m;
m := module()
	export num_flops, execute;
	num_flops := (M,N,num_iterations) -> (M-1)*(N-1)*num_iterations*6.0;

	execute := proc(omega::float, G::array, num_iterations::integer)
		local M, N, omega_over_four, one_minus_omega, Mm1, Nm1, p, i,
			j, Gi, Gim1, Gip1;
		M := op(2, eval(G));
		N := op(2, eval(G[op(1, M)]));

		omega_over_four := omega * 0.25;
		one_minus_omega := 1.0 - omega;

		Mm1 := op(2, M) - 1;
		Nm1 := op(2, N) - 1; 
		for p from 1 to num_iterations do
			for i from op(1, M) + 1 to Mm1  do
				Gi := G[i];
				Gim1 := G[i - 1];
				Gip1 := G[i + 1];
				for j from op(1, N)+1 to Nm1 do
					G[i][j] := omega_over_four  * (Gim1[j] + Gip1[j]  + Gi[j - 1] + Gi[j + 1])  + one_minus_omega * Gi[j];
				end do:
			end do:
		end do:
	end proc:
end module:
return m;
end proc: # }}}

SparseCompRow := proc() # {{{
local m;
m := module()
	export num_flops, matmult;
	num_flops := (N,nz,num_iterations) -> iquo(nz,N)*N*2.0*num_iterations;

	matmult := proc( y::array, val::array, row::array,
		col::array, x::array, NUM_ITERATIONS::float)
		local M, reps, r, s, rowR, rowRp1, i;
		M := op(2, eval(row));

		for reps from 1 to NUM_ITERATIONS do
			for r from op(1, M) to op(2, M) - 1 do
				s := 0.0; 
				rowR := floor(row[r]);
				rowRp1 := floor(row[r + 1]);
				for i from rowR to rowRp1 - 1 do
					s := s + x[floor(col[i])]*val[i];
				end do:
				y[r] := s;
			end do:
		end do:
	end proc:

end module:
return m;
end proc: # }}}

LU := proc() # {{{
local m;
m := module()
	export num_flops, solve, factor;

	num_flops := (N) -> 2.0 * N ^ 3 / 3.0;

	factor := proc(A::array,  pivot::array)
		local N, M, minMN, j, jp, t, i, ab, tA, recp, k, ii, Aj, AiiJ,
			jj, Aii;
		N := op(2, eval(A));
		M := op(2, eval(A[op(1, N)]));

		minMN := min(op(2, M) - op(1, M), op(2, N) - op(1, N));
		for j from op(1, M) to minMN do
			jp := j;

			t := abs(A[j][j]);
			for i from j+1 to op(2, M) do
				ab := abs(A[i][j]);
				if ab > t then 
					jp := i;
					t := ab;
				end if;
			end do;

			pivot[j] := jp;

			if A[jp][j] = 0 then return 1; end if;

			if jp <> j then
				tA := eval(A[j]);
				A[j] := eval(A[jp]);
				A[jp] := eval(tA);
			end if;

			if j < op(2, M) then 
				recp :=  1.0 / A[j][j];
				for k from j+1 to op(2, M) do 
					A[k][j] := A[k][j]*recp;
				end do;
			end if;

			if j < minMN then
				for ii from j+1 to op(2, M) do
					Aii := A[ii];
					Aj := A[j];
					AiiJ := Aii[j];
					for jj from j+1 to op(2, N) do
						 Aii[jj] := Aii[jj] 
							- AiiJ * Aj[jj];
					end do;
				end do;
			end if;
		end do;

		return 0;
	end proc;

	solve := proc(LU::array, pvt::array, b::array)
		local M, N, ii, i, ip, sum, j;
		M := op(2, eval(LU));
		N := op(2, eval(LU[op(1, M)]));
		ii := 0;

		for i from op(1, M) to op(2,M) do
			ip := pvt[i];
			sum := b[ip];
			b[ip] := b[i];
			if ii = 0 then
				for j from ii to i-1 do
					sum := sum - LU[i][j] * b[j];
				end do;
			else 
				if sum = 0.0 then ii := i; end if;
			end if;
			b[i] := sum;
		end do;

		for i from op(2, N) by -1 to op(1, N) do
			sum := b[i];
			for j from i+1 to op(2, N) do
				sum := sum - LU[i][j] * b[j];
			end do;
			b[i] := sum / LU[i][i];
		end do;
	end proc;               

end module:
return m;
end proc: # }}}

kernel := module() # {{{
	export measureFFT, measureMonteCarlo, measureSOR, measureSparseMatmult,
		measureLU, test;
	local RandomMatrix, RandomVector, CopyMatrix, NewVectorCopy,
		normabs, matvec, matvec3;

	measureFFT := proc(N::integer, mintime::float, R::`module`)
		local x, cycles, Q, i, EPS;
		x := RandomVector(2*N, R);
		cycles := 1;
		Q := Stopwatch();

		while true do
			Q:-strt();
			# print("here", cycles);
			for i from 0 to cycles-1 do
				# print("before:", x);
				smFFT():-transform(x); # forward transform
				# print("middle:", x);
				smFFT():-inverse(x);   # backward transform
				# print("after:", x); quit;
			end do;
			# print("there");
			Q:-stp();
			if Q:-rd() >= mintime then break; end if;

			cycles := cycles * 2;
		end do;

		EPS := 1.0e-6;
		if smFFT():-test(x) / N > EPS then
			return 0.0;
		end if;

		return smFFT():-num_flops(N)*cycles / Q:-rd() * 1.0e-6;
	end proc;

	measureSOR := proc(N::integer, min_time::float, R::`module`)
		local Q, G, cycles;
		G := RandomMatrix(N, N, R);

		Q := Stopwatch();
		cycles := 1;
		while true do
			Q:-strt();
			# print("before: G=", G);
			SOR():-execute(1.25, G, cycles);
			# print("after: G=", G); quit;
			Q:-stp();
			if Q:-rd() >= min_time then break; end if;

			cycles := cycles * 2;
		end do:
		return SOR():-num_flops(N, N, cycles) / Q:-rd() * 1.0e-6;
	end proc:

	measureMonteCarlo := proc(min_time, R)
		local Q, cycles, val;
		Q := Stopwatch();

		cycles := 1;
		while true do
			Q:-strt();
			val := MonteCarlo:-integrate(cycles);
			Q:-stp();
			# print(val);
			if Q:-rd() >= min_time then break; end if;

			cycles := cycles * 2;
		end do;
		return MonteCarlo:-num_flops(cycles) / Q:-rd() * 1.0e-6;
	end proc;
	
	measureSparseMatmult := proc(N::integer, nz::integer, min_time::float,
					R::`module`)
		local x, y, nr, anz, val, col, row, r, rowr, step, i, Q, cycles;
		x := RandomVector(N, R);
		y := array(0..N - 1);
		nr := floor(nz / N);
		anz := nr * N;
		val := RandomVector(anz, R);
		col := array(0..anz - 1);
		row := array(0..N);
		row[0] := 0.0;
		for r from 0 to (N - 1) do
			rowr := row[r];
			row[r + 1] := evalhf(rowr + nr);
			step := floor(r / nr);
			if step < 1 then step := 1 end if:
			for i from 0 to (nr - 1) do
				col[floor(rowr + i)] := evalhf(i * step);
			end do:
		end do:
		Q := Stopwatch();
		cycles := 1.;
		while true do
			Q:-strt();
			SparseCompRow():-matmult(y, val, row, 
					eval(col), x, cycles);
			#print(y); quit;
			Q:-stp();
			if Q:-rd() >= min_time then break; end if:
			cycles := cycles * 2;
		end do:
		return SparseCompRow():-num_flops(N, nz, cycles) / Q:-rd() 
			* 1.0e-6;
	end proc:

	measureLU := proc(N::integer, min_time::float, R::`module`)
		local A, lu, pivot, Q, i, cycles, b, x, EPS;

		A := RandomMatrix(N, N,  R);
		lu := array(0..N-1);
		for i from 0 to N-1 do
			lu[i] := array(0..N-1);
		end do;
		pivot := array(0..N-1);

		Q := Stopwatch();

		cycles := 1;
		while true do
			Q:-strt();
			for i from 0 to cycles-1 do
				CopyMatrix(lu, A);
				LU():-factor(lu, pivot);
			end do;
			Q:-stp();
			if Q:-rd() >= min_time then break; end if;

			cycles := cycles * 2;
		end do:

		b := RandomVector(N, R);
		x := NewVectorCopy(b);
		LU():-solve(lu, pivot, x);
		EPS := 1.0e-6;
		if normabs(b, matvec(A,x)) / N > EPS then 
			return 0.0; 
		end if;

		return LU():-num_flops(N) * cycles / Q:-rd() * 1.0e-6;
	end proc:


	RandomMatrix := proc(M::integer, N::integer, R::`module`)
		local A, i, B, j;
		A := array(0..M-1);
        	for i from 0 to (M - 1) do
			B := array(0..N-1);
			for j from 0 to N - 1 do
            			B[j] := R:-nextDouble();
			end do:
			A[i] := eval(B);
		end do:
		return A;
	end proc:
	
	RandomVector := proc(N::integer, R::`module`)
		local A, i;
		A := array(0..N - 1);
		for i from 0 to N - 1 do
			A[i] := R:-nextDouble();
		end do:
		return eval(A);
	end proc:

	CopyMatrix := proc(B::array, A::array)
		local M, N, remainder, i, Bi, Ai, j;
		M := op(2, eval(A));
		N := op(2, eval(A[op(1, M)]));
		remainder := irem (op(2,N)-op(1,N)+1, 4);

		for i from op(1, M) to op(2, M) do
			Bi := B[i];
			Ai := A[i];
			for j from 0 to remainder-1 do
				Bi[j] := Ai[j];
			end do;
			for j from remainder by 4 to op(2, N)-1 do
				Bi[j] := Ai[j];
				Bi[j+1] := Ai[j+1];
				Bi[j+2] := Ai[j+2];
				Bi[j+3] := Ai[j+3];
			end do;
		end do;
	end proc;

	NewVectorCopy := proc(x::array)
		local N, y, i;
		N := op(2, eval(x));

		y := array(N);
		for i from op(1, N) to op(2, N) do
			y[i] := x[i];
		end do;

		return y;
	end proc;

	normabs := proc(x::array, y::array)
		local N, sum, i;
		N := op(2, eval(x));
		sum := 0.0;

		for i from op(1, N) to op(2, N) do
			sum := sum + abs(x[i]-y[i]);
		end do;

		return sum;
	end proc;

	matvec := proc(A::array, x::array)
		local N, y;
		N := op(2, eval(x));
		y := array(N);

		matvec3(A, x, y);

		return eval(y);
	end proc;

	matvec3 := proc(A::array, x::array, y::array)
		local M, N, i, sum, Ai, j;
		M := op(2, eval(A));
		N := op(2, eval(A[op(1, M)]));

		for i from op(1, M) to op(2, M) do
			sum := 0.0;
			Ai := A[i];
			for j from op(1, N) to op(2, N) do
				sum := sum + Ai[j] * x[j];
			end do;

			y[i] := sum;
		end do;
	end proc;

	test := proc()
		local r, m, lu, i, p, b, x, EPS;
		r := Random(1);
		m := RandomMatrix(3, 3, r);
		lu := array(0..3-1);
		for i from 0 to 3-1 do
			lu[i] := array(0..3-1);
		end do;
		CopyMatrix(lu, m);
		p := array(0..3-1);
		print("m=", m);
		LU():-factor(lu, p);
		print("m=", eval(m), "p=", p);
		b := RandomVector(3, r);
		x := NewVectorCopy(b);
		print("b=", b);
		LU():-solve(lu, p, x);
		print("x=", x);
		EPS := 1.0e-5;
		print("matvec", matvec(m, x));
		print("normabs", normabs(b, matvec(m,x)));
		if normabs(b, matvec(m,x)) / 3 > EPS then 
			print("here ");
			return 0.0; 
		end if;
	end proc;
end module: # }}}

commandline := module() # {{{
	export main;
	main := proc(arg) 
		local res, R, min_time, FFT_size, SOR_size, Sparse_size_M, 
			Sparse_size_nz, LU_size;

		min_time := Constants:-RESOLUTION_DEFAULT;

		FFT_size := Constants:-FFT_SIZE;
		SOR_size := Constants:-SOR_SIZE;
		Sparse_size_M := Constants:-SPARSE_SIZE_M;
		Sparse_size_nz := Constants:-SPARSE_SIZE_nz;
		LU_size := Constants:-LU_SIZE;
		#POLY_size := Constants:-POLY_SIZE;
		#RecLU_size := Constants:-RecLU_SIZE;

		# if (args.length > 0) {

			# if (args[0].equalsIgnoreCase("-h")
					# || args[0].equalsIgnoreCase("-help")) {
				# print("Usage: [-large] [minimum_time]");
				# return;
			# }

			# int current_arg = 0;
			# if (args[current_arg].equalsIgnoreCase("-large")) {
				# FFT_size = Constants.LG_FFT_SIZE;
				# SOR_size = Constants.LG_SOR_SIZE;
				# Sparse_size_M = Constants.LG_SPARSE_SIZE_M;
				# Sparse_size_nz = Constants.LG_SPARSE_SIZE_nz;
				# LU_size = Constants.LG_LU_SIZE;
				# POLY_size = Constants.LG_POLY_SIZE;

				# current_arg++;
			# }

			# if (args.length > current_arg)
				# min_time = Double.valueOf(args[current_arg]).doubleValue();
		# }

		res := Array(0..7);
		R := Random(Constants:-RANDOM_SEED);

		print();
		print("SciGMark 1.0 - Maple - specialized");
		print();

		res[1] := kernel:-measureFFT(FFT_size, min_time, R);
		if res[1] = 0.0 then
			print("FFT: ERROR, INVALID NUMERICAL RESULT!");
		else
			print("FFT (", FFT_size, "): ", res[1]);
		end if;

		res[2] := kernel:-measureSOR(SOR_size, min_time, R);
		print("SOR (", SOR_size, "x", SOR_size, "): ", res[2]);
		res[3] := kernel:-measureMonteCarlo(min_time, R);
		print("Monte Carlo : ", res[3]);
		res[4] := kernel:-measureSparseMatmult(Sparse_size_M, 
				Sparse_size_nz, min_time, R);
		print("Sparse matmult (N=", Sparse_size_M, "nz=",
				Sparse_size_nz, "):", res[4]);
		res[5] := kernel:-measureLU(LU_size, min_time, R);
		if res[5] = 0.0 then
			print("LU: ERROR, INVALID NUMERICAL RESULT!");
		else
			print("LU (", LU_size, "x", LU_size, "): ", res[5]);
		end if;
		# res[6] = kernel.measureMultPoly(POLY_size, min_time);
		# print("PolyMult (N=" + POLY_size + "): " + res[6]);
		# res[7] = kernel.measureRecLU(RecLU_size, min_time, R);
		# System.out.print("RecLU (" + RecLU_size + "x" + RecLU_size + "): ");
		# if (res[7] == 0.0)
			# print(" ERROR, INVALID NUMERICAL RESULT!");
		# else
			# print(res[7]);

		# res[0] = (res[1] + res[2] + res[3] + res[4] + res[5] + res[6] + res[7]) / 7;
		res[0] := (res[1] + res[2] + res[3] + res[4] + res[5]) / 5;
		print();
		print("Composite Score: ", res[0]);
	end proc;


end module: # }}}

# {{{ Main program
writeto("maple.txt");commandline:-main();
R := Random(11):
kernel:-measureSOR(Constants:-SOR_SIZE, Constants:-RESOLUTION_DEFAULT, R);
kernel:-measureMonteCarlo(Constants:-RESOLUTION_DEFAULT, R);
kernel:-measureSparseMatmult(3, 5, 1., R);

quit;