Changed handling of duplicates and storage format. Changed interface

to both geins and spasb. New and better strategies.
psblas3-type-indexed
Salvatore Filippone 19 years ago
parent 77106751a4
commit 4b2f930cf6

@ -80,28 +80,31 @@
!
! Queries into spmat%info
!
integer, parameter :: psb_root_=0
integer, parameter :: psb_nztotreq_=1, psb_nzrowreq_=2
integer, parameter :: psb_nzsizereq_=3
!
! Entries and values for spmat%info
!
integer, parameter :: psb_nnz_=1, psb_dupl_=5
integer, parameter :: psb_del_bnd_=6, psb_srtd_=7
integer, parameter :: psb_state_=8, psb_upd_=9
integer, parameter :: psb_upd_pnt_=10, psb_ifasize_=10
integer, parameter :: psb_nnz_=1
integer, parameter :: psb_del_bnd_=7, psb_srtd_=8
integer, parameter :: psb_state_=9
integer, parameter :: psb_upd_pnt_=10
integer, parameter :: psb_dupl_=11, psb_upd_=12
integer, parameter :: psb_ifasize_=16
integer, parameter :: psb_spmat_null_=0, psb_spmat_bld_=1
integer, parameter :: psb_spmat_asb_=2, psb_spmat_upd_=4
integer, parameter :: psb_ireg_flgs_=10, psb_ip2_=0
integer, parameter :: psb_iflag_=2, psb_ichk_=3
integer, parameter :: psb_nnzt_=4, psb_zero_=5,psb_ipc_=6
integer, parameter :: psb_dupl_err_ =1
integer, parameter :: psb_dupl_ovwrt_=2
integer, parameter :: psb_dupl_add_ =3
integer, parameter :: psb_perm_update_=98765
integer, parameter :: psb_srch_update_=98764
integer, parameter :: psb_isrtdcoo_ =98761
integer, parameter :: psb_dupl_ovwrt_ = 0
integer, parameter :: psb_dupl_add_ = 1
integer, parameter :: psb_dupl_err_ = 2
integer, parameter :: psb_dupl_def_ = psb_dupl_ovwrt_
integer, parameter :: psb_upd_dflt_ = 0
integer, parameter :: psb_upd_perm_ = 98765
integer, parameter :: psb_upd_srch_ = 98764
integer, parameter :: psb_isrtdcoo_ = 98761
integer, parameter :: psb_maxjdrows_=8, psb_minjdrows_=4
integer, parameter :: psb_dbleint_=2
!

@ -68,36 +68,37 @@ Contains
if(psb_get_errstatus().ne.0) return
info=0
if (associated(rrax)) then
dim=size(rrax)
If (dim /= len) Then
Allocate(tmp(len),stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
if (.true.) then
do i=1, min(len,dim)
tmp(i)=rrax(i)
end do
else
tmp(1:min(len,dim))=rrax(1:min(len,dim))
end if
deallocate(rrax,stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
rrax=>tmp
end if
else
allocate(rrax(len),stat=info)
if (info /= 0) then
dim=size(rrax)
If (dim /= len) Then
Allocate(tmp(len),stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
end if
if (.true.) then
do i=1, min(len,dim)
tmp(i)=rrax(i)
end do
else
tmp(1:min(len,dim))=rrax(1:min(len,dim))
end if
deallocate(rrax,stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
rrax=>tmp
end if
else
dim = 0
allocate(rrax(len),stat=info)
if (info /= 0) then
err=4000
call psb_errpush(err,name)
goto 9999
end if
endif
if (present(pad)) then
rrax(dim+1:len) = pad

@ -47,7 +47,7 @@ module psb_spmat_type
! describe some chacteristics of sparse matrix
character(len=11) :: descra
! Contains some additional informations on sparse matrix
integer :: infoa(10)
integer :: infoa(psb_ifasize_)
! Contains sparse matrix coefficients
real(kind(1.d0)), pointer :: aspk(:)=>null()
! Contains indeces that describes sparse matrix structure
@ -63,7 +63,7 @@ module psb_spmat_type
! describe some chacteristics of sparse matrix
character(len=11) :: descra
! Contains some additional informations on sparse matrix
integer :: infoa(10)
integer :: infoa(psb_ifasize_)
! Contains sparse matrix coefficients
complex(kind(1.d0)), pointer :: aspk(:)=>null()
! Contains indeces that describes sparse matrix structure
@ -80,6 +80,14 @@ module psb_spmat_type
module procedure psb_dspclone, psb_zspclone
end interface
interface psb_sp_setifld
module procedure psb_dsp_setifld, psb_zsp_setifld
end interface
interface psb_sp_getifld
module procedure psb_dsp_getifld, psb_zsp_getifld
end interface
interface psb_sp_transfer
module procedure psb_dsp_transfer, psb_zsp_transfer
end interface
@ -212,7 +220,7 @@ contains
a%m=max(0,m)
a%k=max(0,k)
call psb_sp_reall(a,nnz,info)
if (debug) write(0,*) 'Check in ALLOCATE ',info,associated(a%pl),associated(a%pr)
a%pl(1)=0
a%pr(1)=0
! set infoa fields
@ -296,6 +304,8 @@ contains
call psb_realloc(max(1,a%m),a%pl,info)
if (info /= 0) return
call psb_realloc(max(1,a%k),a%pr,info)
if (debug) write(0,*) associated(a%ia1),associated(a%ia2),&
& associated(a%aspk),associated(a%pl),associated(a%pr),info
if (info /= 0) return
Return
@ -420,6 +430,63 @@ contains
End Subroutine psb_dsp_transfer
Subroutine psb_dsp_setifld(val,field,a,info)
implicit none
!....Parameters...
Type(psb_dspmat_type), intent(inout) :: A
Integer, intent(in) :: field,val
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
info = 0
!!$ call psb_realloc(psb_ifasize_,a%infoa,info)
if (info == 0) &
& call psb_setifield(val,field,a%infoa,psb_ifasize_,info)
Return
end subroutine psb_dsp_setifld
function psb_dsp_getifld(field,a,info)
implicit none
!....Parameters...
Type(psb_dspmat_type), intent(in) :: A
Integer, intent(in) :: field
Integer :: psb_dsp_getifld
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
integer :: val
info = 0
val = -1
if ((field < 1).or.(field > psb_ifasize_)) then
info = -1
psb_dsp_getifld = val
return
endif
!!$ if (.not.associated(a%infoa)) then
!!$ info = -2
!!$ return
!!$ endif
call psb_getifield(val,field,a%infoa,psb_ifasize_,info)
psb_dsp_getifld = val
Return
end function psb_dsp_getifld
subroutine psb_dsp_free(a,info)
implicit none
!....Parameters...
@ -765,6 +832,66 @@ contains
End Subroutine psb_zsp_transfer
Subroutine psb_zsp_setifld(val,field,a,info)
implicit none
!....Parameters...
Type(psb_zspmat_type), intent(inout) :: A
Integer, intent(in) :: field,val
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
info = 0
!!$ call psb_realloc(psb_ifasize_,a%infoa,info)
if (info == 0) &
& call psb_setifield(val,field,a%infoa,psb_ifasize_,info)
Return
end subroutine psb_zsp_setifld
function psb_zsp_getifld(field,a,info)
implicit none
!....Parameters...
Type(psb_zspmat_type), intent(in) :: A
Integer, intent(in) :: field
Integer :: psb_zsp_getifld
Integer, intent(out) :: info
!locals
logical, parameter :: debug=.false.
integer :: val
info = 0
val = -1
if ((field < 1).or.(field > psb_ifasize_)) then
info = -1
psb_zsp_getifld = val
return
endif
!!$ if (.not.associated(a%infoa)) then
!!$ info = -2
!!$ return
!!$ endif
call psb_getifield(val,field,a%infoa,psb_ifasize_,info)
psb_zsp_getifld = val
Return
end function psb_zsp_getifld
subroutine psb_zsp_free(a,info)
implicit none

@ -293,7 +293,7 @@ Module psb_tools_mod
interface psb_geins
! 2-D double precision version
subroutine psb_dins(m, n, x, ix, jx, blck, desc_a, info,&
& iblck, jblck)
& iblck, jblck,dupl)
use psb_descriptor_type
integer, intent(in) :: m,n
type(psb_desc_type), intent(in) :: desc_a
@ -302,10 +302,11 @@ Module psb_tools_mod
real(kind(1.d0)), intent(in) :: blck(:,:)
integer,intent(out) :: info
integer, optional, intent(in) :: iblck,jblck
integer, optional, intent(in) :: dupl
end subroutine psb_dins
! 2-D double precision square version
subroutine psb_dinsvm(m, x, ix, jx, blck, desc_a,info,&
& iblck)
& iblck,dupl)
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
@ -314,10 +315,11 @@ Module psb_tools_mod
real(kind(1.d0)), intent(in) :: blck(:)
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: dupl
end subroutine psb_dinsvm
! 1-D double precision version
subroutine psb_dinsvv(m, x, ix, blck, desc_a, info,&
& iblck,insflag)
& iblck,insflag,dupl)
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
@ -327,10 +329,11 @@ Module psb_tools_mod
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: insflag
integer, optional, intent(in) :: dupl
end subroutine psb_dinsvv
! 2-D integer version
subroutine psb_iins(m, n, x, ix, jx, blck, desc_a, info,&
& iblck, jblck)
& iblck, jblck,dupl)
use psb_descriptor_type
integer, intent(in) :: m,n
type(psb_desc_type), intent(in) :: desc_a
@ -339,10 +342,11 @@ Module psb_tools_mod
integer, intent(in) :: blck(:,:)
integer,intent(out) :: info
integer, optional, intent(in) :: iblck,jblck
integer, optional, intent(in) :: dupl
end subroutine psb_iins
! 2-D integer square version
subroutine psb_iinsvm(m, x, ix, jx, blck, desc_a,info,&
& iblck)
& iblck,dupl)
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
@ -351,10 +355,11 @@ Module psb_tools_mod
integer, intent(in) :: blck(:)
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: dupl
end subroutine psb_iinsvm
! 1-D integer version
subroutine psb_iinsvv(m, x, ix, blck, desc_a, info,&
& iblck,insflag)
& iblck,insflag,dupl)
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
@ -364,10 +369,11 @@ Module psb_tools_mod
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: insflag
integer, optional, intent(in) :: dupl
end subroutine psb_iinsvv
! 2-D double precision version
subroutine psb_zins(m, n, x, ix, jx, blck, desc_a, info,&
& iblck, jblck)
& iblck, jblck,dupl)
use psb_descriptor_type
integer, intent(in) :: m,n
type(psb_desc_type), intent(in) :: desc_a
@ -376,10 +382,11 @@ Module psb_tools_mod
complex(kind(1.d0)), intent(in) :: blck(:,:)
integer,intent(out) :: info
integer, optional, intent(in) :: iblck,jblck
integer, optional, intent(in) :: dupl
end subroutine psb_zins
! 2-D double precision square version
subroutine psb_zinsvm(m, x, ix, jx, blck, desc_a,info,&
& iblck)
& iblck,dupl)
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
@ -388,10 +395,11 @@ Module psb_tools_mod
complex(kind(1.d0)), intent(in) :: blck(:)
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: dupl
end subroutine psb_zinsvm
! 1-D double precision version
subroutine psb_zinsvv(m, x, ix, blck, desc_a, info,&
& iblck,insflag)
& iblck,insflag,dupl)
use psb_descriptor_type
integer, intent(in) :: m
type(psb_desc_type), intent(in) :: desc_a
@ -401,6 +409,7 @@ Module psb_tools_mod
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: insflag
integer, optional, intent(in) :: dupl
end subroutine psb_zinsvv
end interface
@ -523,23 +532,23 @@ Module psb_tools_mod
end interface
interface psb_spasb
subroutine psb_dspasb(a,desc_a, info, afmt, up, dup)
subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl)
use psb_descriptor_type
use psb_spmat_type
type(psb_dspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
integer,optional, intent(in) :: dup
character, optional, intent(in) :: afmt*5, up
integer,optional, intent(in) :: dupl, upd
character, optional, intent(in) :: afmt*5
end subroutine psb_dspasb
subroutine psb_zspasb(a,desc_a, info, afmt, up, dup)
subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl)
use psb_descriptor_type
use psb_spmat_type
type(psb_zspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
integer,optional, intent(in) :: dup
character, optional, intent(in) :: afmt*5, up
integer,optional, intent(in) :: dupl, upd
character, optional, intent(in) :: afmt*5
end subroutine psb_zspasb
end interface

@ -165,7 +165,8 @@ contains
goto 9999
end if
b%infoa(psb_upd_) = 6
call psb_sp_setifld(psb_dupl_ovwrt_,psb_dupl_,b,info)
call psb_sp_setifld(psb_upd_dflt_,psb_upd_,b,info)
b%fida = 'COO'
b%m=a%m
b%k=a%k

@ -165,7 +165,8 @@ contains
goto 9999
end if
b%infoa(psb_upd_) = 6
call psb_sp_setifld(psb_dupl_ovwrt_,psb_dupl_,b,info)
call psb_sp_setifld(psb_upd_dflt_,psb_upd_,b,info)
b%fida = 'COO'
b%m=a%m
b%k=a%k

@ -457,14 +457,14 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
! check on blacs grid
call blacs_gridinfo(icontxt, nprow, npcol, myrow, mycol)
if (nprow == -1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol /= 1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
endif
ia = 1
@ -477,25 +477,25 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
ib = 1
if (present(doswap)) then
idoswap = doswap
idoswap = doswap
else
idoswap = 1
idoswap = 1
endif
if (present(trans)) then
if ( (toupper(trans).eq.'N').or.(toupper(trans).eq.'T')) then
itrans = toupper(trans)
else if (toupper(trans).eq.'C') then
info = 3020
call psb_errpush(info,name)
goto 9999
else
info = 70
call psb_errpush(info,name)
goto 9999
end if
if ( (toupper(trans).eq.'N').or.(toupper(trans).eq.'T')) then
itrans = toupper(trans)
else if (toupper(trans).eq.'C') then
info = 3020
call psb_errpush(info,name)
goto 9999
else
info = 70
call psb_errpush(info,name)
goto 9999
end if
else
itrans = 'N'
itrans = 'N'
endif
m = desc_a%matrix_data(psb_m_)
@ -509,143 +509,143 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
liwork= 2*ncol
if (a%pr(1) /= 0) liwork = liwork + n * ik
if (a%pl(1) /= 0) liwork = liwork + m * ik
! write(0,*)'---->>>',work(1)
! write(0,*)'---->>>',work(1)
if (present(work)) then
if(size(work).ge.liwork) then
iwork => work
liwork=size(work)
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
if(size(work).ge.liwork) then
iwork => work
liwork=size(work)
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
end if
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
! checking for matrix correctness
call psb_chkmat(m,n,ia,ja,desc_a%matrix_data,info,iia,jja)
if(info.ne.0) then
info=4010
ch_err='psb_chkmat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
info=4010
ch_err='psb_chkmat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (itrans.eq.'N') then
! Matrix is not transposed
if((ja.ne.ix).or.(ia.ne.iy)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
! Matrix is not transposed
if((ja.ne.ix).or.(ia.ne.iy)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(n,ik,size(x),ix,jx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(m,ik,size(y),iy,jy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(n,ik,size(x),ix,jx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(m,ik,size(y),iy,jy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if (idoswap == 0) then
x(nrow+1:ncol)=dzero
else
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& dzero,x,desc_a,iwork,info,data=psb_comm_halo_)
end if
if (idoswap == 0) then
x(nrow+1:ncol)=dzero
else
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& dzero,x,desc_a,iwork,info,data=psb_comm_halo_)
end if
! local Matrix-vector product
call psb_csmm(alpha,a,x(iix:lldx),beta,y(iiy:lldy),info)
! local Matrix-vector product
call psb_csmm(alpha,a,x(iix:lldx),beta,y(iiy:lldy),info)
if(info.ne.0) then
info = 4011
call psb_errpush(info,name)
goto 9999
end if
if(info.ne.0) then
info = 4011
call psb_errpush(info,name)
goto 9999
end if
else
! Matrix is transposed
if((ja.ne.iy).or.(ia.ne.ix)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
! Matrix is transposed
if((ja.ne.iy).or.(ia.ne.ix)) then
! this case is not yet implemented
info = 3030
call psb_errpush(info,name)
goto 9999
end if
if(desc_a%ovrlap_elem(1).ne.-1) then
info = 3070
call psb_errpush(info,name)
goto 9999
end if
if(desc_a%ovrlap_elem(1).ne.-1) then
info = 3070
call psb_errpush(info,name)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(m,ik,size(x),ix,jx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(n,ik,size(y),iy,jy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
! checking for vectors correctness
call psb_chkvect(m,ik,size(x),ix,jx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(n,ik,size(y),iy,jy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
call psb_errpush(info,name)
goto 9999
end if
xp => x(iix:lldx)
yp => y(iiy:lldy)
xp => x(iix:lldx)
yp => y(iiy:lldy)
yp(nrow+1:ncol)=dzero
yp(nrow+1:ncol)=dzero
! local Matrix-vector product
call psb_csmm(alpha,a,xp,beta,yp,info,trans=itrans)
! local Matrix-vector product
call psb_csmm(alpha,a,xp,beta,yp,info,trans=itrans)
if(info.ne.0) then
info = 4010
ch_err='dcsmm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(idoswap /= 0)&
& call psi_swaptran(ior(psb_swap_send_,psb_swap_recv_),&
& done,yp,desc_a,iwork,info)
if(info.ne.0) then
info = 4010
ch_err='PSI_dSwapTran'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(info.ne.0) then
info = 4010
ch_err='dcsmm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(idoswap /= 0)&
& call psi_swaptran(ior(psb_swap_send_,psb_swap_recv_),&
& done,yp,desc_a,iwork,info)
if(info.ne.0) then
info = 4010
ch_err='PSI_dSwapTran'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
if(.not.present(work)) then
deallocate(iwork)
deallocate(iwork)
end if
nullify(iwork)
@ -656,8 +656,8 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(icontxt)
return
call psb_error(icontxt)
return
end if
return
end subroutine psb_dspmv

@ -81,6 +81,7 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
use psi_mod
use psb_check_mod
use psb_error_mod
use psb_string_mod
implicit none
real(kind(1.D0)), intent(in) :: alpha, beta
@ -159,16 +160,16 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
endif
if (present(unitd)) then
lunitd = unitd
lunitd = toupper(unitd)
else
lunitd = 'U'
endif
if (present(trans)) then
if ((trans.eq.'N').or.(trans.eq.'T')&
& .or.(trans.eq.'n').or.(trans.eq.'t')) then
itrans = trans
else if ((trans.eq.'C').or.(trans.eq.'c')) then
itrans = toupper(trans)
if((itrans.eq.'N').or.(itrans.eq.'T')) then
! Ok
else if (itrans.eq.'C') then
info = 3020
call psb_errpush(info,name)
goto 9999
@ -388,6 +389,7 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
use psi_mod
use psb_check_mod
use psb_error_mod
use psb_string_mod
real(kind(1.D0)), intent(in) :: alpha, beta
real(kind(1.d0)), intent(in), target :: x(:)
@ -422,14 +424,14 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
! check on blacs grid
call blacs_gridinfo(icontxt, nprow, npcol, myrow, mycol)
if (nprow == -1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol /= 1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
endif
! just this case right now
@ -448,25 +450,26 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
endif
if (present(unitd)) then
lunitd = unitd
lunitd = toupper(unitd)
else
lunitd = 'U'
endif
if (present(trans)) then
if((trans.eq.'N').or.(trans.eq.'T')) then
itrans = trans
else if (trans.eq.'C') then
info = 3020
call psb_errpush(info,name)
goto 9999
else
info = 70
call psb_errpush(info,name)
goto 9999
end if
itrans = toupper(trans)
if((itrans.eq.'N').or.(itrans.eq.'T')) then
! Ok
else if (itrans.eq.'C') then
info = 3020
call psb_errpush(info,name)
goto 9999
else
info = 70
call psb_errpush(info,name)
goto 9999
end if
else
itrans = 'N'
itrans = 'N'
endif
m = desc_a%matrix_data(psb_m_)
@ -476,9 +479,9 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
lldy = size(y)
if((lldx.lt.ncol).or.(lldy.lt.ncol)) then
info=3010
call psb_errpush(info,name)
goto 9999
info=3010
call psb_errpush(info,name)
goto 9999
end if
! check for presence/size of a work area
@ -486,34 +489,34 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
if (a%pr(1) /= 0) llwork = liwork + m * ik
if (a%pl(1) /= 0) llwork = llwork + m * ik
if (present(work)) then
if(size(work).lt.liwork) then
call psb_realloc(liwork,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
iwork => work
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
if(size(work).lt.liwork) then
call psb_realloc(liwork,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
end if
iwork => work
else
call psb_realloc(liwork,iwork,info)
if(info.ne.0) then
info=4010
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
iwork(1)=0.d0
if(present(d)) then
lld = size(d)
id => d
if (present(d)) then
lld = size(d)
id => d
else
lld=1
allocate(id(1))
id=1.d0
lld=1
allocate(id(1))
id=1.d0
end if
! checking for matrix correctness
@ -522,25 +525,25 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
call psb_chkvect(m,ik,size(x),ix,jx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(m,ik,size(y),iy,jy,desc_a%matrix_data,info,iiy,jjy)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect/mat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
info=4010
ch_err='psb_chkvect/mat'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(ja.ne.ix) then
! this case is not yet implemented
info = 3030
! this case is not yet implemented
info = 3030
end if
if((iix.ne.1).or.(iiy.ne.1)) then
! this case is not yet implemented
info = 3040
! this case is not yet implemented
info = 3040
end if
if(info.ne.0) then
call psb_errpush(info,name)
goto 9999
call psb_errpush(info,name)
goto 9999
end if
! Perform local triangular system solve
@ -549,43 +552,43 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
call psb_cssm(alpha,a,xp,beta,yp,info,unitd=lunitd,d=id,trans=itrans)
if(info.ne.0) then
info = 4010
ch_err='dcssm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
info = 4010
ch_err='dcssm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
! update overlap elements
if(lchoice.gt.0) then
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& done,yp,desc_a,iwork,info)
i=0
! switch on update type
select case (lchoice)
case(psb_square_root_)
do while(desc_a%ovrlap_elem(i).ne.-ione)
y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_)) =&
& y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_))/&
& sqrt(real(desc_a%ovrlap_elem(i+psb_n_dom_ovr_)))
i = i+2
end do
case(psb_avg_)
do while(desc_a%ovrlap_elem(i).ne.-ione)
y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_)) =&
& y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_))/&
& real(desc_a%ovrlap_elem(i+psb_n_dom_ovr_))
i = i+2
end do
case(psb_sum_)
! do nothing
case default
! wrong value for choice argument
info = 70
int_err=(/10,lchoice,0,0,0/)
call psb_errpush(info,name,i_err=int_err)
goto 9999
end select
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& done,yp,desc_a,iwork,info)
i=0
! switch on update type
select case (lchoice)
case(psb_square_root_)
do while(desc_a%ovrlap_elem(i).ne.-ione)
y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_)) =&
& y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_))/&
& sqrt(real(desc_a%ovrlap_elem(i+psb_n_dom_ovr_)))
i = i+2
end do
case(psb_avg_)
do while(desc_a%ovrlap_elem(i).ne.-ione)
y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_)) =&
& y(desc_a%ovrlap_elem(i+psb_ovrlp_elem_))/&
& real(desc_a%ovrlap_elem(i+psb_n_dom_ovr_))
i = i+2
end do
case(psb_sum_)
! do nothing
case default
! wrong value for choice argument
info = 70
int_err=(/10,lchoice,0,0,0/)
call psb_errpush(info,name,i_err=int_err)
goto 9999
end select
end if
if(.not.present(work)) deallocate(iwork)
@ -598,8 +601,8 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(icontxt)
return
call psb_error(icontxt)
return
end if
return
end subroutine psb_dspsv

