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

First round of changes: fixed SYMBMM,NUMBMM to use GETROW.

Browse Source 
THis has to be updated in two ways:
1. the symbmm/numbmm equivalent must be made smarter (perhaps a cache
for rows extracted from B?) 
2. the whole getrow/getblk/clip chain must be changed.
psblas3-type-indexed
Salvatore Filippone 18 years ago
parent 3cd8e0722d
commit f56e369210
13 changed files with 98 additions and 182 deletions
Show Stats Download Patch File Download Diff File

4
base/modules/psb_error_mod.F90
Unescape Escape View File

@ -80,8 +80,8 @@ contains
! also changes error action to "return" ! also changes error action to "return"
subroutine psb_erractionsave(err_act) subroutine psb_erractionsave(err_act)
integer, intent(out) :: err_act integer, intent(out) :: err_act
err_act=err_action err_act = err_action
err_action=act_ret err_action = act_ret
end subroutine psb_erractionsave end subroutine psb_erractionsave

11
base/modules/psb_serial_mod.f90
Unescape Escape View File

@ -136,6 +136,14 @@ module psb_serial_mod
end subroutine psb_zcsmm end subroutine psb_zcsmm
end interface end interface
interface psb_cest
subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, iup, info)
integer, intent(in) :: m,n,nnz,iup
integer, intent(out) :: lia1, lia2, lar, info
character, intent(inout) :: afmt*5
end subroutine psb_cest
end interface
interface psb_fixcoo interface psb_fixcoo
subroutine psb_dfixcoo(a,info,idir) subroutine psb_dfixcoo(a,info,idir)
use psb_spmat_type use psb_spmat_type
@ -417,7 +425,7 @@ module psb_serial_mod
end interface end interface
interface psb_sp_getrow interface psb_sp_getrow
subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw) subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw,bw)
use psb_spmat_type use psb_spmat_type
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
integer, intent(in) :: irw integer, intent(in) :: irw
@ -427,6 +435,7 @@ module psb_serial_mod
integer, intent(in), target, optional :: iren(:) integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw integer, intent(in), optional :: lrw
integer, intent(out) :: info integer, intent(out) :: info
type(psb_dspmat_type), intent(inout), optional, target :: bw
end subroutine psb_dspgetrow end subroutine psb_dspgetrow
subroutine psb_zspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw) subroutine psb_zspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
use psb_spmat_type use psb_spmat_type

10
base/serial/psb_dcsdp.f90
Unescape Escape View File

@ -48,7 +48,7 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl)
use psb_error_mod use psb_error_mod
use psb_spmat_type use psb_spmat_type
use psb_string_mod use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_dcsdp
implicit none implicit none
!....Parameters... !....Parameters...
Type(psb_dspmat_type), intent(in) :: A Type(psb_dspmat_type), intent(in) :: A
@ -70,14 +70,6 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl)
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, iup, info)
integer, intent(in) :: m,n,nnz,iup
integer, intent(out) :: lia1, lia2, lar, info
character, intent(inout) :: afmt*5
end subroutine psb_cest
end interface
name='psb_csdp' name='psb_csdp'
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)

2
base/serial/psb_dcsrws.f90
Unescape Escape View File

@ -42,7 +42,7 @@ subroutine psb_dcsrws(rw,a,info,trans)
integer :: info integer :: info
character, optional :: trans character, optional :: trans
Interface dcsrws Interface
subroutine dcsrws(trans,m,n,fida,descra,a,ia1,ia2,& subroutine dcsrws(trans,m,n,fida,descra,a,ia1,ia2,&
& infoa,rowsum,ierror) & infoa,rowsum,ierror)
integer, intent(in) :: m,n integer, intent(in) :: m,n

28
base/serial/psb_dnumbmm.f90
Unescape Escape View File

