Thank you for these answers.
@ctkelley @cgeoga Having already tried to implement the Gram-Schmidt algorithm in the past, I fear I won’t be able to match the stability and the performance of the library’s QR method. Thus my question was toward using the existing function of LinearAlgebra
to achieve my goal, if possible.
@stevengj Here is a usecase
I want to create an orthogonal matrix distributed on different processors, but this matrix can be very large, so it shouldn’t be generated all at once and then split the columns.
I generate the different partitions on the different processors, and then I orthogonalize the partitions against one another.
function orthonormalize(A::DArray)
for p in workers()
# Compute the local orthonormal partition
@spawnat p A[:l] = Matrix(LinearAlgebra.qr(A[:l]).Q)
# On processors p+1 .. P, orthogonalize against the basis of p
for p2 in workers()[p:length(workers())]
@spawnat p2 begin
A_on_p = @fetchfrom p A[:l]
A[:l] = (I - A_on_p * A_on_p') * A[:l]
end
end
end
end
# main
U = drand((m,m), workers(), [1,4])
orthonormalize(U)
See the line 4, instead I would like to perform the inplace orthonormalization of A[:l]
(the local part of the input matrix).
If I use @spawnat p LinearAlgebra.qr!(A[:l])
on line 4 instead, I don’t understand how to use A[:l]
afterwards to orthogonalize the other partitions, because it contains only the reflectors. I cannot do A_on_p * A_on_p'
for example.