alter procedure dbo.gaussianElimination
(
	@matrix nvarchar(max)
	,@verbose bit = 1
)
as
/*
Example usage:
To solve Ax=B

    |25  5  1|       |106.8 |
A = |64  8  1| , B = |177.21|
    |144 12 1|       |279.21|

exec dbo.gaussianElimination
	N'<matrix>
		<row values="25,5,1,106.8"/>
		<row values="64,8,1,177.21"/>
		<row values="144,12,1,279.21"/>
	</matrix>', @verbose = 1
*/
set nocount on
set xact_abort on

--=================== SETUP ===============================================================================
-- parse the input into a table
declare	@handle int

declare	@vals table
(
	m smallint identity(0,1) primary key
	,vals varchar(max)
)

exec sp_xml_preparedocument @handle output, @matrix

insert	@vals
(
	vals
)
select	vals
from	openxml (@handle, '/matrix/row',2)
	with (vals  varchar(max) '@values')

exec sp_xml_removedocument @handle

-- validate input (basic)
if exists(select * from @vals where left(vals,1) not like '[-+0-9 ]' or right(vals,1) not like '[0-9 ]' or patindex('%[^-+.,0-9 ]%',vals) > 0)
begin
	raiserror('Malformed input in matrix values this row:',16,1)
	select m as row, vals from @vals where left(vals,1) not like '[-+0-9 ]' or right(vals,1) not like '[0-9 ]' or patindex('%[^-+.,0-9 ]%',vals) > 0
	return
end

