Improvemnts to MAT ASB in OpenMP

omp-threadsafe
Salvatore Filippone 2 years ago
parent 776c755112
commit 4d988ea3db

@ -131,7 +131,7 @@ Contains
! ...Local Variables
complex(psb_spk_),allocatable :: tmp(:)
integer(psb_mpk_) :: dim, lb_, lbi,ub_
integer(psb_mpk_) :: dim, lb_, lbi,ub_, i
integer(psb_ipk_) :: err_act,err
character(len=30) :: name
logical, parameter :: debug=.false.
@ -179,7 +179,10 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb_-1+dim+1:lb_-1+len) = pad
!$omp parallel do private(i) shared(dim,len)
do i=lb_-1+dim+1,lb_-1+len
rrax(i) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -204,7 +207,7 @@ Contains
complex(psb_spk_),allocatable :: tmp(:,:)
integer(psb_ipk_) :: err_act,err
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2, i
character(len=30) :: name
name='psb_r_m_c_rk2'
@ -267,8 +270,14 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb1_-1+dim+1:lb1_-1+len1,:) = pad
rrax(lb1_:lb1_-1+dim,lb2_-1+dim2+1:lb2_-1+len2) = pad
!$omp parallel do private(i) shared(lb1_,dim,len1)
do i=lb1_-1+dim+1,lb1_-1+len1
rrax(i,:) = pad
end do
!$omp parallel do private(i) shared(lb1_,dim,len1,lb2_,dim2,len2)
do i=lb1_,lb1_-1+len1
rrax(i,lb2_-1+dim2+1:lb2_-1+len2) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -982,48 +991,7 @@ Contains
goto 9999
end if
!!$ If (len > psb_size(v)) Then
!!$ if (present(newsz)) then
!!$ isz = (max(len+1,newsz))
!!$ else
!!$ if (present(addsz)) then
!!$ isz = len+max(1,addsz)
!!$ else
!!$ isz = max(len+10, int(1.25*len))
!!$ endif
!!$ endif
!!$
!!$ call psb_realloc(isz,v,info,pad=pad)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ call psb_errpush(info,name,a_err='psb_realloc')
!!$ goto 9999
!!$ End If
!!$ end If
If (len > psb_size(v)) Then
#if defined(OPENMP)
!$OMP CRITICAL
if (len > psb_size(v)) then
if (present(newsz)) then
isz = (max(len+1,newsz))
else
if (present(addsz)) then
isz = len+max(1,addsz)
else
isz = max(len+10, int(1.25*len))
endif
endif
call psb_realloc(isz,v,info,pad=pad)
end if
!$OMP END CRITICAL
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
end if
#else
if (present(newsz)) then
isz = (max(len+1,newsz))
else
@ -1040,7 +1008,6 @@ Contains
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
End If
#endif
end If
call psb_erractionrestore(err_act)

@ -131,7 +131,7 @@ Contains
! ...Local Variables
real(psb_dpk_),allocatable :: tmp(:)
integer(psb_mpk_) :: dim, lb_, lbi,ub_
integer(psb_mpk_) :: dim, lb_, lbi,ub_, i
integer(psb_ipk_) :: err_act,err
character(len=30) :: name
logical, parameter :: debug=.false.
@ -179,7 +179,10 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb_-1+dim+1:lb_-1+len) = pad
!$omp parallel do private(i) shared(dim,len)
do i=lb_-1+dim+1,lb_-1+len
rrax(i) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -204,7 +207,7 @@ Contains
real(psb_dpk_),allocatable :: tmp(:,:)
integer(psb_ipk_) :: err_act,err
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2, i
character(len=30) :: name
name='psb_r_m_d_rk2'
@ -267,8 +270,14 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb1_-1+dim+1:lb1_-1+len1,:) = pad
rrax(lb1_:lb1_-1+dim,lb2_-1+dim2+1:lb2_-1+len2) = pad
!$omp parallel do private(i) shared(lb1_,dim,len1)
do i=lb1_-1+dim+1,lb1_-1+len1
rrax(i,:) = pad
end do
!$omp parallel do private(i) shared(lb1_,dim,len1,lb2_,dim2,len2)
do i=lb1_,lb1_-1+len1
rrax(i,lb2_-1+dim2+1:lb2_-1+len2) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -982,48 +991,7 @@ Contains
goto 9999
end if
!!$ If (len > psb_size(v)) Then
!!$ if (present(newsz)) then
!!$ isz = (max(len+1,newsz))
!!$ else
!!$ if (present(addsz)) then
!!$ isz = len+max(1,addsz)
!!$ else
!!$ isz = max(len+10, int(1.25*len))
!!$ endif
!!$ endif
!!$
!!$ call psb_realloc(isz,v,info,pad=pad)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ call psb_errpush(info,name,a_err='psb_realloc')
!!$ goto 9999
!!$ End If
!!$ end If
If (len > psb_size(v)) Then
#if defined(OPENMP)
!$OMP CRITICAL
if (len > psb_size(v)) then
if (present(newsz)) then
isz = (max(len+1,newsz))
else
if (present(addsz)) then
isz = len+max(1,addsz)
else
isz = max(len+10, int(1.25*len))
endif
endif
call psb_realloc(isz,v,info,pad=pad)
end if
!$OMP END CRITICAL
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
end if
#else
if (present(newsz)) then
isz = (max(len+1,newsz))
else
@ -1040,7 +1008,6 @@ Contains
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
End If
#endif
end If
call psb_erractionrestore(err_act)