@ -41,26 +41,13 @@
subroutine psb_dnumbmm(a,b,c) subroutine psb_dnumbmm(a,b,c)
use psb_spmat_type use psb_spmat_type
use psb_serial_mod, psb_protect_name => psb_dnumbmm
implicit none implicit none
type(psb_dspmat_type) :: a,b,c type(psb_dspmat_type) :: a,b,c
real(kind(1.d0)), allocatable :: temp(:) real(kind(1.d0)), allocatable :: temp(:)
integer :: info integer :: info
logical :: csra, csrb
interface psb_sp_getrow
subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
use psb_spmat_type
type(psb_dspmat_type), intent(in) :: a
integer, intent(in) :: irw
integer, intent(out) :: nz
integer, intent(inout) :: ia(:), ja(:)
real(kind(1.d0)), intent(inout) :: val(:)
integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw
integer, intent(out) :: info
end subroutine psb_dspgetrow
end interface
allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info) allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info)
if (info /= 0) then if (info /= 0) then
@ -71,7 +58,10 @@ subroutine psb_dnumbmm(a,b,c)
! Note: we still have to test about possible performance hits. ! Note: we still have to test about possible performance hits.
! !
! !
if (.true.) then csra = (toupper(a%fida(1:3))=='CSR')
csrb = (toupper(b%fida(1:3))=='CSR')
if (csra.and.csrb) then
call numbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,& call numbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,&
& b%ia2,b%ia1,0,b%aspk,& & b%ia2,b%ia1,0,b%aspk,&
& c%ia2,c%ia1,0,c%aspk,temp) & c%ia2,c%ia1,0,c%aspk,temp)
@ -94,7 +84,7 @@ contains
real(kind(1.d0)), allocatable :: aval(:),bval(:) real(kind(1.d0)), allocatable :: aval(:),bval(:)
integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,ii,minlm,minmn,minln integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,ii,minlm,minmn,minln
real(kind(1.d0)) :: ajj real(kind(1.d0)) :: ajj
type(psb_dspmat_type) :: w
n = a%m n = a%m
m = a%k m = a%k
@ -114,7 +104,7 @@ contains
minmn = min(m,n) minmn = min(m,n)
do i = 1,n do i = 1,n
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info) call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info,bw=w)
do jj=1, nazr do jj=1, nazr
j=iacl(jj) j=iacl(jj)
@ -125,7 +115,7 @@ contains
return return
endif endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info) call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info,bw=w)
do k=1,nbzr do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn

51
base/serial/psb_dspgetrow.f90
Unescape Escape View File

