From bd27d1a7bb7abe3f855211e624060a26544ae426 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 18 Apr 2006 09:20:17 +0000 Subject: [PATCH] Fixed bugs in duplicate/regeneration handling. --- src/modules/psb_const.fh | 2 +- src/modules/psb_spmat_type.f90 | 175 +++++++++++++++++++++ src/serial/aux/Makefile | 2 +- src/serial/aux/ibsrch.f | 1 - src/serial/aux/issrch.f | 47 ++++++ src/serial/dp/dcoco.f | 4 +- src/serial/dp/dcocr.f | 34 +++-- src/serial/dp/zcoco.f | 2 +- src/serial/dp/zcocr.f | 2 +- src/serial/psb_dcoins.f90 | 253 +++++++++++++++++++++++++++++-- src/serial/psb_dcsdp.f90 | 68 +++++---- src/serial/psb_zcoins.f90 | 269 ++++++++++++++++++++++++++++++--- src/serial/psb_zcsdp.f90 | 47 +++++- src/tools/psb_dspasb.f90 | 11 +- src/tools/psb_zspasb.f90 | 12 +- 15 files changed, 829 insertions(+), 100 deletions(-) create mode 100644 src/serial/aux/issrch.f diff --git a/src/modules/psb_const.fh b/src/modules/psb_const.fh index e221e8ae..97030853 100644 --- a/src/modules/psb_const.fh +++ b/src/modules/psb_const.fh @@ -102,8 +102,8 @@ 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_upd_perm_ = 98765 integer, parameter :: psb_isrtdcoo_ = 98761 integer, parameter :: psb_maxjdrows_=8, psb_minjdrows_=4 integer, parameter :: psb_dbleint_=2 diff --git a/src/modules/psb_spmat_type.f90 b/src/modules/psb_spmat_type.f90 index 30d5324d..44486dcb 100644 --- a/src/modules/psb_spmat_type.f90 +++ b/src/modules/psb_spmat_type.f90 @@ -92,6 +92,10 @@ module psb_spmat_type module procedure psb_dsp_transfer, psb_zsp_transfer end interface + interface psb_sp_trimsize + module procedure psb_dsp_trimsize, psb_zsp_trimsize + end interface + interface psb_sp_reall module procedure psb_dspreallocate, psb_dspreall3, & & psb_zspreallocate, psb_zspreall3 @@ -110,6 +114,10 @@ module psb_spmat_type module procedure psb_dspreinit, psb_zspreinit end interface + interface psb_sp_sizeof + module procedure psb_dspsizeof, psb_zspsizeof + end interface + contains subroutine psb_nullify_dsp(mat) @@ -323,7 +331,9 @@ contains logical, parameter :: debug=.false. info = 0 + if (debug) write(0,*) 'Before realloc',nd,size(a%aspk) call psb_realloc(nd,a%aspk,info) + if (debug) write(0,*) 'After realloc',nd,size(a%aspk),info if (info /= 0) return call psb_realloc(ni2,a%ia2,info) if (info /= 0) return @@ -453,6 +463,52 @@ contains end subroutine psb_dsp_setifld + subroutine psb_dsp_trimsize(a, i1,i2,ia,info) + use psb_string_mod + implicit none + !....Parameters... + Type(psb_dspmat_type), intent(in) :: A + Integer, intent(out) :: i1, i2, ia, info + + !locals + Integer :: nza,nz1, nz2, nzl, nzr + logical, parameter :: debug=.false. + + info = 0 + if (psb_sp_getifld(psb_upd_,a,info) == psb_upd_perm_) then + info = -1 + i1 = size(a%ia1) + i2 = size(a%ia2) + ia = size(a%aspk) + return + endif + select case(toupper(a%fida)) + case('CSR') + nza = a%ia2(a%m+1)-1 + ia = nza + i1 = nza + i2 = a%m + 1 + case('COO') + nza = a%infoa(psb_nnz_) + i1 = nza + i2 = nza + ia = nza + case('JAD') + ! Feeling lazy today + i1 = size(a%ia1) + i2 = size(a%ia2) + ia = size(a%aspk) + case default + info = -2 + i1 = size(a%ia1) + i2 = size(a%ia2) + ia = size(a%aspk) + end select + + Return + + End Subroutine psb_dsp_trimsize + function psb_dsp_getifld(field,a,info) implicit none !....Parameters... @@ -487,6 +543,42 @@ contains end function psb_dsp_getifld + function psb_dspsizeof(a) + implicit none + !....Parameters... + + Type(psb_dspmat_type), intent(in) :: A + Integer :: psb_dspsizeof + + !locals + logical, parameter :: debug=.false. + integer :: val + + val = 4*size(a%infoa) + + if (associated(a%aspk)) then + val = val + 8 * size(a%aspk) + endif + + if (associated(a%ia1)) then + val = val + 4 * size(a%ia1) + endif + if (associated(a%ia2)) then + val = val + 4 * size(a%ia2) + endif + if (associated(a%pl)) then + val = val + 4 * size(a%pl) + endif + if (associated(a%pr)) then + val = val + 4 * size(a%pr) + endif + + + psb_dspsizeof = val + Return + + end function psb_dspsizeof + subroutine psb_dsp_free(a,info) implicit none !....Parameters... @@ -857,6 +949,52 @@ contains end subroutine psb_zsp_setifld + subroutine psb_zsp_trimsize(a, i1,i2,ia,info) + use psb_string_mod + implicit none + !....Parameters... + Type(psb_zspmat_type), intent(in) :: A + Integer, intent(out) :: i1, i2, ia, info + + !locals + Integer :: nza,nz1, nz2, nzl, nzr + logical, parameter :: debug=.false. + + info = 0 + if (psb_sp_getifld(psb_upd_,a,info) == psb_upd_perm_) then + info = -1 + i1 = size(a%ia1) + i2 = size(a%ia2) + ia = size(a%aspk) + return + endif + select case(toupper(a%fida)) + case('CSR') + nza = a%ia2(a%m+1)-1 + ia = nza + i1 = nza + i2 = a%m + 1 + case('COO') + nza = a%infoa(psb_nnz_) + i1 = nza + i2 = nza + ia = nza + case('JAD') + ! Feeling lazy today + i1 = size(a%ia1) + i2 = size(a%ia2) + ia = size(a%aspk) + case default + info = -2 + i1 = size(a%ia1) + i2 = size(a%ia2) + ia = size(a%aspk) + end select + + Return + + End Subroutine psb_zsp_trimsize + function psb_zsp_getifld(field,a,info) implicit none !....Parameters... @@ -891,6 +1029,43 @@ contains end function psb_zsp_getifld + function psb_zspsizeof(a) + implicit none + !....Parameters... + + Type(psb_zspmat_type), intent(in) :: A + Integer :: psb_zspsizeof + + !locals + logical, parameter :: debug=.false. + integer :: val + + val = 4*size(a%infoa) + + if (associated(a%aspk)) then + val = val + 16 * size(a%aspk) + endif + + if (associated(a%ia1)) then + val = val + 4 * size(a%ia1) + endif + if (associated(a%ia2)) then + val = val + 4 * size(a%ia2) + endif + if (associated(a%pl)) then + val = val + 4 * size(a%pl) + endif + if (associated(a%pr)) then + val = val + 4 * size(a%pr) + endif + + + psb_zspsizeof = val + Return + + end function psb_zspsizeof + + subroutine psb_zsp_free(a,info) diff --git a/src/serial/aux/Makefile b/src/serial/aux/Makefile index 55164e62..48fb997f 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 imsr.o imsrx.o + mrgsrt.o isaperm.o ibsrch.o issrch.o imsr.o imsrx.o OBJS=$(FOBJS) diff --git a/src/serial/aux/ibsrch.f b/src/serial/aux/ibsrch.f index 518259ba..eb676b7b 100644 --- a/src/serial/aux/ibsrch.f +++ b/src/serial/aux/ibsrch.f @@ -50,7 +50,6 @@ C lb = m + 1 end if enddo - return end diff --git a/src/serial/aux/issrch.f b/src/serial/aux/issrch.f new file mode 100644 index 00000000..c5119c15 --- /dev/null +++ b/src/serial/aux/issrch.f @@ -0,0 +1,47 @@ +C +C Parallel Sparse BLAS v2.0 +C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +C Alfredo Buttari University of Rome Tor Vergata +C +C Redistribution and use in source and binary forms, with or without +C modification, are permitted provided that the following conditions +C are met: +C 1. Redistributions of source code must retain the above copyright +C notice, this list of conditions and the following disclaimer. +C 2. Redistributions in binary form must reproduce the above copyright +C notice, this list of conditions, and the following disclaimer in the +C documentation and/or other materials provided with the distribution. +C 3. The name of the PSBLAS group or the names of its contributors may +C not be used to endorse or promote products derived from this +C software without specific written permission. +C +C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +C POSSIBILITY OF SUCH DAMAGE. +C +C + subroutine issrch(ipos,key,n,v) + integer ipos, key, n + integer v(n) + + integer i + + + ipos = -1 + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + return + end + diff --git a/src/serial/dp/dcoco.f b/src/serial/dp/dcoco.f index acfab50e..36c0605a 100644 --- a/src/serial/dp/dcoco.f +++ b/src/serial/dp/dcoco.f @@ -71,14 +71,14 @@ c ierror = 0 call fcpsb_erractionsave(err_act) - call psb_getifield(check_flag,psb_dupl_,info,psb_ifasize_,ierror) + call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror) if (trans.eq.'N') then scale = (unitd.eq.'L') ! meaningless p1(1) = 0 p2(1) = 0 - call psb_getifield(nnz,psb_nnz_,info,psb_ifasize_,ierror) + call psb_getifield(nnz,psb_nnz_,infon,psb_ifasize_,ierror) if (debug) then write(*,*) 'on entry to dcoco: nnz laux ', + nnz,laux,larn,lia1n,lia2n diff --git a/src/serial/dp/dcocr.f b/src/serial/dp/dcocr.f index 55c6e5df..63b60bd7 100644 --- a/src/serial/dp/dcocr.f +++ b/src/serial/dp/dcocr.f @@ -73,7 +73,7 @@ C IERROR = 0 CALL FCPSB_ERRACTIONSAVE(ERR_ACT) - call psb_getifield(check_flag,psb_dupl_,info,psb_ifasize_,ierror) + call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror) IF ((TRANS.EQ.'N').or.(TRANS.EQ.'n')) THEN @@ -230,16 +230,22 @@ C ... Error, there are duplicated elements ... C ... Insert only the last duplicated element ... ARN(ELEM_CSR-1) = AR(ELEM) ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1 + if (debug) write(0,*) 'Duplicated overwrite perm ', + + elem_csr-1,elem 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 + if (debug) write(0,*) 'Duplicated add perm ', + + elem_csr-1,elem END IF ENDIF ELEM = ELEM + 1 ENDDO IAN2(ROW+1) = ELEM_CSR ENDDO + + ELSE c$$$ write(0,*) 'Going for ELSE !!!?' C .... Order with key IA ... @@ -304,12 +310,12 @@ C ... Error, there are duplicated elements ... 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', + if (debug) write(0,*) 'Duplicated overwrite srch', + elem_csr-1,elem 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', + if (debug) write(0,*) 'Duplicated add srch', + elem_csr-1,elem END IF ENDIF @@ -321,11 +327,11 @@ C ... Sum the duplicated element ... if (debug) write(0,*)'Done Rebuild CSR', + ian2(m+1),ia(elem) - if (debug) then - do i=ian2(m+1), nnz - write(0,*) 'Overflow check :',ia(i),ja(i),ar(i) - enddo - endif +c$$$ if (debug) then +c$$$ do i=ian2(m+1), nnz +c$$$ write(0,*) 'Overflow check :',ia(i),ja(i),ar(i) +c$$$ enddo +c$$$ endif ELSE IF (DESCRA(1:1).EQ.'S' .AND. DESCRA(2:2).EQ.'U') THEN @@ -396,19 +402,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 +500,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/zcoco.f b/src/serial/dp/zcoco.f index e40e3b76..f39d2d51 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) - call psb_getifield(check_flag,psb_dupl_,info,psb_ifasize_,ierror) + call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror) if (trans.eq.'N') then scale = (unitd.eq.'L') ! meaningless p1(1) = 0 diff --git a/src/serial/dp/zcocr.f b/src/serial/dp/zcocr.f index a55799a0..ef0d4fd5 100644 --- a/src/serial/dp/zcocr.f +++ b/src/serial/dp/zcocr.f @@ -73,7 +73,7 @@ C IERROR = 0 CALL FCPSB_ERRACTIONSAVE(ERR_ACT) - call psb_getifield(check_flag,psb_dupl_,info,psb_ifasize_,ierror) + call psb_getifield(check_flag,psb_dupl_,infon,psb_ifasize_,ierror) IF ((TRANS.EQ.'N').or.(TRANS.EQ.'n')) THEN diff --git a/src/serial/psb_dcoins.f90 b/src/serial/psb_dcoins.f90 index 00508aaa..2e1d855d 100644 --- a/src/serial/psb_dcoins.f90 +++ b/src/serial/psb_dcoins.f90 @@ -117,7 +117,7 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) goto 9999 endif endif - call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,& + call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,ng,& & imin,imax,jmin,jmax,info) if(info.ne.izero) then info=4010 @@ -144,7 +144,10 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,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,& + nza = a%ia2(ip1+psb_nnz_) + nzl = a%infoa(psb_del_bnd_) + + call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,ng,& & imin,imax,jmin,jmax,nzl,info) if(info.ne.izero) then @@ -163,16 +166,32 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza case (psb_upd_dflt_, psb_upd_srch_) - write(0,*) 'Default & search inner update to be implemented' - info = 2230 - call psb_errpush(info,name) - goto 9999 + + select case(toupper(a%fida)) + case ('CSR') +!!$ write(0,*) 'Calling csr_inner_upd' + call csr_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) +!!$ write(0,*) 'From csr_inner_upd:',info + case ('COO') + call coo_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) + + case default + + info = 2230 + call psb_errpush(info,name) + goto 9999 + end select + case default + info = 2231 call psb_errpush(info,name) goto 9999 end select + case default info = 2232 call psb_errpush(info,name) @@ -192,10 +211,12 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) return contains - subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,gtl,imin,imax,jmin,jmax,nzl,info) + + subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) implicit none - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng integer, intent(in) :: ia(*),ja(*),gtl(*) integer, intent(inout) :: nza real(kind(1.d0)), intent(in) :: val(*) @@ -203,8 +224,6 @@ contains integer, intent(out) :: info integer :: i,ir,ic - info = 0 - if (nza >= nzl) then do i=1, nz nza = nza + 1 @@ -227,11 +246,11 @@ contains end subroutine psb_inner_upd - subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,gtl,& + subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,gtl,ng,& & imin,imax,jmin,jmax,info) implicit none - integer, intent(in) :: nz, imin,imax,jmin,jmax + integer, intent(in) :: nz, imin,imax,jmin,jmax,ng integer, intent(in) :: ia(*),ja(*),gtl(*) integer, intent(inout) :: nza,ia1(*),ia2(*) real(kind(1.d0)), intent(in) :: val(*) @@ -258,5 +277,215 @@ contains end do end subroutine psb_inner_ins + + + subroutine csr_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) + implicit none + + type(psb_dspmat_type), intent(inout) :: a + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng + integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(inout) :: nza + real(kind(1.d0)), intent(in) :: val(*) + integer, intent(out) :: info + integer :: i,ir,ic,check_flag, ilr, ilc, ip, & + & i1,i2,nc,lb,ub,m,nnz,dupl + + info = 0 + + dupl = psb_sp_getifld(psb_dupl_,a,info) + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + i1 = a%ia2(ir) + i2 = a%ia2(ir+1) + nc=i2-i1 + + + if (.true.) then + call issrch(ip,ic,nc,a%ia1(i1:i2-1)) + if (ip>0) then + a%aspk(i1+ip-1) = val(i) + else + write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1) + info = -1 + return + end if + + else +!!$ + ip = -1 + lb = i1 + ub = i2-1 + do + if (lb > ub) exit + m = (lb+ub)/2 +!!$ write(0,*) 'Debug: ',lb,m,ub + if (ic == a%ia1(m)) then + ip = m + lb = ub + 1 + else if (ic < a%ia1(m)) then + ub = m-1 + else + lb = m + 1 + end if + enddo + + if (ip>0) then + a%aspk(ip) = val(i) + else + write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1) + info = -1 + return + end if + + end if + end if + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + i1 = a%ia2(ir) + i2 = a%ia2(ir+1) + nc = i2-i1 + call issrch(ip,ic,nc,a%ia1(i1:i2-1)) + if (ip>0) then + a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) + else + info = -2 + return + end if + end if + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + end subroutine csr_inner_upd + + subroutine coo_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) + implicit none + + type(psb_dspmat_type), intent(inout) :: a + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng + integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(inout) :: nza + real(kind(1.d0)), intent(in) :: val(*) + integer, intent(out) :: info + integer :: i,ir,ic,check_flag, ilr, ilc, ip, & + & i1,i2,nc,lb,ub,m,nnz,dupl,isrt + + info = 0 + + dupl = psb_sp_getifld(psb_dupl_,a,info) + + if (psb_sp_getifld(psb_srtd_,a,info) /= psb_isrtdcoo_) then + info = -4 + return + end if + + ilr = -1 + ilc = -1 + nnz = psb_sp_getifld(psb_nnz_,a,info) + + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if (ir /= ilr) then + call ibsrch(i1,ir,nnz,a%ia1) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia1(i2+1) /= a%ia1(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia1(i1-1) /= a%ia1(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + end if + nc = i2-i1+1 + call issrch(ip,ic,nc,a%ia2(i1:i2)) + if (ip>0) then + a%aspk(i1+ip-1) = val(i) + else + info = -2 + return + end if + end if + end do + case(psb_dupl_add_) + ! Add + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if (ir /= ilr) then + call ibsrch(i1,ir,nnz,a%ia1) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia1(i2+1) /= a%ia1(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia1(i1-1) /= a%ia1(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + end if + nc = i2-i1+1 + call issrch(ip,ic,nc,a%ia2(i1:i2)) + if (ip>0) then + a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) + else + info = -2 + return + end if + end if + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + + end subroutine coo_inner_upd + end subroutine psb_dcoins diff --git a/src/serial/psb_dcsdp.f90 b/src/serial/psb_dcsdp.f90 index af680fe3..7c84bdea 100644 --- a/src/serial/psb_dcsdp.f90 +++ b/src/serial/psb_dcsdp.f90 @@ -64,7 +64,7 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) Integer :: nzr, ntry, ifc_,ierror, ia1_size,& & ia2_size, aspk_size,size_req,n_row,n_col,upd_,dupl_ integer :: ip1, ip2, nnz, iflag, ichk, nnzt,& - & ipc, i, count, err_act, ierrv(5) + & ipc, i, count, err_act, ierrv(5), i1, i2, ia character :: check_,trans_,unitd_, up Integer, Parameter :: maxtry=8 logical, parameter :: debug=.false. @@ -104,12 +104,12 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) check_ = 'N' endif if (present(trans)) then - trans_ = toupper(trans ) + trans_ = toupper(trans) else trans_ = 'N' endif if (present(unitd)) then - unitd_ = toupper(unitd ) + unitd_ = toupper(unitd) else unitd_ = 'U' endif @@ -161,21 +161,37 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) call psb_sp_setifld(dupl,psb_dupl_,b,info) end if - + upd_ = psb_sp_getifld(psb_upd_,b,info) + select case(upd_) + case(psb_upd_dflt_,psb_upd_srch_,psb_upd_perm_) + ! Legal value, do nothing + case default + ! Fix bad value + upd_ = psb_upd_dflt_ + call psb_sp_setifld(upd_,psb_upd_,b,info) + end select + dupl_ = psb_sp_getifld(psb_dupl_,b,info) + select case(dupl_) + case(psb_dupl_ovwrt_,psb_dupl_add_,psb_dupl_err_) + ! Legal value, do nothing + case default + ! Fix bad value + dupl_ = psb_dupl_def_ + call psb_sp_setifld(dupl_,psb_dupl_,b,info) + end select + ! ...matrix conversion... b%m=a%m b%k=a%k call psb_spinfo(psb_nztotreq_,a,size_req,info) if (debug) write(0,*) 'DCSDP : size_req 1:',size_req ! - upd_ = psb_sp_getifld(psb_upd_,b,info) - + n_row=b%m n_col=b%k call psb_cest(b%fida, n_row,n_col,size_req,& & ia1_size, ia2_size, aspk_size, upd_,info) -!!$ write(0,*) 'ESTIMATE : ',ia1_size,ia2_size,aspk_Size,upd_ if (info /= no_err) then info=4010 ch_err='psb_cest' @@ -208,12 +224,6 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) case ('CSR') -!!$ -!!$ ia1_size=a%infoa(psb_nnz_) -!!$ ia2_size=a%m+1 -!!$ aspk_size=a%infoa(psb_nnz_) -!!$ call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info) - call dcrcr(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& & a%ia1, a%ia2, a%infoa, b%pl, b%descra, b%aspk, b%ia1,& & b%ia2, b%infoa, b%pr, size(b%aspk), size(b%ia1),& @@ -231,10 +241,6 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) !...converting to JAD !...output matrix may not be big enough -!!$ ia1_size=a%infoa(psb_nnz_) -!!$ ia2_size=a%m+1 -!!$ aspk_size=a%infoa(psb_nnz_) -!!$ call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info) do call dcrjd(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& @@ -274,9 +280,6 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) case ('COO') -!!$ aspk_size=max(size(a%aspk),a%ia2(a%m+1)) -!!$ call psb_sp_reall(b,aspk_size,info) -!!$ write(0,*) 'From DCSDP90:',b%fida,size(b%aspk),info call dcrco(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& & a%ia1, a%ia2, a%infoa, b%pl, b%descra, b%aspk, b%ia1,& & b%ia2, b%infoa, b%pr, size(b%aspk), size(b%ia1),& @@ -295,8 +298,6 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) case ('CSR') -!!$ aspk_size=max(size(a%aspk),a%ia2(a%m+1)) -!!$ call psb_sp_reall(b,aspk_size,info) call dcocr(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& & a%ia2, a%ia1, a%infoa, b%pl, b%descra, b%aspk, b%ia1,& & b%ia2, b%infoa, b%pr, size(b%aspk), size(b%ia1),& @@ -364,12 +365,8 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) end do - - case ('COO') -!!$ aspk_size=max(size(a%aspk),a%ia2(a%m+1)) -!!$ call psb_sp_reall(b,aspk_size,info) call dcoco(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& & a%ia1, a%ia2, a%infoa, b%pl, b%descra, b%aspk, b%ia1,& & b%ia2, b%infoa, b%pr, size(b%aspk), size(b%ia1),& @@ -383,18 +380,32 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) end select +!!$ write(0,*) 'End of assembly', psb_sp_getifld(psb_upd_,b,info) ,psb_upd_perm_ + if (psb_sp_getifld(psb_upd_,b,info) /= psb_upd_perm_) then +!!$ write(0,*) 'Going for trimsize',size(b%ia1),size(b%ia2),size(b%aspk) + call psb_sp_trimsize(b,i1,i2,ia,info) +!!$ write(0,*) 'From trimsize',i1,i2,ia,info + if (info == 0) call psb_sp_reall(b,i1,i2,ia,info) +!!$ write(0,*) 'From realloc',size(b%ia1),size(b%ia2),size(b%aspk) + endif + else if (check_=='R') then + !...Regenerating matrix + if (psb_sp_getifld(psb_state_,b,info) /= psb_spmat_upd_) then info = 8888 call psb_errpush(info,name) goto 9999 endif + ! + ! dupl_ and upd_ fields should not be changed. + ! select case(psb_sp_getifld(psb_upd_,b,info)) case(psb_upd_perm_) - + if (debug) write(0,*) 'Regeneration with 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_) @@ -484,7 +495,8 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) endif case(psb_upd_dflt_,psb_upd_srch_) - ! Nothing to be done + ! Nothing to be done here. + if (debug) write(0,*) 'Going through on regeneration with psb_upd_srch_' case default ! Wrong value info = 8888 diff --git a/src/serial/psb_zcoins.f90 b/src/serial/psb_zcoins.f90 index 57dfd166..d852673d 100644 --- a/src/serial/psb_zcoins.f90 +++ b/src/serial/psb_zcoins.f90 @@ -28,7 +28,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: psbdcoins.f90 +! File: psbzcoins.f90 ! Subroutine: ! Parameters: subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) @@ -50,7 +50,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), iupd - logical, parameter :: debug=.true. + logical, parameter :: debug=.false. character(len=20) :: name, ch_err name='psb_zcoins' @@ -91,6 +91,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) spstate = psb_sp_getifld(psb_state_,a,info) select case(spstate) + case(psb_spmat_bld_) if ((ufida /= 'COO').and.(ufida/='COI')) then info = 134 @@ -116,7 +117,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) goto 9999 endif endif - call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,& + call psb_inner_ins(nz,ia,ja,val,nza,a%ia1,a%ia2,a%aspk,gtl,ng,& & imin,imax,jmin,jmax,info) if(info.ne.izero) then info=4010 @@ -129,10 +130,10 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) 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 + if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then + call psb_sp_setifld(nza,psb_del_bnd_,a,info) endif - a%infoa(psb_nnz_) = nza + call psb_sp_setifld(nza,psb_nnz_,a,info) case(psb_spmat_upd_) @@ -143,8 +144,12 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,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,& + nza = a%ia2(ip1+psb_nnz_) + nzl = a%infoa(psb_del_bnd_) + + call psb_inner_upd(nz,ia,ja,val,nza,a%aspk,gtl,ng,& & imin,imax,jmin,jmax,nzl,info) + if(info.ne.izero) then info=4010 ch_err='psb_inner_upd' @@ -156,19 +161,37 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) 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 + 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 + + select case(toupper(a%fida)) + case ('CSR') +!!$ write(0,*) 'Calling csr_inner_upd' + call csr_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) +!!$ write(0,*) 'From csr_inner_upd:',info + case ('COO') + call coo_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) + + case default + + info = 2230 + call psb_errpush(info,name) + goto 9999 + end select + case default + info = 2231 call psb_errpush(info,name) goto 9999 end select + case default info = 2232 call psb_errpush(info,name) @@ -188,10 +211,12 @@ subroutine psb_zcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) return contains - subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,gtl,imin,imax,jmin,jmax,nzl,info) + + subroutine psb_inner_upd(nz,ia,ja,val,nza,aspk,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) implicit none - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng integer, intent(in) :: ia(*),ja(*),gtl(*) integer, intent(inout) :: nza complex(kind(1.d0)), intent(in) :: val(*) @@ -199,8 +224,6 @@ contains integer, intent(out) :: info integer :: i,ir,ic - info = 0 - if (nza >= nzl) then do i=1, nz nza = nza + 1 @@ -223,11 +246,11 @@ contains end subroutine psb_inner_upd - subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,gtl,& + subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,gtl,ng,& & imin,imax,jmin,jmax,info) implicit none - integer, intent(in) :: nz, imin,imax,jmin,jmax + integer, intent(in) :: nz, imin,imax,jmin,jmax,ng integer, intent(in) :: ia(*),ja(*),gtl(*) integer, intent(inout) :: nza,ia1(*),ia2(*) complex(kind(1.d0)), intent(in) :: val(*) @@ -254,5 +277,215 @@ contains end do end subroutine psb_inner_ins + + + subroutine csr_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) + implicit none + + type(psb_zspmat_type), intent(inout) :: a + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng + integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(inout) :: nza + complex(kind(1.d0)), intent(in) :: val(*) + integer, intent(out) :: info + integer :: i,ir,ic,check_flag, ilr, ilc, ip, & + & i1,i2,nc,lb,ub,m,nnz,dupl + + info = 0 + + dupl = psb_sp_getifld(psb_dupl_,a,info) + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + i1 = a%ia2(ir) + i2 = a%ia2(ir+1) + nc=i2-i1 + + + if (.true.) then + call issrch(ip,ic,nc,a%ia1(i1:i2-1)) + if (ip>0) then + a%aspk(i1+ip-1) = val(i) + else + write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1) + info = -1 + return + end if + + else +!!$ + ip = -1 + lb = i1 + ub = i2-1 + do + if (lb > ub) exit + m = (lb+ub)/2 +!!$ write(0,*) 'Debug: ',lb,m,ub + if (ic == a%ia1(m)) then + ip = m + lb = ub + 1 + else if (ic < a%ia1(m)) then + ub = m-1 + else + lb = m + 1 + end if + enddo + + if (ip>0) then + a%aspk(ip) = val(i) + else + write(0,*)'Was searching ',ic,' in: ',i1,i2,' : ',a%ia1(i1:i2-1) + info = -1 + return + end if + + end if + end if + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + i1 = a%ia2(ir) + i2 = a%ia2(ir+1) + nc = i2-i1 + call issrch(ip,ic,nc,a%ia1(i1:i2-1)) + if (ip>0) then + a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) + else + info = -2 + return + end if + end if + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + end subroutine csr_inner_upd + + subroutine coo_inner_upd(nz,ia,ja,val,nza,a,gtl,ng,& + & imin,imax,jmin,jmax,nzl,info) + implicit none + + type(psb_zspmat_type), intent(inout) :: a + integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl,ng + integer, intent(in) :: ia(*),ja(*),gtl(*) + integer, intent(inout) :: nza + complex(kind(1.d0)), intent(in) :: val(*) + integer, intent(out) :: info + integer :: i,ir,ic,check_flag, ilr, ilc, ip, & + & i1,i2,nc,lb,ub,m,nnz,dupl,isrt + + info = 0 + + dupl = psb_sp_getifld(psb_dupl_,a,info) + + if (psb_sp_getifld(psb_srtd_,a,info) /= psb_isrtdcoo_) then + info = -4 + return + end if + + ilr = -1 + ilc = -1 + nnz = psb_sp_getifld(psb_nnz_,a,info) + + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if (ir /= ilr) then + call ibsrch(i1,ir,nnz,a%ia1) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia1(i2+1) /= a%ia1(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia1(i1-1) /= a%ia1(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + end if + nc = i2-i1+1 + call issrch(ip,ic,nc,a%ia2(i1:i2)) + if (ip>0) then + a%aspk(i1+ip-1) = val(i) + else + info = -2 + return + end if + end if + end do + case(psb_dupl_add_) + ! Add + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if (ir /= ilr) then + call ibsrch(i1,ir,nnz,a%ia1) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia1(i2+1) /= a%ia1(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia1(i1-1) /= a%ia1(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + end if + nc = i2-i1+1 + call issrch(ip,ic,nc,a%ia2(i1:i2)) + if (ip>0) then + a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) + else + info = -2 + return + end if + end if + end do + + case default + info = -3 + write(0,*) 'Duplicate handling: ',dupl + end select + + end subroutine coo_inner_upd + end subroutine psb_zcoins diff --git a/src/serial/psb_zcsdp.f90 b/src/serial/psb_zcsdp.f90 index a112488e..14dec36a 100644 --- a/src/serial/psb_zcsdp.f90 +++ b/src/serial/psb_zcsdp.f90 @@ -64,7 +64,7 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) Integer :: nzr, ntry, ifc_,ierror, ia1_size,& & ia2_size, aspk_size,size_req,n_row,n_col,upd_,dupl_ integer :: ip1, ip2, nnz, iflag, ichk, nnzt,& - & ipc, i, count, err_act, ierrv(5) + & ipc, i, count, err_act, ierrv(5), i1, i2, ia character :: check_,trans_,unitd_, up Integer, Parameter :: maxtry=8 logical, parameter :: debug=.false. @@ -104,12 +104,12 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) check_ = 'N' endif if (present(trans)) then - trans_ = toupper(trans ) + trans_ = toupper(trans) else trans_ = 'N' endif if (present(unitd)) then - unitd_ = toupper(unitd ) + unitd_ = toupper(unitd) else unitd_ = 'U' endif @@ -120,6 +120,7 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) allocate(work(max(size(a%ia1),size(b%ia1),& & size(a%ia2),size(b%ia2))+max(a%m,b%m)+1000),stat=info) endif + if (info /= 0) then info=2040 call psb_errpush(info,name) @@ -157,15 +158,32 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) call psb_sp_setifld(dupl,psb_dupl_,b,info) end if - + upd_ = psb_sp_getifld(psb_upd_,b,info) + select case(upd_) + case(psb_upd_dflt_,psb_upd_srch_,psb_upd_perm_) + ! Legal value, do nothing + case default + ! Fix bad value + upd_ = psb_upd_dflt_ + call psb_sp_setifld(upd_,psb_upd_,b,info) + end select + dupl_ = psb_sp_getifld(psb_dupl_,b,info) + select case(dupl_) + case(psb_dupl_ovwrt_,psb_dupl_add_,psb_dupl_err_) + ! Legal value, do nothing + case default + ! Fix bad value + dupl_ = psb_dupl_def_ + call psb_sp_setifld(dupl_,psb_dupl_,b,info) + end select + ! ...matrix conversion... b%m=a%m b%k=a%k call psb_spinfo(psb_nztotreq_,a,size_req,info) if (debug) write(0,*) 'DCSDP : size_req 1:',size_req ! - upd_ = psb_sp_getifld(psb_upd_,b,info) - + n_row=b%m n_col=b%k call psb_cest(b%fida, n_row,n_col,size_req,& @@ -364,6 +382,15 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) end select +!!$ write(0,*) 'End of assembly', psb_sp_getifld(psb_upd_,b,info) ,psb_upd_perm_ + if (psb_sp_getifld(psb_upd_,b,info) /= psb_upd_perm_) then +!!$ write(0,*) 'Going for trimsize',size(b%ia1),size(b%ia2),size(b%aspk) + call psb_sp_trimsize(b,i1,i2,ia,info) +!!$ write(0,*) 'From trimsize',i1,i2,ia,info + if (info == 0) call psb_sp_reall(b,i1,i2,ia,info) +!!$ write(0,*) 'From realloc',size(b%ia1),size(b%ia2),size(b%aspk) + endif + else if (check_=='R') then !...Regenerating matrix if (psb_sp_getifld(psb_state_,b,info) /= psb_spmat_upd_) then @@ -372,10 +399,13 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) goto 9999 endif + ! + ! dupl_ and upd_ fields should not be changed. + ! select case(psb_sp_getifld(psb_upd_,b,info)) case(psb_upd_perm_) - + if (debug) write(0,*) 'Regeneration with 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_) @@ -465,7 +495,8 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) endif case(psb_upd_dflt_,psb_upd_srch_) - ! Nothing to be done + ! Nothing to be done here. + if (debug) write(0,*) 'Going through on regeneration with psb_upd_srch_' case default ! Wrong value info = 8888 diff --git a/src/tools/psb_dspasb.f90 b/src/tools/psb_dspasb.f90 index 37d7cea5..373d5081 100644 --- a/src/tools/psb_dspasb.f90 +++ b/src/tools/psb_dspasb.f90 @@ -125,13 +125,12 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl) 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 + select case(dupl) + case (psb_dupl_ovwrt_,psb_dupl_add_,psb_dupl_err_) dupl_ = dupl - endif + case default + dupl_ = psb_dupl_def_ + end select else dupl_ = psb_dupl_def_ endif diff --git a/src/tools/psb_zspasb.f90 b/src/tools/psb_zspasb.f90 index 4bdb4505..d50d8670 100644 --- a/src/tools/psb_zspasb.f90 +++ b/src/tools/psb_zspasb.f90 @@ -125,18 +125,16 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl) 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 + select case(dupl) + case (psb_dupl_ovwrt_,psb_dupl_add_,psb_dupl_err_) dupl_ = dupl - endif + case default + dupl_ = psb_dupl_def_ + end select else dupl_ = psb_dupl_def_ endif - a%m = n_row a%k = n_col