@ -11,7 +11,8 @@ FOBJS = psb_cest.o psb_dcoins.o psb_dcsdp.o psb_dcsmm.o psb_dcsmv.o \
psb_zcsnmi.o psb_zcsrws.o psb_zcssm.o psb_zcssv.o psb_zcsdp.o\
psb_zfixcoo.o psb_zipcoo2csr.o psb_zipcsr2coo.o psb_zipcoo2csc.o \
psb_zcoins.o psb_zcsprt.o psb_zneigh.o psb_ztransp.o \
psb_zrwextd.o psb_zsymbmm.o psb_znumbmm.o psb_zspinfo.o psb_zspscal.o
psb_zrwextd.o psb_zsymbmm.o psb_znumbmm.o psb_zspinfo.o psb_zspscal.o\
psb_getifield.o psb_setifield.o
INCDIRS = -I ../../lib -I .

@ -4,7 +4,7 @@ include ../../../Make.inc
#
FOBJS = isr.o isrx.o \
mrgsrt.o isaperm.o ibsrch.o
mrgsrt.o isaperm.o ibsrch.o imsr.o imsrx.o
OBJS=$(FOBJS)

@ -30,7 +30,8 @@
!!$
! File: imsr.f90
! Subroutine:
! Parameters:subroutine imsr(n,x)
! Parameters:
subroutine imsr(n,x)
integer :: n
integer :: x(n)

@ -30,7 +30,8 @@
!!$
! File: imsrx.f90
! Subroutine:
! Parameters:subroutine imsrx(n,x,indx)
! Parameters:
subroutine imsrx(n,x,indx)
integer :: n
integer :: x(n)
integer :: indx(n)

@ -71,13 +71,14 @@ c
ierror = 0
call fcpsb_erractionsave(err_act)
check_flag=ibits(info(psb_upd_),1,2)
call psb_getifield(check_flag,psb_dupl_,info,psb_ifasize_,ierror)
if (trans.eq.'N') then
scale = (unitd.eq.'L') ! meaningless
p1(1) = 0
p2(1) = 0
nnz = info(psb_nnz_)
call psb_getifield(nnz,psb_nnz_,info,psb_ifasize_,ierror)
if (debug) then
write(*,*) 'on entry to dcoco: nnz laux ',
+ nnz,laux,larn,lia1n,lia2n
@ -171,16 +172,16 @@ c ... insert remaining element ...
do elem_in = 2, nnz
if ((ia1n(elem_in).eq.ia1n(elem_out)).and.
+ (ia2n(elem_in).eq.ia2n(elem_out))) then
if (check_flag.eq.1) then
if (check_flag.eq.psb_dupl_err_) then
c ... error, there are duplicated elements ...
ierror = 130
call fcpsb_errpush(ierror,name,int_val)
goto 9999
else if (check_flag.eq.2) then
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the first duplicated element ...
ia2n(ip2+aux(ipx+elem_in-1)-1) = elem_out
else if (check_flag.eq.3) then
c ... sum the duplicated element ...
else if (check_flag.eq.psb_dupl_add_) then
c ... add the duplicated element ...
arn(elem_out) = arn(elem_out) + arn(elem_in)
ia2n(ip2+aux(ipx+elem_in-1)-1) = elem_out
end if
@ -219,15 +220,15 @@ c ... insert remaining element ...
do elem_in = 2, nnz
if ((ia1n(elem_in).eq.ia1n(elem_out)).and.
+ (ia2n(elem_in).eq.ia2n(elem_out))) then
if (check_flag.eq.1) then
if (check_flag.eq.psb_dupl_err_) then
c ... error, there are duplicated elements ...
ierror = 130
call fcpsb_errpush(ierror,name,int_val)
goto 9999
else if (check_flag.eq.2) then
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the first duplicated element ...
else if (check_flag.eq.3) then
c ... sum the duplicated element ...
else if (check_flag.eq.psb_dupl_add_) then
c ... add the duplicated element ...
arn(elem_out) = arn(elem_out) + arn(elem_in)
end if
else