@ -131,7 +131,7 @@ Contains
! ...Local Variables
integer(psb_epk_),allocatable :: tmp(:)
integer(psb_mpk_) :: dim, lb_, lbi,ub_
integer(psb_mpk_) :: dim, lb_, lbi,ub_, i
integer(psb_ipk_) :: err_act,err
character(len=30) :: name
logical, parameter :: debug=.false.
@ -179,7 +179,10 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb_-1+dim+1:lb_-1+len) = pad
!$omp parallel do private(i) shared(dim,len)
do i=lb_-1+dim+1,lb_-1+len
rrax(i) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -204,7 +207,7 @@ Contains
integer(psb_epk_),allocatable :: tmp(:,:)
integer(psb_ipk_) :: err_act,err
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2, i
character(len=30) :: name
name='psb_r_m_e_rk2'
@ -267,8 +270,14 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb1_-1+dim+1:lb1_-1+len1,:) = pad
rrax(lb1_:lb1_-1+dim,lb2_-1+dim2+1:lb2_-1+len2) = pad
!$omp parallel do private(i) shared(lb1_,dim,len1)
do i=lb1_-1+dim+1,lb1_-1+len1
rrax(i,:) = pad
end do
!$omp parallel do private(i) shared(lb1_,dim,len1,lb2_,dim2,len2)
do i=lb1_,lb1_-1+len1
rrax(i,lb2_-1+dim2+1:lb2_-1+len2) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -982,48 +991,7 @@ Contains
goto 9999
end if
!!$ If (len > psb_size(v)) Then
!!$ if (present(newsz)) then
!!$ isz = (max(len+1,newsz))
!!$ else
!!$ if (present(addsz)) then
!!$ isz = len+max(1,addsz)
!!$ else
!!$ isz = max(len+10, int(1.25*len))
!!$ endif
!!$ endif
!!$
!!$ call psb_realloc(isz,v,info,pad=pad)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ call psb_errpush(info,name,a_err='psb_realloc')
!!$ goto 9999
!!$ End If
!!$ end If
If (len > psb_size(v)) Then
#if defined(OPENMP)
!$OMP CRITICAL
if (len > psb_size(v)) then
if (present(newsz)) then
isz = (max(len+1,newsz))
else
if (present(addsz)) then
isz = len+max(1,addsz)
else
isz = max(len+10, int(1.25*len))
endif
endif
call psb_realloc(isz,v,info,pad=pad)
end if
!$OMP END CRITICAL
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
end if
#else
if (present(newsz)) then
isz = (max(len+1,newsz))
else
@ -1040,7 +1008,6 @@ Contains
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
End If
#endif
end If
call psb_erractionrestore(err_act)

@ -131,7 +131,7 @@ Contains
! ...Local Variables
integer(psb_i2pk_),allocatable :: tmp(:)
integer(psb_mpk_) :: dim, lb_, lbi,ub_
integer(psb_mpk_) :: dim, lb_, lbi,ub_, i
integer(psb_ipk_) :: err_act,err
character(len=30) :: name
logical, parameter :: debug=.false.
@ -179,7 +179,10 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb_-1+dim+1:lb_-1+len) = pad
!$omp parallel do private(i) shared(dim,len)
do i=lb_-1+dim+1,lb_-1+len
rrax(i) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -204,7 +207,7 @@ Contains
integer(psb_i2pk_),allocatable :: tmp(:,:)
integer(psb_ipk_) :: err_act,err
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2, i
character(len=30) :: name
name='psb_r_m_i2_rk2'
@ -267,8 +270,14 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb1_-1+dim+1:lb1_-1+len1,:) = pad
rrax(lb1_:lb1_-1+dim,lb2_-1+dim2+1:lb2_-1+len2) = pad
!$omp parallel do private(i) shared(lb1_,dim,len1)
do i=lb1_-1+dim+1,lb1_-1+len1
rrax(i,:) = pad
end do
!$omp parallel do private(i) shared(lb1_,dim,len1,lb2_,dim2,len2)
do i=lb1_,lb1_-1+len1
rrax(i,lb2_-1+dim2+1:lb2_-1+len2) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -982,48 +991,7 @@ Contains
goto 9999
end if
!!$ If (len > psb_size(v)) Then
!!$ if (present(newsz)) then
!!$ isz = (max(len+1,newsz))
!!$ else
!!$ if (present(addsz)) then
!!$ isz = len+max(1,addsz)
!!$ else
!!$ isz = max(len+10, int(1.25*len))
!!$ endif
!!$ endif
!!$
!!$ call psb_realloc(isz,v,info,pad=pad)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ call psb_errpush(info,name,a_err='psb_realloc')
!!$ goto 9999
!!$ End If
!!$ end If
If (len > psb_size(v)) Then
#if defined(OPENMP)
!$OMP CRITICAL
if (len > psb_size(v)) then
if (present(newsz)) then
isz = (max(len+1,newsz))
else
if (present(addsz)) then
isz = len+max(1,addsz)
else
isz = max(len+10, int(1.25*len))
endif
endif
call psb_realloc(isz,v,info,pad=pad)
end if
!$OMP END CRITICAL
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
end if
#else
if (present(newsz)) then
isz = (max(len+1,newsz))
else
@ -1040,7 +1008,6 @@ Contains
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
End If
#endif
end If
call psb_erractionrestore(err_act)

