fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fixed bug in matrix regeneration stuff.

Browse Source
psblas3-type-indexed
Salvatore Filippone 21 years ago
parent e3741fa0cf
commit 4adaeaf5d5
5 changed files with 292 additions and 225 deletions
Show Stats Download Patch File Download Diff File

11
src/serial/dp/dcocr.f
Unescape Escape View File

@ -74,7 +74,7 @@ C
CALL FCPSB_ERRACTIONSAVE(ERR_ACT) CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
CHECK_FLAG=IBITS(INFO(PSB_UPD_),1,2) CHECK_FLAG=IBITS(INFO(PSB_UPD_),1,2)
c$$$ write(0,*) 'DCOCR FLAG ',info(upd_),check_flag c$$$ write(0,*) 'DCOCR FLAG ',info(psb_upd_),check_flag
IF (TRANS.EQ.'N') THEN IF (TRANS.EQ.'N') THEN
SCALE = (UNITD.EQ.'L') ! meaningless SCALE = (UNITD.EQ.'L') ! meaningless
@ -126,6 +126,9 @@ c$$$ do k=1,nnz
c$$$ write(*,*) k,ia(k),ja(k),ar(k) c$$$ write(*,*) k,ia(k),ja(k),ar(k)
c$$$ enddo c$$$ enddo
c$$$ endif c$$$ endif
c$$$ write(0,*) 'DCOCR Sizes ',lian2,((m+1)+nnz+psb_ireg_flgs_+1),
c$$$ + m+1,nnz,psb_ireg_flgs_,
c$$$ + laux,2*(2+nnz)
if ((lian2.ge.((m+1)+nnz+psb_ireg_flgs_+1)) if ((lian2.ge.((m+1)+nnz+psb_ireg_flgs_+1))
+ .and.(laux.ge.2*(2+nnz))) then + .and.(laux.ge.2*(2+nnz))) then
C C
@ -145,7 +148,10 @@ c
ian2(ip1+psb_nnz_) = 0 ian2(ip1+psb_nnz_) = 0
ian2(ip1+psb_ichk_) = nnz+check_flag ian2(ip1+psb_ichk_) = nnz+check_flag
c$$$ write(0,*)'DCOCR m,ip1,ip2,nnz',m, c$$$ write(0,*)'DCOCR Check: ',ip2,ian2(ip1+psb_iflag_),
c$$$ + ian2(ip1+psb_nnzt_), ian2(ip1+psb_nnz_),
c$$$ + ian2(ip1+psb_ichk_)
c$$$ + ip1,ip2,nnz,ian2(ip1+nnzt_) c$$$ + ip1,ip2,nnz,ian2(ip1+nnzt_)
if (debug) write(0,*) 'Build check :',ian2(ip1+psb_nnzt_) if (debug) write(0,*) 'Build check :',ian2(ip1+psb_nnzt_)
@ -235,6 +241,7 @@ C ... Sum the duplicated element ...
IAN2(ROW+1) = ELEM_CSR IAN2(ROW+1) = ELEM_CSR
ENDDO ENDDO
ELSE ELSE
c$$$ write(0,*) 'Going for ELSE !!!?'
C .... Order with key IA ... C .... Order with key IA ...
CALL MRGSRT(NNZ,IA,AUX,IRET) CALL MRGSRT(NNZ,IA,AUX,IRET)
IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX) IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX)

8
src/serial/psb_cest.f90
Unescape Escape View File

@ -28,14 +28,14 @@
!!$ POSSIBILITY OF SUCH DAMAGE. !!$ POSSIBILITY OF SUCH DAMAGE.
!!$ !!$
!!$ !!$
subroutine psb_cest(afmt, nnz, lia1, lia2, lar, up, info) subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info)
use psb_error_mod use psb_error_mod
use psb_const_mod use psb_const_mod
implicit none implicit none
! .. scalar arguments .. ! .. scalar arguments ..
integer :: nnz, lia1, lia2, lar, info integer :: m,n,nnz, lia1, lia2, lar, info
character :: up character :: up
! .. array arguments.. ! .. array arguments..
character(len=5) :: afmt character(len=5) :: afmt
@ -52,7 +52,7 @@ subroutine psb_cest(afmt, nnz, lia1, lia2, lar, up, info)
if ((up.eq.'y').or.(up.eq.'Y')) then if ((up.eq.'y').or.(up.eq.'Y')) then
if (afmt.eq.'JAD') then if (afmt.eq.'JAD') then
lia1 = 2*(nnz + nnz/5) +1000 lia1 = 2*(nnz + nnz/5) +1000
lia2 = 2*(nnz + nnz/5) +1000 lia2 = 2*(nnz + nnz/5) +1000 +m
lar = nnz + nnz/5 lar = nnz + nnz/5
else if (afmt.eq.'COO') then else if (afmt.eq.'COO') then
lia1 = nnz lia1 = nnz
@ -60,7 +60,7 @@ subroutine psb_cest(afmt, nnz, lia1, lia2, lar, up, info)
lar = nnz lar = nnz
else if(afmt.eq.'CSR') then else if(afmt.eq.'CSR') then
lia1 = nnz lia1 = nnz
lia2 = 2*nnz + 1000 lia2 = 2*nnz + 1000 + m + 1
lar = nnz lar = nnz
else else
info = 136 info = 136

