This is somewhat better:
function ModGS2(A::Array{Float64,2})
delta = 1e-16
(m,n) = size(A)
for j = 1: n
q = slice(A, :,j)
q[:] = q / (norm(q)+delta)
A[:,j] = q
for k = j+1: n
a = slice(A, :,k)
r = dot(q, a)
for i =1:m; A[i,k] = a[i] - r*q[i]; end
end
end
return nothing
end
function test()
A = rand(1000, 500);
B = copy(A);
@time ModGS(A);
@time ModGS2(B);
@show any(A .!= B)
end
julia> test()
0.585496 seconds (875.25 k allocations: 2.831 GB, 14.72% gc time)
0.249571 seconds (376.25 k allocations: 17.231 MB)
any(A .!= B) = false
false
On Thursday, June 9, 2016 at 6:28:27 PM UTC-4, [email protected] wrote:
>
> Hi Julia Users
>
> I write a small function to do Gram-Schmidt orthogonalization:
>
> function ModGS(A::Array{Float64,2})
> delta = 1e-16
> (m,n) = size(A)
> for j = 1: n
> q = A[:,j]
> q = q / (norm(q)+delta)
> A[:,j] = q
> for k = j+1: n
> a = A[:,k]
> r = dot(q, a)
> A[:,k] = a - r*q
> end
> end
> return nothing
> end
>
> A = rand(1000, 500);
> @time ModGS(A);
>
> 0.696324 seconds (875.25 k allocations: 2.831 GB, 15.78% gc time)
>
>
> it cause huge amount of memory allocation, Does anybody know the reason? How
> to improve it?
>
>
> I also use memory track tool, this is the output:
>
> - function ModGS(A::Array{Float64, 2})
>
> 65440991 delta = 1e-16
>
> 0 (m,n) = size(A)
>
> 0 for j = 1: n
>
> 4056000 q = A[:,j]
>
> 4040000 q = q / (norm(q)+delta)
>
> 0 A[:,j] = q
>
> 0 for k = j+1: n
>
> 1011972000 a = A[:,k]
>
> 0 r = dot(q, a)
>
> 2019952000 A[:,k] = a - r*q
>
> - end
>
> - end
>
> 0 return nothing
>
> - end
>
> -
>
> -
>
> - A = randn(1000, 500);
>
> - ModGS(A)
>
> -
>
>
> Thanks
>
>
>
>
>
