Commit 57897def authored by Nick R. Papior's avatar Nick R. Papior
Browse files

enh: cleanup in TS, phases, ts-dq

This brings a much needed clean-up of comments and
array indices in the TS approach.
Basically many of the comments were wrong based on
a transposed view of the matrix + phases.

Additionally we had the dm_update routine doing
the transposing of the matrix calculations.
This was very awkward. This is only for the EQ
with

 DM_eq \propto (Gf - Gf^\dagger)

Now we do it a little easier and do it in the
respective add_DM routines, this makes
dm_update much easier to comprehend and it
should also be faster since searching the sparse
elements should be slower than indexing the transpose
array element in add_DM.
This means we have a problem for ts_dq in the MUMPS@k
solver. The difference in minor but may result in
slower dEf convergence.

Cleaned up all phase calculations, now they are
all based on dot_product(k, sc_off) which is easier
to comprehend.

Cleaned up MUMPS scattering matrix calculation. It is
easier to understand now.

For the self-energy calculations we now rely on Siesta
to provide a Hermitian matrix. I.e. we do not calculate
H = (H + H^A) / 2. Generally Siesta behaves very well and
we should trust this. This may introduce slight noises
in the calculations. We could, if large noises occur make
this feature optional. Secondly, the transfer matrix
is not symmetrized, so one could argue that this needs to
be done in any case.
SPEEDUP.

For Hk we also remove this symmetrization and shift the
Fermi-level in one go. This removes a complete loop of
the sparse matrix elements in Hk.
SPEEDUP.

ts_dq calculations now have more is_nan checks and also
prints out the charge accummulated at Ef.
The interpolation from the file also tries to decide
whether the change in Ef co-incides with the sign of dq.
If not the file interpolation will be turned off for
that iteration. Generally a tight dq tolerance requires
a smaller TS.dQ.Factor.

All phase changes have been applied to the tbtrans code.
This also encouraged a small clean-up in the orb_current
calculation, merged some if-statements.