@ -73,8 +73,8 @@ C
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
CHECK_FLAG=IBITS(INFO(PSB_UPD_),1,2)
c$$$ write(0,*) 'DCOCR FLAG ',info(psb_upd_),check_flag
call psb_getifield(check_flag,psb_dupl_,info,psb_ifasize_,ierror)
IF ((TRANS.EQ.'N').or.(TRANS.EQ.'n')) THEN
SCALE = (UNITD.EQ.'L') ! meaningless
@ -221,16 +221,16 @@ C ... Insert other element of row ...
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr
ELEM_CSR = ELEM_CSR+1
ELSE
IF (CHECK_FLAG.EQ.1) THEN
IF (CHECK_FLAG.EQ.psb_dupl_err_) THEN
C ... Error, there are duplicated elements ...
IERROR = 130
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ELSE IF (CHECK_FLAG.EQ.2) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_ovwrt_) THEN
C ... Insert only the last duplicated element ...
ARN(ELEM_CSR-1) = AR(ELEM)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1
ELSE IF (CHECK_FLAG.EQ.3) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_add_) THEN
C ... Sum the duplicated element ...
ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1
@ -296,17 +296,17 @@ C ... Insert other element of row ...
ARN(ELEM_CSR) = AR(ELEM)
ELEM_CSR = ELEM_CSR+1
ELSE
IF (CHECK_FLAG.EQ.1) THEN
IF (CHECK_FLAG.EQ.psb_dupl_err_) THEN
C ... Error, there are duplicated elements ...
IERROR = 130
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ELSE IF (CHECK_FLAG.EQ.2) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_ovwrt_) THEN
C ... Insert only the last duplicated element ...
ARN(ELEM_CSR-1) = AR(ELEM)
if (debug) write(0,*) 'Duplicated overwrite',
+ elem_csr-1,elem
ELSE IF (CHECK_FLAG.EQ.3) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_add_) THEN
C ... Sum the duplicated element ...
ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM)
if (debug) write(0,*) 'Duplicated add',

@ -59,7 +59,9 @@ c .. Local Arrays ..
POINT_AR = 1
POINT_JA = 0
CHECK_FLAG=IBITS(INFON(PSB_UPD_),1,2)
call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror)
IF ((LARN.LT.POINT_AR).OR.(LKA.LT.POINT_AR)) THEN
IERROR = 60

@ -71,7 +71,7 @@ c
ierror = 0
call fcpsb_erractionsave(err_act)
check_flag=ibits(info(psb_upd_),1,2)
call psb_getifield(check_flag,psb_dupl_,info,psb_ifasize_,ierror)
if (trans.eq.'N') then
scale = (unitd.eq.'L') ! meaningless
p1(1) = 0
@ -172,15 +172,15 @@ c ... insert remaining element ...
do elem_in = 2, nnz
if ((ia1n(elem_in).eq.ia1n(elem_out)).and.
+ (ia2n(elem_in).eq.ia2n(elem_out))) then
if (check_flag.eq.1) then
if (check_flag.eq.psb_dupl_err_) then
c ... error, there are duplicated elements ...
ierror = 130
call fcpsb_errpush(ierror,name,int_val)
goto 9999
else if (check_flag.eq.2) then
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the first duplicated element ...
ia2n(ip2+aux(ipx+elem_in-1)-1) = elem_out
else if (check_flag.eq.3) then
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_out) = arn(elem_out) + arn(elem_in)
ia2n(ip2+aux(ipx+elem_in-1)-1) = elem_out
@ -220,14 +220,14 @@ c ... insert remaining element ...
do elem_in = 2, nnz
if ((ia1n(elem_in).eq.ia1n(elem_out)).and.
+ (ia2n(elem_in).eq.ia2n(elem_out))) then
if (check_flag.eq.1) then
if (check_flag.eq.psb_dupl_err_) then
c ... error, there are duplicated elements ...
ierror = 130
call fcpsb_errpush(ierror,name,int_val)
goto 9999
else if (check_flag.eq.2) then
else if (check_flag.eq.psb_dupl_ovwrt_) then
c ... insert only the first duplicated element ...
else if (check_flag.eq.3) then
else if (check_flag.eq.psb_dupl_add_) then
c ... sum the duplicated element ...
arn(elem_out) = arn(elem_out) + arn(elem_in)
end if

@ -73,8 +73,8 @@ C
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
CHECK_FLAG=IBITS(INFO(PSB_UPD_),1,2)
c$$$ write(0,*) 'ZCOCR FLAG ',info(psb_upd_),check_flag
call psb_getifield(check_flag,psb_dupl_,info,psb_ifasize_,ierror)
IF ((TRANS.EQ.'N').or.(TRANS.EQ.'n')) THEN
SCALE = (UNITD.EQ.'L') ! meaningless
@ -221,16 +221,16 @@ C ... Insert other element of row ...
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr
ELEM_CSR = ELEM_CSR+1
ELSE
IF (CHECK_FLAG.EQ.1) THEN
IF (CHECK_FLAG.EQ.psb_dupl_err_) THEN
C ... Error, there are duplicated elements ...
IERROR = 130
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ELSE IF (CHECK_FLAG.EQ.2) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_ovwrt_) THEN
C ... Insert only the last duplicated element ...
ARN(ELEM_CSR-1) = AR(ELEM)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1
ELSE IF (CHECK_FLAG.EQ.3) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_add_) THEN
C ... Sum the duplicated element ...
ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM)
ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1
@ -296,17 +296,17 @@ C ... Insert other element of row ...
ARN(ELEM_CSR) = AR(ELEM)
ELEM_CSR = ELEM_CSR+1
ELSE
IF (CHECK_FLAG.EQ.1) THEN
IF (CHECK_FLAG.EQ.psb_dupl_err_) THEN
C ... Error, there are duplicated elements ...
IERROR = 130
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ELSE IF (CHECK_FLAG.EQ.2) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_ovwrt_) THEN
C ... Insert only the last duplicated element ...
ARN(ELEM_CSR-1) = AR(ELEM)
if (debug) write(0,*) 'Duplicated overwrite',
+ elem_csr-1,elem
ELSE IF (CHECK_FLAG.EQ.3) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_add_) THEN
C ... Sum the duplicated element ...
ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM)
if (debug) write(0,*) 'Duplicated add',
@ -396,19 +396,19 @@ C ... Insert other element of row ...
ELEM_CSR = ELEM_CSR+1
ENDIF
ELSE
IF (CHECK_FLAG.EQ.1) THEN
IF (CHECK_FLAG.EQ.psb_dupl_err_) THEN
C ... Error, there are duplicated elements ...
IERROR = 130
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ELSE IF (CHECK_FLAG.EQ.2) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_ovwrt_) THEN
C ... Insert only the last duplicated element ...
IF(JA(ELEM).GT.IA(ELEM)) THEN
ARN(ELEM_CSR-1) = AR(ELEM)
ENDIF
if (debug) write(0,*) 'Duplicated overwrite',
+ elem_csr-1,elem
ELSE IF (CHECK_FLAG.EQ.3) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_add_) THEN
C ... Sum the duplicated element ...
IF(JA(ELEM).GT.IA(ELEM)) THEN
ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM)
@ -494,19 +494,19 @@ C ... Insert other element of row ...
ELEM_CSR = ELEM_CSR+1
ENDIF
ELSE
IF (CHECK_FLAG.EQ.1) THEN
IF (CHECK_FLAG.EQ.psb_dupl_err_) THEN
C ... Error, there are duplicated elements ...
IERROR = 130
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ELSE IF (CHECK_FLAG.EQ.2) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_ovwrt_) THEN
C ... Insert only the last duplicated element ...
IF(JA(ELEM).LT.IA(ELEM)) THEN
ARN(ELEM_CSR-1) = AR(ELEM)
ENDIF
if (debug) write(0,*) 'Duplicated overwrite',
+ elem_csr-1,elem
ELSE IF (CHECK_FLAG.EQ.3) THEN
ELSE IF (CHECK_FLAG.EQ.psb_dupl_add_) THEN
C ... Sum the duplicated element ...
IF(JA(ELEM).LT.IA(ELEM)) THEN
ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM)

@ -59,7 +59,7 @@ c .. Local Arrays ..
POINT_AR = 1
POINT_JA = 0
CHECK_FLAG=IBITS(INFON(PSB_UPD_),1,2)
call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror)
IF ((LARN.LT.POINT_AR).OR.(LKA.LT.POINT_AR)) THEN
IERROR = 60