@ -39,7 +39,7 @@
!* format. * !* format. *
!* * !* *
!***************************************************************************** !*****************************************************************************
subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw) subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw,bw)
use psb_spmat_type use psb_spmat_type
use psb_string_mod use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_dspgetrow use psb_serial_mod, psb_protect_name => psb_dspgetrow
@ -52,19 +52,22 @@ subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
integer, intent(in), target, optional :: iren(:) integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw integer, intent(in), optional :: lrw
integer, intent(out) :: info integer, intent(out) :: info
type(psb_dspmat_type), intent(inout), optional, target :: bw
integer :: lrw_, ierr(5), err_act integer :: lrw_, ierr(5), err_act
type(psb_dspmat_type) :: b type(psb_dspmat_type), target :: b
type(psb_dspmat_type), pointer :: b_
integer, pointer :: iren_(:) integer, pointer :: iren_(:)
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
name='psb_sp_getrow' name = 'psb_sp_getrow'
info = 0 info = 0
call psb_erractionsave(err_act) !!$ call psb_erractionsave(err_act)
call psb_set_erraction(0) !!$ call psb_set_erraction(0)
call psb_nullify_sp(b)
if (present(lrw)) then if (present(lrw)) then
lrw_ = lrw lrw_ = lrw
@ -75,13 +78,26 @@ subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
write(0,*) 'SPGETROW input error: fixing lrw',irw,lrw_ write(0,*) 'SPGETROW input error: fixing lrw',irw,lrw_
lrw_ = irw lrw_ = irw
end if end if
call psb_sp_all(lrw_-irw+1,lrw_-irw+1,b,info) if (present(bw)) then
b_ => bw
else
b_ => b
end if
call psb_nullify_sp(b_)
if (.not.(allocated(b_%aspk).and.allocated(b_%ia1).and.&
& allocated(b_%ia2))) then
write(0,*) 'First allocation for B in SPGETROW'
call psb_sp_all(lrw_-irw+1,lrw_-irw+1,b_,info)
end if
if (present(iren)) then if (present(iren)) then
call psb_sp_getblk(irw,a,b,info,iren=iren,lrw=lrw_) call psb_sp_getblk(irw,a,b_,info,iren=iren,lrw=lrw_)
else else
call psb_sp_getblk(irw,a,b,info,lrw=lrw_) call psb_sp_getblk(irw,a,b_,info,lrw=lrw_)
end if end if
if (info /= 0) then if (info /= 0) then
info=136 info=136
ch_err=a%fida(1:3) ch_err=a%fida(1:3)
@ -89,17 +105,17 @@ subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
goto 9999 goto 9999
end if end if
if (toupper(b%fida) /= 'COO') then if (toupper(b_%fida) /= 'COO') then
info=4010 info=4010
ch_err=a%fida(1:3) ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
endif endif
nz = b%infoa(psb_nnz_) nz = b_%infoa(psb_nnz_)
if (size(ia)>= nz) then if (size(ia)>= nz) then
ia(1:nz) = b%ia1(1:nz) ia(1:nz) = b_%ia1(1:nz)
else else
info = 135 info = 135
ierr(1) = 4 ierr(1) = 4
@ -109,7 +125,7 @@ subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
endif endif
if (size(ja)>= nz) then if (size(ja)>= nz) then
ja(1:nz) = b%ia2(1:nz) ja(1:nz) = b_%ia2(1:nz)
else else
info = 135 info = 135
ierr(1) = 5 ierr(1) = 5
@ -119,7 +135,7 @@ subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
endif endif
if (size(val)>= nz) then if (size(val)>= nz) then
val(1:nz) = b%aspk(1:nz) val(1:nz) = b_%aspk(1:nz)
else else
info = 135 info = 135
ierr(1) = 6 ierr(1) = 6
@ -129,13 +145,14 @@ subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
endif endif
call psb_sp_free(b,info) !!$ call psb_sp_free(b,info)
call psb_erractionrestore(err_act) !!$ call psb_erractionrestore(err_act)
return return
9999 continue 9999 continue
call psb_erractionrestore(err_act) !!$ call psb_erractionrestore(err_act)
call psb_erractionsave(err_act)
if (err_act.eq.psb_act_abort_) then if (err_act.eq.psb_act_abort_) then
call psb_error() call psb_error()
return return

15
base/serial/psb_dspgtblk.f90
Unescape Escape View File

