diff --git a/doc/hpl.j b/doc/hpl.j index bc57f43663a52..d2ee5322d4a8e 100644 --- a/doc/hpl.j +++ b/doc/hpl.j @@ -4,38 +4,38 @@ function hpl (A::Matrix, b::Vector, blocksize::Number) -global A; -n = length(A); -A = [A b]; +global A +n = length(A) +A = [A b] -B_rows = 0:blocksize:n; -B_rows(end) = n; -B_cols = [B_rows n+1]; -nB = length(B_rows); -depend = cell(nB, nB); +B_rows = 0:blocksize:n +B_rows(end) = n +B_cols = [B_rows n+1] +nB = length(B_rows) +depend = cell(nB, nB) # Add a ghost row of dependencies to boostrap the computation for j=1:nB - depend{1,j} = true; + depend{1,j} = true end for i=1:(nB-1) # Threads for panel factorizations - I = (B_rows[i]+1):B_rows[i+1]; - [depend{i+1,i}, panel_p] = spawn(panel_factor, I, depend{i,i}); + I = (B_rows[i]+1):B_rows[i+1] + [depend{i+1,i}, panel_p] = spawn(panel_factor, I, depend{i,i}) # Threads for trailing updates for j=(i+1):nB - J = (B_cols[j]+1):B_cols[j+1]; - depend{i+1,j} = spawn(trailing_update, I, J, panel_p, depend{i+1,i},depend{i,j}); + J = (B_cols[j]+1):B_cols[j+1] + depend{i+1,j} = spawn(trailing_update, I, J, panel_p, depend{i+1,i},depend{i,j}) end end # Completion of the last diagonal block signals termination -wait(depend{nB, nB}); +wait(depend{nB, nB}) # Solve the triangular system -x = triu(A[1:n,1:n]) \ A[:,n+1]; +x = triu(A[1:n,1:n]) \ A[:,n+1] return x @@ -48,19 +48,19 @@ function panel_factor(I, col_dep) # Panel factorization global A -n = size (A, 1); +n = size (A, 1) # Enforce dependencies -wait(col_dep); +wait(col_dep) # Factorize a panel -K = I[1]:n; -[A[K,I], panel_p] = lu(A[K,I],'vector','economy'); +K = I[1]:n +[A[K,I], panel_p] = lu(A[K,I],'vector','economy') # Panel permutation -panel_p = K(panel_p); +panel_p = K(panel_p) -done = true; +done = true return (done, panel_p) @@ -73,25 +73,25 @@ function trailing_update(I, J, panel_p, row_dep, col_dep) # Trailing submatrix update global A -n = size (A, 1); +n = size (A, 1) # Enforce dependencies -wait(row_dep, col_dep); +wait(row_dep, col_dep) # Apply permutation from pivoting -K = (I[end]+1):n; -A(I[1]:n, J) = A(panel_p, J); +K = (I[end]+1):n +A(I[1]:n, J) = A(panel_p, J) # Compute blocks of U -L = tril(A[I,I],-1) + eye(length(I)); -A[I, J] = L \ A[I, J]; +L = tril(A[I,I],-1) + eye(length(I)) +A[I, J] = L \ A[I, J] # Trailing submatrix update if !isempty(K) - A[K,J] = A[K,J] - A[K,I]*A[I,J]; + A[K,J] = A[K,J] - A[K,I]*A[I,J] end -done = true; +done = true return done