mirror of
https://github.com/sfilippone/psblas3.git
synced 2026-10-16 10:52:08 +00:00
Fixed bugs in duplicate/regeneration handling.
This commit is contained in:
@@ -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
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -50,7 +50,6 @@ C
|
||||
lb = m + 1
|
||||
end if
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
+20
-14
@@ -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)
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
+241
-12
@@ -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
|
||||
|
||||
|
||||
+40
-28
@@ -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
|
||||
|
||||
+251
-18
@@ -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
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user