@ -61,7 +61,7 @@ subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
name='psb_spgtblk' name='psb_spgtblk'
info = 0 info = 0
call psb_erractionsave(err_act) !!$ call psb_erractionsave(err_act)
irw_ = irw irw_ = irw
if (present(lrw)) then if (present(lrw)) then
@ -79,7 +79,7 @@ subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
append_=.false. append_=.false.
endif endif
if (present(iren)) then if (present(iren)) then
iren_=>iren iren_ => iren
else else
iren_ => null() iren_ => null()
end if end if
@ -110,11 +110,12 @@ subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
end if end if
call psb_erractionrestore(err_act) !!$ call psb_erractionrestore(err_act)
return return
9999 continue 9999 continue
call psb_erractionrestore(err_act) !!$ call psb_erractionrestore(err_act)
call psb_erractionsave(err_act)
if (err_act.eq.psb_act_abort_) then if (err_act.eq.psb_act_abort_) then
call psb_error() call psb_error()
return return
@ -156,6 +157,7 @@ contains
end do end do
if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then
!!$ write(0,*) 'Realloc in gtblk 1',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
call psb_sp_reall(b,nzb+nz,iret) call psb_sp_reall(b,nzb+nz,iret)
endif endif
@ -197,6 +199,7 @@ contains
nz = a%ia2(idx+nr) - a%ia2(idx) nz = a%ia2(idx+nr) - a%ia2(idx)
if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then
!!$ write(0,*) 'Realloc in gtblk 2',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
call psb_sp_reall(b,nzb+nz,iret) call psb_sp_reall(b,nzb+nz,iret)
endif endif
b%fida='COO' b%fida='COO'
@ -310,6 +313,7 @@ contains
! Now do the copy. ! Now do the copy.
nz = jp - ip +1 nz = jp - ip +1
if (size(b%ia1) < nzb+nz) then if (size(b%ia1) < nzb+nz) then
!!$ write(0,*) 'Realloc in gtblk 3',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
call psb_sp_reall(b,nzb+nz,iret) call psb_sp_reall(b,nzb+nz,iret)
endif endif
b%fida='COO' b%fida='COO'
@ -335,6 +339,7 @@ contains
nz = (nza*(lrw-irw+1))/max(a%m,1) nz = (nza*(lrw-irw+1))/max(a%m,1)
if (size(b%ia1) < nzb+nz) then if (size(b%ia1) < nzb+nz) then
!!$ write(0,*) 'Realloc in gtblk 4',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
call psb_sp_reall(b,nzb+nz,iret) call psb_sp_reall(b,nzb+nz,iret)
endif endif
@ -345,6 +350,7 @@ contains
k = k + 1 k = k + 1
if (k > nz) then if (k > nz) then
nz = k nz = k
!!$ write(0,*) 'Realloc in gtblk 5',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
call psb_sp_reall(b,nzb+nz,iret) call psb_sp_reall(b,nzb+nz,iret)
end if end if
b%aspk(nzb+k) = a%aspk(i) b%aspk(nzb+k) = a%aspk(i)
@ -359,6 +365,7 @@ contains
k = k + 1 k = k + 1
if (k > nz) then if (k > nz) then
nz = k nz = k
!!$ write(0,*) 'Realloc in gtblk 6',size(b%ia1),size(b%ia2),size(b%aspk),nzb,nz
call psb_sp_reall(b,nzb+nz,iret) call psb_sp_reall(b,nzb+nz,iret)
end if end if
b%aspk(nzb+k) = a%aspk(i) b%aspk(nzb+k) = a%aspk(i)

21
base/serial/psb_dspgtdiag.f90
Unescape Escape View File

@ -34,9 +34,6 @@
!***************************************************************************** !*****************************************************************************
!* * !* *
!* Takes a specified row from matrix A and copies into matrix B (possibly *
!* appending to B). Output is always COO. Input might be anything, once *
!* we get to actually write the code..... *
!* * !* *
!***************************************************************************** !*****************************************************************************
subroutine psb_dspgtdiag(a,d,info) subroutine psb_dspgtdiag(a,d,info)
@ -45,25 +42,13 @@ subroutine psb_dspgtdiag(a,d,info)
use psb_spmat_type use psb_spmat_type
use psb_error_mod use psb_error_mod
use psb_const_mod use psb_const_mod
use psb_serial_mod, psb_protect_name => psb_dspgtdiag
implicit none implicit none
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
real(kind(1.d0)), intent(inout) :: d(:) real(kind(1.d0)), intent(inout) :: d(:)
integer, intent(out) :: info integer, intent(out) :: info
interface psb_spgtblk
subroutine psb_dspgtblk(irw,a,b,info,append,iren,lrw)
use psb_spmat_type
type(psb_dspmat_type), intent(in) :: a
integer, intent(in) :: irw
type(psb_dspmat_type), intent(inout) :: b
logical, intent(in), optional :: append
integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw
integer, intent(out) :: info
end subroutine psb_dspgtblk
end interface
type(psb_dspmat_type) :: tmpa type(psb_dspmat_type) :: tmpa
integer :: i,j,k,nr, nz, err_act, ii, rng, irb, nrb integer :: i,j,k,nr, nz, err_act, ii, rng, irb, nrb
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -102,10 +87,10 @@ subroutine psb_dspgtdiag(a,d,info)
write(0,*)'in spgtdiag' write(0,*)'in spgtdiag'
do i=1, rng, nrb do i=1, rng, nrb
irb=min(i+nrb-1,rng) irb=min(i+nrb-1,rng)
call psb_spgtblk(i,a,tmpa,info,lrw=irb) call psb_sp_getblk(i,a,tmpa,info,lrw=irb)
if(info.ne.0) then if(info.ne.0) then
info=4010 info=4010
ch_err='psb_spgtblk' ch_err='psb_sp_getblk'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if