if (select count(distinct len(vals)-len(replace(vals,',',''))) from @vals) <> 1
begin
	raiserror('The matrix is not symmetrical, different rows have different number of values',16,1)
	select len(vals)-len(replace(vals,',','')) as [#values], vals as inputMatrix from @vals order by m
	return
end

if @verbose = 1
	select vals as inputMatrix from @vals order by m

-- create our matrix table, and insert the provided matrix
create table #matrix
(
	m smallint	-- the row # 0 based ix
	,n smallint	-- the col # 0 based ix
	,value float
	,primary key(m,n)
)

declare @m smallint
	,@n smallint
	,@valuelist varchar(max)
	,@ix int
	,@value float
	,@i smallint
	,@j smallint
	,@maxi smallint
	,@k smallint
	,@u smallint

select	@m = 0
	,@n = 0

-- parse the values into our matrix table
while 1=1
begin
	set @ix = 0
	set @n = 0
	set @valuelist = (select vals from @vals where m = @m)
	if @valuelist is null break
	while 1=1
	begin
		select	@ix = charindex(',',@valuelist)
			,@value = case when @ix = 0 then @valuelist else left(@valuelist,@ix-1) end
			,@valuelist = substring(@valuelist,@ix+1,len(@valuelist))
		insert	#matrix(m,n,value)
		values	(@m,@n,@value)
		if @ix = 0 break
		set @n = @n + 1
	end
	set @m = @m+1
end

select	@m = (select 1+max(m) from #matrix)
	,@n = (select 1+max(n) from #matrix)
--=================== END SETUP ===============================================================================


--=================== GAUSSIAN ELIMINATION ===============================================================================
-- Adapted from http://en.wikipedia.org/wiki/Gaussian_elimination
select	@i = 0 ,@j = 0

while 1=1
begin
	-- find pivot in column @j starting in row @i
	set @maxi = (select top 1 m from #matrix where m >= @i and n = @j order by abs(value) desc)
	if (select value from #matrix where n = @j and m = @maxi) <> 0
	begin
		-- swap rows i and maxi, but do not change the value of i
		-- Now A[i,j] will contain the old value of A[maxi,j]
		update	#matrix
		set	m = case when m = @i then @maxi else @i end
		where	m in(@i,@maxi)
			and @i <> @maxi

		-- divide each entry in row i by A[i,j]
		-- Now A[i,j] will have the value 1. important step, this property is used in the back-elimination
		update	#matrix
		set	value = value / (select value from #matrix where m = @i and n = @j)
		where	m = @i

		-- subtract A[u,j] * row i from all rows below @i (pivot row)
		update	mv1
		set	value = mv1.value - ((select value from #matrix where m = mv1.m and n = @j) * mv2.value)
		from	#matrix mv1
			join #matrix mv2
				on mv1.n = mv2.n
		where	mv1.m > @i
			and mv2.m = @i
	end
	select	@i = @i+1 ,@j = @j+1
	if @i = @m or @j = @n-1
		break
end
-- back substitution
select	@i = @m-1, @j = @n-2
while 1=1
begin
	-- the back substitution presupposes a row-echelon matrix with 1 as the leading value in each row
	if (select count(*) from #matrix where m = @i and n < @n-1 and value <> 0) = 1
	begin
		set	@value = (select value from #matrix where m = @i and n = @n-1)
		update	m
		set	value = case when n = @j then 0 else value - @value * (select value from #matrix m2 where n = @j and m.m = m2.m) end
		from	#matrix m
		where	n in(@j,@n-1) and m < @i
	end
	else
	begin	
		set @j=@j+1 -- stay on same column 
	end

	select	@i = @i-1, @j = @j-1
	if @i < 0 or @j < 0
		break
end		
--=================== END GAUSSIAN ELIMINATION ===============================================================================

--=================== DISPLAY RESULTS ===============================================================================
declare @res varchar(max)
	,@plus varchar(1)
declare @restab table(result varchar(max))
set	@i = 0
while @i < @m
begin
	select	@res = ''
		,@j = 0
	while @j < @n - 1
	begin
		select	@value = (select value from #matrix where m = @i and n = @j)
			,@plus = ''
		if @value <> 0
		begin
			if sign(@value) = 1 and @res <> ''
				set @plus = '+'
			if abs(@value) = 1
				set @res = @res + @plus + 'x' + ltrim(@j) + ' '
			else
				set @res = @res + @plus + case when @value-1.0*floor(@value) = 0.0 then ltrim(cast(@value as int)) else ltrim(str(@value,14,6)) end + '*x' + ltrim(@j) + ' '
		end
		set @j = @j + 1
	end
	if @res <> ''
	begin
		set @res = @res + ' = ' + (select case when value-1.0*floor(value) = 0.0 then ltrim(cast(value as int)) else ltrim(str(value,14,6)) end from #matrix where m = @i and n = @j)
		insert @restab values(@res)
	end
	else
	begin
		-- check unsolvable situation 0 = x
		if (select value from #matrix where m = @i and n = @j) <> 0
		begin
			set @res = '0 = ' + (select case when value-1.0*floor(value) = 0.0 then ltrim(cast(value as int)) else ltrim(str(value,14,6)) end from #matrix where m = @i and n = @j)
			insert @restab values(@res+' -- unsolvable')
		end
	end
	set @i = @i + 1
end
select result from @restab order by result

if @verbose = 1 and not exists(select * from @restab where result like '%unsolvable')
begin
	-- this is to set up testing scenario
	set @n = 0
	declare @declaration nvarchar(max); set @declaration = 'set nocount on'
	while @n < (select max(n) from #matrix)
	begin
		set @declaration = @declaration + char(10) + 'declare @x' + ltrim(@n) + ' float'
		if @n = (select max(n)-1 from #matrix)
			update @vals set vals = stuff(vals,charindex(',',vals),1,' * @x' + ltrim(@n) + ') --= ' )
		else
			update @vals set vals = stuff(vals,charindex(',',vals),1,' * @x' + ltrim(@n) + ') + (' )
		set @n = @n + 1
	end
	-- display test scenario
	select @declaration as testingScenario
	union all
	select replace(result,'x','set @x') from @restab where len(result)-len(replace(result,'x','')) = 1
	union all
	select 'select (' + vals from @vals
end

if @verbose = 1 --or @verbose = 0
begin
	print 'Resultant row echelon matrix after solving the augmented input matrix' + char(10) + replicate('=',120)
	-- display the #matrix
	declare	@pvt1 varchar(max)
		,@pvt2 varchar(max)
		,@sql nvarchar(max)
	select	@n = 1
		,@pvt1 = '[0] as n0'
		,@pvt2 = '[0]'
	
	while @n < (select max(n)+1 from #matrix)
	begin
		select	@pvt1 = @pvt1 + ',[' + ltrim(@n) + '] as n' + ltrim(@n)
			,@pvt2 = @pvt2 + ',[' + ltrim(@n) + ']'
			,@n = @n + 1
	end
	
	select	@sql = replace(replace(
	'select		''m'' + ltrim(m) as m,>pivotselect<
	from		(
				select	n,m,value
				from	#matrix
			) matrix
	pivot		(
				max(value) for n in(>pivotlist<)
			) as pvt
	order by	m'
	,'>pivotselect<',@pvt1),'>pivotlist<',@pvt2)
	exec sp_executesql @sql
end
--=================== END DISPLAY RESULTS ===============================================================================


-- FINISH
-- cleanup
drop table #matrix
return 0

quote:

---

Originally posted by jezemine

but why?

my instinct is this kind of thing is better done in compiled code with library such as LAPACK or LAPACK++

---