@ -131,7 +131,7 @@ Contains
! ...Local Variables
integer(psb_mpk_),allocatable :: tmp(:)
integer(psb_mpk_) :: dim, lb_, lbi,ub_
integer(psb_mpk_) :: dim, lb_, lbi,ub_, i
integer(psb_ipk_) :: err_act,err
character(len=30) :: name
logical, parameter :: debug=.false.
@ -179,7 +179,10 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb_-1+dim+1:lb_-1+len) = pad
!$omp parallel do private(i) shared(dim,len)
do i=lb_-1+dim+1,lb_-1+len
rrax(i) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -204,7 +207,7 @@ Contains
integer(psb_mpk_),allocatable :: tmp(:,:)
integer(psb_ipk_) :: err_act,err
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2, i
character(len=30) :: name
name='psb_r_m_m_rk2'
@ -267,8 +270,14 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb1_-1+dim+1:lb1_-1+len1,:) = pad
rrax(lb1_:lb1_-1+dim,lb2_-1+dim2+1:lb2_-1+len2) = pad
!$omp parallel do private(i) shared(lb1_,dim,len1)
do i=lb1_-1+dim+1,lb1_-1+len1
rrax(i,:) = pad
end do
!$omp parallel do private(i) shared(lb1_,dim,len1,lb2_,dim2,len2)
do i=lb1_,lb1_-1+len1
rrax(i,lb2_-1+dim2+1:lb2_-1+len2) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -982,48 +991,7 @@ Contains
goto 9999
end if
!!$ If (len > psb_size(v)) Then
!!$ if (present(newsz)) then
!!$ isz = (max(len+1,newsz))
!!$ else
!!$ if (present(addsz)) then
!!$ isz = len+max(1,addsz)
!!$ else
!!$ isz = max(len+10, int(1.25*len))
!!$ endif
!!$ endif
!!$
!!$ call psb_realloc(isz,v,info,pad=pad)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ call psb_errpush(info,name,a_err='psb_realloc')
!!$ goto 9999
!!$ End If
!!$ end If
If (len > psb_size(v)) Then
#if defined(OPENMP)
!$OMP CRITICAL
if (len > psb_size(v)) then
if (present(newsz)) then
isz = (max(len+1,newsz))
else
if (present(addsz)) then
isz = len+max(1,addsz)
else
isz = max(len+10, int(1.25*len))
endif
endif
call psb_realloc(isz,v,info,pad=pad)
end if
!$OMP END CRITICAL
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
end if
#else
if (present(newsz)) then
isz = (max(len+1,newsz))
else
@ -1040,7 +1008,6 @@ Contains
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
End If
#endif
end If
call psb_erractionrestore(err_act)

@ -131,7 +131,7 @@ Contains
! ...Local Variables
real(psb_spk_),allocatable :: tmp(:)
integer(psb_mpk_) :: dim, lb_, lbi,ub_
integer(psb_mpk_) :: dim, lb_, lbi,ub_, i
integer(psb_ipk_) :: err_act,err
character(len=30) :: name
logical, parameter :: debug=.false.
@ -179,7 +179,10 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb_-1+dim+1:lb_-1+len) = pad
!$omp parallel do private(i) shared(dim,len)
do i=lb_-1+dim+1,lb_-1+len
rrax(i) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -204,7 +207,7 @@ Contains
real(psb_spk_),allocatable :: tmp(:,:)
integer(psb_ipk_) :: err_act,err
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2, i
character(len=30) :: name
name='psb_r_m_s_rk2'
@ -267,8 +270,14 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb1_-1+dim+1:lb1_-1+len1,:) = pad
rrax(lb1_:lb1_-1+dim,lb2_-1+dim2+1:lb2_-1+len2) = pad
!$omp parallel do private(i) shared(lb1_,dim,len1)
do i=lb1_-1+dim+1,lb1_-1+len1
rrax(i,:) = pad
end do
!$omp parallel do private(i) shared(lb1_,dim,len1,lb2_,dim2,len2)
do i=lb1_,lb1_-1+len1
rrax(i,lb2_-1+dim2+1:lb2_-1+len2) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -982,48 +991,7 @@ Contains
goto 9999
end if
!!$ If (len > psb_size(v)) Then
!!$ if (present(newsz)) then
!!$ isz = (max(len+1,newsz))
!!$ else
!!$ if (present(addsz)) then
!!$ isz = len+max(1,addsz)
!!$ else
!!$ isz = max(len+10, int(1.25*len))
!!$ endif
!!$ endif
!!$
!!$ call psb_realloc(isz,v,info,pad=pad)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ call psb_errpush(info,name,a_err='psb_realloc')
!!$ goto 9999
!!$ End If
!!$ end If
If (len > psb_size(v)) Then
#if defined(OPENMP)
!$OMP CRITICAL
if (len > psb_size(v)) then
if (present(newsz)) then
isz = (max(len+1,newsz))
else
if (present(addsz)) then
isz = len+max(1,addsz)
else
isz = max(len+10, int(1.25*len))
endif
endif
call psb_realloc(isz,v,info,pad=pad)
end if
!$OMP END CRITICAL
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
end if
#else
if (present(newsz)) then
isz = (max(len+1,newsz))
else
@ -1040,7 +1008,6 @@ Contains
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
End If
#endif
end If
call psb_erractionrestore(err_act)