46
base/serial/psb_dsymbmm.f90
Unescape Escape View File

@ -41,7 +41,7 @@
subroutine psb_dsymbmm(a,b,c,info) subroutine psb_dsymbmm(a,b,c,info)
use psb_spmat_type use psb_spmat_type
use psb_string_mod use psb_string_mod
use psb_serial_mod, only : psb_msort use psb_serial_mod, psb_protect_name => psb_dsymbmm
implicit none implicit none
type(psb_dspmat_type) :: a,b,c type(psb_dspmat_type) :: a,b,c
@ -56,41 +56,14 @@ subroutine psb_dsymbmm(a,b,c,info)
integer, allocatable :: ic(:),jc(:) integer, allocatable :: ic(:),jc(:)
end subroutine symbmm end subroutine symbmm
end interface end interface
interface psb_sp_getrow character(len=20) :: name, ch_err
subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw) integer :: err_act
use psb_spmat_type logical :: csra, csrb
type(psb_dspmat_type), intent(in) :: a
integer, intent(in) :: irw
integer, intent(out) :: nz
integer, intent(inout) :: ia(:), ja(:)
real(kind(1.d0)), intent(inout) :: val(:)
integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw
integer, intent(out) :: info
end subroutine psb_dspgetrow
end interface
character(len=20) :: name, ch_err
integer :: err_act
name='psb_symbmm' name='psb_symbmm'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select case(toupper(a%fida(1:3))) csra = (toupper(a%fida(1:3))=='CSR')
case ('CSR') csrb = (toupper(b%fida(1:3))=='CSR')
case default
info=135
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
select case(toupper(b%fida(1:3)))
case ('CSR')
case default
info=136
ch_err=b%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
if (b%m /= a%k) then if (b%m /= a%k) then
write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k
@ -105,7 +78,7 @@ subroutine psb_dsymbmm(a,b,c,info)
! Note: we need to test whether there is a performance impact ! Note: we need to test whether there is a performance impact
! in not using the original Douglas & Bank code. ! in not using the original Douglas & Bank code.
! !
if (.true.) then if (csra.and.csrb) then
call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,& call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,&
& b%ia2,b%ia1,0,& & b%ia2,b%ia1,0,&
& c%ia2,c%ia1,0,itemp) & c%ia2,c%ia1,0,itemp)
@ -140,6 +113,7 @@ contains
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
real(kind(1.d0)), allocatable :: aval(:),bval(:) real(kind(1.d0)), allocatable :: aval(:),bval(:)
integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,ii,minlm,minmn integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,ii,minlm,minmn
type(psb_dspmat_type) :: w
n = a%m n = a%m
@ -169,7 +143,7 @@ contains
main: do i=1,n main: do i=1,n
istart=-1 istart=-1
length=0 length=0
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info) call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info,bw=w)
do jj=1, nazr do jj=1, nazr
j=iacl(jj) j=iacl(jj)
@ -179,7 +153,7 @@ contains
info = 1 info = 1
return return
endif endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info) call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info,bw=w)
do k=1,nbzr do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn

10
base/serial/psb_zcsdp.f90
Unescape Escape View File