@ -49,8 +49,8 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
character(len=5) :: ufida
integer :: i,j,ir,ic,nr,nc, ng, nza, isza,spstate, nnz,&
& ip1, nzl, err_act, int_err(5)
logical, parameter :: debug=.true.
& ip1, nzl, err_act, int_err(5), iupd
logical, parameter :: debug=.false.
character(len=20) :: name, ch_err
name='psb_dcoins'
@ -59,113 +59,124 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
info = 0
if (nz <= 0) then
info = 10
int_err(1)=1
call psb_errpush(info,name,i_err=int_err)
goto 9999
info = 10
int_err(1)=1
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
if (size(ia) < nz) then
info = 35
int_err(1)=2
call psb_errpush(info,name,i_err=int_err)
goto 9999
info = 35
int_err(1)=2
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
if (size(ja) < nz) then
info = 35
int_err(1)=3
call psb_errpush(info,name,i_err=int_err)
goto 9999
info = 35
int_err(1)=3
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
if (size(val) < nz) then
info = 35
int_err(1)=4
call psb_errpush(info,name,i_err=int_err)
goto 9999
info = 35
int_err(1)=4
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
!!$ ufida = toupper(a%fida)
call touppers(a%fida,ufida)
ng = size(gtl)
spstate = a%infoa(psb_state_)
spstate = psb_sp_getifld(psb_state_,a,info)
select case(spstate)
case(psb_spmat_bld_)
if ((ufida /= 'COO').and.(ufida/='COI')) then
info = 134
ch_err(1:3)=ufida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_spinfo(psb_nztotreq_,a,nza,info)
call psb_spinfo(psb_nzsizereq_,a,isza,info)
if(info.ne.izero) then
info=4010
ch_err='psb_spinfo'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if ((nza+nz)>isza) then
call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info)
if(info.ne.izero) then
info=4010
ch_err='psb_sp_reall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
endif
call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,&
& imin,imax,jmin,jmax,info)
if(info.ne.izero) then
if ((ufida /= 'COO').and.(ufida/='COI')) then
info = 134
ch_err(1:3)=ufida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_spinfo(psb_nztotreq_,a,nza,info)
call psb_spinfo(psb_nzsizereq_,a,isza,info)
if(info.ne.izero) then
info=4010
ch_err='psb_spinfo'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if ((nza+nz)>isza) then
call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info)
if(info.ne.izero) then
info=4010
ch_err='psb_inner_ins'
ch_err='psb_sp_reall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debug) then
if ((nza - a%infoa(psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: insert discarded items '
end if
end if
if ((nza - a%infoa(psb_nnz_)) /= nz) then
a%infoa(psb_del_bnd_) = nza
endif
a%infoa(psb_nnz_) = nza
endif
endif
call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,&
& imin,imax,jmin,jmax,info)
if(info.ne.izero) then
info=4010
ch_err='psb_inner_ins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debug) then
if ((nza - a%infoa(psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: insert discarded items '
end if
end if
if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then
call psb_sp_setifld(nza,psb_del_bnd_,a,info)
endif
call psb_sp_setifld(nza,psb_nnz_,a,info)
case(psb_spmat_upd_)
if (ibits(a%infoa(psb_upd_),2,1).eq.1) then
ip1 = a%infoa(psb_upd_pnt_)
nza = a%ia2(ip1+psb_nnz_)
nzl = a%infoa(psb_del_bnd_)
call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,&
& imin,imax,jmin,jmax,nzl,info)
if(info.ne.izero) then
info=4010
ch_err='psb_inner_upd'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
iupd = psb_sp_getifld(psb_upd_,a,info)
select case (iupd)
case (psb_upd_perm_)
ip1 = psb_sp_getifld(psb_upd_pnt_,a,info)
nzl = psb_sp_getifld(psb_del_bnd_,a,info)
nza = a%ia2(ip1+psb_nnz_)
call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,&
& imin,imax,jmin,jmax,nzl,info)
if(info.ne.izero) then
info=4010
ch_err='psb_inner_upd'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debug) then
if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: update discarded items '
end if
end if
a%ia2(ip1+psb_nnz_) = nza
a%ia2(ip1+psb_nnz_) = nza
else
info = 2231
call psb_errpush(info,name)
goto 9999
endif
if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza
case (psb_upd_dflt_, psb_upd_srch_)
write(0,*) 'Default & search inner update to be implemented'
info = 2230
call psb_errpush(info,name)
goto 9999
case default
info = 2231
call psb_errpush(info,name)
goto 9999
end select
case default
info = 2232
call psb_errpush(info,name)
goto 9999
info = 2232
call psb_errpush(info,name)
goto 9999
end select
return
@ -175,8 +186,8 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error()
return
call psb_error()
return
end if
return

@ -47,6 +47,7 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
use psb_const_mod
use psb_error_mod
use psb_spmat_type
use psb_string_mod
implicit none
!....Parameters...
@ -68,13 +69,14 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
Integer, Parameter :: maxtry=8
logical, parameter :: debug=.false.
character(len=20) :: name, ch_err
interface psb_cest
subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info)
integer, intent(in) :: m,n,nnz
integer, intent(out) :: lia1, lia2, lar, info
character, intent(inout) :: afmt*5
character, intent(in) :: up
end subroutine psb_cest
subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info)
integer, intent(in) :: m,n,nnz
integer, intent(out) :: lia1, lia2, lar, info
character, intent(inout) :: afmt*5
character, intent(in) :: up
end subroutine psb_cest
end interface
interface psb_spinfo
@ -98,17 +100,17 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
ifc_ = 1
endif
if (present(check)) then
check_ = check
check_ = toupper(check)
else
check_ = 'N'
endif
if (present(trans)) then
trans_ = trans
trans_ = toupper(trans )
else
trans_ = 'N'
endif
if (present(unitd)) then
unitd_ = unitd
unitd_ = toupper(unitd )
else
unitd_ = 'U'
endif
@ -132,7 +134,7 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
endif
if((check_=='Y').or.(check_=='C')) then
if(a%fida(1:3)=='CSR') then
if(toupper(a%fida(1:3))=='CSR') then
call dcsrck(trans,a%m,a%k,a%descra,a%aspk,a%ia1,a%ia2,work,size(work),info)
if(info /= 0) then
info=4010
@ -153,13 +155,14 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
b%k=a%k
call psb_spinfo(psb_nztotreq_,a,size_req,info)
if (debug) write(0,*) 'DCSDP : size_req 1:',size_req
!! PULL IUP FROM INFOA FIELD
iup = iand(b%infoa(psb_upd_),4)
if (iup > 0) then
!
iup = psb_sp_getifld(psb_upd_,b,info)
if (iup == psb_upd_perm_) then
up = 'Y'
else
up = 'N'
endif
n_row=b%m
n_col=b%k
call psb_cest(b%fida, n_row,n_col,size_req,&
@ -190,11 +193,11 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
b%pr(:) = 0
select case (a%fida(1:3))
select case (toupper(a%fida(1:3)))
case ('CSR')
select case (b%fida(1:3))
select case (toupper(b%fida(1:3)))
case ('CSR')
@ -281,7 +284,7 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
case ('COO','COI')
select case (b%fida(1:3))
select case (toupper(b%fida(1:3)))
case ('CSR')
@ -375,99 +378,116 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
else if (check_=='R') then
!...Regenerating matrix
if (b%infoa(psb_state_) /= psb_spmat_upd_) then
info = 8888
call psb_errpush(info,name)
goto 9999
endif
if (ibits(b%infoa(psb_upd_),2,1).eq.0) then
!
! Nothing to be done......
!
if (psb_sp_getifld(psb_state_,b,info) /= psb_spmat_upd_) then
info = 8888
call psb_errpush(info,name)
goto 9999
endif
select case(psb_sp_getifld(psb_upd_,b,info))
if (b%fida(1:3)/='JAD') then
ip1 = b%infoa(psb_upd_pnt_)
ip2 = b%ia2(ip1+psb_ip2_)
nnz = b%ia2(ip1+psb_nnz_)
iflag = b%ia2(ip1+psb_iflag_)
ichk = b%ia2(ip1+psb_ichk_)
nnzt = b%ia2(ip1+psb_nnzt_)
if (debug) write(*,*) 'Regeneration start: ',&
& b%infoa(psb_upd_),psb_perm_update_,nnz,nnzt ,iflag,info
if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then
info = 8889
write(*,*) 'Regeneration start error: ',&
& b%infoa(psb_upd_),psb_perm_update_,nnz,nnzt ,iflag,ichk
call psb_errpush(info,name)
goto 9999
endif
do i= 1, nnz
work(i) = 0.d0
enddo
if (iflag.eq.2) then
do i=1, nnz
work(b%ia2(ip2+i-1)) = b%aspk(i)
case(psb_upd_perm_)
if (toupper(b%fida(1:3))/='JAD') then
ip1 = psb_sp_getifld(psb_upd_pnt_,b,info)
ip2 = b%ia2(ip1+psb_ip2_)
nnz = b%ia2(ip1+psb_nnz_)
iflag = b%ia2(ip1+psb_iflag_)
ichk = b%ia2(ip1+psb_ichk_)
nnzt = b%ia2(ip1+psb_nnzt_)
if (debug) write(*,*) 'Regeneration start: ',&
& b%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,info
if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then
info = 8889
write(*,*) 'Regeneration start error: ',&
& b%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk
call psb_errpush(info,name)
goto 9999
endif
do i= 1, nnz
work(i) = 0.d0
enddo
else if (iflag.eq.3) then
select case(iflag)
case(psb_dupl_ovwrt_,psb_dupl_err_)
do i=1, nnz
work(b%ia2(ip2+i-1)) = b%aspk(i)
enddo
case(psb_dupl_add_)
do i=1, nnz
work(b%ia2(ip2+i-1)) = b%aspk(i) + work(b%ia2(ip2+i-1))
enddo
case default
info = 8887
call psb_errpush(info,name)
goto 9999
end select
do i=1, nnz
work(b%ia2(ip2+i-1)) = b%aspk(i) + work(b%ia2(ip2+i-1))
b%aspk(i) = work(i)
enddo
endif
do i=1, nnz
b%aspk(i) = work(i)
enddo
else if (b%fida(1:3) == 'JAD') then
ip1 = b%infoa(psb_upd_pnt_)
ip2 = b%ia1(ip1+psb_ip2_)
count = b%ia1(ip1+psb_zero_)
ipc = b%ia1(ip1+psb_ipc_)
nnz = b%ia1(ip1+psb_nnz_)
iflag = b%ia1(ip1+psb_iflag_)
ichk = b%ia1(ip1+psb_ichk_)
nnzt = b%ia1(ip1+psb_nnzt_)
if (debug) write(*,*) 'Regeneration start: ',&
& b%infoa(psb_upd_),psb_perm_update_,nnz,nnzt,count, &
& iflag,info
if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then
info = 10
write(*,*) 'Regeneration start error: ',&
& b%infoa(psb_upd_),psb_perm_update_,nnz,nnzt ,iflag,ichk
call psb_errpush(info,name)
goto 9999
endif
do i= 1, nnz+count
work(i) = 0.d0
enddo
if (iflag.eq.2) then
do i=1, nnz
work(b%ia1(ip2+i-1)) = b%aspk(i)
else if (toupper(b%fida(1:3)) == 'JAD') then
ip1 = psb_sp_getifld(psb_upd_pnt_,b,info)
ip2 = b%ia1(ip1+psb_ip2_)
count = b%ia1(ip1+psb_zero_)
ipc = b%ia1(ip1+psb_ipc_)
nnz = b%ia1(ip1+psb_nnz_)
iflag = b%ia1(ip1+psb_iflag_)
ichk = b%ia1(ip1+psb_ichk_)
nnzt = b%ia1(ip1+psb_nnzt_)
if (debug) write(*,*) 'Regeneration start: ',&
& b%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt,count, &
& iflag,info
if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then
info = 10
write(*,*) 'Regeneration start error: ',&
& b%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk
call psb_errpush(info,name)
goto 9999
endif
do i= 1, nnz+count
work(i) = 0.d0
enddo
else if (iflag.eq.3) then
do i=1, nnz
work(b%ia1(ip2+i-1)) = b%aspk(i) + work(b%ia1(ip2+i-1))
select case(iflag)
case(psb_dupl_ovwrt_,psb_dupl_err_)
do i=1, nnz
work(b%ia1(ip2+i-1)) = b%aspk(i)
enddo
case(psb_dupl_add_)
do i=1, nnz
work(b%ia1(ip2+i-1)) = b%aspk(i) + work(b%ia1(ip2+i-1))
enddo
case default
info = 8887
call psb_errpush(info,name)
goto 9999
end select
do i=1, nnz+count
b%aspk(i) = work(i)
enddo
do i=1, count
b%aspk(b%ia1(ipc+i-1)) = 0.d0
end do
endif
do i=1, nnz+count
b%aspk(i) = work(i)
enddo
do i=1, count
b%aspk(b%ia1(ipc+i-1)) = 0.d0
end do
endif
case(psb_upd_dflt_,psb_upd_srch_)
! Nothing to be done
case default
! Wrong value
info = 8888
call psb_errpush(info,name)
goto 9999
end select
end if
b%infoa(psb_state_) = psb_spmat_asb_
call psb_sp_setifld(psb_spmat_asb_,psb_state_,b,info)
call psb_erractionrestore(err_act)
return

@ -0,0 +1,14 @@
subroutine psb_getifield(val,field,info,isize,ierr)
integer :: val,field,isize,ierr
integer :: info(*)
ierr = 0
val = -1
if ((field < 1).or.(field > isize)) then
ierr = -1
return
endif
val = info(field)
return
end subroutine psb_getifield

@ -0,0 +1,14 @@
subroutine psb_setifield(val,field,info,isize,ierr)
integer :: val,field,isize,ierr
integer :: info(*)
ierr = 0
if ((field < 1).or.(field > isize)) then
ierr = -1
return
endif
info(field) = val
return
end subroutine psb_setifield

@ -49,7 +49,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
character(len=5) :: ufida
integer :: i,j,ir,ic,nr,nc, ng, nza, isza,spstate, nnz,&
& ip1, nzl, err_act, int_err(5)
& ip1, nzl, err_act, int_err(5), iupd
logical, parameter :: debug=.true.
character(len=20) :: name, ch_err
@ -59,113 +59,120 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
info = 0
if (nz <= 0) then
info = 10
int_err(1)=1
call psb_errpush(info,name,i_err=int_err)
goto 9999
info = 10
int_err(1)=1
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
if (size(ia) < nz) then
info = 35
int_err(1)=2
call psb_errpush(info,name,i_err=int_err)
goto 9999
info = 35
int_err(1)=2
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
if (size(ja) < nz) then
info = 35
int_err(1)=3
call psb_errpush(info,name,i_err=int_err)
goto 9999
info = 35
int_err(1)=3
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
if (size(val) < nz) then
info = 35
int_err(1)=4
call psb_errpush(info,name,i_err=int_err)
goto 9999
info = 35
int_err(1)=4
call psb_errpush(info,name,i_err=int_err)
goto 9999
end if
!!$ ufida = toupper(a%fida)
call touppers(a%fida,ufida)
ng = size(gtl)
spstate = a%infoa(psb_state_)
spstate = psb_sp_getifld(psb_state_,a,info)
select case(spstate)
case(psb_spmat_bld_)
if ((ufida /= 'COO').and.(ufida/='COI')) then
info = 134
ch_err(1:3)=ufida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_spinfo(psb_nztotreq_,a,nza,info)
call psb_spinfo(psb_nzsizereq_,a,isza,info)
if(info.ne.izero) then
info=4010
ch_err='psb_spinfo'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if ((nza+nz)>isza) then
call psb_sp_reall(a,nza+nz,info)
if(info.ne.izero) then
info=4010
ch_err='psb_sp_reall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
endif
call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,&
& imin,imax,jmin,jmax,info)
if(info.ne.izero) then
if ((ufida /= 'COO').and.(ufida/='COI')) then
info = 134
ch_err(1:3)=ufida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_spinfo(psb_nztotreq_,a,nza,info)
call psb_spinfo(psb_nzsizereq_,a,isza,info)
if(info.ne.izero) then
info=4010
ch_err='psb_spinfo'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if ((nza+nz)>isza) then
call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info)
if(info.ne.izero) then
info=4010
ch_err='psb_inner_ins'
ch_err='psb_sp_reall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debug) then
if ((nza - a%infoa(psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: insert discarded items '
end if
end if
if ((nza - a%infoa(psb_nnz_)) /= nz) then
a%infoa(psb_del_bnd_) = nza
endif
a%infoa(psb_nnz_) = nza
endif
endif
call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,&
& imin,imax,jmin,jmax,info)
if(info.ne.izero) then
info=4010
ch_err='psb_inner_ins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debug) then
if ((nza - a%infoa(psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: insert discarded items '
end if
end if
if ((nza - a%infoa(psb_nnz_)) /= nz) then
a%infoa(psb_del_bnd_) = nza
endif
a%infoa(psb_nnz_) = nza
case(psb_spmat_upd_)
if (ibits(a%infoa(psb_upd_),2,1).eq.1) then
ip1 = a%infoa(psb_upd_pnt_)
nza = a%ia2(ip1+psb_nnz_)
nzl = a%infoa(psb_del_bnd_)
call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,&
& imin,imax,jmin,jmax,nzl,info)
if(info.ne.izero) then
info=4010
ch_err='psb_inner_upd'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
iupd = psb_sp_getifld(psb_upd_,a,info)
select case (iupd)
case (psb_upd_perm_)
ip1 = psb_sp_getifld(psb_upd_pnt_,a,info)
nzl = psb_sp_getifld(psb_del_bnd_,a,info)
nza = a%ia2(ip1+psb_nnz_)
call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,&
& imin,imax,jmin,jmax,nzl,info)
if(info.ne.izero) then
info=4010
ch_err='psb_inner_upd'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debug) then
if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: update discarded items '
end if
end if
a%ia2(ip1+psb_nnz_) = nza
else
info = 2231
call psb_errpush(info,name)
goto 9999
endif
a%ia2(ip1+psb_nnz_) = nza
case (psb_upd_dflt_, psb_upd_srch_)
write(0,*) 'Default & search inner update to be implemented'
info = 2230
call psb_errpush(info,name)
goto 9999
case default
info = 2231
call psb_errpush(info,name)
goto 9999
end select
case default
info = 2232
call psb_errpush(info,name)
goto 9999
info = 2232
call psb_errpush(info,name)
goto 9999
end select
return
@ -175,8 +182,8 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error()
return
call psb_error()
return
end if
return
@ -197,7 +204,7 @@ contains
if (nza >= nzl) then
do i=1, nz
nza = nza + 1
a%aspk(nza) = val(i)
aspk(nza) = val(i)
end do
else
do i=1, nz
@ -208,7 +215,7 @@ contains
ic = gtl(ic)
if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then
nza = nza + 1
a%aspk(nza) = val(i)
aspk(nza) = val(i)
end if
end if
end do
@ -239,9 +246,9 @@ contains
ic = gtl(ic)
if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then
nza = nza + 1
a%ia1(nza) = ir
a%ia2(nza) = ic
a%aspk(nza) = val(i)
ia1(nza) = ir
ia2(nza) = ic
aspk(nza) = val(i)
end if
end if
end do

@ -47,6 +47,7 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd)
use psb_const_mod
use psb_error_mod
use psb_spmat_type
use psb_string_mod
implicit none
!....Parameters...
@ -70,12 +71,12 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd)
character(len=20) :: name, ch_err
interface psb_cest
subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info)
integer, intent(in) :: m,n,nnz
integer, intent(out) :: lia1, lia2, lar, info
character, intent(inout) :: afmt*5
character, intent(in) :: up
end subroutine psb_cest
subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info)
integer, intent(in) :: m,n,nnz
integer, intent(out) :: lia1, lia2, lar, info
character, intent(inout) :: afmt*5
character, intent(in) :: up
end subroutine psb_cest
end interface
interface psb_spinfo
@ -99,17 +100,17 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd)
ifc_ = 1
endif
if (present(check)) then
check_ = check
check_ = toupper(check)
else
check_ = 'N'
endif
if (present(trans)) then
trans_ = trans
trans_ = toupper(trans )
else
trans_ = 'N'
endif
if (present(unitd)) then
unitd_ = unitd
unitd_ = toupper(unitd )
else
unitd_ = 'U'
endif
@ -133,7 +134,7 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd)
endif
if((check_=='Y').or.(check_=='C')) then
if(a%fida(1:3)=='CSR') then
if(toupper(a%fida(1:3))=='CSR') then
call zcsrck(trans,a%m,a%k,a%descra,a%aspk,a%ia1,a%ia2,work,size(work),info)
if(info /= 0) then
info=4010
@ -154,13 +155,14 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd)
b%k=a%k
call psb_spinfo(psb_nztotreq_,a,size_req,info)
if (debug) write(0,*) 'DCSDP : size_req 1:',size_req
!! PULL IUP FROM INFOA FIELD
iup = iand(b%infoa(psb_upd_),4)
if (iup > 0) then
!
iup = psb_sp_getifld(psb_upd_,b,info)
if (iup == psb_upd_perm_) then
up = 'Y'
else
up = 'N'
endif
n_row=b%m
n_col=b%k
call psb_cest(b%fida, n_row,n_col,size_req,&
@ -191,11 +193,11 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd)
b%pr(:) = 0
select case (a%fida(1:3))
select case (toupper(a%fida(1:3)))
case ('CSR')
select case (b%fida(1:3))
select case (toupper(b%fida(1:3)))
case ('CSR')
@ -271,7 +273,7 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd)
case ('COO','COI')
select case (b%fida(1:3))
select case (toupper(b%fida(1:3)))
case ('CSR')
@ -361,99 +363,116 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd)
else if (check_=='R') then
!...Regenerating matrix
if (b%infoa(psb_state_) /= psb_spmat_upd_) then
info = 8888
call psb_errpush(info,name)
goto 9999
endif
if (ibits(b%infoa(psb_upd_),2,1).eq.0) then
!
! Nothing to be done......
!
if (psb_sp_getifld(psb_state_,b,info) /= psb_spmat_upd_) then
info = 8888
call psb_errpush(info,name)
goto 9999
endif
select case(psb_sp_getifld(psb_upd_,b,info))
if (b%fida(1:3)/='JAD') then
ip1 = b%infoa(psb_upd_pnt_)
ip2 = b%ia2(ip1+psb_ip2_)
nnz = b%ia2(ip1+psb_nnz_)
iflag = b%ia2(ip1+psb_iflag_)
ichk = b%ia2(ip1+psb_ichk_)
nnzt = b%ia2(ip1+psb_nnzt_)
if (debug) write(*,*) 'Regeneration start: ',&
& b%infoa(psb_upd_),psb_perm_update_,nnz,nnzt ,iflag,info
if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then
info = 8889
write(*,*) 'Regeneration start error: ',&
& b%infoa(psb_upd_),psb_perm_update_,nnz,nnzt ,iflag,ichk
call psb_errpush(info,name)
goto 9999
endif
do i= 1, nnz
work(i) = 0.d0
enddo
if (iflag.eq.2) then
do i=1, nnz
work(b%ia2(ip2+i-1)) = b%aspk(i)
case(psb_upd_perm_)
if (toupper(b%fida(1:3))/='JAD') then
ip1 = psb_sp_getifld(psb_upd_pnt_,b,info)
ip2 = b%ia2(ip1+psb_ip2_)
nnz = b%ia2(ip1+psb_nnz_)
iflag = b%ia2(ip1+psb_iflag_)
ichk = b%ia2(ip1+psb_ichk_)
nnzt = b%ia2(ip1+psb_nnzt_)
if (debug) write(*,*) 'Regeneration start: ',&
& b%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,info
if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then
info = 8889
write(*,*) 'Regeneration start error: ',&
& b%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk
call psb_errpush(info,name)
goto 9999
endif
do i= 1, nnz
work(i) = 0.d0
enddo
else if (iflag.eq.3) then
select case(iflag)
case(psb_dupl_ovwrt_,psb_dupl_err_)
do i=1, nnz
work(b%ia2(ip2+i-1)) = b%aspk(i)
enddo
case(psb_dupl_add_)
do i=1, nnz
work(b%ia2(ip2+i-1)) = b%aspk(i) + work(b%ia2(ip2+i-1))
enddo
case default
info = 8887
call psb_errpush(info,name)
goto 9999
end select
do i=1, nnz
work(b%ia2(ip2+i-1)) = b%aspk(i) + work(b%ia2(ip2+i-1))
b%aspk(i) = work(i)
enddo
endif
do i=1, nnz
b%aspk(i) = work(i)
enddo
else if (b%fida(1:3) == 'JAD') then
ip1 = b%infoa(psb_upd_pnt_)
ip2 = b%ia1(ip1+psb_ip2_)
count = b%ia1(ip1+psb_zero_)
ipc = b%ia1(ip1+psb_ipc_)
nnz = b%ia1(ip1+psb_nnz_)
iflag = b%ia1(ip1+psb_iflag_)
ichk = b%ia1(ip1+psb_ichk_)
nnzt = b%ia1(ip1+psb_nnzt_)
if (debug) write(*,*) 'Regeneration start: ',&
& b%infoa(psb_upd_),psb_perm_update_,nnz,nnzt,count, &
& iflag,info
if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then
info = 10
write(*,*) 'Regeneration start error: ',&
& b%infoa(psb_upd_),psb_perm_update_,nnz,nnzt ,iflag,ichk
call psb_errpush(info,name)
goto 9999
endif
do i= 1, nnz+count
work(i) = 0.d0
enddo
if (iflag.eq.2) then
do i=1, nnz
work(b%ia1(ip2+i-1)) = b%aspk(i)
else if (toupper(b%fida(1:3)) == 'JAD') then
ip1 = psb_sp_getifld(psb_upd_pnt_,b,info)
ip2 = b%ia1(ip1+psb_ip2_)
count = b%ia1(ip1+psb_zero_)
ipc = b%ia1(ip1+psb_ipc_)
nnz = b%ia1(ip1+psb_nnz_)
iflag = b%ia1(ip1+psb_iflag_)
ichk = b%ia1(ip1+psb_ichk_)
nnzt = b%ia1(ip1+psb_nnzt_)
if (debug) write(*,*) 'Regeneration start: ',&
& b%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt,count, &
& iflag,info
if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then
info = 10
write(*,*) 'Regeneration start error: ',&
& b%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk
call psb_errpush(info,name)
goto 9999
endif
do i= 1, nnz+count
work(i) = 0.d0
enddo
else if (iflag.eq.3) then
do i=1, nnz
work(b%ia1(ip2+i-1)) = b%aspk(i) + work(b%ia1(ip2+i-1))
select case(iflag)
case(psb_dupl_ovwrt_,psb_dupl_err_)
do i=1, nnz
work(b%ia1(ip2+i-1)) = b%aspk(i)
enddo
case(psb_dupl_add_)
do i=1, nnz
work(b%ia1(ip2+i-1)) = b%aspk(i) + work(b%ia1(ip2+i-1))
enddo
case default
info = 8887
call psb_errpush(info,name)
goto 9999
end select
do i=1, nnz+count
b%aspk(i) = work(i)
enddo
do i=1, count
b%aspk(b%ia1(ipc+i-1)) = 0.d0
end do
endif
do i=1, nnz+count
b%aspk(i) = work(i)
enddo
do i=1, count
b%aspk(b%ia1(ipc+i-1)) = 0.d0
end do
endif
case(psb_upd_dflt_,psb_upd_srch_)
! Nothing to be done
case default
! Wrong value
info = 8888
call psb_errpush(info,name)
goto 9999
end select
end if
b%infoa(psb_state_) = psb_spmat_asb_
call psb_sp_setifld(psb_spmat_asb_,psb_state_,b,info)
call psb_erractionrestore(err_act)
return