10
src/serial/psb_dcoins.f90
Unescape Escape View File

@ -149,11 +149,11 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info)
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
endif endif
!!$ if (debug) then if (debug) then
!!$ if ((nza - a%ia2(ip1+nnz_)) /= nz) then if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then
!!$ write(0,*) 'PSB_COINS: update discarded items ' write(0,*) 'PSB_COINS: update discarded items '
!!$ end if end if
!!$ end if end if
a%ia2(ip1+psb_nnz_) = nza a%ia2(ip1+psb_nnz_) = nza
else else

91
src/serial/psb_dcsdp.f90
Unescape Escape View File

@ -61,13 +61,31 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
real(kind(1.d0)), allocatable :: work(:) real(kind(1.d0)), allocatable :: work(:)
type(psb_dspmat_type) :: temp_a type(psb_dspmat_type) :: temp_a
Integer :: nzr, ntry, ifc_,ierror, ia1_size,& Integer :: nzr, ntry, ifc_,ierror, ia1_size,&
& ia2_size, aspk_size & ia2_size, aspk_size,size_req,n_row,n_col,iup
integer :: ip1, ip2, nnz, iflag, ichk, nnzt,& integer :: ip1, ip2, nnz, iflag, ichk, nnzt,&
& ipc, i, count, err_act, ierrv(5) & ipc, i, count, err_act, ierrv(5)
character :: check_,trans_,unitd_ character :: check_,trans_,unitd_, up
Integer, Parameter :: maxtry=8 Integer, Parameter :: maxtry=8
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
character(len=20) :: name, ch_err 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
end interface
interface psb_spinfo
subroutine psb_dspinfo(ireq,a,ires,info,iaux)
use psb_spmat_type
type(psb_dspmat_type), intent(in) :: a
integer, intent(in) :: ireq
integer, intent(out) :: ires, info
integer, intent(in), optional :: iaux
end subroutine psb_dspinfo
end interface
name='psb_dcsdp' name='psb_dcsdp'
info = 0 info = 0
@ -133,6 +151,45 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
! ...matrix conversion... ! ...matrix conversion...
b%m=a%m b%m=a%m
b%k=a%k 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
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,&
& ia1_size, ia2_size, aspk_size, up,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
if ((size(b%aspk) < aspk_size) .or.&
&(size(b%ia1) < ia1_size) .or.&
&(size(b%ia2) < ia2_size) .or.&
&(size(b%pl) < b%m) .or.&
&(size(b%pr) < b%k )) then
call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info)
endif
if (info /= no_err) then
info=4010
ch_err='psb_sp_reall'
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 (a%fida(1:3))
case ('CSR') case ('CSR')
@ -141,11 +198,11 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
case ('CSR') case ('CSR')
!!$
ia1_size=a%infoa(psb_nnz_) !!$ ia1_size=a%infoa(psb_nnz_)
ia2_size=a%m+1 !!$ ia2_size=a%m+1
aspk_size=a%infoa(psb_nnz_) !!$ aspk_size=a%infoa(psb_nnz_)
call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info) !!$ 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,& 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,& & a%ia1, a%ia2, a%infoa, b%pl, b%descra, b%aspk, b%ia1,&
@ -164,10 +221,10 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
!...converting to JAD !...converting to JAD
!...output matrix may not be big enough !...output matrix may not be big enough
ia1_size=a%infoa(psb_nnz_) !!$ ia1_size=a%infoa(psb_nnz_)
ia2_size=a%m+1 !!$ ia2_size=a%m+1
aspk_size=a%infoa(psb_nnz_) !!$ aspk_size=a%infoa(psb_nnz_)
call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info) !!$ call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info)
do do
call dcrjd(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& call dcrjd(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,&
@ -207,8 +264,8 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
case ('COO') case ('COO')
aspk_size=max(size(a%aspk),a%ia2(a%m+1)) !!$ aspk_size=max(size(a%aspk),a%ia2(a%m+1))
call psb_sp_reall(b,aspk_size,info) !!$ call psb_sp_reall(b,aspk_size,info)
!!$ write(0,*) 'From DCSDP90:',b%fida,size(b%aspk),info !!$ write(0,*) 'From DCSDP90:',b%fida,size(b%aspk),info
call dcrco(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& 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,& & a%ia1, a%ia2, a%infoa, b%pl, b%descra, b%aspk, b%ia1,&
@ -228,8 +285,8 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
case ('CSR') case ('CSR')
aspk_size=max(size(a%aspk),a%ia2(a%m+1)) !!$ aspk_size=max(size(a%aspk),a%ia2(a%m+1))
call psb_sp_reall(b,aspk_size,info) !!$ call psb_sp_reall(b,aspk_size,info)
call dcocr(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& 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,& & 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),& & b%ia2, b%infoa, b%pr, size(b%aspk), size(b%ia1),&
@ -301,8 +358,8 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd)
case ('COO') case ('COO')
aspk_size=max(size(a%aspk),a%ia2(a%m+1)) !!$ aspk_size=max(size(a%aspk),a%ia2(a%m+1))
call psb_sp_reall(b,aspk_size,info) !!$ call psb_sp_reall(b,aspk_size,info)
call dcoco(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& 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,& & 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),& & b%ia2, b%infoa, b%pr, size(b%aspk), size(b%ia1),&

49
src/tools/psb_dspasb.f90
Unescape Escape View File

@ -54,10 +54,11 @@ subroutine psb_dspasb(a,desc_a, info, afmt, up, dup)
implicit none implicit none
interface psb_cest interface psb_cest
subroutine psb_cest(afmt, nnz, lia1, lia2, lar, up, info) subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info)
integer, intent(in) :: nnz integer, intent(in) :: m,n,nnz
integer, intent(out) :: lia1, lia2, lar, info integer, intent(out) :: lia1, lia2, lar, info
character, intent(in) :: afmt*5, up character, intent(inout) :: afmt*5
character, intent(in) :: up
end subroutine psb_cest end subroutine psb_cest
end interface end interface
@ -179,26 +180,28 @@ subroutine psb_dspasb(a,desc_a, info, afmt, up, dup)
! work area requested must be fixed to ! work area requested must be fixed to
! No of Grid'd processes and NNZ+2 ! No of Grid'd processes and NNZ+2
! !
size_req = max(a%infoa(psb_nnz_),1)+3 !!$ size_req = max(a%infoa(psb_nnz_),1)+3
if (debug) write(0,*) 'DCSDP : size_req 1:',size_req !!$ if (debug) write(0,*) 'DCSDP : size_req 1:',size_req
call psb_cest(a%fida, size_req, ia1_size, ia2_size, aspk_size, iup,info) !!$ call psb_cest(a%fida, n_row,n_col,size_req,&
if (info /= no_err) then !!$ & ia1_size, ia2_size, aspk_size, iup,info)
info=4010 !!$ write(0,*) 'ESTIMATE : ',ia1_size,ia2_size,aspk_Size,iup
ch_err='psb_cest' !!$ if (info /= no_err) then
call psb_errpush(info,name,a_err=ch_err) !!$ info=4010
goto 9999 !!$ ch_err='psb_cest'
endif !!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
call psb_sp_reall(a,ia1_size,ia2_size,aspk_size,info) !!$ endif
if (info /= no_err) then !!$
info=4010 !!$ call psb_sp_reall(a,ia1_size,ia2_size,aspk_size,info)
ch_err='psb_sp_reall' !!$ if (info /= no_err) then
call psb_errpush(info,name,a_err=ch_err) !!$ info=4010
goto 9999 !!$ ch_err='psb_sp_reall'
endif !!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
a%pl(:) = 0 !!$ endif
a%pr(:) = 0 !!$
!!$ a%pl(:) = 0
!!$ a%pr(:) = 0
if (debugwrt) then if (debugwrt) then
iout = 30+myrow iout = 30+myrow

Write Preview
Loading…
Cancel
Save