@ -48,7 +48,7 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl)
use psb_error_mod use psb_error_mod
use psb_spmat_type use psb_spmat_type
use psb_string_mod use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_zcsdp
implicit none implicit none
!....Parameters... !....Parameters...
Type(psb_zspmat_type), intent(in) :: A Type(psb_zspmat_type), intent(in) :: A
@ -70,14 +70,6 @@ subroutine psb_zcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl)
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, iup, info)
integer, intent(in) :: m,n,nnz,iup
integer, intent(out) :: lia1, lia2, lar, info
character, intent(inout) :: afmt*5
end subroutine psb_cest
end interface
name='psb_csdp' name='psb_csdp'
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)

24
base/serial/psb_znumbmm.f90
Unescape Escape View File

@ -28,7 +28,7 @@
!!$ POSSIBILITY OF SUCH DAMAGE. !!$ POSSIBILITY OF SUCH DAMAGE.
!!$ !!$
!!$ !!$
! File: psb_znumbmm.f90 ! File: psb_dnumbmm.f90
! Subroutine: ! Subroutine:
! Parameters: ! Parameters:
! !
@ -41,26 +41,13 @@
subroutine psb_znumbmm(a,b,c) subroutine psb_znumbmm(a,b,c)
use psb_spmat_type use psb_spmat_type
use psb_serial_mod, psb_protect_name => psb_znumbmm
implicit none implicit none
type(psb_zspmat_type) :: a,b,c type(psb_zspmat_type) :: a,b,c
complex(kind(1.d0)), allocatable :: temp(:) complex(kind(1.d0)), allocatable :: temp(:)
integer :: info integer :: info
logical :: csra, csrb
interface psb_sp_getrow
subroutine psb_zspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
use psb_spmat_type
type(psb_zspmat_type), intent(in) :: a
integer, intent(in) :: irw
integer, intent(out) :: nz
integer, intent(inout) :: ia(:), ja(:)
complex(kind(1.d0)), intent(inout) :: val(:)
integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw
integer, intent(out) :: info
end subroutine psb_zspgetrow
end interface
allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info) allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info)
if (info /= 0) then if (info /= 0) then
@ -71,7 +58,10 @@ subroutine psb_znumbmm(a,b,c)
! Note: we still have to test about possible performance hits. ! Note: we still have to test about possible performance hits.
! !
! !
if (.true.) then csra = (toupper(a%fida(1:3))=='CSR')
csrb = (toupper(b%fida(1:3))=='CSR')
if (csra.and.csrb) then
call znumbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,& call znumbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,&
& b%ia2,b%ia1,0,b%aspk,& & b%ia2,b%ia1,0,b%aspk,&
& c%ia2,c%ia1,0,c%aspk,temp) & c%ia2,c%ia1,0,c%aspk,temp)

19
base/serial/psb_zspgtdiag.f90
Unescape Escape View File

@ -34,9 +34,6 @@
!***************************************************************************** !*****************************************************************************
!* * !* *
!* Takes a specified row from matrix A and copies into matrix B (possibly *
!* appending to B). Output is always COO. Input might be anything, once *
!* we get to actually write the code..... *
!* * !* *
!***************************************************************************** !*****************************************************************************
subroutine psb_zspgtdiag(a,d,info) subroutine psb_zspgtdiag(a,d,info)
@ -45,25 +42,13 @@ subroutine psb_zspgtdiag(a,d,info)
use psb_spmat_type use psb_spmat_type
use psb_error_mod use psb_error_mod
use psb_const_mod use psb_const_mod
use psb_serial_mod, psb_protect_name => psb_zspgtdiag
implicit none implicit none
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in) :: a
complex(kind(1.d0)), intent(inout) :: d(:) complex(kind(1.d0)), intent(inout) :: d(:)
integer, intent(out) :: info integer, intent(out) :: info
interface psb_spgtblk
subroutine psb_zspgtblk(irw,a,b,info,append,iren,lrw)
use psb_spmat_type
type(psb_zspmat_type), intent(in) :: a
integer, intent(in) :: irw
type(psb_zspmat_type), intent(inout) :: b
logical, intent(in), optional :: append
integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw
integer, intent(out) :: info
end subroutine psb_zspgtblk
end interface
type(psb_zspmat_type) :: tmpa type(psb_zspmat_type) :: tmpa
integer :: i,j,k,nr, nz, err_act, ii, rng, irb, nrb integer :: i,j,k,nr, nz, err_act, ii, rng, irb, nrb
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -102,7 +87,7 @@ subroutine psb_zspgtdiag(a,d,info)
write(0,*)'in spgtdiag' write(0,*)'in spgtdiag'
do i=1, rng, nrb do i=1, rng, nrb
irb=min(i+nrb-1,rng) irb=min(i+nrb-1,rng)
call psb_spgtblk(i,a,tmpa,info,lrw=irb) call psb_sp_getblk(i,a,tmpa,info,lrw=irb)
if(info.ne.0) then if(info.ne.0) then
info=4010 info=4010
ch_err='psb_spgtblk' ch_err='psb_spgtblk'

