From 4b2f930cf6b4eadcbd6d03786e798389c5b88c75 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 14 Apr 2006 11:56:05 +0000 Subject: [PATCH] Changed handling of duplicates and storage format. Changed interface to both geins and spasb. New and better strategies. --- src/modules/psb_const.fh | 25 +- src/modules/psb_realloc_mod.f90 | 55 ++-- src/modules/psb_spmat_type.f90 | 133 ++++++++- src/modules/psb_tools_mod.f90 | 39 ++- src/prec/psb_dbldaggrmat.f90 | 3 +- src/prec/psb_zbldaggrmat.f90 | 3 +- src/psblas/psb_dspmm.f90 | 256 ++++++++--------- src/psblas/psb_dspsm.f90 | 189 ++++++------ src/serial/Makefile | 3 +- src/serial/aux/Makefile | 2 +- src/serial/{ => aux}/imsr.f90 | 3 +- src/serial/{ => aux}/imsrx.f90 | 3 +- src/serial/dp/dcoco.f | 21 +- src/serial/dp/dcocr.f | 16 +- src/serial/dp/dgindex.f | 4 +- src/serial/dp/zcoco.f | 14 +- src/serial/dp/zcocr.f | 28 +- src/serial/dp/zgindex.f | 2 +- src/serial/psb_dcoins.f90 | 177 ++++++------ src/serial/psb_dcsdp.f90 | 212 +++++++------- src/serial/psb_getifield.f90 | 14 + src/serial/psb_setifield.f90 | 14 + src/serial/psb_zcoins.f90 | 181 ++++++------ src/serial/psb_zcsdp.f90 | 211 +++++++------- src/tools/psb_dins.f90 | 210 ++++++++++---- src/tools/psb_dspasb.f90 | 317 +++++++++----------- src/tools/psb_dsprn.f90 | 57 ++-- src/tools/psb_iins.f90 | 415 +++++++++++++++----------- src/tools/psb_zins.f90 | 495 +++++++++++++++++++------------- src/tools/psb_zspasb.f90 | 209 ++++++-------- test/Fileread/mat_dist.f90 | 12 +- test/pargen/ppde90.f90 | 2 +- 32 files changed, 1878 insertions(+), 1447 deletions(-) rename src/serial/{ => aux}/imsr.f90 (98%) rename src/serial/{ => aux}/imsrx.f90 (98%) create mode 100644 src/serial/psb_getifield.f90 create mode 100644 src/serial/psb_setifield.f90 diff --git a/src/modules/psb_const.fh b/src/modules/psb_const.fh index 96502cc7..e221e8ae 100644 --- a/src/modules/psb_const.fh +++ b/src/modules/psb_const.fh @@ -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 ! diff --git a/src/modules/psb_realloc_mod.f90 b/src/modules/psb_realloc_mod.f90 index 9ed0a2c6..84ffd52b 100644 --- a/src/modules/psb_realloc_mod.f90 +++ b/src/modules/psb_realloc_mod.f90 @@ -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 diff --git a/src/modules/psb_spmat_type.f90 b/src/modules/psb_spmat_type.f90 index a92fa66c..30d5324d 100644 --- a/src/modules/psb_spmat_type.f90 +++ b/src/modules/psb_spmat_type.f90 @@ -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 diff --git a/src/modules/psb_tools_mod.f90 b/src/modules/psb_tools_mod.f90 index 086dda36..bbfde9fd 100644 --- a/src/modules/psb_tools_mod.f90 +++ b/src/modules/psb_tools_mod.f90 @@ -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 diff --git a/src/prec/psb_dbldaggrmat.f90 b/src/prec/psb_dbldaggrmat.f90 index 0e77fdb0..f8671321 100644 --- a/src/prec/psb_dbldaggrmat.f90 +++ b/src/prec/psb_dbldaggrmat.f90 @@ -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 diff --git a/src/prec/psb_zbldaggrmat.f90 b/src/prec/psb_zbldaggrmat.f90 index cd8b16c0..02f3810b 100644 --- a/src/prec/psb_zbldaggrmat.f90 +++ b/src/prec/psb_zbldaggrmat.f90 @@ -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 diff --git a/src/psblas/psb_dspmm.f90 b/src/psblas/psb_dspmm.f90 index a517c55c..14cfc57e 100644 --- a/src/psblas/psb_dspmm.f90 +++ b/src/psblas/psb_dspmm.f90 @@ -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 diff --git a/src/psblas/psb_dspsm.f90 b/src/psblas/psb_dspsm.f90 index 666a5a35..ce62b20b 100644 --- a/src/psblas/psb_dspsm.f90 +++ b/src/psblas/psb_dspsm.f90 @@ -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 diff --git a/src/serial/Makefile b/src/serial/Makefile index d2449540..51fb2e6e 100644 --- a/src/serial/Makefile +++ b/src/serial/Makefile @@ -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 . diff --git a/src/serial/aux/Makefile b/src/serial/aux/Makefile index fe02038b..55164e62 100644 --- a/src/serial/aux/Makefile +++ b/src/serial/aux/Makefile @@ -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) diff --git a/src/serial/imsr.f90 b/src/serial/aux/imsr.f90 similarity index 98% rename from src/serial/imsr.f90 rename to src/serial/aux/imsr.f90 index e0754fc3..dda3301e 100644 --- a/src/serial/imsr.f90 +++ b/src/serial/aux/imsr.f90 @@ -30,7 +30,8 @@ !!$ ! File: imsr.f90 ! Subroutine: - ! Parameters:subroutine imsr(n,x) + ! Parameters: +subroutine imsr(n,x) integer :: n integer :: x(n) diff --git a/src/serial/imsrx.f90 b/src/serial/aux/imsrx.f90 similarity index 98% rename from src/serial/imsrx.f90 rename to src/serial/aux/imsrx.f90 index 59ed5741..ce4725b3 100644 --- a/src/serial/imsrx.f90 +++ b/src/serial/aux/imsrx.f90 @@ -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) diff --git a/src/serial/dp/dcoco.f b/src/serial/dp/dcoco.f index 3f89cf51..acfab50e 100644 --- a/src/serial/dp/dcoco.f +++ b/src/serial/dp/dcoco.f @@ -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 diff --git a/src/serial/dp/dcocr.f b/src/serial/dp/dcocr.f index f2987404..55c6e5df 100644 --- a/src/serial/dp/dcocr.f +++ b/src/serial/dp/dcocr.f @@ -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', diff --git a/src/serial/dp/dgindex.f b/src/serial/dp/dgindex.f index 2b5baed3..9b8a96da 100644 --- a/src/serial/dp/dgindex.f +++ b/src/serial/dp/dgindex.f @@ -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 diff --git a/src/serial/dp/zcoco.f b/src/serial/dp/zcoco.f index a94b0033..e40e3b76 100644 --- a/src/serial/dp/zcoco.f +++ b/src/serial/dp/zcoco.f @@ -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 diff --git a/src/serial/dp/zcocr.f b/src/serial/dp/zcocr.f index 6b9a6884..a55799a0 100644 --- a/src/serial/dp/zcocr.f +++ b/src/serial/dp/zcocr.f @@ -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) diff --git a/src/serial/dp/zgindex.f b/src/serial/dp/zgindex.f index 92ddaf93..f5c8f9a7 100644 --- a/src/serial/dp/zgindex.f +++ b/src/serial/dp/zgindex.f @@ -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 diff --git a/src/serial/psb_dcoins.f90 b/src/serial/psb_dcoins.f90 index f54b20df..00508aaa 100644 --- a/src/serial/psb_dcoins.f90 +++ b/src/serial/psb_dcoins.f90 @@ -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 + + if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza - a%ia2(ip1+psb_nnz_) = nza - else - info = 2231 - call psb_errpush(info,name) - goto 9999 - endif + 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 @@ -191,7 +202,7 @@ contains real(kind(1.d0)), intent(inout) :: aspk(*) integer, intent(out) :: info integer :: i,ir,ic - + info = 0 if (nza >= nzl) then @@ -213,7 +224,7 @@ contains end if end do end if - + end subroutine psb_inner_upd subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,gtl,& diff --git a/src/serial/psb_dcsdp.f90 b/src/serial/psb_dcsdp.f90 index e4bbb522..26d0f843 100644 --- a/src/serial/psb_dcsdp.f90 +++ b/src/serial/psb_dcsdp.f90 @@ -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,& @@ -185,16 +188,16 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd) call psb_errpush(info,name,a_err=ch_err) goto 9999 endif - + b%pl(:) = 0 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 - do i=1, nnz - work(b%ia2(ip2+i-1)) = b%aspk(i) + work(b%ia2(ip2+i-1)) + 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 + 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 diff --git a/src/serial/psb_getifield.f90 b/src/serial/psb_getifield.f90 new file mode 100644 index 00000000..021537ec --- /dev/null +++ b/src/serial/psb_getifield.f90 @@ -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 diff --git a/src/serial/psb_setifield.f90 b/src/serial/psb_setifield.f90 new file mode 100644 index 00000000..76647bfd --- /dev/null +++ b/src/serial/psb_setifield.f90 @@ -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 diff --git a/src/serial/psb_zcoins.f90 b/src/serial/psb_zcoins.f90 index 91c5f66e..57dfd166 100644 --- a/src/serial/psb_zcoins.f90 +++ b/src/serial/psb_zcoins.f90 @@ -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 @@ -191,13 +198,13 @@ contains complex(kind(1.d0)), intent(inout) :: aspk(*) integer, intent(out) :: info integer :: i,ir,ic - + info = 0 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,12 +215,12 @@ 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 end if - + end subroutine psb_inner_upd subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,gtl,& @@ -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 diff --git a/src/serial/psb_zcsdp.f90 b/src/serial/psb_zcsdp.f90 index 6e32c29d..94a669c2 100644 --- a/src/serial/psb_zcsdp.f90 +++ b/src/serial/psb_zcsdp.f90 @@ -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,& @@ -186,16 +188,16 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd) call psb_errpush(info,name,a_err=ch_err) goto 9999 endif - + b%pl(:) = 0 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 - do i=1, nnz - work(b%ia2(ip2+i-1)) = b%aspk(i) + work(b%ia2(ip2+i-1)) + 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 + 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 diff --git a/src/tools/psb_dins.f90 b/src/tools/psb_dins.f90 index 9b52f6c5..9f60ab86 100644 --- a/src/tools/psb_dins.f90 +++ b/src/tools/psb_dins.f90 @@ -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 diff --git a/src/tools/psb_dspasb.f90 b/src/tools/psb_dspasb.f90 index 71fb0c88..783e8b77 100644 --- a/src/tools/psb_dspasb.f90 +++ b/src/tools/psb_dspasb.f90 @@ -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 diff --git a/src/tools/psb_dsprn.f90 b/src/tools/psb_dsprn.f90 index 916400dd..6895fbf2 100644 --- a/src/tools/psb_dsprn.f90 +++ b/src/tools/psb_dsprn.f90 @@ -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 diff --git a/src/tools/psb_iins.f90 b/src/tools/psb_iins.f90 index 151c1875..629753a4 100644 --- a/src/tools/psb_iins.f90 +++ b/src/tools/psb_iins.f90 @@ -45,13 +45,13 @@ ! 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 use psb_error_mod implicit none - + !....parameters... integer, intent(in) :: m,n type(psb_desc_type), intent(in) :: desc_a @@ -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 @@ -73,110 +74,138 @@ subroutine psb_iins(m, n, x, ix, jx, blck, desc_a, info,& call psb_erractionsave(err_act) name = 'psb_iins' - + 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,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 + endif + if (present(dupl)) then + dupl_ = dupl else - jblock = 1 + 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,11 +213,11 @@ 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 - + end subroutine psb_iins @@ -238,23 +267,23 @@ 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 use psb_error_mod implicit none - + ! m rows number of submatrix belonging to blck to be inserted - + ! iblck first row of submatrix belonging to blck to be inserted - + ! ix x global-row corresponding to position at which blck submatrix ! must be inserted - + ! jx x global-col corresponding to position at which blck submatrix ! must be inserted - + !....parameters... integer, intent(in) :: m type(psb_desc_type), intent(in) :: desc_a @@ -263,42 +292,71 @@ 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 info=0 name = 'psb_iinsvm' call psb_erractionsave(err_act) - - + + 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 - 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,11 +364,11 @@ 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 - + end subroutine psb_iinsvm @@ -356,20 +414,20 @@ 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 use psb_error_mod implicit none - + ! m rows number of submatrix belonging to blck to be inserted - + ! iblck first row of submatrix belonging to blck to be inserted - + ! ix x global-row corresponding to position at which blck submatrix ! must be inserted - + !....parameters... integer, intent(in) :: m type(psb_desc_type), intent(in) :: desc_a @@ -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,34 +506,62 @@ 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 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 diff --git a/src/tools/psb_zins.f90 b/src/tools/psb_zins.f90 index 144a4dca..c05d7245 100644 --- a/src/tools/psb_zins.f90 +++ b/src/tools/psb_zins.f90 @@ -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 diff --git a/src/tools/psb_zspasb.f90 b/src/tools/psb_zspasb.f90 index dd788f9b..80233f5e 100644 --- a/src/tools/psb_zspasb.f90 +++ b/src/tools/psb_zspasb.f90 @@ -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) diff --git a/test/Fileread/mat_dist.f90 b/test/Fileread/mat_dist.f90 index b5967484..4205fcd7 100644 --- a/test/Fileread/mat_dist.f90 +++ b/test/Fileread/mat_dist.f90 @@ -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 diff --git a/test/pargen/ppde90.f90 b/test/pargen/ppde90.f90 index 1898dbb7..cce97926 100644 --- a/test/pargen/ppde90.f90 +++ b/test/pargen/ppde90.f90 @@ -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