@ -131,7 +131,7 @@ Contains
! ...Local Variables
complex(psb_dpk_),allocatable :: tmp(:)
integer(psb_mpk_) :: dim, lb_, lbi,ub_
integer(psb_mpk_) :: dim, lb_, lbi,ub_, i
integer(psb_ipk_) :: err_act,err
character(len=30) :: name
logical, parameter :: debug=.false.
@ -179,7 +179,10 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb_-1+dim+1:lb_-1+len) = pad
!$omp parallel do private(i) shared(dim,len)
do i=lb_-1+dim+1,lb_-1+len
rrax(i) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -204,7 +207,7 @@ Contains
complex(psb_dpk_),allocatable :: tmp(:,:)
integer(psb_ipk_) :: err_act,err
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2
integer(psb_mpk_) :: dim,dim2,lb1_, lb2_, ub1_, ub2_,lbi1, lbi2, i
character(len=30) :: name
name='psb_r_m_z_rk2'
@ -267,8 +270,14 @@ Contains
end if
endif
if (present(pad)) then
rrax(lb1_-1+dim+1:lb1_-1+len1,:) = pad
rrax(lb1_:lb1_-1+dim,lb2_-1+dim2+1:lb2_-1+len2) = pad
!$omp parallel do private(i) shared(lb1_,dim,len1)
do i=lb1_-1+dim+1,lb1_-1+len1
rrax(i,:) = pad
end do
!$omp parallel do private(i) shared(lb1_,dim,len1,lb2_,dim2,len2)
do i=lb1_,lb1_-1+len1
rrax(i,lb2_-1+dim2+1:lb2_-1+len2) = pad
end do
endif
call psb_erractionrestore(err_act)
return
@ -982,48 +991,7 @@ Contains
goto 9999
end if
!!$ If (len > psb_size(v)) Then
!!$ if (present(newsz)) then
!!$ isz = (max(len+1,newsz))
!!$ else
!!$ if (present(addsz)) then
!!$ isz = len+max(1,addsz)
!!$ else
!!$ isz = max(len+10, int(1.25*len))
!!$ endif
!!$ endif
!!$
!!$ call psb_realloc(isz,v,info,pad=pad)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ call psb_errpush(info,name,a_err='psb_realloc')
!!$ goto 9999
!!$ End If
!!$ end If
If (len > psb_size(v)) Then
#if defined(OPENMP)
!$OMP CRITICAL
if (len > psb_size(v)) then
if (present(newsz)) then
isz = (max(len+1,newsz))
else
if (present(addsz)) then
isz = len+max(1,addsz)
else
isz = max(len+10, int(1.25*len))
endif
endif
call psb_realloc(isz,v,info,pad=pad)
end if
!$OMP END CRITICAL
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
end if
#else
if (present(newsz)) then
isz = (max(len+1,newsz))
else
@ -1040,7 +1008,6 @@ Contains
call psb_errpush(info,name,a_err='psb_realloc')
goto 9999
End If
#endif
end If
call psb_erractionrestore(err_act)

@ -1865,26 +1865,29 @@ module psb_c_base_mat_mod
integer(psb_ipk_), intent(in), optional :: idir
end subroutine psb_c_fix_coo_inner
end interface
interface
subroutine psb_c_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,info)
subroutine psb_c_fix_coo_inner_colmajor(nr,nc,nzin,dupl,&
& ia,ja,val,iaux,nzout,info)
import
integer(psb_ipk_), intent(in) :: nr,nc,nzin,dupl
integer(psb_ipk_), intent(inout) :: ia(:), ja(:), iaux(:)
complex(psb_spk_), intent(inout) :: val(:)
integer(psb_ipk_), intent(out) :: nzout
integer(psb_ipk_), intent(out) :: info
end subroutine psb_c_fix_coo_inner_rowmajor
end subroutine psb_c_fix_coo_inner_colmajor
end interface
interface
subroutine psb_c_fix_coo_inner_colmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,info)
subroutine psb_c_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,&
& ia,ja,val,iaux,nzout,info)
import
integer(psb_ipk_), intent(in) :: nr,nc,nzin,dupl
integer(psb_ipk_), intent(inout) :: ia(:), ja(:), iaux(:)
complex(psb_spk_), intent(inout) :: val(:)
integer(psb_ipk_), intent(out) :: nzout
integer(psb_ipk_), intent(out) :: info
end subroutine psb_c_fix_coo_inner_colmajor
end subroutine psb_c_fix_coo_inner_rowmajor
end interface
!

@ -1865,26 +1865,29 @@ module psb_d_base_mat_mod
integer(psb_ipk_), intent(in), optional :: idir
end subroutine psb_d_fix_coo_inner
end interface
interface
subroutine psb_d_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,info)
subroutine psb_d_fix_coo_inner_colmajor(nr,nc,nzin,dupl,&
& ia,ja,val,iaux,nzout,info)
import
integer(psb_ipk_), intent(in) :: nr,nc,nzin,dupl
integer(psb_ipk_), intent(inout) :: ia(:), ja(:), iaux(:)
real(psb_dpk_), intent(inout) :: val(:)
integer(psb_ipk_), intent(out) :: nzout
integer(psb_ipk_), intent(out) :: info
end subroutine psb_d_fix_coo_inner_rowmajor
end subroutine psb_d_fix_coo_inner_colmajor
end interface
interface
subroutine psb_d_fix_coo_inner_colmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,info)
subroutine psb_d_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,&
& ia,ja,val,iaux,nzout,info)
import
integer(psb_ipk_), intent(in) :: nr,nc,nzin,dupl
integer(psb_ipk_), intent(inout) :: ia(:), ja(:), iaux(:)
real(psb_dpk_), intent(inout) :: val(:)
integer(psb_ipk_), intent(out) :: nzout
integer(psb_ipk_), intent(out) :: info
end subroutine psb_d_fix_coo_inner_colmajor
end subroutine psb_d_fix_coo_inner_rowmajor
end interface
!