37
base/serial/psb_zsymbmm.f90
Unescape Escape View File

@ -41,6 +41,7 @@
subroutine psb_zsymbmm(a,b,c,info) subroutine psb_zsymbmm(a,b,c,info)
use psb_spmat_type use psb_spmat_type
use psb_string_mod use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_zsymbmm
implicit none implicit none
type(psb_zspmat_type) :: a,b,c type(psb_zspmat_type) :: a,b,c
@ -56,40 +57,14 @@ subroutine psb_zsymbmm(a,b,c,info)
end subroutine symbmm end subroutine symbmm
end interface end interface
interface psb_sp_getrow
subroutine psb_zspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
use psb_spmat_type
type(psb_zspmat_type), intent(in) :: a
integer, intent(in) :: irw
integer, intent(out) :: nz
integer, intent(inout) :: ia(:), ja(:)
complex(kind(1.d0)), intent(inout) :: val(:)
integer, intent(in), target, optional :: iren(:)
integer, intent(in), optional :: lrw
integer, intent(out) :: info
end subroutine psb_zspgetrow
end interface
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
integer :: err_act integer :: err_act
logical :: csra, csrb
name='psb_symbmm' name='psb_symbmm'
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select case(toupper(a%fida(1:3))) csra = (toupper(a%fida(1:3))=='CSR')
case ('CSR') csrb = (toupper(b%fida(1:3))=='CSR')
case default
info=135
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
select case(toupper(b%fida(1:3)))
case ('CSR')
case default
info=136
ch_err=b%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
if (b%m /= a%k) then if (b%m /= a%k) then
write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k
@ -100,12 +75,11 @@ subroutine psb_zsymbmm(a,b,c,info)
endif endif
nze = max(a%m+1,2*a%m) nze = max(a%m+1,2*a%m)
call psb_sp_reall(c,nze,info) call psb_sp_reall(c,nze,info)
! !
! Note: we need to test whether there is a performance impact ! Note: we need to test whether there is a performance impact
! in not using the original Douglas & Bank code. ! in not using the original Douglas & Bank code.
! !
if (.true.) then if (csra.and.csrb) then
call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,& call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,&
& b%ia2,b%ia1,0,& & b%ia2,b%ia1,0,&
& c%ia2,c%ia1,0,itemp) & c%ia2,c%ia1,0,itemp)
@ -113,6 +87,7 @@ subroutine psb_zsymbmm(a,b,c,info)
call inner_symbmm(a,b,c,itemp,info) call inner_symbmm(a,b,c,itemp,info)
endif endif
call psb_realloc(size(c%ia1),c%aspk,info) call psb_realloc(size(c%ia1),c%aspk,info)
c%pl(1) = 0 c%pl(1) = 0
c%pr(1) = 0 c%pr(1) = 0
c%m=a%m c%m=a%m

Write Preview
Loading…
Cancel
Save