--MatrixOperations.d

face IntegerMatrixOperations is

	function IntegerMatrixMultiply( n:integer, A:array[1..n,1..n] of integer,
			B:array[1..n,1..n] of integer ) return array[1..n,1..n] of integer 

	function IntegerMatrixTranspose(s:integer,B:array [1..s,1..s] of integer) 
		return array [1..s,1..s] of integer

	function gauss_jordan( size:integer, A :  array[1..size,1..size] of floating, b :  array[1..size] of floating)
	   return array[1..size] of floating		



end face	--of MatrixOperations

----------------------------------------------------------------

body imo: IntegerMatrixOperations is

	function IntegerMatrixMultiply( n:integer, A:array[1..n,1..n] of integer,
			B:array[1..n,1..n] of integer ) return array[1..n,1..n] of integer 
	is
	declare variable C :array[1..n,1..n] of integer
	begin
	forall r:integer in 1..n 
		begin
		forall c:integer in 1..n
			declare variable dotProductResult : shared integer := 0;
			begin
			forall i:integer in 1..n
				declare variable aTerm dummy:integer := 0;
				begin
				aTerm := A[r,i]*B[i,c]
				;
				fetchadd(dotProductResult,aTerm,dummy)
				end	--of forall i
			;
			C[r,c] := dotProductResult
			end	--of forall c
		end	--of forall r
	;
	{true}
	return C
	{true}
	end	--of IntegerMatrixMultiply

	function IntegerMatrixTranspose(s:integer,B:array [1..s,1..s] of integer) 
		return array [1..s,1..s] of integer
	is
	declare
		variable C : spread array [1..s,1..s] of integer
	begin
		forall i : integer in 1..s
			begin
			forall j : integer in 1..s
				begin
				C[i,j]:= B[j,i]
				end
			end				
		;
		return C
	end	--of IntegerMatrix Transpose

function gauss_jordan( size:integer, A :  array[1..size,1..size] of floating, b :  array[1..size] of floating)
   return array[1..size] of floating
is
declare
  variable Ap : spread array[1..size,1..size] of floating := A;	
  variable x : spread array[1..size] of floating := b;
 variable j : integer := 1;
begin
{pre:A: all m:integer in 1.. size are A[m,m] <> 0}	--non-zero diagonal
  do
    (j <= size) =>
	begin
	  forall i:integer in 1..size
	  	declare
	  		variable sizePlusOne:integer := size+1;
	  	begin
	  	forall k:integer in j..sizePlusOne
	  		begin
			if
			   i = j => skip
			[]
			   not(i = j) =>
				    if
					(k > size) => x[i] := x[i] - (Ap[i,j]/Ap[j,j])*x[j]
				    []
					not(k > size) =>
						Ap[i,k] := Ap[i,k] - (Ap[i,j]/Ap[j,j])*Ap[j,k]
				    fi
			fi
			end	--of forall k
		end	--of forall i
	  ;
	  j := j+1
	end	--of (j <= size) =>
  od
  ;
  forall i:integer in 1..size 
  	begin
   	x[i] := x[i]/Ap[i,i]
   	end
  ;
  return x
end --of Gauss-Jordan


procedure matrix_invert( size : integer, A : array[1..size,1..size] of floating,
					singular : boolean) 
is
declare
  variable C : spread array[1..size,1..size] of floating := 0;
  variable D : spread array[1..size,1..size] of floating := A; 
  variable i : shared integer:=0;
  variable j : integer:=1;
begin
  singular == false;
  forall i:integer in 1..size 
  	begin 
  	C[i,i] := 1
  	end
--  {defn:I:C = I}
  ;
  {inv:: C = A-D}
  {bound:: size-j >= 0}
  do
    (j <= size) and not(singular) =>
	begin
	--normalize diagonal element
	  if
	    (D[j,j] <> 0) =>
		forall i:integer in 1..size 
			begin
		  	C[j,i]  := C[j,i]/D[j,j] 
		  	&
		  	D[j,i] := D[j,i]/D[j,j]
		  	end
	  []
	    (D[j,j] = 0) =>
	    --pivot matrix by adding a lower row with a non-zero element
		declare
		  variable k : integer:=j+1;
		begin
		  do
		    (k <= size) and (D[k,j] = 0) => k := k+1
		  od
		  ;
		  if
		    (k > size) => singular == true
		  []
		    not(k > size) =>
			begin
			  --pivot found in row k, add to row j
			  forall i:integer in 1..size 
			  	begin
			   	C[j,i] := C[j,i]+C[k,i] 
			   	&
			   	D[j,i] := D[j,i] + D[j,i]
			   	end
			   ;
			  --now normalize
			  forall i:integer in 1..size 
			  	begin
			   	C[j,i]  := C[j,i]/D[j,j] 
			   	&
			   	D[j,i] := D[j,i]/D[j,j]
			   	end
			end
		  fi
		end --of pivoting to gen non-zero diagonal element
	  fi	--of normalizing diagonal element
	--now subtract that row (scaled by D[i,j]) from all the others
	;
	forall i:integer in 1..size 
		begin
		  if
		    (i = j) => skip	--the skip command does nothing
		  []
		    (i <> j) =>
			declare
			 	variable k : shared integer := 0;
			begin	--D[i,j] is the multiple of row j to subtract
			  forall k:integer in 1..size 
			  	begin
			  	C[i,k] := C[i,k] - C[i,k]/D[i,j] 
			  	&
			  	D[i,k]  :=D[i,k] - D[i,k]/D[i,j]
			  	end
			end
		  fi
		 end
	  ;
	  j := j+1
  end
  od
--  {C = A-I AND D=I AND (j>size) OR singular}
  ;
  --copy inverse back to A if not singular
  if
    (singular) => skip
  []
    not(singular) => A := C
  fi
end --of matrix inversion

	
end body 	--body imo: IntegerMatrix Operations