@ -1865,26 +1865,29 @@ module psb_s_base_mat_mod
integer(psb_ipk_), intent(in), optional :: idir
end subroutine psb_s_fix_coo_inner
end interface
interface
subroutine psb_s_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,info)
subroutine psb_s_fix_coo_inner_colmajor(nr,nc,nzin,dupl,&
& ia,ja,val,iaux,nzout,info)
import
integer(psb_ipk_), intent(in) :: nr,nc,nzin,dupl
integer(psb_ipk_), intent(inout) :: ia(:), ja(:), iaux(:)
real(psb_spk_), intent(inout) :: val(:)
integer(psb_ipk_), intent(out) :: nzout
integer(psb_ipk_), intent(out) :: info
end subroutine psb_s_fix_coo_inner_rowmajor
end subroutine psb_s_fix_coo_inner_colmajor
end interface
interface
subroutine psb_s_fix_coo_inner_colmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,info)
subroutine psb_s_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,&
& ia,ja,val,iaux,nzout,info)
import
integer(psb_ipk_), intent(in) :: nr,nc,nzin,dupl
integer(psb_ipk_), intent(inout) :: ia(:), ja(:), iaux(:)
real(psb_spk_), intent(inout) :: val(:)
integer(psb_ipk_), intent(out) :: nzout
integer(psb_ipk_), intent(out) :: info
end subroutine psb_s_fix_coo_inner_colmajor
end subroutine psb_s_fix_coo_inner_rowmajor
end interface
!

@ -1865,26 +1865,29 @@ module psb_z_base_mat_mod
integer(psb_ipk_), intent(in), optional :: idir
end subroutine psb_z_fix_coo_inner
end interface
interface
subroutine psb_z_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,info)
subroutine psb_z_fix_coo_inner_colmajor(nr,nc,nzin,dupl,&
& ia,ja,val,iaux,nzout,info)
import
integer(psb_ipk_), intent(in) :: nr,nc,nzin,dupl
integer(psb_ipk_), intent(inout) :: ia(:), ja(:), iaux(:)
complex(psb_dpk_), intent(inout) :: val(:)
integer(psb_ipk_), intent(out) :: nzout
integer(psb_ipk_), intent(out) :: info
end subroutine psb_z_fix_coo_inner_rowmajor
end subroutine psb_z_fix_coo_inner_colmajor
end interface
interface
subroutine psb_z_fix_coo_inner_colmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,info)
subroutine psb_z_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,&
& ia,ja,val,iaux,nzout,info)
import
integer(psb_ipk_), intent(in) :: nr,nc,nzin,dupl
integer(psb_ipk_), intent(inout) :: ia(:), ja(:), iaux(:)
complex(psb_dpk_), intent(inout) :: val(:)
integer(psb_ipk_), intent(out) :: nzout
integer(psb_ipk_), intent(out) :: info
end subroutine psb_z_fix_coo_inner_colmajor
end subroutine psb_z_fix_coo_inner_rowmajor
end interface
!

@ -728,7 +728,7 @@ subroutine psb_c_coo_print(iout,a,iv,head,ivr,ivc)
character(len=80) :: frmt
integer(psb_ipk_) :: i,j, ni, nr, nc, nz
write(iout,'(a)') '%%MatrixMarket matrix coordinate complex general'
write(iout,'(a)') '%%MatrixMarket matrix coordinate complex general'
if (present(head)) write(iout,'(a,a)') '% ',head
write(iout,'(a)') '%'
write(iout,'(a,a)') '% COO'
@ -3172,9 +3172,9 @@ subroutine psb_c_cp_coo_from_coo(a,b,info)
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, nz
a%ia(i) = b%ia(i)
a%ja(i) = b%ja(i)
a%val(i) = b%val(i)
a%ia(i) = b%ia(i)
a%ja(i) = b%ja(i)
a%val(i) = b%val(i)
end do
end block
#else
@ -5284,13 +5284,13 @@ function psb_lc_coo_maxval(a) result(res)
if (allocated(a%val)) then
nnz = min(nnz,size(a%val))
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
do i=1, nnz
res = max(res,abs(a%val(i)))
end do
end block
block
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, nnz
res = max(res,abs(a%val(i)))
end do
end block
#else
res = maxval(abs(a%val(1:nnz)))
#endif
@ -5351,13 +5351,13 @@ function psb_lc_coo_csnmi(a) result(res)
vt(i) = vt(i) + abs(a%val(j))
end do
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
do i=1, m
res = max(res,abs(vt(i)))
end do
end block
block
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, m
res = max(res,abs(vt(i)))
end do
end block
#else
res = maxval(vt(1:m))
#endif
@ -5403,7 +5403,7 @@ function psb_lc_coo_csnm1(a) result(res)
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
!$omp parallel do private(i)
do i=1, n
res = max(res,abs(vt(i)))
end do
@ -6856,7 +6856,7 @@ subroutine psb_lc_coo_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info)
if (nz < 0) then
info = psb_err_iarg_neg_
3 call psb_errpush(info,name,i_err=(/1_psb_ipk_/))
3 call psb_errpush(info,name,i_err=(/1_psb_ipk_/))
goto 9999
end if
if (size(ia) < nz) then

