For statements involving operations on matrices or vectors, Julia almost always creates data items on the heap. For example, q = A[:,j] allocates a new heap space for A[:,j], and possibly more heap space for q. Even the 'sub' operation creates a heap item.
For this reason, if you are operating on matrices or vectors, especially small ones, and you want high performance, it is best to write matrix operations that occur inside loops as scalar rather than vector operations. (In other words, the Gram-Schmidt program needs to be rewritten with more nested loops in order to get high performance.) There are several macro packages available such as Einsum.jl or Devectorize that allow you to write code in the syntactic style of matrix/vector operations while the macro 'under the hood' generates the necessary for-loops. -- Steve Vavasis 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 > > > > >