@ -45,7 +45,7 @@
! iblck - integer(optional). First row of submatrix belonging to blck to be inserted.
! jblck - integer(optional). First col of submatrix belonging to blck to be inserted.
subroutine psb_dins(m, n, x, ix, jx, blck, desc_a, info,&
& iblck, jblck)
& iblck, jblck,dupl)
!....insert dense submatrix to dense matrix .....
use psb_descriptor_type
use psb_const_mod
@ -61,12 +61,13 @@ subroutine psb_dins(m, n, x, ix, jx, blck, desc_a, info,&
real(kind(1.d0)), intent(in) :: blck(:,:)
integer,intent(out) :: info
integer, optional, intent(in) :: iblck,jblck
integer, optional, intent(in) :: dupl
!locals.....
integer :: icontxt,i,loc_row,glob_row,row,k,err_act,&
& nprocs,mode, loc_cols,col,iblock, jblock, mglob, int_err(5), err
integer :: nprow,npcol, me ,mypcol
integer :: nprow,npcol, me ,mypcol,dupl_
character :: temp_descra*11,temp_fida*5
character(len=20) :: name, char_err
@ -161,22 +162,52 @@ subroutine psb_dins(m, n, x, ix, jx, blck, desc_a, info,&
else
jblock = 1
endif
if (present(dupl)) then
dupl_ = dupl
else
dupl_ = psb_dupl_ovwrt_
endif
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
do col = 1, n
x(loc_row,jx+col-1) = x(loc_row,jx+col-1) + blck(iblock+i-1,jblock+col-1)
enddo
end if
enddo
select case(dupl_)
case(psb_dupl_ovwrt_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
do col = 1, n
x(loc_row,jx+col-1) = blck(iblock+i-1,jblock+col-1)
enddo
end if
enddo
case(psb_dupl_add_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
do col = 1, n
x(loc_row,jx+col-1) = x(loc_row,jx+col-1) + blck(iblock+i-1,jblock+col-1)
enddo
end if
enddo
case default
info = 321
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
@ -240,7 +271,7 @@ end subroutine psb_dins
! info - integer. Eventually returns an error code
! iblck - integer(optional). First row of submatrix belonging to blck to be inserted.
subroutine psb_dinsvm(m, x, ix, jx, blck, desc_a,info,&
& iblck)
& iblck,dupl)
!....insert dense submatrix to dense matrix .....
use psb_descriptor_type
use psb_const_mod
@ -265,10 +296,11 @@ subroutine psb_dinsvm(m, x, ix, jx, blck, desc_a,info,&
real(kind(1.d0)), intent(in) :: blck(:)
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: dupl
!locals.....
integer :: icontxt,i,loc_row,glob_row,loc_cols,mglob,err_act, int_err(5),err
integer :: nprow,npcol, me ,mypcol, iblock
integer :: nprow,npcol, me ,mypcol, iblock, dupl_
character(len=20) :: name, char_err
if(psb_get_errstatus().ne.0) return
@ -352,21 +384,48 @@ subroutine psb_dinsvm(m, x, ix, jx, blck, desc_a,info,&
iblock = 1
endif
do i = 1, m
!loop over all blck's rows
if (present(dupl)) then
dupl_ = dupl
else
dupl_ = psb_dupl_ovwrt_
endif
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
select case(dupl_)
case(psb_dupl_ovwrt_)
do i = 1, m
!loop over all blck's rows
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row,jx) = x(loc_row,jx) + blck(iblock+i-1)
end if
enddo
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row,jx) = blck(iblock+i-1)
end if
enddo
case(psb_dupl_add_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row,jx) = x(loc_row,jx) + blck(iblock+i-1)
end if
enddo
case default
info = 321
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
@ -427,7 +486,7 @@ end subroutine psb_dinsvm
! iblck - integer(optional). First row of submatrix belonging to blck to be inserted.
! insflag - integer(optional). ???
subroutine psb_dinsvv(m, x, ix, blck, desc_a, info,&
& iblck,insflag)
& iblck,insflag,dupl)
!....insert dense submatrix to dense matrix .....
use psb_descriptor_type
use psb_const_mod
@ -450,11 +509,12 @@ subroutine psb_dinsvv(m, x, ix, blck, desc_a, info,&
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: insflag
integer, optional, intent(in) :: dupl
!locals.....
integer :: icontxt,i,loc_row,glob_row,row,k,&
& loc_rows,loc_cols,iblock, liflag,mglob,err_act, int_err(5), err
integer :: nprow,npcol, me ,mypcol
integer :: nprow,npcol, me ,mypcol,dupl_
character(len=20) :: name, char_err
if(psb_get_errstatus().ne.0) return
@ -528,32 +588,70 @@ subroutine psb_dinsvv(m, x, ix, blck, desc_a, info,&
else
liflag = psb_upd_glbnum_
end if
if (liflag == psb_upd_glbnum_) then
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row) = x(loc_row) + blck(iblock+i-1)
end if
enddo
else if (liflag == psb_upd_locnum_) then
k = min(ix+m-1,loc_rows)
do i=ix,k
x(i) = x(i) + blck(i-ix+1)
enddo
if (present(dupl)) then
dupl_ = dupl
else
info=-1
dupl_ = psb_dupl_ovwrt_
endif
select case(dupl_)
case(psb_dupl_ovwrt_)
if (liflag == psb_upd_glbnum_) then
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row) = blck(iblock+i-1)
end if
enddo
else if (liflag == psb_upd_locnum_) then
k = min(ix+m-1,loc_rows)
do i=ix,k
x(i) = blck(i-ix+1)
enddo
else
info=-1
call psb_errpush(info,name)
goto 9999
endif
case(psb_dupl_add_)
if (liflag == psb_upd_glbnum_) then
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row) = x(loc_row) + blck(iblock+i-1)
end if
enddo
else if (liflag == psb_upd_locnum_) then
k = min(ix+m-1,loc_rows)
do i=ix,k
x(i) = x(i) + blck(i-ix+1)
enddo
else
info=-1
call psb_errpush(info,name)
goto 9999
endif
case default
info = 321
call psb_errpush(info,name)
goto 9999
endif
end select
call psb_erractionrestore(err_act)
return

@ -42,7 +42,7 @@
! up - character(optional). ???
! dup - integer(optional). ???
!
subroutine psb_dspasb(a,desc_a, info, afmt, up, dup)
subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl)
use psb_descriptor_type
use psb_spmat_type
@ -50,33 +50,25 @@ subroutine psb_dspasb(a,desc_a, info, afmt, up, dup)
use psb_const_mod
use psi_mod
use psb_error_mod
use psb_string_mod
implicit none
interface psb_cest
subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info)
integer, intent(in) :: m,n,nnz
integer, intent(out) :: lia1, lia2, lar, info
character, intent(inout) :: afmt*5
character, intent(in) :: up
end subroutine psb_cest
end interface
!...Parameters....
type(psb_dspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
integer,optional, intent(in) :: dup
character, optional, intent(in) :: afmt*5, up
integer,optional, intent(in) :: dupl, upd
character, optional, intent(in) :: afmt*5
!....Locals....
integer :: int_err(5)
type(psb_dspmat_type) :: atemp
real(kind(1.d0)) :: real_err(5)
integer :: ia1_size,ia2_size,aspk_size,m,i,err,&
& nprow,npcol,myrow,mycol ,size_req,idup,n_col,iout, err_act
integer :: dscstate, spstate, nr,k,j, iupdup
& nprow,npcol,myrow,mycol ,size_req,n_col,iout, err_act
integer :: dscstate, spstate, nr,k,j
integer :: upd_, dupl_
integer :: icontxt,temp(2),isize(2),n_row
character :: iup
logical, parameter :: debug=.false., debugwrt=.false.
character(len=20) :: name, ch_err
@ -93,21 +85,21 @@ subroutine psb_dspasb(a,desc_a, info, afmt, up, dup)
! check on BLACS grid
call blacs_gridinfo(icontxt, nprow, npcol, myrow, mycol)
if (nprow.eq.-1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol /= 1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name)
goto 9999
endif
if (.not.psb_is_asb_dec(dscstate)) then
info = 600
int_err(1) = dscstate
call psb_errpush(info,name)
goto 9999
info = 600
int_err(1) = dscstate
call psb_errpush(info,name)
goto 9999
endif
if (debug) Write (*, *) ' Begin matrix assembly...'
@ -116,165 +108,128 @@ subroutine psb_dspasb(a,desc_a, info, afmt, up, dup)
spstate = a%infoa(psb_state_)
if (spstate == psb_spmat_bld_) then
!
! First case: we come from a fresh build.
!
n_row = desc_a%matrix_data(psb_n_row_)
n_col = desc_a%matrix_data(psb_n_col_)
!
! Second step: handle the local matrix part.
!
iupdup = 0
if (present(up)) then
if(up.eq.'Y') then
iupdup = 4
iup = up
else if (up /= 'N') then
write(0,*)'Wrong value for update input in ASB...'
write(0,*)'Changing to default'
iup = 'N'
else
iup = 'N'
endif
else
iup = 'N'
endif
if (present(dup)) then
if((dup.lt.1).or.(dup.gt.3)) then
write(0,*)'Wrong value for duplicate input in ASB...'
write(0,*)'Changing to default'
idup = 1
else
idup = dup
endif
else
idup = 1
endif
iupdup = ieor(iupdup,idup)
a%infoa(psb_upd_)=iupdup
if (debug) write(0,*)'in ASB',psb_upd_,iupdup
a%m = n_row
a%k = n_col
call psb_sp_clone(a,atemp,info)
if(info /= no_err) then
info=4010
ch_err='psb_sp_clone'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
! convert to user requested format after the temp copy
end if
if (present(afmt)) then
a%fida = afmt
else
a%fida = '???'
endif
!
! work area requested must be fixed to
! No of Grid'd processes and NNZ+2
!
!!$ size_req = max(a%infoa(psb_nnz_),1)+3
!!$ if (debug) write(0,*) 'DCSDP : size_req 1:',size_req
!!$ call psb_cest(a%fida, n_row,n_col,size_req,&
!!$ & ia1_size, ia2_size, aspk_size, iup,info)
!!$ write(0,*) 'ESTIMATE : ',ia1_size,ia2_size,aspk_Size,iup
!!$ if (info /= no_err) then
!!$ info=4010
!!$ ch_err='psb_cest'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$ endif
!!$
!!$ call psb_sp_reall(a,ia1_size,ia2_size,aspk_size,info)
!!$ if (info /= no_err) then
!!$ info=4010
!!$ ch_err='psb_sp_reall'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$ endif
!!$
!!$ a%pl(:) = 0
!!$ a%pr(:) = 0
if (debugwrt) then
iout = 30+myrow
open(iout)
call psb_csprt(iout,atemp,head='Input mat')
close(iout)
endif
! Do the real conversion into the requested storage formatmode
! result is put in A
call psb_csdp(atemp,a,info,ifc=2)
IF (debug) WRITE (*, *) myrow,' ASB: From DCSDP',info,' ',A%FIDA
if (info /= no_err) then
info=4010
ch_err='psb_csdp'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debugwrt) then
iout = 60+myrow
open(iout)
call psb_csprt(iout,a,head='Output mat')
close(iout)
endif
call psb_sp_free(atemp,info)
!
! First case: we come from a fresh build.
!
n_row = desc_a%matrix_data(psb_n_row_)
n_col = desc_a%matrix_data(psb_n_col_)
!
! Second step: handle the local matrix part.
!
if (present(upd)) then
upd_=upd
else
upd_ = psb_upd_dflt_
endif
if (present(dupl)) then
if((dupl < psb_dupl_ovwrt_).or.(dupl > psb_dupl_err_)) then
write(0,*)'Wrong value for duplicate input in ASB...'
write(0,*)'Changing to default'
dupl_ = psb_dupl_def_
else
dupl_ = dupl
endif
else
dupl_ = psb_dupl_def_
endif
call psb_sp_setifld(upd_,psb_upd_,a,info)
call psb_sp_setifld(dupl_,psb_dupl_,a,info)
a%m = n_row
a%k = n_col
call psb_sp_clone(a,atemp,info)
if(info /= no_err) then
info=4010
ch_err='psb_sp_clone'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
! convert to user requested format after the temp copy
end if
if (present(afmt)) then
a%fida = afmt
else
a%fida = '???'
endif
if (debugwrt) then
iout = 30+myrow
open(iout)
call psb_csprt(iout,atemp,head='Input mat')
close(iout)
endif
! Do the real conversion into the requested storage format
! result is put in A
call psb_csdp(atemp,a,info,ifc=2)
IF (debug) WRITE (*, *) myrow,' ASB: From DCSDP',info,' ',A%FIDA
if (info /= no_err) then
info=4010
ch_err='psb_csdp'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debugwrt) then
iout = 60+myrow
open(iout)
call psb_csprt(iout,a,head='Output mat')
close(iout)
endif
call psb_sp_free(atemp,info)
else if (spstate == psb_spmat_upd_) then
!
! Second case: we come from an update loop.
!
! Right now, almost nothing to be done, but this
! may change in the future
! as we revise the implementation of the update routine.
call psb_sp_all(atemp,1,info)
atemp%m=a%m
atemp%k=a%k
! check on allocation
if (info /= no_err) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
call psb_csdp(atemp,a,info,check='R')
! check on error retuned by dcsdp
if (info /= no_err) then
info = 4010
ch_err='psb_csdp90'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_sp_free(atemp,info)
if (info /= no_err) then
info = 4010
ch_err='sp_free'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
!
! Second case: we come from an update loop.
!
! Right now, almost nothing to be done, but this
! may change in the future
! as we revise the implementation of the update routine.
call psb_sp_all(atemp,1,info)
atemp%m=a%m
atemp%k=a%k
! check on allocation
if (info /= no_err) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
call psb_csdp(atemp,a,info,check='R')
! check on error retuned by dcsdp
if (info /= no_err) then
info = 4010
ch_err='psb_csdp90'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_sp_free(atemp,info)
if (info /= no_err) then
info = 4010
ch_err='sp_free'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
else
info = 600
call psb_errpush(info,name)
goto 9999
if (debug) write(0,*) 'Sparse matrix state:',spstate,psb_spmat_bld_,psb_spmat_upd_
info = 600
call psb_errpush(info,name)
goto 9999
if (debug) write(0,*) 'Sparse matrix state:',spstate,psb_spmat_bld_,psb_spmat_upd_
endif
@ -284,8 +239,8 @@ subroutine psb_dspasb(a,desc_a, info, afmt, up, dup)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(icontxt)
return
call psb_error(icontxt)
return
end if
return