@ -1328,8 +1328,8 @@ function psb_c_csr_csnmi(a) result(res)
res = szero
if (a%is_dev()) call a%sync()
!$omp parallel do private(i,acc) reduction(max: res)
do i = 1, a%get_nrows()
!$omp parallel do private(i,j,acc) reduction(max: res)
do i = 1, a%get_nrows()
acc = szero
do j=a%irp(i),a%irp(i+1)-1
acc = acc + abs(a%val(j))
@ -1562,8 +1562,12 @@ subroutine psb_c_csr_get_diag(a,d,info)
if (a%is_unit()) then
d(1:mnm) = cone
!$omp parallel do private(i)
do i=1, mnm
d(i) = cone
end do
else
!$omp parallel do private(i,j,k)
do i=1, mnm
d(i) = czero
do k=a%irp(i),a%irp(i+1)-1
@ -1574,6 +1578,7 @@ subroutine psb_c_csr_get_diag(a,d,info)
enddo
end do
end if
!$omp parallel do private(i)
do i=mnm+1,size(d)
d(i) = czero
end do
@ -1629,6 +1634,7 @@ subroutine psb_c_csr_scal(d,a,info,side)
goto 9999
end if
!$omp parallel do private(i,j)
do i=1, m
do j = a%irp(i), a%irp(i+1) -1
a%val(j) = a%val(j) * d(i)
@ -1643,6 +1649,7 @@ subroutine psb_c_csr_scal(d,a,info,side)
goto 9999
end if
!$omp parallel do private(i,j)
do i=1,a%get_nzeros()
j = a%ja(i)
a%val(i) = a%val(i) * d(j)
@ -1681,6 +1688,7 @@ subroutine psb_c_csr_scals(d,a,info)
call a%make_nonunit()
end if
!$omp parallel do private(i)
do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d
enddo

@ -728,7 +728,7 @@ subroutine psb_d_coo_print(iout,a,iv,head,ivr,ivc)
character(len=80) :: frmt
integer(psb_ipk_) :: i,j, ni, nr, nc, nz
write(iout,'(a)') '%%MatrixMarket matrix coordinate real general'
write(iout,'(a)') '%%MatrixMarket matrix coordinate real general'
if (present(head)) write(iout,'(a,a)') '% ',head
write(iout,'(a)') '%'
write(iout,'(a,a)') '% COO'
@ -3172,9 +3172,9 @@ subroutine psb_d_cp_coo_from_coo(a,b,info)
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, nz
a%ia(i) = b%ia(i)
a%ja(i) = b%ja(i)
a%val(i) = b%val(i)
a%ia(i) = b%ia(i)
a%ja(i) = b%ja(i)
a%val(i) = b%val(i)
end do
end block
#else
@ -5284,13 +5284,13 @@ function psb_ld_coo_maxval(a) result(res)
if (allocated(a%val)) then
nnz = min(nnz,size(a%val))
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
do i=1, nnz
res = max(res,abs(a%val(i)))
end do
end block
block
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, nnz
res = max(res,abs(a%val(i)))
end do
end block
#else
res = maxval(abs(a%val(1:nnz)))
#endif
@ -5351,13 +5351,13 @@ function psb_ld_coo_csnmi(a) result(res)
vt(i) = vt(i) + abs(a%val(j))
end do
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
do i=1, m
res = max(res,abs(vt(i)))
end do
end block
block
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, m
res = max(res,abs(vt(i)))
end do
end block
#else
res = maxval(vt(1:m))
#endif
@ -5403,7 +5403,7 @@ function psb_ld_coo_csnm1(a) result(res)
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
!$omp parallel do private(i)
do i=1, n
res = max(res,abs(vt(i)))
end do
@ -6856,7 +6856,7 @@ subroutine psb_ld_coo_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info)
if (nz < 0) then
info = psb_err_iarg_neg_
3 call psb_errpush(info,name,i_err=(/1_psb_ipk_/))
3 call psb_errpush(info,name,i_err=(/1_psb_ipk_/))
goto 9999
end if
if (size(ia) < nz) then

@ -1328,8 +1328,8 @@ function psb_d_csr_csnmi(a) result(res)
res = dzero
if (a%is_dev()) call a%sync()
!$omp parallel do private(i,acc) reduction(max: res)
do i = 1, a%get_nrows()
!$omp parallel do private(i,j,acc) reduction(max: res)
do i = 1, a%get_nrows()
acc = dzero
do j=a%irp(i),a%irp(i+1)-1
acc = acc + abs(a%val(j))
@ -1562,8 +1562,12 @@ subroutine psb_d_csr_get_diag(a,d,info)
if (a%is_unit()) then
d(1:mnm) = done
!$omp parallel do private(i)
do i=1, mnm
d(i) = done
end do
else
!$omp parallel do private(i,j,k)
do i=1, mnm
d(i) = dzero
do k=a%irp(i),a%irp(i+1)-1
@ -1574,6 +1578,7 @@ subroutine psb_d_csr_get_diag(a,d,info)
enddo
end do
end if
!$omp parallel do private(i)
do i=mnm+1,size(d)
d(i) = dzero
end do
@ -1629,6 +1634,7 @@ subroutine psb_d_csr_scal(d,a,info,side)
goto 9999
end if
!$omp parallel do private(i,j)
do i=1, m
do j = a%irp(i), a%irp(i+1) -1
a%val(j) = a%val(j) * d(i)
@ -1643,6 +1649,7 @@ subroutine psb_d_csr_scal(d,a,info,side)
goto 9999
end if
!$omp parallel do private(i,j)
do i=1,a%get_nzeros()
j = a%ja(i)
a%val(i) = a%val(i) * d(j)
@ -1681,6 +1688,7 @@ subroutine psb_d_csr_scals(d,a,info)
call a%make_nonunit()
end if
!$omp parallel do private(i)
do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d
enddo