I ran all BTD, LAPACK and MUMPS calculations
for gamma and k-points. All give the same results.
parent b01d92f4
......@@ -50,6 +50,7 @@ contains
H, S, DM, EDM, Ef, &
Qtot)
use, intrinsic :: ieee_arithmetic, only: ieee_is_nan
use units, only : eV
use alloc, only : re_alloc, de_alloc
......@@ -335,8 +336,7 @@ contains
if ( i_F < 2 ) then
call ts_dq%calculate_dEf(Qtot, &
sp_dist, nspin, n_nzs, DM, S, dEf)
call ts_dq%calculate_dEf(Qtot, sp_dist, nspin, n_nzs, DM, S, dEf)
! In case the charge *is* converged
! Then we further reduce dEf since
! it causes SCF fluctuations
......@@ -355,18 +355,28 @@ contains
! In case we have accumulated 2 or more points
dEf = Ef
call interp_spline(i_F,Q_Ef(1:i_F,1),Q_Ef(1:i_F,2),Qtot,Ef)
dEf = Ef - dEf
if ( ieee_is_nan(Ef) ) then
! There was an error in the interpolation
! So we'll stick with the current one and continue
! It probably means we are close anyway.
Ef = dEf
converged = .true.
! Truncate to the maximum allowed change in Fermi-level
call ts_dq_truncate(Q_Ef(i_F,2), &
TS_DQ_FERMI_MAX, Ef, truncated=truncated)
converged = converged .and. (.not. truncated)
else
dEf = Ef - dEf
if ( IONode ) then
write(*,'(a,es11.4)') 'ts-dq-iscf: cubic spline dq = ', &
Q_Ef(i_F,1) - qtot
! Truncate to the maximum allowed change in Fermi-level
call ts_dq_truncate(Q_Ef(i_F,2), &
TS_DQ_FERMI_MAX, Ef, truncated=truncated)
converged = converged .and. (.not. truncated)
if ( IONode ) then
write(*,'(a,es11.4)') 'ts-dq-iscf: cubic spline dq = ', &
Q_Ef(i_F,1) - qtot
end if
end if
end if
! Signal only doing this for the first run
......@@ -441,11 +451,16 @@ contains
!
! Guess-stimate the actual Fermi-shift
! typically will the above be "too" little
! the above will typically be "too" little.
! So we interpolate between all previous
! estimations for this geometry...
call ts_dq_Fermi_file(Ef, dEf)
Ef = Ef + dEf
call ts_dq_Fermi_file(Q_Ef(1,2), dEf, dq=Q_Ef(1,1) - Qtot)
! If the fermi-file extrapolation does not work
! we get back a dEf == 0.
! In this case we rely on the simpler dQ/dQ@Ef
if ( abs(dEf) > 1.e-11_dp ) then
Ef = Q_Ef(1,2) + dEf
end if
end if
......
This diff is collapsed.
......@@ -1537,8 +1537,13 @@ contains
! The algorithm outside should take care of the
! nullification of the k-point in the semi-infinite direction
do i = 0 , n_s - 1
ph(i) = exp(cmplx(0._dp, sum(kq*sc_off(:,i)),kind=dp) )
! While the Siesta convention matches xijg = R(j) - R(i)
! and this means ij -> (i,j) we cannot simply take the negative phase here.
! This is because the transfer matrix (H01) is not symmetric
! and therefore we *must* setup the Hamiltonian in the
! correct order (not relying on symmetries)
do is = 0 , n_s - 1
ph(is) = exp(cmplx(0._dp, dot_product(kq,sc_off(:,is))))
end do
! Initialize arrays
......@@ -1547,61 +1552,31 @@ contains
Hk_T(:,:) = z_0
Sk_T(:,:) = z_0
!$OMP parallel default(shared), private(i,j,io,jo,ind,is)
!$OMP parallel default(shared), private(io,jo,ind,is)
! We will not have any data-race condition here
!$OMP do
do io = 1 , no
! Create 00
do j = 1 , ncol00(io)
ind = l_ptr00(io) + j
do ind = l_ptr00(io) + 1 , l_ptr00(io) + ncol00(io)
jo = ucorb(l_col00(ind),no)
is = (l_col00(ind)-1) / no
Hk(io,jo) = Hk(io,jo) + H00(ind,ispin) * ph(is)
Sk(io,jo) = Sk(io,jo) + S00(ind) * ph(is)
enddo
Hk(io,jo) = Hk(io,jo) + (H00(ind,ispin) - Ef * S00(ind)) * ph(is)
Sk(io,jo) = Sk(io,jo) + S00(ind) * ph(is)
end do
! Create 01
do j = 1 , ncol01(io)
ind = l_ptr01(io) + j
do ind = l_ptr01(io) + 1 , l_ptr01(io) + ncol01(io)
jo = ucorb(l_col01(ind),no)
is = (l_col01(ind)-1) / no
Hk_T(io,jo) = Hk_T(io,jo) + H01(ind,ispin) * ph(is)
Sk_T(io,jo) = Sk_T(io,jo) + S01(ind) * ph(is)
Hk_T(io,jo) = Hk_T(io,jo) + (H01(ind,ispin) - Ef * S01(ind)) * ph(is)
Sk_T(io,jo) = Sk_T(io,jo) + S01(ind) * ph(is)
end do
end do
!$OMP end do
! Symmetrize 00 and make EF the energy-zero
! We will not have any data-race condition here
!$OMP do
do io = 1 , no
do jo = 1 , io - 1
Sk(io,jo) = 0.5_dp*( Sk(io,jo) + conjg(Sk(jo,io)) )
Sk(jo,io) = conjg(Sk(io,jo))
Hk(io,jo) = 0.5_dp*( Hk(io,jo) + conjg(Hk(jo,io)) ) - &
Ef * Sk(io,jo)
Hk(jo,io) = conjg(Hk(io,jo))
! Transfer matrix is not symmetric, so do not symmetrize
Hk_T(jo,io) = Hk_T(jo,io) - Ef * Sk_T(jo,io)
Hk_T(io,jo) = Hk_T(io,jo) - Ef * Sk_T(io,jo)
end do
Sk(io,io) = real(Sk(io,io),dp)
Hk(io,io) = real(Hk(io,io),dp) - Ef * Sk(io,io)
! Transfer matrix
Hk_T(io,io) = Hk_T(io,io) - Ef * Sk_T(io,io)
end do
!$OMP end do nowait
......
......@@ -1786,7 +1786,7 @@ contains
if ( this%inf_dir == INF_NEGATIVE ) then
tm(this%t_dir) = -1
else if ( this%inf_dir == INF_POSITIVE ) then
tm(this%t_dir) = 1
tm(this%t_dir) = 1
else
call die('Electrode direction not recognized')
end if
......
......@@ -280,8 +280,7 @@ contains
tsup_sp_uc, Calc_Forces)
! Include spin factor and 1/\pi
kw = 1._dp / Pi
if ( nspin == 1 ) kw = kw * 2._dp
kw = 2._dp / (Pi * nspin)
#ifdef TRANSIESTA_TIMING
call timer('TS_HS',1)
......@@ -629,10 +628,10 @@ contains
! Arrays needed for looping the sparsity
type(Sparsity), pointer :: s
integer, pointer :: l_ncol(:), l_ptr(:), l_col(:)
integer, pointer :: s_ncol(:), s_ptr(:), s_col(:)
integer, pointer :: s_ncol(:), s_ptr(:), s_col(:), sp_col(:)
real(dp), pointer :: D(:,:), E(:,:), Sg(:)
integer :: io, ind, nr
integer :: s_ptr_begin, s_ptr_end, sin
integer :: sp, sind
integer :: iu, ju, i1, i2
logical :: lis_eq, hasEDM, calc_q
......@@ -665,28 +664,29 @@ contains
if ( calc_q .and. hasEDM ) then
!$OMP parallel do default(shared), &
!$OMP&private(io,iu,ind,ju,s_ptr_begin,s_ptr_end,sin)
!$OMP& private(io,iu,ind,ju,sp,sp_col,sind), &
!$OMP& reduction(-:q)
do io = 1 , nr
! Quickly go past the buffer atoms...
if ( l_ncol(io) /= 0 ) then
s_ptr_begin = s_ptr(io) + 1
s_ptr_end = s_ptr(io) + s_ncol(io)
sp = s_ptr(io)
sp_col => s_col(sp+1:sp+s_ncol(io))
! The update region equivalent GF part
iu = io - orb_offset(io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ju = l_col(ind) - orb_offset(l_col(ind)) - &
offset(N_Elec,Elecs,l_col(ind))
ju = l_col(ind) - orb_offset(l_col(ind)) - offset(N_Elec,Elecs,l_col(ind))
! Search for overlap index
! spS is transposed, so we have to conjugate the
! S value, then we may take the imaginary part.
sin = s_ptr_begin - 1 + SFIND(s_col(s_ptr_begin:s_ptr_end), l_col(ind))
sind = sp + SFIND(sp_col, l_col(ind))
if ( sp < sind ) then
! Gf = -Gf^\dagger @ \Gamma, no need to worry about Gf.S
q = q - aimag(Gf(iu,ju) * Sg(sind))
end if
if ( sin >= s_ptr_begin ) q = q - aimag(GF(iu,ju) * Sg(sin))
D(ind,i1) = D(ind,i1) - aimag( GF(iu,ju) * DMfact )
E(ind,i2) = E(ind,i2) - aimag( GF(iu,ju) * EDMfact )
......@@ -703,8 +703,7 @@ contains
if ( l_ncol(io) /= 0 ) then
iu = io - orb_offset(io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ju = l_col(ind) - orb_offset(l_col(ind)) - &
offset(N_Elec,Elecs,l_col(ind))
ju = l_col(ind) - orb_offset(l_col(ind)) - offset(N_Elec,Elecs,l_col(ind))
D(ind,i1) = D(ind,i1) - aimag( GF(iu,ju) * DMfact )
E(ind,i2) = E(ind,i2) - aimag( GF(iu,ju) * EDMfact )
end do
......@@ -715,17 +714,17 @@ contains
else if ( calc_q ) then
!$OMP parallel do default(shared), &
!$OMP&private(io,iu,ind,ju,s_ptr_begin,s_ptr_end,sin)
!$OMP& private(io,iu,sp,sp_col,ind,ju,sind), &
!$OMP& reduction(-:q)
do io = 1 , nr
if ( l_ncol(io) /= 0 ) then
s_ptr_begin = s_ptr(io) + 1
s_ptr_end = s_ptr(io) + s_ncol(io)
sp = s_ptr(io)
sp_col => s_col(sp+1:sp+s_ncol(io))
iu = io - orb_offset(io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ju = l_col(ind) - orb_offset(l_col(ind)) - &
offset(N_Elec,Elecs,l_col(ind))
sin = s_ptr_begin - 1 + SFIND(s_col(s_ptr_begin:s_ptr_end), l_col(ind))
if ( sin >= s_ptr_begin ) q = q - aimag(GF(iu,ju) * Sg(sin))
ju = l_col(ind) - orb_offset(l_col(ind)) - offset(N_Elec,Elecs,l_col(ind))
sind = sp + SFIND(sp_col, l_col(ind))
if ( sp < sind ) q = q - aimag(GF(iu,ju) * Sg(sind))
D(ind,i1) = D(ind,i1) - aimag( GF(iu,ju) * DMfact )
end do
end if
......@@ -740,8 +739,7 @@ contains
if ( l_ncol(io) /= 0 ) then
iu = io - orb_offset(io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ju = l_col(ind) - orb_offset(l_col(ind)) - &
offset(N_Elec,Elecs,l_col(ind))
ju = l_col(ind) - orb_offset(l_col(ind)) - offset(N_Elec,Elecs,l_col(ind))
D(ind,i1) = D(ind,i1) - aimag( GF(iu,ju) * DMfact )
end do
end if
......@@ -785,9 +783,6 @@ contains
end if
end if
! For ts_dq we should not multiply by 2 since we don't do G + G^\dagger for Gamma-only
! this is because G is ensured symmetric for Gamma-point and thus it is not needed.
contains
pure function offset(N_Elec,Elecs,io)
......@@ -818,7 +813,7 @@ contains
! Remark that we need the left buffer orbitals
! to calculate the actual orbital of the sparse matrices...
integer, intent(in) :: no_u
complex(dp), intent(out) :: GFinv(no_u**2)
complex(dp), intent(out) :: GFinv(no_u,no_u)
integer, intent(in) :: N_Elec
type(Elec), intent(in) :: Elecs(N_Elec)
! The Hamiltonian and overlap sparse matrices
......@@ -829,7 +824,7 @@ contains
type(Sparsity), pointer :: sp
integer, pointer :: l_ncol(:), l_ptr(:), l_col(:)
real(dp), pointer :: H(:), S(:)
integer :: io, iu, ind, ioff, nr
integer :: io, iu, ind, nr
if ( cE%fake ) return
......@@ -847,9 +842,9 @@ contains
nrows_g=nr)
! Initialize
GFinv(:) = cmplx(0._dp,0._dp, dp)
GFinv(:,:) = cmplx(0._dp,0._dp, dp)
!$OMP parallel default(shared), private(io,ioff,iu,ind)
!$OMP parallel default(shared), private(io,iu,ind)
! We will only loop in the central region
! We have constructed the sparse array to only contain
......@@ -859,16 +854,11 @@ contains
if ( l_ncol(io) /= 0 ) then
ioff = orb_offset(io)
iu = (io - ioff - 1) * no_u
iu = io - orb_offset(io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ioff = orb_offset(l_col(ind))
! Notice that we transpose S and H back here
! See symmetrize_HS_Gamma (H is hermitian)
GFinv(iu+l_col(ind)-ioff) = Z * S(ind) - H(ind)
GFinv(l_col(ind)-orb_offset(l_col(ind)),iu) = Z * S(ind) - H(ind)
end do
......
......@@ -698,11 +698,11 @@ contains
! Arrays needed for looping the sparsity
type(Sparsity), pointer :: s
integer, pointer :: l_ncol(:), l_ptr(:), l_col(:)
integer, pointer :: s_ncol(:), s_ptr(:), s_col(:)
integer, pointer :: s_ncol(:), s_ptr(:), s_col(:), sp_col(:)
complex(dp), pointer :: D(:,:), E(:,:), Sk(:)
integer :: io, ind, nr
integer :: s_ptr_begin, s_ptr_end, sin
integer :: iu, ju, i1, i2
integer :: sp, sind
integer :: iu, iuT, ju, juT, i1, i2
logical :: lis_eq, hasEDM, calc_q
lis_eq = .true.
......@@ -733,30 +733,35 @@ contains
if ( calc_q .and. hasEDM ) then
!$OMP parallel do default(shared), private(io,iu,ind,ju,s_ptr_begin,s_ptr_end,sin)
!$OMP parallel do default(shared), &
!$OMP& private(io,sp,sp_col,iu,iuT,ind,ju,juT,sind), &
!$OMP& reduction(-:q)
do io = 1 , nr
! Quickly go past the buffer atoms...
if ( l_ncol(io) /= 0 ) then
s_ptr_begin = s_ptr(io) + 1
s_ptr_end = s_ptr(io) + s_ncol(io)
sp = s_ptr(io)
sp_col => s_col(sp+1:sp+s_ncol(io))
! The update region equivalent GF part
iu = io - orb_offset(io)
iuT = iu - offset(N_Elec,Elecs,io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ju = l_col(ind) - orb_offset(l_col(ind)) - &
offset(N_Elec,Elecs,l_col(ind))
! Search for overlap index
! spS is transposed, so we have to conjugate the
! S value, then we may take the imaginary part.
sin = s_ptr_begin - 1 + SFIND(s_col(s_ptr_begin:s_ptr_end), l_col(ind))
if ( sin >= s_ptr_begin ) q = q - aimag(GF(iu,ju) * Sk(sin))
D(ind,i1) = D(ind,i1) - Gf(iu,ju) * DMfact
E(ind,i2) = E(ind,i2) - Gf(iu,ju) * EDMfact
juT = l_col(ind) - orb_offset(l_col(ind))
ju = juT - offset(N_Elec,Elecs,l_col(ind))
! Search for overlap index, note it is transposed (- phase)
! So this is S(ju,iu)
sind = sp + SFIND(sp_col, l_col(ind))
if ( sp < sind ) then
! For charge calculation it is Tr[(G-G^\dagger).S] i.e. the matrix
! multiplication
q = q - aimag((Gf(iu,ju) - conjg(Gf(juT,iuT))) * Sk(sind))
end if
D(ind,i1) = D(ind,i1) + Gf(iu,ju) * DMfact - conjg(Gf(juT,iuT) * DMfact)
E(ind,i2) = E(ind,i2) + Gf(iu,ju) * EDMfact - conjg(Gf(juT,iuT) * EDMfact)
end do
......@@ -766,15 +771,16 @@ contains
else if ( hasEDM ) then
!$OMP parallel do default(shared), private(io,iu,ind,ju)
!$OMP parallel do default(shared), private(io,iu,iuT,ind,ju,juT)
do io = 1 , nr
if ( l_ncol(io) /= 0 ) then
iu = io - orb_offset(io)
iuT = iu - offset(N_Elec,Elecs,io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ju = l_col(ind) - orb_offset(l_col(ind)) - &
offset(N_Elec,Elecs,l_col(ind))
D(ind,i1) = D(ind,i1) - Gf(iu,ju) * DMfact
E(ind,i2) = E(ind,i2) - Gf(iu,ju) * EDMfact
juT = l_col(ind) - orb_offset(l_col(ind))
ju = juT - offset(N_Elec,Elecs,l_col(ind))
D(ind,i1) = D(ind,i1) + Gf(iu,ju) * DMfact - conjg(Gf(juT,iuT) * DMfact)
E(ind,i2) = E(ind,i2) + Gf(iu,ju) * EDMfact - conjg(Gf(juT,iuT) * EDMfact)
end do
end if
end do
......@@ -782,33 +788,39 @@ contains
else if ( calc_q ) then
!$OMP parallel do default(shared), private(io,iu,ind,ju,s_ptr_begin,s_ptr_end,sin)
!$OMP parallel do default(shared), &
!$OMP& private(io,sp,sp_col,iu,iuT,ind,ju,juT,sind), &
!$OMP& reduction(-:q)
do io = 1 , nr
if ( l_ncol(io) /= 0 ) then
s_ptr_begin = s_ptr(io) + 1
s_ptr_end = s_ptr(io) + s_ncol(io)
sp = s_ptr(io)
sp_col => s_col(sp+1:sp+s_ncol(io))
iu = io - orb_offset(io)
iuT = iu - offset(N_Elec,Elecs,io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ju = l_col(ind) - orb_offset(l_col(ind)) - &
offset(N_Elec,Elecs,l_col(ind))
sin = s_ptr_begin - 1 + SFIND(s_col(s_ptr_begin:s_ptr_end), l_col(ind))
if ( sin >= s_ptr_begin ) q = q - aimag(GF(iu,ju) * Sk(sin))
D(ind,i1) = D(ind,i1) - Gf(iu,ju) * DMfact
juT = l_col(ind) - orb_offset(l_col(ind))
ju = juT - offset(N_Elec,Elecs,l_col(ind))
sind = sp + SFIND(sp_col, l_col(ind))
if ( sp < sind ) then
q = q - aimag((Gf(iu,ju) - conjg(Gf(juT,iuT))) * Sk(sind))
end if
D(ind,i1) = D(ind,i1) + Gf(iu,ju) * DMfact - conjg(Gf(juT,iuT) * DMfact)
end do
end if
end do
!$OMP end parallel do
else
!$OMP parallel do default(shared), private(io,iu,ind,ju)
!$OMP parallel do default(shared), private(io,iu,iuT,ind,ju,juT)
do io = 1 , nr
if ( l_ncol(io) /= 0 ) then
iu = io - orb_offset(io)
iuT = iu - offset(N_Elec,Elecs,io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ju = l_col(ind) - orb_offset(l_col(ind)) - &
offset(N_Elec,Elecs,l_col(ind))
D(ind,i1) = D(ind,i1) - Gf(iu,ju) * DMfact
juT = l_col(ind) - orb_offset(l_col(ind))
ju = juT - offset(N_Elec,Elecs,l_col(ind))
D(ind,i1) = D(ind,i1) + Gf(iu,ju) * DMfact - conjg(Gf(juT,iuT) * DMfact)
end do
end if
end do
......@@ -818,6 +830,10 @@ contains
else ! lis_eq
#ifndef TS_NOCHECKS
if ( no1 /= no2 ) call die("Error in matrix dimensions")
#endif
if ( hasEDM ) then
!$OMP parallel do default(shared), private(io,iu,ind,ju)
......@@ -850,8 +866,6 @@ contains
end if
end if
if ( calc_q ) q = q * 2._dp
contains
function offset(N_Elec,Elecs,io)
......@@ -891,7 +905,7 @@ contains
type(Sparsity), pointer :: sp
integer, pointer :: l_ncol(:), l_ptr(:), l_col(:)
complex(dp), pointer :: H(:), S(:)
integer :: io, iu, ind, ioff, nr
integer :: io, iu, ind, nr
if ( cE%fake ) return
......@@ -911,7 +925,7 @@ contains
! Initialize
GFinv(:,:) = cmplx(0._dp,0._dp,dp)
!$OMP parallel default(shared), private(io,ioff,iu,ind)
!$OMP parallel default(shared), private(io,iu,ind)
! We will only loop in the central region
! We have constructed the sparse array to only contain
......@@ -925,11 +939,8 @@ contains
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
ioff = orb_offset(l_col(ind))
! Notice that we transpose S and H back here
! See symmetrize_HS_Gamma (H is hermitian)
GFinv(l_col(ind)-ioff,iu) = Z * S(ind) - H(ind)
! Transpose to match phase convention in ts_sparse_helper (- phase)
GFinv(l_col(ind)-orb_offset(l_col(ind)),iu) = Z * S(ind) - H(ind)
end do
end if
......@@ -937,7 +948,7 @@ contains
!$OMP end do
do io = 1 , N_Elec
call insert_Self_Energies(no_u, Gfinv, Elecs(io))
call insert_Self_Energies(no_u, Gfinv, Elecs(io))
end do
!$OMP end parallel
......
......@@ -220,9 +220,9 @@ contains
ioff = io - orb_offset(io)
do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io)
mum%IRN(ind) = l_col(ind) - orb_offset(l_col(ind))
mum%JCN(ind) = ioff
mum%IRN(ind) = ioff
mum%JCN(ind) = l_col(ind) - orb_offset(l_col(ind))
end do
......@@ -316,16 +316,14 @@ contains
! TODO, this requires that the sparsity pattern is symmetric
! Which it always is!
mum%IRHS_PTR(:) = 1
!$OMP parallel do default(shared), private(io,j,ind)
!$OMP parallel do default(shared), private(io,ind)
do io = 1 , nr