82 lines
2.2 KiB
Text
82 lines
2.2 KiB
Text
-- demo/rosettacode/QRdecomposition.exw
|
|
function vtranspose(sequence v)
|
|
-- transpose a vector of length m into an mx1 matrix,
|
|
-- eg {1,2,3} -> {{1},{2},{3}}
|
|
for i=1 to length(v) do v[i] = {v[i]} end for
|
|
return v
|
|
end function
|
|
|
|
function mat_col(sequence a, integer col)
|
|
sequence res = repeat(0,length(a))
|
|
for i=col to length(a) do
|
|
res[i] = a[i,col]
|
|
end for
|
|
return res
|
|
end function
|
|
|
|
function mat_norm(sequence a)
|
|
atom res = 0
|
|
for i=1 to length(a) do
|
|
res += a[i]*a[i]
|
|
end for
|
|
res = sqrt(res)
|
|
return res
|
|
end function
|
|
|
|
function mat_ident(integer n)
|
|
sequence res = repeat(repeat(0,n),n)
|
|
for i=1 to n do
|
|
res[i,i] = 1
|
|
end for
|
|
return res
|
|
end function
|
|
|
|
function QRHouseholder(sequence a)
|
|
integer columns = length(a[1]),
|
|
rows = length(a),
|
|
m = max(columns,rows),
|
|
n = min(rows,columns)
|
|
sequence q, I = mat_ident(m), Q = I, u, v
|
|
|
|
--
|
|
-- Programming note: The code of this main loop was not as easily
|
|
-- written as the first glance might suggest. Explicitly setting
|
|
-- to 0 any a[i,j] [etc] that should be 0 but have inadvertently
|
|
-- gotten set to +/-1e-15 or thereabouts may be advisable. The
|
|
-- commented-out code was retrieved from a backup and should be
|
|
-- treated as an example and not be trusted (iirc, it made no
|
|
-- difference to the test cases used, so I deleted it, and then
|
|
-- had second thoughts a few days later).
|
|
--
|
|
for j=1 to min(m-1,n) do
|
|
u = mat_col(a,j)
|
|
u[j] -= mat_norm(u)
|
|
v = sq_div(u,mat_norm(u))
|
|
q = sq_sub(I,sq_mul(2,matrix_mul(vtranspose(v),{v})))
|
|
a = matrix_mul(q,a)
|
|
-- for row=j+1 to length(a) do
|
|
-- a[row][j] = 0
|
|
-- end for
|
|
Q = matrix_mul(Q,q)
|
|
end for
|
|
|
|
-- Get the upper triangular matrix R.
|
|
sequence R = repeat(repeat(0,n),m)
|
|
for i=1 to n do -- (logically 1 to m(>=n), but no need)
|
|
for j=i to n do
|
|
R[i,j] = a[i,j]
|
|
end for
|
|
end for
|
|
|
|
return {Q,R}
|
|
end function
|
|
|
|
sequence a = {{12, -51, 4},
|
|
{ 6, 167, -68},
|
|
{-4, 24, -41}},
|
|
{q,r} = QRHouseholder(a)
|
|
|
|
?"A" pp(a,{pp_Nest,1})
|
|
?"Q" pp(q,{pp_Nest,1})
|
|
?"R" pp(r,{pp_Nest,1})
|
|
?"Q * R" pp(matrix_mul(q,r),{pp_Nest,1})
|