@ -728,7 +728,7 @@ subroutine psb_s_coo_print(iout,a,iv,head,ivr,ivc)
character(len=80) :: frmt
integer(psb_ipk_) :: i,j, ni, nr, nc, nz
write(iout,'(a)') '%%MatrixMarket matrix coordinate real general'
write(iout,'(a)') '%%MatrixMarket matrix coordinate real general'
if (present(head)) write(iout,'(a,a)') '% ',head
write(iout,'(a)') '%'
write(iout,'(a,a)') '% COO'
@ -3172,9 +3172,9 @@ subroutine psb_s_cp_coo_from_coo(a,b,info)
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, nz
a%ia(i) = b%ia(i)
a%ja(i) = b%ja(i)
a%val(i) = b%val(i)
a%ia(i) = b%ia(i)
a%ja(i) = b%ja(i)
a%val(i) = b%val(i)
end do
end block
#else
@ -5284,13 +5284,13 @@ function psb_ls_coo_maxval(a) result(res)
if (allocated(a%val)) then
nnz = min(nnz,size(a%val))
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
do i=1, nnz
res = max(res,abs(a%val(i)))
end do
end block
block
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, nnz
res = max(res,abs(a%val(i)))
end do
end block
#else
res = maxval(abs(a%val(1:nnz)))
#endif
@ -5351,13 +5351,13 @@ function psb_ls_coo_csnmi(a) result(res)
vt(i) = vt(i) + abs(a%val(j))
end do
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
do i=1, m
res = max(res,abs(vt(i)))
end do
end block
block
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, m
res = max(res,abs(vt(i)))
end do
end block
#else
res = maxval(vt(1:m))
#endif
@ -5403,7 +5403,7 @@ function psb_ls_coo_csnm1(a) result(res)
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
!$omp parallel do private(i)
do i=1, n
res = max(res,abs(vt(i)))
end do
@ -6856,7 +6856,7 @@ subroutine psb_ls_coo_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info)
if (nz < 0) then
info = psb_err_iarg_neg_
3 call psb_errpush(info,name,i_err=(/1_psb_ipk_/))
3 call psb_errpush(info,name,i_err=(/1_psb_ipk_/))
goto 9999
end if
if (size(ia) < nz) then

@ -1328,8 +1328,8 @@ function psb_s_csr_csnmi(a) result(res)
res = szero
if (a%is_dev()) call a%sync()
!$omp parallel do private(i,acc) reduction(max: res)
do i = 1, a%get_nrows()
!$omp parallel do private(i,j,acc) reduction(max: res)
do i = 1, a%get_nrows()
acc = szero
do j=a%irp(i),a%irp(i+1)-1
acc = acc + abs(a%val(j))
@ -1562,8 +1562,12 @@ subroutine psb_s_csr_get_diag(a,d,info)
if (a%is_unit()) then
d(1:mnm) = sone
!$omp parallel do private(i)
do i=1, mnm
d(i) = sone
end do
else
!$omp parallel do private(i,j,k)
do i=1, mnm
d(i) = szero
do k=a%irp(i),a%irp(i+1)-1
@ -1574,6 +1578,7 @@ subroutine psb_s_csr_get_diag(a,d,info)
enddo
end do
end if
!$omp parallel do private(i)
do i=mnm+1,size(d)
d(i) = szero
end do
@ -1629,6 +1634,7 @@ subroutine psb_s_csr_scal(d,a,info,side)
goto 9999
end if
!$omp parallel do private(i,j)
do i=1, m
do j = a%irp(i), a%irp(i+1) -1
a%val(j) = a%val(j) * d(i)
@ -1643,6 +1649,7 @@ subroutine psb_s_csr_scal(d,a,info,side)
goto 9999
end if
!$omp parallel do private(i,j)
do i=1,a%get_nzeros()
j = a%ja(i)
a%val(i) = a%val(i) * d(j)
@ -1681,6 +1688,7 @@ subroutine psb_s_csr_scals(d,a,info)
call a%make_nonunit()
end if
!$omp parallel do private(i)
do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d
enddo

@ -728,7 +728,7 @@ subroutine psb_z_coo_print(iout,a,iv,head,ivr,ivc)
character(len=80) :: frmt
integer(psb_ipk_) :: i,j, ni, nr, nc, nz
write(iout,'(a)') '%%MatrixMarket matrix coordinate complex general'
write(iout,'(a)') '%%MatrixMarket matrix coordinate complex general'
if (present(head)) write(iout,'(a,a)') '% ',head
write(iout,'(a)') '%'
write(iout,'(a,a)') '% COO'
@ -3172,9 +3172,9 @@ subroutine psb_z_cp_coo_from_coo(a,b,info)
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, nz
a%ia(i) = b%ia(i)
a%ja(i) = b%ja(i)
a%val(i) = b%val(i)
a%ia(i) = b%ia(i)
a%ja(i) = b%ja(i)
a%val(i) = b%val(i)
end do
end block
#else
@ -5284,13 +5284,13 @@ function psb_lz_coo_maxval(a) result(res)
if (allocated(a%val)) then
nnz = min(nnz,size(a%val))
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
do i=1, nnz
res = max(res,abs(a%val(i)))
end do
end block
block
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, nnz
res = max(res,abs(a%val(i)))
end do
end block
#else
res = maxval(abs(a%val(1:nnz)))
#endif
@ -5351,13 +5351,13 @@ function psb_lz_coo_csnmi(a) result(res)
vt(i) = vt(i) + abs(a%val(j))
end do
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
do i=1, m
res = max(res,abs(vt(i)))
end do
end block
block
integer(psb_ipk_) :: i
!$omp parallel do private(i)
do i=1, m
res = max(res,abs(vt(i)))
end do
end block
#else
res = maxval(vt(1:m))
#endif
@ -5403,7 +5403,7 @@ function psb_lz_coo_csnm1(a) result(res)
#if defined(OPENMP)
block
integer(psb_ipk_) :: i
!$omp parallel do private(i) reduction(max: res)
!$omp parallel do private(i)
do i=1, n
res = max(res,abs(vt(i)))
end do
@ -6856,7 +6856,7 @@ subroutine psb_lz_coo_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info)
if (nz < 0) then
info = psb_err_iarg_neg_
3 call psb_errpush(info,name,i_err=(/1_psb_ipk_/))
3 call psb_errpush(info,name,i_err=(/1_psb_ipk_/))
goto 9999
end if
if (size(ia) < nz) then