@ -73,40 +73,43 @@ Subroutine psb_dsprn(a, desc_a,info)
! ....verify blacs grid correctness..
if (npcol.ne.1) then
info = 2030
call psb_errpush(info,name)
goto 9999
info = 2030
call psb_errpush(info,name)
goto 9999
endif
if (debug) write(*,*) 'got through igamx2d '
if (.not.psb_is_asb_dec(desc_a%matrix_data(psb_dec_type_))) then
info=590
call psb_errpush(info,name)
goto 9999
info=590
call psb_errpush(info,name)
goto 9999
endif
if (a%infoa(psb_state_) == psb_spmat_asb_) then
a%aspk(:) = dzero
if (ibits(a%infoa(psb_upd_),2,1)==1) then
if(a%fida(1:3).eq.'JAD') then
a%ia1(a%infoa(psb_upd_pnt_)+psb_nnz_) = 0
else
a%ia2(a%infoa(psb_upd_pnt_)+psb_nnz_) = 0
endif
endif
a%infoa(psb_state_) = psb_spmat_upd_
else if (a%infoa(psb_state_) == psb_spmat_bld_) then
! in this case do nothing. this allows sprn to be called
! right after allocate, with spins doing the right thing.
! hopefully :-)
else if (a%infoa(psb_state_) == psb_spmat_upd_) then
else
info=591
call psb_errpush(info,name)
endif
select case(psb_sp_getifld(psb_state_,a,info))
case(psb_spmat_asb_)
a%aspk(:) = dzero
if (psb_sp_getifld(psb_upd_,a,info)==psb_upd_perm_) then
if(a%fida(1:3).eq.'JAD') then
a%ia1(a%infoa(psb_upd_pnt_)+psb_nnz_) = 0
else
a%ia2(a%infoa(psb_upd_pnt_)+psb_nnz_) = 0
endif
endif
a%infoa(psb_state_) = psb_spmat_upd_
case(psb_spmat_bld_)
! in this case do nothing. this allows sprn to be called
! right after allocate, with spins doing the right thing.
! hopefully :-)
case( psb_spmat_upd_)
case default
info=591
call psb_errpush(info,name)
end select
if (info /= 0) goto 9999

