diff --git a/devel/svdwf/random_svd.irp.f b/devel/svdwf/random_svd.irp.f
index 9146372..9ac9119 100644
--- a/devel/svdwf/random_svd.irp.f
+++ b/devel/svdwf/random_svd.irp.f
@@ -13,8 +13,7 @@ subroutine run
   include 'constants.include.F'
   double precision, allocatable :: U(:,:), V(:,:), D(:), A(:,:)
   integer :: i, j, k, l, q, r, m, n, iter
-  double precision,allocatable    :: Z(:,:), P(:,:), Yt(:,:), UYt(:,:)
-  double precision :: r1,r2
+  double precision,allocatable    :: Z(:,:), P(:,:), Yt(:,:), UYt(:,:), r1(:,:)
 
   m = n_det_alpha_unique
   n = n_det_beta_unique
@@ -25,52 +24,75 @@ subroutine run
 
   ! Z(m,r) = A(m,n).P(n,r)                                                                    
   Z(:,:) = 0.d0
+  allocate(P(n,r))
+
+  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,r1)
+  allocate(r1(N_det,2))
+  !$OMP DO 
   do l=1,r
+   call random_number(r1)
+   r1(:,1) = dsqrt(-2.d0*dlog(r1(:,1)))
+   r1(:,1) = r1(:,1) * dcos(dtwo_pi*r1(:,2))
    do k=1,N_det
      i = psi_bilinear_matrix_rows(k)
      j = psi_bilinear_matrix_columns(k)
-     call random_number(r1)
-     call random_number(r2)
-     r1 = dsqrt(-2.d0*dlog(r1))
-     r2 = dtwo_pi*r2
-     Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1*dcos(r2)
+     Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1(k,1) 
    enddo
   enddo
+  !$OMP END DO
+  deallocate(r1)
+  !$OMP END PARALLEL
 
   ! Power iterations
-  allocate(P(n,r))
   do iter=1,20
+   print *, 'Power iteration ', iter, '/', 20
+
    ! P(n,r) = At(n,m).Z(m,r)
-   P(:,:) = 0.d0
+   !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
+   !$OMP DO
    do l=1,r
+    P(:,l) = 0.d0
     do k=1,N_det
       i = psi_bilinear_matrix_rows(k)
       j = psi_bilinear_matrix_columns(k)
       P(j,l) = P(j,l) + psi_bilinear_matrix_values(k,1) * Z(i,l)
     enddo
    enddo
-   Z(:,:) = 0.d0
+   !$OMP END DO
+
+   !$OMP BARRIER
+
+   !$OMP DO
    do l=1,r
+    Z(:,l) = 0.d0
     do k=1,N_det
       i = psi_bilinear_matrix_rows(k)
       j = psi_bilinear_matrix_columns(k)
       Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * P(j,l)
     enddo
    enddo
+   !$OMP END DO
+
+   !$OMP END PARALLEL
+
    ! Compute QR
    call ortho_qr(Z,size(Z,1),m,r)
   enddo
 
   ! Y(r,n) = Zt(r,m).A(m,n)
   allocate(Yt(n,r))
-  Yt(:,:) = 0.d0
+  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
+  !$OMP DO
   do l=1,r
+    Yt(:,l) = 0.d0
     do k=1,N_det
      i = psi_bilinear_matrix_rows(k)
      j = psi_bilinear_matrix_columns(k)
      Yt(j,l) = Yt(j,l) + Z(i,l) * psi_bilinear_matrix_values(k,1)
    enddo
   enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
 
   allocate(D(r),V(n,r), UYt(r,r))
   call svd(Yt,size(Yt,1),V,size(V,1),D,UYt,size(UYt,1),n,r)