@ -1328,8 +1328,8 @@ function psb_z_csr_csnmi(a) result(res)
res = dzero
if (a%is_dev()) call a%sync()
!$omp parallel do private(i,acc) reduction(max: res)
do i = 1, a%get_nrows()
!$omp parallel do private(i,j,acc) reduction(max: res)
do i = 1, a%get_nrows()
acc = dzero
do j=a%irp(i),a%irp(i+1)-1
acc = acc + abs(a%val(j))
@ -1562,8 +1562,12 @@ subroutine psb_z_csr_get_diag(a,d,info)
if (a%is_unit()) then
d(1:mnm) = zone
!$omp parallel do private(i)
do i=1, mnm
d(i) = zone
end do
else
!$omp parallel do private(i,j,k)
do i=1, mnm
d(i) = zzero
do k=a%irp(i),a%irp(i+1)-1
@ -1574,6 +1578,7 @@ subroutine psb_z_csr_get_diag(a,d,info)
enddo
end do
end if
!$omp parallel do private(i)
do i=mnm+1,size(d)
d(i) = zzero
end do
@ -1629,6 +1634,7 @@ subroutine psb_z_csr_scal(d,a,info,side)
goto 9999
end if
!$omp parallel do private(i,j)
do i=1, m
do j = a%irp(i), a%irp(i+1) -1
a%val(j) = a%val(j) * d(i)
@ -1643,6 +1649,7 @@ subroutine psb_z_csr_scal(d,a,info,side)
goto 9999
end if
!$omp parallel do private(i,j)
do i=1,a%get_nzeros()
j = a%ja(i)
a%val(i) = a%val(i) * d(j)
@ -1681,6 +1688,7 @@ subroutine psb_z_csr_scals(d,a,info)
call a%make_nonunit()
end if
!$omp parallel do private(i)
do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d
enddo

@ -569,112 +569,112 @@ contains
call psb_cdasb(desc_a,info,mold=imold)
tcdasb = psb_wtime()-t1
if (.false.) then
!
! Add extra rows to test remote build.
!
block
integer(psb_ipk_) :: ks, i
ks = desc_a%get_local_cols()-desc_a%get_local_rows()
if (ks > 0) ks = max(1,ks / 10)
mysz = nlr+ks
call psb_realloc(mysz,myidx,info)
do i=nlr+1, mysz
myidx(i) = i
end do
call desc_a%l2gv1(myidx(nlr+1:mysz),info)
!write(0,*) iam,' Check on extra nodes ',nlr,mysz,':',myidx(nlr+1:mysz)
do ii= nlr+1, mysz, nb
ib = min(nb,mysz-ii+1)
icoeff = 1
do k=1,ib
i=ii+k-1
! local matrix pointer
glob_row=myidx(i)
! compute gridpoint coordinates
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
! x, y, z coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
z = (iz-1)*deltah
zt(k) = f_(x,y,z)
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
if (ix == 1) then
zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y-1,z)
val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
if (iy == 1) then
zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z-1)
val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
if (iz == 1) then
zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z)
val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
& + c(x,y,z)
call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
!
! Add extra rows
!
block
integer(psb_ipk_) :: ks, i
ks = desc_a%get_local_cols()-desc_a%get_local_rows()
if (ks > 0) ks = max(1,ks / 10)
mysz = nlr+ks
call psb_realloc(mysz,myidx,info)
do i=nlr+1, mysz
myidx(i) = i
end do
call desc_a%l2gv1(myidx(nlr+1:mysz),info)
!write(0,*) iam,' Check on extra nodes ',nlr,mysz,':',myidx(nlr+1:mysz)
do ii= nlr+1, mysz, nb
ib = min(nb,mysz-ii+1)
icoeff = 1
do k=1,ib
i=ii+k-1
! local matrix pointer
glob_row=myidx(i)
! compute gridpoint coordinates
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
! x, y, z coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
z = (iz-1)*deltah
zt(k) = f_(x,y,z)
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
if (ix == 1) then
zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
! term depending on (x,y,z+1)
val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
if (iz == idim) then
zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y+1,z)
val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
if (iy == idim) then
zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x+1,y,z)
val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
if (ix==idim) then
zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
endif
! term depending on (x,y-1,z)
val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
if (iy == 1) then
zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z-1)
val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
if (iz == 1) then
zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z)
val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
& + c(x,y,z)
call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
! term depending on (x,y,z+1)
val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
if (iz == idim) then
zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y+1,z)
val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
if (iy == idim) then
zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x+1,y,z)
val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
if (ix==idim) then
zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
end do
call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
if(info /= psb_success_) exit
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
if(info /= psb_success_) exit
zt(:)=dzero
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
if(info /= psb_success_) exit
end do
call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
if(info /= psb_success_) exit
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
if(info /= psb_success_) exit
zt(:)=dzero
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
if(info /= psb_success_) exit
end do
end block
end block
end if
call psb_barrier(ctxt)
t1 = psb_wtime()

Loading…
Cancel
Save