@ -45,7 +45,7 @@
! iblck - integer(optional). First row of submatrix belonging to blck to be inserted.
! jblck - integer(optional). First col of submatrix belonging to blck to be inserted.
subroutine psb_iins(m, n, x, ix, jx, blck, desc_a, info,&
& iblck, jblck)
& iblck, jblck,dupl)
!....insert dense submatrix to dense matrix .....
use psb_descriptor_type
use psb_const_mod
@ -60,11 +60,12 @@ subroutine psb_iins(m, n, x, ix, jx, blck, desc_a, info,&
integer, intent(in) :: blck(:,:)
integer, intent(out) :: info
integer, optional, intent(in) :: iblck,jblck
integer, optional, intent(in) :: dupl
!locals.....
integer :: icontxt,i,loc_row,glob_row,&
& loc_cols,col,iblock, jblock, mglob
& loc_cols,col,iblock, jblock, mglob,dupl_
integer :: nprow,npcol, myrow ,mycol, int_err(5),err_act
character(len=20) :: name, ch_err
@ -75,9 +76,9 @@ subroutine psb_iins(m, n, x, ix, jx, blck, desc_a, info,&
if ((.not.associated(desc_a%matrix_data))) then
info=3110
call psb_errpush(info, name)
return
info=3110
call psb_errpush(info, name)
return
end if
icontxt=desc_a%matrix_data(psb_ctxt_)
@ -85,98 +86,126 @@ subroutine psb_iins(m, n, x, ix, jx, blck, desc_a, info,&
! check on blacs grid
call blacs_gridinfo(icontxt, nprow, npcol, myrow, mycol)
if (nprow.eq.-1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol.ne.1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
endif
if (.not.associated(desc_a%glob_to_loc)) then
info=3110
call psb_errpush(info,name,int_err)
goto 9999
info=3110
call psb_errpush(info,name,int_err)
goto 9999
end if
!... check parameters....
if (m.lt.0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
else if (n.lt.0) then
info = 10
int_err(1) = 2
int_err(2) = n
call psb_errpush(info,name,int_err)
goto 9999
info = 10
int_err(1) = 2
int_err(2) = n
call psb_errpush(info,name,int_err)
goto 9999
else if (ix.lt.1) then
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
else if (jx.lt.1) then
info = 20
int_err(1) = 7
int_err(2) = jx
call psb_errpush(info,name,int_err)
goto 9999
info = 20
int_err(1) = 7
int_err(2) = jx
call psb_errpush(info,name,int_err)
goto 9999
else if (.not.psb_is_ok_dec(desc_a%matrix_data(psb_dec_type_))) then
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
else if (size(x, dim=1).lt.desc_a%matrix_data(psb_n_row_)) then
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
else if (size(x, dim=2).lt.n) then
! check if dimension of x is greater than dimension of submatrix
! to insert
info = 320
int_err(1) = 2
int_err(2) = size(x, dim=2)
int_err(3) = n
call psb_errpush(info,name,int_err)
goto 9999
! check if dimension of x is greater than dimension of submatrix
! to insert
info = 320
int_err(1) = 2
int_err(2) = size(x, dim=2)
int_err(3) = n
call psb_errpush(info,name,int_err)
goto 9999
endif
loc_cols = desc_a%matrix_data(psb_n_col_)
mglob = desc_a%matrix_data(psb_m_)
if (present(iblck)) then
iblock = iblck
iblock = iblck
else
iblock = 1
iblock = 1
endif
if (present(jblck)) then
jblock = jblck
jblock = jblck
else
jblock = 1
jblock = 1
endif
if (present(dupl)) then
dupl_ = dupl
else
dupl_ = psb_dupl_ovwrt_
endif
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
do col = 1, n
x(loc_row,jx+col-1) = x(loc_row,jx+col-1) + blck(iblock+i-1,jblock+col-1)
enddo
end if
enddo
select case(dupl_)
case(psb_dupl_ovwrt_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
do col = 1, n
x(loc_row,jx+col-1) = blck(iblock+i-1,jblock+col-1)
enddo
end if
enddo
case(psb_dupl_add_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
do col = 1, n
x(loc_row,jx+col-1) = x(loc_row,jx+col-1) + blck(iblock+i-1,jblock+col-1)
enddo
end if
enddo
case default
info = 321
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
@ -184,8 +213,8 @@ subroutine psb_iins(m, n, x, ix, jx, blck, desc_a, info,&
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(icontxt)
return
call psb_error(icontxt)
return
end if
return
@ -238,7 +267,7 @@ end subroutine psb_iins
! info - integer. Eventually returns an error code
! iblck - integer(optional). First row of submatrix belonging to blck to be inserted.
subroutine psb_iinsvm(m, x, ix, jx, blck, desc_a, info,&
& iblck)
& iblck,dupl)
!....insert dense submatrix to dense matrix .....
use psb_descriptor_type
use psb_const_mod
@ -263,11 +292,12 @@ subroutine psb_iinsvm(m, x, ix, jx, blck, desc_a, info,&
integer, intent(in) :: blck(:)
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: dupl
!locals.....
integer :: icontxt,i,loc_row,glob_row,&
& loc_cols,iblock, jblock,mglob, err_act, int_err(5)
integer :: nprow,npcol, myrow ,mycol
integer :: nprow,npcol, myrow ,mycol,dupl_
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
@ -280,25 +310,53 @@ subroutine psb_iinsvm(m, x, ix, jx, blck, desc_a, info,&
mglob = desc_a%matrix_data(psb_m_)
if (present(iblck)) then
iblock = iblck
iblock = iblck
else
iblock = 1
iblock = 1
endif
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
if (present(dupl)) then
dupl_ = dupl
else
dupl_ = psb_dupl_ovwrt_
endif
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row,jx) = x(loc_row,jx) + blck(iblock+i-1)
end if
enddo
select case(dupl_)
case(psb_dupl_ovwrt_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row,jx) = blck(iblock+i-1)
end if
enddo
case(psb_dupl_add_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row,jx) = x(loc_row,jx) + blck(iblock+i-1)
end if
enddo
case default
info = 321
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
@ -306,8 +364,8 @@ subroutine psb_iinsvm(m, x, ix, jx, blck, desc_a, info,&
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(icontxt)
return
call psb_error(icontxt)
return
end if
return
@ -356,7 +414,7 @@ end subroutine psb_iinsvm
! info - integer. Eventually returns an error code
! iblck - integer(optional). First row of submatrix belonging to blck to be inserted.
subroutine psb_iinsvv(m, x, ix, blck, desc_a, info,&
& iblck)
& iblck,dupl)
!....insert dense submatrix to dense matrix .....
use psb_descriptor_type
use psb_const_mod
@ -378,11 +436,12 @@ subroutine psb_iinsvv(m, x, ix, blck, desc_a, info,&
integer, intent(in) :: blck(:)
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: dupl
!locals.....
integer :: icontxt,i,loc_row,glob_row,k,&
& loc_rows,loc_cols,col,iblock, jblock, mglob, err_act, int_err(5)
integer :: nprow,npcol, myrow ,mycol
integer :: nprow,npcol, myrow ,mycol,dupl_
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
@ -391,55 +450,55 @@ subroutine psb_iinsvv(m, x, ix, blck, desc_a, info,&
call psb_erractionsave(err_act)
if ((.not.associated(desc_a%matrix_data))) then
info=3110
call psb_errpush(info,name)
return
info=3110
call psb_errpush(info,name)
return
end if
icontxt=desc_a%matrix_data(psb_ctxt_)
! check on blacs grid
call blacs_gridinfo(icontxt, nprow, npcol, myrow, mycol)
if (nprow.eq.-1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol.ne.1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
endif
if (.not.associated(desc_a%glob_to_loc)) then
info=3110
call psb_errpush(info,name)
goto 9999
info=3110
call psb_errpush(info,name)
goto 9999
end if
!... check parameters....
if (m.lt.0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
else if (ix.lt.1) then
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
else if (.not.psb_is_ok_dec(desc_a%matrix_data(psb_dec_type_))) then
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
else if (size(x, dim=1).lt.desc_a%matrix_data(psb_n_row_)) then
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
endif
loc_rows=desc_a%matrix_data(psb_n_row_)
@ -447,25 +506,53 @@ subroutine psb_iinsvv(m, x, ix, blck, desc_a, info,&
mglob = desc_a%matrix_data(psb_m_)
if (present(iblck)) then
iblock = iblck
iblock = iblck
else
iblock = 1
iblock = 1
endif
if (present(dupl)) then
dupl_ = dupl
else
dupl_ = psb_dupl_ovwrt_
endif
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
select case(dupl_)
case(psb_dupl_ovwrt_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row) = blck(iblock+i-1)
end if
enddo
case(psb_dupl_add_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row) = x(loc_row) + blck(iblock+i-1)
end if
enddo
end if
enddo
case default
info = 321
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
@ -473,8 +560,8 @@ subroutine psb_iinsvv(m, x, ix, blck, desc_a, info,&
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error(icontxt)
return
call psb_error(icontxt)
return
end if
return

@ -45,7 +45,7 @@
! iblck - integer(optional). First row of submatrix belonging to blck to be inserted.
! jblck - integer(optional). First col of submatrix belonging to blck to be inserted.
subroutine psb_zins(m, n, x, ix, jx, blck, desc_a, info,&
& iblck, jblck)
& iblck, jblck,dupl)
!....insert dense submatrix to dense matrix .....
use psb_descriptor_type
use psb_const_mod
@ -61,12 +61,13 @@ subroutine psb_zins(m, n, x, ix, jx, blck, desc_a, info,&
complex(kind(1.d0)), intent(in) :: blck(:,:)
integer,intent(out) :: info
integer, optional, intent(in) :: iblck,jblck
integer, optional, intent(in) :: dupl
!locals.....
integer :: icontxt,i,loc_row,glob_row,row,k,err_act,&
& nprocs,mode, loc_cols,col,iblock, jblock, mglob, int_err(5), err
integer :: nprow,npcol, me ,mypcol
integer :: nprow,npcol, me ,mypcol,dupl_
character :: temp_descra*11,temp_fida*5
character(len=20) :: name, char_err
@ -76,14 +77,14 @@ subroutine psb_zins(m, n, x, ix, jx, blck, desc_a, info,&
name = 'psb_zins'
if (.not.associated(desc_a%glob_to_loc)) then
info=3110
call psb_errpush(info,name)
return
info=3110
call psb_errpush(info,name)
return
end if
if ((.not.associated(desc_a%matrix_data))) then
int_err(1)=3110
call psb_errpush(info,name)
return
int_err(1)=3110
call psb_errpush(info,name)
return
end if
icontxt=desc_a%matrix_data(psb_ctxt_)
@ -91,92 +92,120 @@ subroutine psb_zins(m, n, x, ix, jx, blck, desc_a, info,&
! check on blacs grid
call blacs_gridinfo(icontxt, nprow, npcol, me, mypcol)
if (nprow.eq.-1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol.ne.1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
endif
!... check parameters....
if (m.lt.0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
else if (n.lt.0) then
info = 10
int_err(1) = 2
int_err(2) = n
call psb_errpush(info,name,int_err)
goto 9999
info = 10
int_err(1) = 2
int_err(2) = n
call psb_errpush(info,name,int_err)
goto 9999
else if (ix.lt.1) then
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
else if (jx.lt.1) then
info = 20
int_err(1) = 7
int_err(2) = jx
call psb_errpush(info,name,int_err)
goto 9999
info = 20
int_err(1) = 7
int_err(2) = jx
call psb_errpush(info,name,int_err)
goto 9999
else if (.not.psb_is_ok_dec(desc_a%matrix_data(psb_dec_type_))) then
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
else if (size(x, dim=1).lt.desc_a%matrix_data(psb_n_row_)) then
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
else if (size(x, dim=2).lt.n) then
! check if dimension of x is greater than dimension of submatrix
! to insert
info = 320
int_err(1) = 2
int_err(2) = size(x, dim=2)
int_err(3) = n
call psb_errpush(info,name,int_err)
goto 9999
! check if dimension of x is greater than dimension of submatrix
! to insert
info = 320
int_err(1) = 2
int_err(2) = size(x, dim=2)
int_err(3) = n
call psb_errpush(info,name,int_err)
goto 9999
endif
loc_cols = desc_a%matrix_data(psb_n_col_)
mglob = desc_a%matrix_data(psb_m_)
if (present(iblck)) then
iblock = iblck
iblock = iblck
else
iblock = 1
iblock = 1
endif
if (present(jblck)) then
jblock = jblck
jblock = jblck
else
jblock = 1
jblock = 1
endif
if (present(dupl)) then
dupl_ = dupl
else
dupl_ = psb_dupl_ovwrt_
endif
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
do col = 1, n
x(loc_row,jx+col-1) = x(loc_row,jx+col-1) + blck(iblock+i-1,jblock+col-1)
enddo
end if
enddo
select case(dupl_)
case(psb_dupl_ovwrt_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
do col = 1, n
x(loc_row,jx+col-1) = blck(iblock+i-1,jblock+col-1)
enddo
end if
enddo
case(psb_dupl_add_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
do col = 1, n
x(loc_row,jx+col-1) = x(loc_row,jx+col-1) + blck(iblock+i-1,jblock+col-1)
enddo
end if
enddo
case default
info = 321
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
@ -185,9 +214,9 @@ subroutine psb_zins(m, n, x, ix, jx, blck, desc_a, info,&
call psb_erractionrestore(err_act)
if (err_act.eq.act_ret) then
return
return
else
call psb_error(icontxt)
call psb_error(icontxt)
end if
return
@ -240,7 +269,7 @@ end subroutine psb_zins
! info - integer. Eventually returns an error code
! iblck - integer(optional). First row of submatrix belonging to blck to be inserted.
subroutine psb_zinsvm(m, x, ix, jx, blck, desc_a,info,&
& iblck)
& iblck,dupl)
!....insert dense submatrix to dense matrix .....
use psb_descriptor_type
use psb_const_mod
@ -265,10 +294,11 @@ subroutine psb_zinsvm(m, x, ix, jx, blck, desc_a,info,&
complex(kind(1.d0)), intent(in) :: blck(:)
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: dupl
!locals.....
integer :: icontxt,i,loc_row,glob_row,loc_cols,mglob,err_act, int_err(5),err
integer :: nprow,npcol, me ,mypcol, iblock
integer :: nprow,npcol, me ,mypcol, iblock,dupl_
character(len=20) :: name, char_err
if(psb_get_errstatus().ne.0) return
@ -277,14 +307,14 @@ subroutine psb_zinsvm(m, x, ix, jx, blck, desc_a,info,&
name = 'psb_zinsvm'
if (.not.associated(desc_a%glob_to_loc)) then
info=3110
call psb_errpush(info,name)
return
info=3110
call psb_errpush(info,name)
return
end if
if ((.not.associated(desc_a%matrix_data))) then
int_err(1)=3110
call psb_errpush(info,name)
return
int_err(1)=3110
call psb_errpush(info,name)
return
end if
icontxt=desc_a%matrix_data(psb_ctxt_)
@ -292,80 +322,107 @@ subroutine psb_zinsvm(m, x, ix, jx, blck, desc_a,info,&
! check on blacs grid
call blacs_gridinfo(icontxt, nprow, npcol, me, mypcol)
if (nprow.eq.-1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol.ne.1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
endif
!... check parameters....
if (m.lt.0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
else if (ix.lt.1) then
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
else if (jx.lt.1) then
info = 20
int_err(1) = 7
int_err(2) = jx
call psb_errpush(info,name,int_err)
goto 9999
info = 20
int_err(1) = 7
int_err(2) = jx
call psb_errpush(info,name,int_err)
goto 9999
else if (.not.psb_is_ok_dec(desc_a%matrix_data(psb_dec_type_))) then
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
else if (size(x, dim=1).lt.desc_a%matrix_data(psb_n_row_)) then
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
else if (size(x, dim=2).lt.1) then
! check if dimension of x is greater than dimension of submatrix
! to insert
info = 320
int_err(1) = 2
int_err(2) = size(x, dim=2)
int_err(3) = 1
call psb_errpush(info,name,int_err)
goto 9999
! check if dimension of x is greater than dimension of submatrix
! to insert
info = 320
int_err(1) = 2
int_err(2) = size(x, dim=2)
int_err(3) = 1
call psb_errpush(info,name,int_err)
goto 9999
endif
loc_cols=desc_a%matrix_data(psb_n_col_)
mglob = desc_a%matrix_data(psb_m_)
if (present(iblck)) then
iblock = iblck
iblock = iblck
else
iblock = 1
endif
if (present(dupl)) then
dupl_ = dupl
else
iblock = 1
dupl_ = psb_dupl_ovwrt_
endif
do i = 1, m
!loop over all blck's rows
select case(dupl_)
case(psb_dupl_ovwrt_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row,jx) = blck(iblock+i-1)
end if
enddo
case(psb_dupl_add_)
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row,jx) = x(loc_row,jx) + blck(iblock+i-1)
end if
enddo
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row,jx) = x(loc_row,jx) + blck(iblock+i-1)
end if
enddo
case default
info = 321
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
@ -374,9 +431,9 @@ subroutine psb_zinsvm(m, x, ix, jx, blck, desc_a,info,&
call psb_erractionrestore(err_act)
if (err_act.eq.act_ret) then
return
return
else
call psb_error(icontxt)
call psb_error(icontxt)
end if
return
@ -427,7 +484,7 @@ end subroutine psb_zinsvm
! iblck - integer(optional). First row of submatrix belonging to blck to be inserted.
! insflag - integer(optional). ???
subroutine psb_zinsvv(m, x, ix, blck, desc_a, info,&
& iblck,insflag)
& iblck,insflag,dupl)
!....insert dense submatrix to dense matrix .....
use psb_descriptor_type
use psb_const_mod
@ -450,11 +507,12 @@ subroutine psb_zinsvv(m, x, ix, blck, desc_a, info,&
integer, intent(out) :: info
integer, optional, intent(in) :: iblck
integer, optional, intent(in) :: insflag
integer, optional, intent(in) :: dupl
!locals.....
integer :: icontxt,i,loc_row,glob_row,row,k,&
& loc_rows,loc_cols,iblock, liflag,mglob,err_act, int_err(5), err
integer :: nprow,npcol, me ,mypcol
integer :: nprow,npcol, me ,mypcol,dupl_
character(len=20) :: name, char_err
if(psb_get_errstatus().ne.0) return
@ -463,14 +521,14 @@ subroutine psb_zinsvv(m, x, ix, blck, desc_a, info,&
name = 'psb_zinsvv'
if (.not.associated(desc_a%glob_to_loc)) then
info=3110
call psb_errpush(info,name)
return
info=3110
call psb_errpush(info,name)
return
end if
if ((.not.associated(desc_a%matrix_data))) then
int_err(1)=3110
call psb_errpush(info,name)
return
int_err(1)=3110
call psb_errpush(info,name)
return
end if
icontxt=desc_a%matrix_data(psb_ctxt_)
@ -478,40 +536,40 @@ subroutine psb_zinsvv(m, x, ix, blck, desc_a, info,&
! check on blacs grid
call blacs_gridinfo(icontxt, nprow, npcol, me, mypcol)
if (nprow.eq.-1) then
info = 2010
call psb_errpush(info,name)
goto 9999
info = 2010
call psb_errpush(info,name)
goto 9999
else if (npcol.ne.1) then
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
info = 2030
int_err(1) = npcol
call psb_errpush(info,name,int_err)
goto 9999
endif
!... check parameters....
if (m.lt.0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,int_err)
goto 9999
else if (ix.lt.1) then
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
info = 20
int_err(1) = 6
int_err(2) = ix
call psb_errpush(info,name,int_err)
goto 9999
else if (.not.psb_is_ok_dec(desc_a%matrix_data(psb_dec_type_))) then
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
info = 3110
int_err(1) = desc_a%matrix_data(psb_dec_type_)
call psb_errpush(info,name,int_err)
goto 9999
else if (size(x, dim=1).lt.desc_a%matrix_data(psb_n_row_)) then
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
info = 310
int_err(1) = 5
int_err(2) = 4
call psb_errpush(info,name,int_err)
goto 9999
endif
loc_rows=desc_a%matrix_data(psb_n_row_)
@ -519,41 +577,80 @@ subroutine psb_zinsvv(m, x, ix, blck, desc_a, info,&
mglob = desc_a%matrix_data(psb_m_)
if (present(iblck)) then
iblock = iblck
iblock = iblck
else
iblock = 1
iblock = 1
endif
if (present(insflag)) then
liflag = insflag
liflag = insflag
else
liflag = psb_upd_glbnum_
liflag = psb_upd_glbnum_
end if
if (present(dupl)) then
dupl_ = dupl
else
dupl_ = psb_dupl_ovwrt_
endif
if (liflag == psb_upd_glbnum_) then
do i = 1, m
!loop over all blck's rows
select case(dupl_)
case(psb_dupl_ovwrt_)
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
if (liflag == psb_upd_glbnum_) then
do i = 1, m
!loop over all blck's rows
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row) = x(loc_row) + blck(iblock+i-1)
end if
enddo
else if (liflag == psb_upd_locnum_) then
k = min(ix+m-1,loc_rows)
do i=ix,k
x(i) = x(i) + blck(i-ix+1)
enddo
else
info=-1
call psb_errpush(info,name)
goto 9999
endif
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row) = blck(iblock+i-1)
end if
enddo
else if (liflag == psb_upd_locnum_) then
k = min(ix+m-1,loc_rows)
do i=ix,k
x(i) = blck(i-ix+1)
enddo
else
info=-1
call psb_errpush(info,name)
goto 9999
endif
case(psb_dupl_add_)
if (liflag == psb_upd_glbnum_) then
do i = 1, m
!loop over all blck's rows
! row actual block row
glob_row=ix+i-1
if (glob_row > mglob) exit
loc_row=desc_a%glob_to_loc(glob_row)
if (loc_row.ge.1) then
! this row belongs to me
! copy i-th row of block blck in x
x(loc_row) = x(loc_row) + blck(iblock+i-1)
end if
enddo
else if (liflag == psb_upd_locnum_) then
k = min(ix+m-1,loc_rows)
do i=ix,k
x(i) = x(i) + blck(i-ix+1)
enddo
else
info=-1
call psb_errpush(info,name)
goto 9999
endif
case default
info = 321
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
@ -562,9 +659,9 @@ subroutine psb_zinsvv(m, x, ix, blck, desc_a, info,&
call psb_erractionrestore(err_act)
if (err_act.eq.act_ret) then
return
return
else
call psb_error(icontxt)
call psb_error(icontxt)
end if
return

@ -42,7 +42,7 @@
! up - character(optional). ???
! dup - integer(optional). ???
!
subroutine psb_zspasb(a,desc_a, info, afmt, up, dup)
subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl)
use psb_descriptor_type
use psb_spmat_type
@ -50,33 +50,25 @@ subroutine psb_zspasb(a,desc_a, info, afmt, up, dup)
use psb_const_mod
use psi_mod
use psb_error_mod
use psb_string_mod
implicit none
interface psb_cest
subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info)
integer, intent(in) :: m,n,nnz
integer, intent(out) :: lia1, lia2, lar, info
character, intent(inout) :: afmt*5
character, intent(in) :: up
end subroutine psb_cest
end interface
!...Parameters....
type(psb_zspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
integer,optional, intent(in) :: dup
character, optional, intent(in) :: afmt*5, up
integer,optional, intent(in) :: dupl, upd
character, optional, intent(in) :: afmt*5
!....Locals....
integer :: int_err(5)
type(psb_zspmat_type) :: atemp
real(kind(1.d0)) :: real_err(5)
integer :: ia1_size,ia2_size,aspk_size,m,i,err,&
& nprow,npcol,myrow,mycol ,size_req,idup,n_col,iout, err_act
integer :: dscstate, spstate, nr,k,j, iupdup
& nprow,npcol,myrow,mycol ,size_req,n_col,iout, err_act
integer :: dscstate, spstate, nr,k,j
integer :: upd_, dupl_
integer :: icontxt,temp(2),isize(2),n_row
character :: iup
logical, parameter :: debug=.false., debugwrt=.false.
character(len=20) :: name, ch_err
@ -116,118 +108,81 @@ subroutine psb_zspasb(a,desc_a, info, afmt, up, dup)
spstate = a%infoa(psb_state_)
if (spstate == psb_spmat_bld_) then
!
! First case: we come from a fresh build.
!
n_row = desc_a%matrix_data(psb_n_row_)
n_col = desc_a%matrix_data(psb_n_col_)
!
! Second step: handle the local matrix part.
!
iupdup = 0
if (present(up)) then
if(up.eq.'Y') then
iupdup = 4
iup = up
else if (up /= 'N') then
write(0,*)'Wrong value for update input in ASB...'
write(0,*)'Changing to default'
iup = 'N'
else
iup = 'N'
endif
else
iup = 'N'
endif
if (present(dup)) then
if((dup.lt.1).or.(dup.gt.3)) then
write(0,*)'Wrong value for duplicate input in ASB...'
write(0,*)'Changing to default'
idup = 1
else
idup = dup
endif
else
idup = 1
endif
iupdup = ieor(iupdup,idup)
a%infoa(psb_upd_)=iupdup
if (debug) write(0,*)'in ASB',psb_upd_,iupdup
a%m = n_row
a%k = n_col
call psb_sp_clone(a,atemp,info)
if(info /= no_err) then
info=4010
ch_err='psb_sp_clone'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
! convert to user requested format after the temp copy
end if
if (present(afmt)) then
a%fida = afmt
else
a%fida = '???'
endif
!
! work area requested must be fixed to
! No of Grid'd processes and NNZ+2
!
!!$ size_req = max(a%infoa(psb_nnz_),1)+3
!!$ if (debug) write(0,*) 'DCSDP : size_req 1:',size_req
!!$ call psb_cest(a%fida, n_row,n_col,size_req,&
!!$ & ia1_size, ia2_size, aspk_size, iup,info)
!!$ write(0,*) 'ESTIMATE : ',ia1_size,ia2_size,aspk_Size,iup
!!$ if (info /= no_err) then
!!$ info=4010
!!$ ch_err='psb_cest'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$ endif
!!$
!!$ call psb_sp_reall(a,ia1_size,ia2_size,aspk_size,info)
!!$ if (info /= no_err) then
!!$ info=4010
!!$ ch_err='psb_sp_reall'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$ endif
!!$
!!$ a%pl(:) = 0
!!$ a%pr(:) = 0
if (debugwrt) then
iout = 30+myrow
open(iout)
call psb_csprt(iout,atemp,head='Input mat')
close(iout)
endif
! Do the real conversion into the requested storage formatmode
! result is put in A
call psb_csdp(atemp,a,info,ifc=2)
IF (debug) WRITE (*, *) myrow,' ASB: From DCSDP',info,' ',A%FIDA
if (info /= no_err) then
info=4010
ch_err='psb_csdp'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debugwrt) then
iout = 60+myrow
open(iout)
call psb_csprt(iout,a,head='Output mat')
close(iout)
endif
!
! First case: we come from a fresh build.
!
n_row = desc_a%matrix_data(psb_n_row_)
n_col = desc_a%matrix_data(psb_n_col_)
!
! Second step: handle the local matrix part.
!
if (present(upd)) then
upd_=upd
else
upd_ = psb_upd_dflt_
endif
if (present(dupl)) then
if((dupl < psb_dupl_ovwrt_).or.(dupl > psb_dupl_err_)) then
write(0,*)'Wrong value for duplicate input in ASB...'
write(0,*)'Changing to default'
dupl_ = psb_dupl_def_
else
dupl_ = dupl
endif
else
dupl_ = psb_dupl_def_
endif
call psb_sp_setifld(upd_,psb_upd_,a,info)
call psb_sp_setifld(dupl_,psb_dupl_,a,info)
a%m = n_row
a%k = n_col
call psb_sp_clone(a,atemp,info)
if(info /= no_err) then
info=4010
ch_err='psb_sp_clone'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
! convert to user requested format after the temp copy
end if
if (present(afmt)) then
a%fida = afmt
else
a%fida = '???'
endif
if (debugwrt) then
iout = 30+myrow
open(iout)
call psb_csprt(iout,atemp,head='Input mat')
close(iout)
endif
! Do the real conversion into the requested storage format
! result is put in A
call psb_csdp(atemp,a,info,ifc=2)
IF (debug) WRITE (*, *) myrow,' ASB: From DCSDP',info,' ',A%FIDA
if (info /= no_err) then
info=4010
ch_err='psb_csdp'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
if (debugwrt) then
iout = 60+myrow
open(iout)
call psb_csprt(iout,a,head='Output mat')
close(iout)
endif
call psb_sp_free(atemp,info)

@ -420,7 +420,7 @@ contains
call blacs_barrier(icontxt,'all')
t2 = mpi_wtime()
call psb_spasb(a,desc_a,info,dup=1,afmt=afmt)
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
t3 = mpi_wtime()
if(info/=0)then
info=4010
@ -437,7 +437,7 @@ contains
else
call psb_spasb(a,desc_a,info,afmt=afmt,dup=1)
call psb_spasb(a,desc_a,info,afmt=afmt,dupl=psb_dupl_err_)
if(info/=0)then
info=4010
ch_err='psspasb'
@ -780,7 +780,7 @@ contains
call blacs_barrier(icontxt,'all')
t2 = mpi_wtime()
call psb_spasb(a,desc_a,info,dup=1,afmt=afmt)
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
t3 = mpi_wtime()
if(info/=0)then
info=4010
@ -1203,7 +1203,7 @@ contains
call blacs_barrier(icontxt,'all')
t2 = mpi_wtime()
call psb_spasb(a,desc_a,info,dup=1,afmt=afmt)
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
t3 = mpi_wtime()
if(info/=0)then
info=4010
@ -1220,7 +1220,7 @@ contains
else
call psb_spasb(a,desc_a,info,afmt=afmt,dup=1)
call psb_spasb(a,desc_a,info,afmt=afmt,dupl=psb_dupl_err_)
if(info/=0)then
info=4010
ch_err='psspasb'
@ -1563,7 +1563,7 @@ contains
call blacs_barrier(icontxt,'all')
t2 = mpi_wtime()
call psb_spasb(a,desc_a,info,dup=1,afmt=afmt)
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
t3 = mpi_wtime()
if(info/=0)then
info=4010

@ -667,7 +667,7 @@ contains
t1 = mpi_wtime()
call psb_cdasb(desc_a,info)
call psb_spasb(a,desc_a,info,dup=1,afmt=afmt)
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
call blacs_barrier(icontxt,'ALL')
tasb = mpi_wtime()-t1
if(info.ne.0) then

Loading…
Cancel
Save