fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-20 12:52:14 +00:00
Code Issues Packages Projects Releases Wiki Activity

psblas3:

Browse Source 
base/modules/psb_serial_mod.f90
 base/serial/Makefile
 base/serial/psb_cnumbmm.f90
 base/serial/psb_crwextd.f90
 base/serial/psb_csymbmm.f90
 base/serial/psb_dsymbmm.f90
 base/serial/psb_snumbmm.f90
 base/serial/psb_srwextd.f90
 base/serial/psb_ssymbmm.f90
 base/serial/psb_znumbmm.f90
 base/serial/psb_zrwextd.f90
 base/serial/psb_zsymbmm.f90

Fixed interface/implementation for rwextd and smmp for S/C/Z. Now
ready for more work on MLD.
This commit is contained in:
Salvatore Filippone
2009-11-09 13:06:09 +00:00
parent d8cf0466bf
commit c31e742f34
12 changed files with 1661 additions and 719 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/modules/psb_serial_mod.f90
+133
View File
@@ -43,6 +43,20 @@ module psb_serial_mod
use psb_mat_mod
interface psb_symbmm
subroutine psb_ssymbmm(a,b,c,info)
use psb_mat_mod
implicit none
type(psb_s_sparse_mat), intent(in) :: a,b
type(psb_s_sparse_mat), intent(out) :: c
integer, intent(out) :: info
end subroutine psb_ssymbmm
subroutine psb_sbase_symbmm(a,b,c,info)
use psb_mat_mod
implicit none
class(psb_s_base_sparse_mat), intent(in) :: a,b
type(psb_s_csr_sparse_mat), intent(out) :: c
integer, intent(out) :: info
end subroutine psb_sbase_symbmm
subroutine psb_dsymbmm(a,b,c,info)
use psb_mat_mod
implicit none
@@ -57,8 +71,49 @@ module psb_serial_mod
type(psb_d_csr_sparse_mat), intent(out) :: c
integer, intent(out) :: info
end subroutine psb_dbase_symbmm
subroutine psb_csymbmm(a,b,c,info)
use psb_mat_mod
implicit none
type(psb_c_sparse_mat), intent(in) :: a,b
type(psb_c_sparse_mat), intent(out) :: c
integer, intent(out) :: info
end subroutine psb_csymbmm
subroutine psb_cbase_symbmm(a,b,c,info)
use psb_mat_mod
implicit none
class(psb_c_base_sparse_mat), intent(in) :: a,b
type(psb_c_csr_sparse_mat), intent(out) :: c
integer, intent(out) :: info
end subroutine psb_cbase_symbmm
subroutine psb_zsymbmm(a,b,c,info)
use psb_mat_mod
implicit none
type(psb_z_sparse_mat), intent(in) :: a,b
type(psb_z_sparse_mat), intent(out) :: c
integer, intent(out) :: info
end subroutine psb_zsymbmm
subroutine psb_zbase_symbmm(a,b,c,info)
use psb_mat_mod
implicit none
class(psb_z_base_sparse_mat), intent(in) :: a,b
type(psb_z_csr_sparse_mat), intent(out) :: c
integer, intent(out) :: info
end subroutine psb_zbase_symbmm
end interface
interface psb_numbmm
subroutine psb_snumbmm(a,b,c)
use psb_mat_mod
implicit none
type(psb_s_sparse_mat), intent(in) :: a,b
type(psb_s_sparse_mat), intent(inout) :: c
end subroutine psb_snumbmm
subroutine psb_sbase_numbmm(a,b,c)
use psb_mat_mod
implicit none
class(psb_s_base_sparse_mat), intent(in) :: a,b
type(psb_s_csr_sparse_mat), intent(inout) :: c
end subroutine psb_sbase_numbmm
subroutine psb_dnumbmm(a,b,c)
use psb_mat_mod
implicit none
@@ -71,9 +126,51 @@ module psb_serial_mod
class(psb_d_base_sparse_mat), intent(in) :: a,b
type(psb_d_csr_sparse_mat), intent(inout) :: c
end subroutine psb_dbase_numbmm
subroutine psb_cnumbmm(a,b,c)
use psb_mat_mod
implicit none
type(psb_c_sparse_mat), intent(in) :: a,b
type(psb_c_sparse_mat), intent(inout) :: c
end subroutine psb_cnumbmm
subroutine psb_cbase_numbmm(a,b,c)
use psb_mat_mod
implicit none
class(psb_c_base_sparse_mat), intent(in) :: a,b
type(psb_c_csr_sparse_mat), intent(inout) :: c
end subroutine psb_cbase_numbmm
subroutine psb_znumbmm(a,b,c)
use psb_mat_mod
implicit none
type(psb_z_sparse_mat), intent(in) :: a,b
type(psb_z_sparse_mat), intent(inout) :: c
end subroutine psb_znumbmm
subroutine psb_zbase_numbmm(a,b,c)
use psb_mat_mod
implicit none
class(psb_z_base_sparse_mat), intent(in) :: a,b
type(psb_z_csr_sparse_mat), intent(inout) :: c
end subroutine psb_zbase_numbmm
end interface
interface psb_rwextd
subroutine psb_srwextd(nr,a,info,b,rowscale)
use psb_mat_mod
implicit none
integer, intent(in) :: nr
type(psb_s_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
type(psb_s_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
end subroutine psb_srwextd
subroutine psb_sbase_rwextd(nr,a,info,b,rowscale)
use psb_mat_mod
implicit none
integer, intent(in) :: nr
class(psb_s_base_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
class(psb_s_base_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
end subroutine psb_sbase_rwextd
subroutine psb_drwextd(nr,a,info,b,rowscale)
use psb_mat_mod
implicit none
@@ -92,6 +189,42 @@ module psb_serial_mod
class(psb_d_base_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
end subroutine psb_dbase_rwextd
subroutine psb_crwextd(nr,a,info,b,rowscale)
use psb_mat_mod
implicit none
integer, intent(in) :: nr
type(psb_c_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
type(psb_c_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
end subroutine psb_crwextd
subroutine psb_cbase_rwextd(nr,a,info,b,rowscale)
use psb_mat_mod
implicit none
integer, intent(in) :: nr
class(psb_c_base_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
class(psb_c_base_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
end subroutine psb_cbase_rwextd
subroutine psb_zrwextd(nr,a,info,b,rowscale)
use psb_mat_mod
implicit none
integer, intent(in) :: nr
type(psb_z_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
type(psb_z_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
end subroutine psb_zrwextd
subroutine psb_zbase_rwextd(nr,a,info,b,rowscale)
use psb_mat_mod
implicit none
integer, intent(in) :: nr
class(psb_z_base_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
class(psb_z_base_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
end subroutine psb_zbase_rwextd
end interface
base/serial/Makefile
+4 -2
View File
@@ -1,8 +1,10 @@
include ../../Make.inc
FOBJS = psb_lsame.o psb_dsymbmm.o psb_dnumbmm.o \
psb_drwextd.o
FOBJS = psb_lsame.o \
psb_ssymbmm.o psb_dsymbmm.o psb_csymbmm.o psb_zsymbmm.o \
psb_snumbmm.o psb_dnumbmm.o psb_cnumbmm.o psb_znumbmm.o \
psb_srwextd.o psb_drwextd.o psb_crwextd.o psb_zrwextd.o
# psb_regen_mod.o psb_lsame.o psb_zspgetrow.o\
# psb_zcsmm.o psb_zcsmv.o psb_zspgtdiag.o psb_zspgtblk.o\
# psb_zcsnmi.o psb_zcsrws.o psb_zcssm.o psb_zcssv.o psb_zspcnv.o\
base/serial/psb_cnumbmm.f90
+151 -42
View File
@@ -29,7 +29,7 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
! File: psb_snumbmm.f90
! File: psb_cnumbmm.f90
! Subroutine:
! Arguments:
!
@@ -41,41 +41,146 @@
!
subroutine psb_cnumbmm(a,b,c)
use psb_spmat_type
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_cnumbmm
implicit none
implicit none
type(psb_cspmat_type) :: a,b,c
complex(psb_spk_), allocatable :: temp(:)
integer :: info
logical :: csra, csrb
allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info)
if (info /= 0) then
return
type(psb_c_sparse_mat), intent(in) :: a,b
type(psb_c_sparse_mat), intent(inout) :: c
integer :: info
integer :: err_act
character(len=*), parameter :: name='psb_numbmm'
call psb_erractionsave(err_act)
info = 0
if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
call psb_realloc(size(c%ia1),c%aspk,info)
select type(aa=>c%a)
type is (psb_c_csr_sparse_mat)
call psb_numbmm(a%a,b%a,aa)
class default
info = 1121
call psb_errpush(info,name)
goto 9999
end select
call c%set_asb()
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_cnumbmm
subroutine psb_cbase_numbmm(a,b,c)
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_cbase_numbmm
implicit none
class(psb_c_base_sparse_mat), intent(in) :: a,b
type(psb_c_csr_sparse_mat), intent(inout) :: c
integer, allocatable :: itemp(:)
integer :: nze, ma,na,mb,nb
character(len=20) :: name
complex(psb_spk_), allocatable :: temp(:)
integer :: info
integer :: err_act
name='psb_numbmm'
call psb_erractionsave(err_act)
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
if ( mb /= na ) then
write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
endif
allocate(temp(max(ma,na,mb,nb)),stat=info)
if (info /= 0) then
info = 4000
call psb_Errpush(info,name)
goto 9999
endif
!
! Note: we still have to test about possible performance hits.
!
!
csra = (psb_toupper(a%fida(1:3))=='CSR')
csrb = (psb_toupper(b%fida(1:3))=='CSR')
if (csra.and.csrb) then
call cnumbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,&
& b%ia2,b%ia1,0,b%aspk,&
& c%ia2,c%ia1,0,c%aspk,temp)
else
call inner_numbmm(a,b,c,temp,info)
call psb_ensure_size(size(c%ja),c%val,info)
select type(a)
type is (psb_c_csr_sparse_mat)
select type(b)
type is (psb_c_csr_sparse_mat)
call csr_numbmm(a,b,c,temp,info)
class default
call gen_numbmm(a,b,c,temp,info)
end select
class default
call gen_numbmm(a,b,c,temp,info)
end select
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call c%set_asb()
deallocate(temp)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
contains
subroutine inner_numbmm(a,b,c,temp,info)
type(psb_cspmat_type) :: a,b,c
subroutine csr_numbmm(a,b,c,temp,info)
type(psb_c_csr_sparse_mat), intent(in) :: a,b
type(psb_c_csr_sparse_mat), intent(inout) :: c
complex(psb_spk_) :: temp(:)
integer, intent(out) :: info
integer :: nze, ma,na,mb,nb
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
call cnumbmm(ma,na,nb,a%irp,a%ja,0,a%val,&
& b%irp,b%ja,0,b%val,&
& c%irp,c%ja,0,c%val,temp)
end subroutine csr_numbmm
subroutine gen_numbmm(a,b,c,temp,info)
class(psb_c_base_sparse_mat), intent(in) :: a,b
type(psb_c_csr_sparse_mat), intent(inout) :: c
integer :: info
complex(psb_spk_) :: temp(:)
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
@@ -83,56 +188,60 @@ contains
integer :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln
complex(psb_spk_) :: ajj
n = a%m
m = a%k
l = b%k
n = a%get_nrows()
m = a%get_ncols()
l = b%get_ncols()
maxlmn = max(l,m,n)
allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
& aval(maxlmn),bval(maxlmn), stat=info)
if (info /= 0) then
if (info /= 0) then
info = 4000
return
endif
do i = 1,maxlmn
temp(i) = czero
temp(i) = dzero
end do
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
do i = 1,n
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
call a%csget(i,i,nazr,iarw,iacl,aval,info)
do jj=1, nazr
j=iacl(jj)
ajj = aval(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' NUMBMM: Problem with A ',i,jj,j,m
info = 1
return
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
info = 2
return
else
temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k)
endif
enddo
end do
do j = c%ia2(i),c%ia2(i+1)-1
if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn)) then
write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn
do j = c%irp(i),c%irp(i+1)-1
if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then
write(0,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn
info = 3
return
else
c%aspk(j) = temp(c%ia1(j))
temp(c%ia1(j)) = czero
c%val(j) = temp(c%ja(j))
temp(c%ja(j)) = dzero
endif
end do
end do
end subroutine gen_numbmm
end subroutine psb_cbase_numbmm
end subroutine inner_numbmm
end subroutine psb_cnumbmm
base/serial/psb_crwextd.f90
+183 -100
View File
@@ -39,126 +39,54 @@
!
!
subroutine psb_crwextd(nr,a,info,b,rowscale)
use psb_spmat_type
use psb_error_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_crwextd
implicit none
! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
integer, intent(in) :: nr
type(psb_cspmat_type), intent(inout) :: a
integer,intent(out) :: info
type(psb_cspmat_type), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer, intent(in) :: nr
type(psb_c_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
type(psb_c_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer :: i,j,ja,jb,err_act,nza,nzb
character(len=20) :: name, ch_err
type(psb_c_coo_sparse_mat) :: actmp
logical rowscale_
name='psb_crwextd'
info = 0
call psb_erractionsave(err_act)
if (present(rowscale)) then
rowscale_ = rowscale
else
rowscale_ = .true.
end if
if (nr > a%m) then
if (psb_toupper(a%fida) == 'CSR') then
call psb_realloc(nr+1,a%ia2,info)
if (nr > a%get_nrows()) then
select type(aa=> a%a)
type is (psb_c_csr_sparse_mat)
if (present(b)) then
nzb = psb_sp_get_nnzeros(b)
call psb_realloc(size(a%ia1)+nzb,a%ia1,info)
call psb_realloc(size(a%aspk)+nzb,a%aspk,info)
if (psb_toupper(b%fida)=='CSR') then
do i=1, min(nr-a%m,b%m)
a%ia2(a%m+i+1) = a%ia2(a%m+i) + b%ia2(i+1) - b%ia2(i)
ja = a%ia2(a%m+i)
jb = b%ia2(i)
do
if (jb >= b%ia2(i+1)) exit
a%aspk(ja) = b%aspk(jb)
a%ia1(ja) = b%ia1(jb)
ja = ja + 1
jb = jb + 1
end do
end do
do j=i,nr-a%m
a%ia2(a%m+i+1) = a%ia2(a%m+i)
end do
else
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
endif
call psb_rwextd(nr,aa,info,b%a,rowscale)
else
do i=a%m+2,nr+1
a%ia2(i) = a%ia2(i-1)
end do
call psb_rwextd(nr,aa,info,rowscale=rowscale)
end if
a%m = nr
a%k = max(a%k,b%k)
else if (psb_toupper(a%fida) == 'COO') then
type is (psb_c_coo_sparse_mat)
if (present(b)) then
nza = psb_sp_get_nnzeros(a)
nzb = psb_sp_get_nnzeros(b)
call psb_sp_reall(a,nza+nzb,info)
if (psb_toupper(b%fida)=='COO') then
if (rowscale_) then
do j=1,nzb
if ((a%m + b%ia1(j)) <= nr) then
nza = nza + 1
a%ia1(nza) = a%m + b%ia1(j)
a%ia2(nza) = b%ia2(j)
a%aspk(nza) = b%aspk(j)
end if
enddo
else
do j=1,nzb
if ((b%ia1(j)) <= nr) then
nza = nza + 1
a%ia1(nza) = b%ia1(j)
a%ia2(nza) = b%ia2(j)
a%aspk(nza) = b%aspk(j)
endif
enddo
endif
a%infoa(psb_nnz_) = nza
else if(psb_toupper(b%fida)=='CSR') then
do i=1, min(nr-a%m,b%m)
do
jb = b%ia2(i)
if (jb >= b%ia2(i+1)) exit
nza = nza + 1
a%aspk(nza) = b%aspk(jb)
a%ia1(nza) = a%m + i
a%ia2(nza) = b%ia1(jb)
jb = jb + 1
end do
end do
a%infoa(psb_nnz_) = nza
call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale)
else
call psb_rwextd(nr,aa,info,rowscale=rowscale)
end if
class default
call aa%mv_to_coo(actmp,info)
if (info == 0) then
if (present(b)) then
call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale)
else
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
endif
endif
a%m = nr
a%k = max(a%k,b%k)
else if (a%fida == 'JAD') then
info=135
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
else
info=136
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_rwextd(nr,actmp,info,rowscale=rowscale)
end if
end if
if (info == 0) call aa%mv_from_coo(actmp,info)
end select
end if
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
@@ -172,3 +100,158 @@ subroutine psb_crwextd(nr,a,info,b,rowscale)
return
end subroutine psb_crwextd
subroutine psb_cbase_rwextd(nr,a,info,b,rowscale)
use psb_error_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_cbase_rwextd
implicit none
! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
integer, intent(in) :: nr
class(psb_c_base_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
class(psb_c_base_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer :: i,j,ja,jb,err_act,nza,nzb, ma, mb, na, nb
character(len=20) :: name, ch_err
logical rowscale_
name='psb_cbase_rwextd'
info = 0
call psb_erractionsave(err_act)
if (present(rowscale)) then
rowscale_ = rowscale
else
rowscale_ = .true.
end if
ma = a%get_nrows()
na = a%get_ncols()
select type(a)
type is (psb_c_csr_sparse_mat)
call psb_ensure_size(nr+1,a%irp,info)
if (present(b)) then
mb = b%get_nrows()
nb = b%get_ncols()
nzb = b%get_nzeros()
select type (b)
type is (psb_c_csr_sparse_mat)
call psb_ensure_size(size(a%ja)+nzb,a%ja,info)
call psb_ensure_size(size(a%val)+nzb,a%val,info)
do i=1, min(nr-ma,mb)
a%irp(ma+i+1) = a%irp(ma+i) + b%irp(i+1) - b%irp(i)
ja = a%irp(ma+i)
jb = b%irp(i)
do
if (jb >= b%irp(i+1)) exit
a%val(ja) = b%val(jb)
a%ja(ja) = b%ja(jb)
ja = ja + 1
jb = jb + 1
end do
end do
do j=i,nr-ma
a%irp(ma+i+1) = a%irp(ma+i)
end do
class default
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
end select
call a%set_ncols(max(na,nb))
else
do i=ma+2,nr+1
a%irp(i) = a%irp(i-1)
end do
end if
call a%set_nrows(nr)
type is (psb_c_coo_sparse_mat)
nza = a%get_nzeros()
if (present(b)) then
mb = b%get_nrows()
nb = b%get_ncols()
nzb = b%get_nzeros()
call a%reallocate(nza+nzb)
select type(b)
type is (psb_c_coo_sparse_mat)
if (rowscale_) then
do j=1,nzb
if ((ma + b%ia(j)) <= nr) then
nza = nza + 1
a%ia(nza) = ma + b%ia(j)
a%ja(nza) = b%ja(j)
a%val(nza) = b%val(j)
end if
enddo
else
do j=1,nzb
if ((ma + b%ia(j)) <= nr) then
nza = nza + 1
a%ia(nza) = b%ia(j)
a%ja(nza) = b%ja(j)
a%val(nza) = b%val(j)
end if
enddo
endif
call a%set_nzeros(nza)
type is (psb_c_csr_sparse_mat)
do i=1, min(nr-ma,mb)
do
jb = b%irp(i)
if (jb >= b%irp(i+1)) exit
nza = nza + 1
a%val(nza) = b%val(jb)
a%ia(nza) = ma + i
a%ja(nza) = b%ja(jb)
jb = jb + 1
end do
end do
call a%set_nzeros(nza)
class default
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
end select
call a%set_ncols(max(na,nb))
endif
call a%set_nrows(nr)
class default
info = 135
ch_err=a%get_fmt()
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_cbase_rwextd
base/serial/psb_csymbmm.f90
+177 -95
View File
@@ -40,62 +40,107 @@
!
subroutine psb_csymbmm(a,b,c,info)
use psb_spmat_type
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_csymbmm
implicit none
type(psb_cspmat_type) :: a,b,c
type(psb_c_sparse_mat), intent(in) :: a,b
type(psb_c_sparse_mat), intent(out) :: c
integer, intent(out) :: info
type(psb_c_csr_sparse_mat), allocatable :: ccsr
integer :: err_act
character(len=*), parameter :: name='psb_symbmm'
call psb_erractionsave(err_act)
info = 0
if ((a%is_null()) .or.(b%is_null())) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
allocate(ccsr, stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call psb_symbmm(a%a,b%a,ccsr,info)
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call move_alloc(ccsr,c%a)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_csymbmm
subroutine psb_cbase_symbmm(a,b,c,info)
use psb_mat_mod
use psb_serial_mod, psb_protect_name => psb_cbase_symbmm
implicit none
class(psb_c_base_sparse_mat), intent(in) :: a,b
type(psb_c_csr_sparse_mat), intent(out) :: c
integer, intent(out) :: info
integer, allocatable :: itemp(:)
integer :: nze,info
interface
subroutine symbmm (n, m, l, ia, ja, diaga, &
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
character(len=20) :: name
integer :: err_act
logical :: csra, csrb
integer :: nze, ma,na,mb,nb
character(len=20) :: name
integer :: err_act
name='psb_symbmm'
call psb_erractionsave(err_act)
info = 0
csra = (psb_toupper(a%fida(1:3))=='CSR')
csrb = (psb_toupper(b%fida(1:3))=='CSR')
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
if (b%m /= a%k) then
write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k
if ( mb /= na ) then
write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
endif
allocate(itemp(max(a%m,a%k,b%m,b%k)),stat=info)
allocate(itemp(max(ma,na,mb,nb)),stat=info)
if (info /= 0) then
return
info = 4000
call psb_Errpush(info,name)
goto 9999
endif
nze = max(a%m+1,2*a%m)
call psb_sp_reall(c,nze,info)
!
! Note: we need to test whether there is a performance impact
! in not using the original Douglas & Bank code.
!
if (csra.and.csrb) then
call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,&
& b%ia2,b%ia1,0,&
& c%ia2,c%ia1,0,itemp)
else
call inner_symbmm(a,b,c,itemp,info)
endif
call psb_realloc(size(c%ia1),c%aspk,info)
select type(a)
type is (psb_c_csr_sparse_mat)
select type(b)
type is (psb_c_csr_sparse_mat)
call csr_symbmm(a,b,c,itemp,info)
class default
call gen_symbmm(a,b,c,itemp,info)
end select
class default
call gen_symbmm(a,b,c,itemp,info)
end select
c%pl(1) = 0
c%pr(1) = 0
c%m=a%m
c%k=b%k
c%fida='CSR'
c%descra='GUN'
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call psb_realloc(size(c%ja),c%val,info)
deallocate(itemp)
call psb_erractionrestore(err_act)
return
@@ -108,82 +153,119 @@ subroutine psb_csymbmm(a,b,c,info)
return
contains
subroutine inner_symbmm(a,b,c,index,info)
type(psb_cspmat_type) :: a,b,c
subroutine csr_symbmm(a,b,c,itemp,info)
type(psb_c_csr_sparse_mat), intent(in) :: a,b
type(psb_c_csr_sparse_mat), intent(out) :: c
integer :: itemp(:)
integer, intent(out) :: info
interface
subroutine symbmm (n, m, l, ia, ja, diaga, &
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
integer :: nze, ma,na,mb,nb
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
nze = max(ma+1,2*ma)
call c%allocate(ma,nb,nze)
call symbmm(ma,na,nb,a%irp,a%ja,0,&
& b%irp,b%ja,0,&
& c%irp,c%ja,0,itemp)
end subroutine csr_symbmm
subroutine gen_symbmm(a,b,c,index,info)
class(psb_c_base_sparse_mat), intent(in) :: a,b
type(psb_c_csr_sparse_mat), intent(out) :: c
integer :: index(:),info
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
complex(psb_spk_), allocatable :: aval(:),bval(:)
integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
integer :: nze, ma,na,mb,nb
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
n = a%m
m = a%k
l = b%k
nze = max(ma+1,2*ma)
call c%allocate(ma,nb,nze)
n = ma
m = na
l = nb
maxlmn = max(l,m,n)
allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
& aval(maxlmn),bval(maxlmn), stat=info)
& stat=info)
if (info /= 0) then
info = 4000
return
endif
if (size(c%ia2) < n+1) then
call psb_realloc(n+1,c%ia2,info)
endif
do i=1,maxlmn
index(i)=0
end do
c%ia2(1)=1
minlm = min(l,m)
minmn = min(m,n)
c%irp(1)=1
minlm = min(l,m)
minmn = min(m,n)
main: do i=1,n
istart=-1
length=0
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
do jj=1, nazr
j=iacl(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
else
if(index(ibcl(k)) == 0) then
index(ibcl(k))=istart
istart=ibcl(k)
length=length+1
endif
endif
end do
end do
main: do i=1,n
istart=-1
length=0
call a%csget(i,i,nazr,iarw,iacl,info)
do jj=1, nazr
c%ia2(i+1)=c%ia2(i)+length
if (c%ia2(i+1) > size(c%ia1)) then
if (n > (2*i)) then
nze = max(c%ia2(i+1), c%ia2(i)*((n+i-1)/i))
j=iacl(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
info = 1
return
endif
call b%csget(j,j,nbzr,ibrw,ibcl,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
info=2
return
else
nze = max(c%ia2(i+1), nint((dble(c%ia2(i))*(dble(n)/i))) )
endif
call psb_realloc(nze,c%ia1,info)
end if
do j= c%ia2(i),c%ia2(i+1)-1
c%ia1(j)=istart
istart=index(istart)
index(c%ia1(j))=0
if(index(ibcl(k)) == 0) then
index(ibcl(k))=istart
istart=ibcl(k)
length=length+1
endif
endif
end do
call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1))
index(i) = 0
end do main
end do
end subroutine inner_symbmm
c%irp(i+1)=c%irp(i)+length
end subroutine psb_csymbmm
if (c%irp(i+1) > size(c%ja)) then
if (n > (2*i)) then
nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i))
else
nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) )
endif
call psb_realloc(nze,c%ja,info)
end if
do j= c%irp(i),c%irp(i+1)-1
c%ja(j)=istart
istart=index(istart)
index(c%ja(j))=0
end do
call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1))
index(i) = 0
end do main
end subroutine gen_symbmm
end subroutine psb_cbase_symbmm
base/serial/psb_dsymbmm.f90
+3 -4
View File
@@ -187,7 +187,6 @@ contains
type(psb_d_csr_sparse_mat), intent(out) :: c
integer :: index(:),info
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
real(psb_dpk_), allocatable :: aval(:),bval(:)
integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
integer :: nze, ma,na,mb,nb
@@ -205,7 +204,7 @@ contains
maxlmn = max(l,m,n)
allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
& aval(maxlmn),bval(maxlmn), stat=info)
& stat=info)
if (info /= 0) then
info = 4000
return
@@ -222,7 +221,7 @@ contains
main: do i=1,n
istart=-1
length=0
call a%csget(i,i,nazr,iarw,iacl,aval,info)
call a%csget(i,i,nazr,iarw,iacl,info)
do jj=1, nazr
j=iacl(jj)
@@ -232,7 +231,7 @@ contains
info = 1
return
endif
call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
call b%csget(j,j,nbzr,ibrw,ibcl,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
base/serial/psb_snumbmm.f90
+144 -43
View File
@@ -41,72 +41,173 @@
!
subroutine psb_snumbmm(a,b,c)
use psb_spmat_type
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_snumbmm
implicit none
implicit none
type(psb_sspmat_type) :: a,b,c
real(psb_spk_), allocatable :: temp(:)
integer :: info
logical :: csra, csrb
allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info)
if (info /= 0) then
return
type(psb_s_sparse_mat), intent(in) :: a,b
type(psb_s_sparse_mat), intent(inout) :: c
integer :: info
integer :: err_act
character(len=*), parameter :: name='psb_numbmm'
call psb_erractionsave(err_act)
info = 0
if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
call psb_realloc(size(c%ia1),c%aspk,info)
select type(aa=>c%a)
type is (psb_s_csr_sparse_mat)
call psb_numbmm(a%a,b%a,aa)
class default
info = 1121
call psb_errpush(info,name)
goto 9999
end select
call c%set_asb()
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_snumbmm
subroutine psb_sbase_numbmm(a,b,c)
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_sbase_numbmm
implicit none
class(psb_s_base_sparse_mat), intent(in) :: a,b
type(psb_s_csr_sparse_mat), intent(inout) :: c
integer, allocatable :: itemp(:)
integer :: nze, ma,na,mb,nb
character(len=20) :: name
real(psb_spk_), allocatable :: temp(:)
integer :: info
integer :: err_act
name='psb_numbmm'
call psb_erractionsave(err_act)
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
if ( mb /= na ) then
write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
endif
allocate(temp(max(ma,na,mb,nb)),stat=info)
if (info /= 0) then
info = 4000
call psb_Errpush(info,name)
goto 9999
endif
!
! Note: we still have to test about possible performance hits.
!
!
csra = (psb_toupper(a%fida(1:3))=='CSR')
csrb = (psb_toupper(b%fida(1:3))=='CSR')
if (csra.and.csrb) then
call snumbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,&
& b%ia2,b%ia1,0,b%aspk,&
& c%ia2,c%ia1,0,c%aspk,temp)
else
call inner_numbmm(a,b,c,temp,info)
if (info /= 0) then
write(0,*) 'Error ',info,' from inner numbmm'
end if
call psb_ensure_size(size(c%ja),c%val,info)
select type(a)
type is (psb_s_csr_sparse_mat)
select type(b)
type is (psb_s_csr_sparse_mat)
call csr_numbmm(a,b,c,temp,info)
class default
call gen_numbmm(a,b,c,temp,info)
end select
class default
call gen_numbmm(a,b,c,temp,info)
end select
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call psb_sp_setifld(psb_spmat_asb_,psb_state_,c,info)
call c%set_asb()
deallocate(temp)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
contains
subroutine inner_numbmm(a,b,c,temp,info)
type(psb_sspmat_type) :: a,b,c
subroutine csr_numbmm(a,b,c,temp,info)
type(psb_s_csr_sparse_mat), intent(in) :: a,b
type(psb_s_csr_sparse_mat), intent(inout) :: c
real(psb_spk_) :: temp(:)
integer, intent(out) :: info
integer :: nze, ma,na,mb,nb
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
call snumbmm(ma,na,nb,a%irp,a%ja,0,a%val,&
& b%irp,b%ja,0,b%val,&
& c%irp,c%ja,0,c%val,temp)
end subroutine csr_numbmm
subroutine gen_numbmm(a,b,c,temp,info)
class(psb_s_base_sparse_mat), intent(in) :: a,b
type(psb_s_csr_sparse_mat), intent(inout) :: c
integer :: info
real(psb_spk_) :: temp(:)
real(psb_spk_) :: temp(:)
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
real(psb_spk_), allocatable :: aval(:),bval(:)
integer :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln
real(psb_spk_) :: ajj
n = a%m
m = a%k
l = b%k
n = a%get_nrows()
m = a%get_ncols()
l = b%get_ncols()
maxlmn = max(l,m,n)
allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
& aval(maxlmn),bval(maxlmn), stat=info)
if (info /= 0) then
if (info /= 0) then
info = 4000
return
endif
do i = 1,maxlmn
temp(i) = szero
temp(i) = dzero
end do
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
do i = 1,n
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
call a%csget(i,i,nazr,iarw,iacl,aval,info)
do jj=1, nazr
j=iacl(jj)
ajj = aval(jj)
@@ -116,7 +217,7 @@ contains
return
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
@@ -127,19 +228,19 @@ contains
endif
enddo
end do
do j = c%ia2(i),c%ia2(i+1)-1
if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn)) then
write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn
do j = c%irp(i),c%irp(i+1)-1
if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then
write(0,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn
info = 3
return
else
c%aspk(j) = temp(c%ia1(j))
temp(c%ia1(j)) = szero
c%val(j) = temp(c%ja(j))
temp(c%ja(j)) = dzero
endif
end do
end do
end subroutine gen_numbmm
end subroutine inner_numbmm
end subroutine psb_snumbmm
end subroutine psb_sbase_numbmm
base/serial/psb_srwextd.f90
+182 -100
View File
@@ -39,127 +39,54 @@
!
!
subroutine psb_srwextd(nr,a,info,b,rowscale)
use psb_spmat_type
use psb_error_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_srwextd
implicit none
! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
integer, intent(in) :: nr
type(psb_sspmat_type), intent(inout) :: a
integer,intent(out) :: info
type(psb_sspmat_type), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer, intent(in) :: nr
type(psb_s_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
type(psb_s_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer :: i,j,ja,jb,err_act,nza,nzb
character(len=20) :: name, ch_err
type(psb_s_coo_sparse_mat) :: actmp
logical rowscale_
name='psb_srwextd'
info = 0
call psb_erractionsave(err_act)
if (present(rowscale)) then
rowscale_ = rowscale
else
rowscale_ = .true.
end if
if (nr > a%m) then
if (psb_toupper(a%fida) == 'CSR') then
call psb_ensure_size(nr+1,a%ia2,info)
if (nr > a%get_nrows()) then
select type(aa=> a%a)
type is (psb_s_csr_sparse_mat)
if (present(b)) then
nzb = psb_sp_get_nnzeros(b)
call psb_ensure_size(size(a%ia1)+nzb,a%ia1,info)
call psb_ensure_size(size(a%aspk)+nzb,a%aspk,info)
if (psb_toupper(b%fida)=='CSR') then
do i=1, min(nr-a%m,b%m)
a%ia2(a%m+i+1) = a%ia2(a%m+i) + b%ia2(i+1) - b%ia2(i)
ja = a%ia2(a%m+i)
jb = b%ia2(i)
do
if (jb >= b%ia2(i+1)) exit
a%aspk(ja) = b%aspk(jb)
a%ia1(ja) = b%ia1(jb)
ja = ja + 1
jb = jb + 1
end do
end do
do j=i,nr-a%m
a%ia2(a%m+i+1) = a%ia2(a%m+i)
end do
else
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
endif
call psb_rwextd(nr,aa,info,b%a,rowscale)
else
do i=a%m+2,nr+1
a%ia2(i) = a%ia2(i-1)
end do
call psb_rwextd(nr,aa,info,rowscale=rowscale)
end if
a%m = nr
a%k = max(a%k,b%k)
else if (psb_toupper(a%fida) == 'COO') then
type is (psb_s_coo_sparse_mat)
if (present(b)) then
nza = psb_sp_get_nnzeros(a)
nzb = psb_sp_get_nnzeros(b)
call psb_sp_reall(a,nza+nzb,info)
if (psb_toupper(b%fida)=='COO') then
if (rowscale_) then
do j=1,nzb
if ((a%m + b%ia1(j)) <= nr) then
nza = nza + 1
a%ia1(nza) = a%m + b%ia1(j)
a%ia2(nza) = b%ia2(j)
a%aspk(nza) = b%aspk(j)
end if
enddo
else
do j=1,nzb
if ((b%ia1(j)) <= nr) then
nza = nza + 1
a%ia1(nza) = b%ia1(j)
a%ia2(nza) = b%ia2(j)
a%aspk(nza) = b%aspk(j)
endif
enddo
endif
a%infoa(psb_nnz_) = nza
else if(psb_toupper(b%fida)=='CSR') then
do i=1, min(nr-a%m,b%m)
do
jb = b%ia2(i)
if (jb >= b%ia2(i+1)) exit
nza = nza + 1
a%aspk(nza) = b%aspk(jb)
a%ia1(nza) = a%m + i
a%ia2(nza) = b%ia1(jb)
jb = jb + 1
end do
end do
a%infoa(psb_nnz_) = nza
call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale)
else
call psb_rwextd(nr,aa,info,rowscale=rowscale)
end if
class default
call aa%mv_to_coo(actmp,info)
if (info == 0) then
if (present(b)) then
call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale)
else
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
endif
endif
a%m = nr
a%k = max(a%k,b%k)
else if (a%fida == 'JAD') then
info=135
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
else
info=136
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_rwextd(nr,actmp,info,rowscale=rowscale)
end if
end if
if (info == 0) call aa%mv_from_coo(actmp,info)
end select
end if
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
@@ -173,3 +100,158 @@ subroutine psb_srwextd(nr,a,info,b,rowscale)
return
end subroutine psb_srwextd
subroutine psb_sbase_rwextd(nr,a,info,b,rowscale)
use psb_error_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_sbase_rwextd
implicit none
! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
integer, intent(in) :: nr
class(psb_s_base_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
class(psb_s_base_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer :: i,j,ja,jb,err_act,nza,nzb, ma, mb, na, nb
character(len=20) :: name, ch_err
logical rowscale_
name='psb_sbase_rwextd'
info = 0
call psb_erractionsave(err_act)
if (present(rowscale)) then
rowscale_ = rowscale
else
rowscale_ = .true.
end if
ma = a%get_nrows()
na = a%get_ncols()
select type(a)
type is (psb_s_csr_sparse_mat)
call psb_ensure_size(nr+1,a%irp,info)
if (present(b)) then
mb = b%get_nrows()
nb = b%get_ncols()
nzb = b%get_nzeros()
select type (b)
type is (psb_s_csr_sparse_mat)
call psb_ensure_size(size(a%ja)+nzb,a%ja,info)
call psb_ensure_size(size(a%val)+nzb,a%val,info)
do i=1, min(nr-ma,mb)
a%irp(ma+i+1) = a%irp(ma+i) + b%irp(i+1) - b%irp(i)
ja = a%irp(ma+i)
jb = b%irp(i)
do
if (jb >= b%irp(i+1)) exit
a%val(ja) = b%val(jb)
a%ja(ja) = b%ja(jb)
ja = ja + 1
jb = jb + 1
end do
end do
do j=i,nr-ma
a%irp(ma+i+1) = a%irp(ma+i)
end do
class default
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
end select
call a%set_ncols(max(na,nb))
else
do i=ma+2,nr+1
a%irp(i) = a%irp(i-1)
end do
end if
call a%set_nrows(nr)
type is (psb_s_coo_sparse_mat)
nza = a%get_nzeros()
if (present(b)) then
mb = b%get_nrows()
nb = b%get_ncols()
nzb = b%get_nzeros()
call a%reallocate(nza+nzb)
select type(b)
type is (psb_s_coo_sparse_mat)
if (rowscale_) then
do j=1,nzb
if ((ma + b%ia(j)) <= nr) then
nza = nza + 1
a%ia(nza) = ma + b%ia(j)
a%ja(nza) = b%ja(j)
a%val(nza) = b%val(j)
end if
enddo
else
do j=1,nzb
if ((ma + b%ia(j)) <= nr) then
nza = nza + 1
a%ia(nza) = b%ia(j)
a%ja(nza) = b%ja(j)
a%val(nza) = b%val(j)
end if
enddo
endif
call a%set_nzeros(nza)
type is (psb_s_csr_sparse_mat)
do i=1, min(nr-ma,mb)
do
jb = b%irp(i)
if (jb >= b%irp(i+1)) exit
nza = nza + 1
a%val(nza) = b%val(jb)
a%ia(nza) = ma + i
a%ja(nza) = b%ja(jb)
jb = jb + 1
end do
end do
call a%set_nzeros(nza)
class default
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
end select
call a%set_ncols(max(na,nb))
endif
call a%set_nrows(nr)
class default
info = 135
ch_err=a%get_fmt()
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_sbase_rwextd
base/serial/psb_ssymbmm.f90
+174 -97
View File
@@ -40,62 +40,107 @@
!
subroutine psb_ssymbmm(a,b,c,info)
use psb_spmat_type
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_ssymbmm
implicit none
type(psb_sspmat_type) :: a,b,c
integer, allocatable :: itemp(:)
integer :: nze,info
type(psb_s_sparse_mat), intent(in) :: a,b
type(psb_s_sparse_mat), intent(out) :: c
integer, intent(out) :: info
type(psb_s_csr_sparse_mat), allocatable :: ccsr
integer :: err_act
character(len=*), parameter :: name='psb_symbmm'
call psb_erractionsave(err_act)
info = 0
interface
subroutine symbmm (n, m, l, ia, ja, diaga, &
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
if ((a%is_null()) .or.(b%is_null())) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
allocate(ccsr, stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call psb_symbmm(a%a,b%a,ccsr,info)
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call move_alloc(ccsr,c%a)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_ssymbmm
subroutine psb_sbase_symbmm(a,b,c,info)
use psb_mat_mod
use psb_serial_mod, psb_protect_name => psb_sbase_symbmm
implicit none
class(psb_s_base_sparse_mat), intent(in) :: a,b
type(psb_s_csr_sparse_mat), intent(out) :: c
integer, intent(out) :: info
integer, allocatable :: itemp(:)
integer :: nze, ma,na,mb,nb
character(len=20) :: name
integer :: err_act
logical :: csra, csrb
name='psb_symbmm'
call psb_erractionsave(err_act)
info = 0
csra = (psb_toupper(a%fida(1:3))=='CSR')
csrb = (psb_toupper(b%fida(1:3))=='CSR')
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
if (b%m /= a%k) then
write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k
if ( mb /= na ) then
write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
endif
allocate(itemp(max(a%m,a%k,b%m,b%k)),stat=info)
allocate(itemp(max(ma,na,mb,nb)),stat=info)
if (info /= 0) then
return
info = 4000
call psb_Errpush(info,name)
goto 9999
endif
nze = max(a%m+1,2*a%m)
call psb_sp_reall(c,nze,info)
!
! Note: we need to test whether there is a performance impact
! in not using the original Douglas & Bank code.
!
if (csra.and.csrb) then
call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,&
& b%ia2,b%ia1,0,&
& c%ia2,c%ia1,0,itemp)
else
call inner_symbmm(a,b,c,itemp,info)
endif
call psb_realloc(size(c%ia1),c%aspk,info)
select type(a)
type is (psb_s_csr_sparse_mat)
select type(b)
type is (psb_s_csr_sparse_mat)
call csr_symbmm(a,b,c,itemp,info)
class default
call gen_symbmm(a,b,c,itemp,info)
end select
class default
call gen_symbmm(a,b,c,itemp,info)
end select
c%pl(1) = 0
c%pr(1) = 0
c%m=a%m
c%k=b%k
c%fida='CSR'
c%descra='GUN'
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call psb_realloc(size(c%ja),c%val,info)
deallocate(itemp)
call psb_erractionrestore(err_act)
return
@@ -109,86 +154,118 @@ subroutine psb_ssymbmm(a,b,c,info)
contains
subroutine inner_symbmm(a,b,c,index,info)
type(psb_sspmat_type) :: a,b,c
subroutine csr_symbmm(a,b,c,itemp,info)
type(psb_s_csr_sparse_mat), intent(in) :: a,b
type(psb_s_csr_sparse_mat), intent(out) :: c
integer :: itemp(:)
integer, intent(out) :: info
interface
subroutine symbmm (n, m, l, ia, ja, diaga, &
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
integer :: nze, ma,na,mb,nb
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
nze = max(ma+1,2*ma)
call c%allocate(ma,nb,nze)
call symbmm(ma,na,nb,a%irp,a%ja,0,&
& b%irp,b%ja,0,&
& c%irp,c%ja,0,itemp)
end subroutine csr_symbmm
subroutine gen_symbmm(a,b,c,index,info)
class(psb_s_base_sparse_mat), intent(in) :: a,b
type(psb_s_csr_sparse_mat), intent(out) :: c
integer :: index(:),info
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
real(psb_spk_), allocatable :: aval(:),bval(:)
integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
integer :: nze, ma,na,mb,nb
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
n = a%m
m = a%k
l = b%k
nze = max(ma+1,2*ma)
call c%allocate(ma,nb,nze)
n = ma
m = na
l = nb
maxlmn = max(l,m,n)
allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
& aval(maxlmn),bval(maxlmn), stat=info)
& stat=info)
if (info /= 0) then
info = 4000
return
endif
if (size(c%ia2) < n+1) then
call psb_realloc(n+1,c%ia2,info)
endif
do i=1,maxlmn
index(i)=0
end do
c%ia2(1)=1
minlm = min(l,m)
minmn = min(m,n)
c%irp(1)=1
minlm = min(l,m)
minmn = min(m,n)
main: do i=1,n
istart=-1
length=0
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
do jj=1, nazr
j=iacl(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
info = 1
main: do i=1,n
istart=-1
length=0
call a%csget(i,i,nazr,iarw,iacl,info)
do jj=1, nazr
j=iacl(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
info = 1
return
endif
call b%csget(j,j,nbzr,ibrw,ibcl,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
info=2
return
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
info=2
return
else
if(index(ibcl(k)) == 0) then
index(ibcl(k))=istart
istart=ibcl(k)
length=length+1
endif
endif
end do
end do
c%ia2(i+1)=c%ia2(i)+length
if (c%ia2(i+1) > size(c%ia1)) then
if (n > (2*i)) then
nze = max(c%ia2(i+1), c%ia2(i)*((n+i-1)/i))
else
nze = max(c%ia2(i+1), nint((dble(c%ia2(i))*(dble(n)/i))) )
endif
call psb_realloc(nze,c%ia1,info)
end if
do j= c%ia2(i),c%ia2(i+1)-1
c%ia1(j)=istart
istart=index(istart)
index(c%ia1(j))=0
if(index(ibcl(k)) == 0) then
index(ibcl(k))=istart
istart=ibcl(k)
length=length+1
endif
endif
end do
call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1))
index(i) = 0
end do main
end do
end subroutine inner_symbmm
c%irp(i+1)=c%irp(i)+length
end subroutine psb_ssymbmm
if (c%irp(i+1) > size(c%ja)) then
if (n > (2*i)) then
nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i))
else
nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) )
endif
call psb_realloc(nze,c%ja,info)
end if
do j= c%irp(i),c%irp(i+1)-1
c%ja(j)=istart
istart=index(istart)
index(c%ja(j))=0
end do
call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1))
index(i) = 0
end do main
end subroutine gen_symbmm
end subroutine psb_sbase_symbmm
base/serial/psb_znumbmm.f90
+150 -41
View File
@@ -29,7 +29,7 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
! File: psb_dnumbmm.f90
! File: psb_znumbmm.f90
! Subroutine:
! Arguments:
!
@@ -41,41 +41,146 @@
!
subroutine psb_znumbmm(a,b,c)
use psb_spmat_type
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_znumbmm
implicit none
implicit none
type(psb_zspmat_type) :: a,b,c
complex(psb_dpk_), allocatable :: temp(:)
integer :: info
logical :: csra, csrb
allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info)
if (info /= 0) then
return
type(psb_z_sparse_mat), intent(in) :: a,b
type(psb_z_sparse_mat), intent(inout) :: c
integer :: info
integer :: err_act
character(len=*), parameter :: name='psb_numbmm'
call psb_erractionsave(err_act)
info = 0
if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
call psb_realloc(size(c%ia1),c%aspk,info)
select type(aa=>c%a)
type is (psb_z_csr_sparse_mat)
call psb_numbmm(a%a,b%a,aa)
class default
info = 1121
call psb_errpush(info,name)
goto 9999
end select
call c%set_asb()
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_znumbmm
subroutine psb_zbase_numbmm(a,b,c)
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_zbase_numbmm
implicit none
class(psb_z_base_sparse_mat), intent(in) :: a,b
type(psb_z_csr_sparse_mat), intent(inout) :: c
integer, allocatable :: itemp(:)
integer :: nze, ma,na,mb,nb
character(len=20) :: name
complex(psb_dpk_), allocatable :: temp(:)
integer :: info
integer :: err_act
name='psb_numbmm'
call psb_erractionsave(err_act)
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
if ( mb /= na ) then
write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
endif
allocate(temp(max(ma,na,mb,nb)),stat=info)
if (info /= 0) then
info = 4000
call psb_Errpush(info,name)
goto 9999
endif
!
! Note: we still have to test about possible performance hits.
!
!
csra = (psb_toupper(a%fida(1:3))=='CSR')
csrb = (psb_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,&
& b%ia2,b%ia1,0,b%aspk,&
& c%ia2,c%ia1,0,c%aspk,temp)
else
call inner_numbmm(a,b,c,temp,info)
call psb_ensure_size(size(c%ja),c%val,info)
select type(a)
type is (psb_z_csr_sparse_mat)
select type(b)
type is (psb_z_csr_sparse_mat)
call csr_numbmm(a,b,c,temp,info)
class default
call gen_numbmm(a,b,c,temp,info)
end select
class default
call gen_numbmm(a,b,c,temp,info)
end select
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call c%set_asb()
deallocate(temp)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
contains
subroutine inner_numbmm(a,b,c,temp,info)
type(psb_zspmat_type) :: a,b,c
subroutine csr_numbmm(a,b,c,temp,info)
type(psb_z_csr_sparse_mat), intent(in) :: a,b
type(psb_z_csr_sparse_mat), intent(inout) :: c
complex(psb_dpk_) :: temp(:)
integer, intent(out) :: info
integer :: nze, ma,na,mb,nb
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
call cnumbmm(ma,na,nb,a%irp,a%ja,0,a%val,&
& b%irp,b%ja,0,b%val,&
& c%irp,c%ja,0,c%val,temp)
end subroutine csr_numbmm
subroutine gen_numbmm(a,b,c,temp,info)
class(psb_z_base_sparse_mat), intent(in) :: a,b
type(psb_z_csr_sparse_mat), intent(inout) :: c
integer :: info
complex(psb_dpk_) :: temp(:)
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
@@ -83,14 +188,14 @@ contains
integer :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln
complex(psb_dpk_) :: ajj
n = a%m
m = a%k
l = b%k
n = a%get_nrows()
m = a%get_ncols()
l = b%get_ncols()
maxlmn = max(l,m,n)
allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
& aval(maxlmn),bval(maxlmn), stat=info)
if (info /= 0) then
if (info /= 0) then
info = 4000
return
endif
@@ -102,37 +207,41 @@ contains
minmn = min(m,n)
do i = 1,n
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
call a%csget(i,i,nazr,iarw,iacl,aval,info)
do jj=1, nazr
j=iacl(jj)
ajj = aval(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' NUMBMM: Problem with A ',i,jj,j,m
info = 1
return
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
info = 2
return
else
temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k)
endif
enddo
end do
do j = c%ia2(i),c%ia2(i+1)-1
if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn)) then
write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn
do j = c%irp(i),c%irp(i+1)-1
if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then
write(0,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn
info = 3
return
else
c%aspk(j) = temp(c%ia1(j))
temp(c%ia1(j)) = dzero
c%val(j) = temp(c%ja(j))
temp(c%ja(j)) = dzero
endif
end do
end do
end subroutine gen_numbmm
end subroutine psb_zbase_numbmm
end subroutine inner_numbmm
end subroutine psb_znumbmm
base/serial/psb_zrwextd.f90
+183 -100
View File
@@ -39,126 +39,54 @@
!
!
subroutine psb_zrwextd(nr,a,info,b,rowscale)
use psb_spmat_type
use psb_error_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_zrwextd
implicit none
! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
integer, intent(in) :: nr
type(psb_zspmat_type), intent(inout) :: a
integer,intent(out) :: info
type(psb_zspmat_type), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer, intent(in) :: nr
type(psb_z_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
type(psb_z_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer :: i,j,ja,jb,err_act,nza,nzb
character(len=20) :: name, ch_err
type(psb_z_coo_sparse_mat) :: actmp
logical rowscale_
name='psb_zrwextd'
info = 0
call psb_erractionsave(err_act)
if (present(rowscale)) then
rowscale_ = rowscale
else
rowscale_ = .true.
end if
if (nr > a%m) then
if (psb_toupper(a%fida) == 'CSR') then
call psb_realloc(nr+1,a%ia2,info)
if (nr > a%get_nrows()) then
select type(aa=> a%a)
type is (psb_z_csr_sparse_mat)
if (present(b)) then
nzb = psb_sp_get_nnzeros(b)
call psb_realloc(size(a%ia1)+nzb,a%ia1,info)
call psb_realloc(size(a%aspk)+nzb,a%aspk,info)
if (psb_toupper(b%fida)=='CSR') then
do i=1, min(nr-a%m,b%m)
a%ia2(a%m+i+1) = a%ia2(a%m+i) + b%ia2(i+1) - b%ia2(i)
ja = a%ia2(a%m+i)
jb = b%ia2(i)
do
if (jb >= b%ia2(i+1)) exit
a%aspk(ja) = b%aspk(jb)
a%ia1(ja) = b%ia1(jb)
ja = ja + 1
jb = jb + 1
end do
end do
do j=i,nr-a%m
a%ia2(a%m+i+1) = a%ia2(a%m+i)
end do
else
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
endif
call psb_rwextd(nr,aa,info,b%a,rowscale)
else
do i=a%m+2,nr+1
a%ia2(i) = a%ia2(i-1)
end do
call psb_rwextd(nr,aa,info,rowscale=rowscale)
end if
a%m = nr
a%k = max(a%k,b%k)
else if (psb_toupper(a%fida) == 'COO') then
type is (psb_z_coo_sparse_mat)
if (present(b)) then
nza = psb_sp_get_nnzeros(a)
nzb = psb_sp_get_nnzeros(b)
call psb_sp_reall(a,nza+nzb,info)
if (psb_toupper(b%fida)=='COO') then
if (rowscale_) then
do j=1,nzb
if ((a%m + b%ia1(j)) <= nr) then
nza = nza + 1
a%ia1(nza) = a%m + b%ia1(j)
a%ia2(nza) = b%ia2(j)
a%aspk(nza) = b%aspk(j)
end if
enddo
else
do j=1,nzb
if ((b%ia1(j)) <= nr) then
nza = nza + 1
a%ia1(nza) = b%ia1(j)
a%ia2(nza) = b%ia2(j)
a%aspk(nza) = b%aspk(j)
endif
enddo
endif
a%infoa(psb_nnz_) = nza
else if(psb_toupper(b%fida)=='CSR') then
do i=1, min(nr-a%m,b%m)
do
jb = b%ia2(i)
if (jb >= b%ia2(i+1)) exit
nza = nza + 1
a%aspk(nza) = b%aspk(jb)
a%ia1(nza) = a%m + i
a%ia2(nza) = b%ia1(jb)
jb = jb + 1
end do
end do
a%infoa(psb_nnz_) = nza
call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale)
else
call psb_rwextd(nr,aa,info,rowscale=rowscale)
end if
class default
call aa%mv_to_coo(actmp,info)
if (info == 0) then
if (present(b)) then
call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale)
else
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
endif
endif
a%m = nr
a%k = max(a%k,b%k)
else if (a%fida == 'JAD') then
info=135
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
else
info=136
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_rwextd(nr,actmp,info,rowscale=rowscale)
end if
end if
if (info == 0) call aa%mv_from_coo(actmp,info)
end select
end if
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
@@ -172,3 +100,158 @@ subroutine psb_zrwextd(nr,a,info,b,rowscale)
return
end subroutine psb_zrwextd
subroutine psb_zbase_rwextd(nr,a,info,b,rowscale)
use psb_error_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_zbase_rwextd
implicit none
! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
integer, intent(in) :: nr
class(psb_z_base_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
class(psb_z_base_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer :: i,j,ja,jb,err_act,nza,nzb, ma, mb, na, nb
character(len=20) :: name, ch_err
logical rowscale_
name='psb_zbase_rwextd'
info = 0
call psb_erractionsave(err_act)
if (present(rowscale)) then
rowscale_ = rowscale
else
rowscale_ = .true.
end if
ma = a%get_nrows()
na = a%get_ncols()
select type(a)
type is (psb_z_csr_sparse_mat)
call psb_ensure_size(nr+1,a%irp,info)
if (present(b)) then
mb = b%get_nrows()
nb = b%get_ncols()
nzb = b%get_nzeros()
select type (b)
type is (psb_z_csr_sparse_mat)
call psb_ensure_size(size(a%ja)+nzb,a%ja,info)
call psb_ensure_size(size(a%val)+nzb,a%val,info)
do i=1, min(nr-ma,mb)
a%irp(ma+i+1) = a%irp(ma+i) + b%irp(i+1) - b%irp(i)
ja = a%irp(ma+i)
jb = b%irp(i)
do
if (jb >= b%irp(i+1)) exit
a%val(ja) = b%val(jb)
a%ja(ja) = b%ja(jb)
ja = ja + 1
jb = jb + 1
end do
end do
do j=i,nr-ma
a%irp(ma+i+1) = a%irp(ma+i)
end do
class default
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
end select
call a%set_ncols(max(na,nb))
else
do i=ma+2,nr+1
a%irp(i) = a%irp(i-1)
end do
end if
call a%set_nrows(nr)
type is (psb_z_coo_sparse_mat)
nza = a%get_nzeros()
if (present(b)) then
mb = b%get_nrows()
nb = b%get_ncols()
nzb = b%get_nzeros()
call a%reallocate(nza+nzb)
select type(b)
type is (psb_z_coo_sparse_mat)
if (rowscale_) then
do j=1,nzb
if ((ma + b%ia(j)) <= nr) then
nza = nza + 1
a%ia(nza) = ma + b%ia(j)
a%ja(nza) = b%ja(j)
a%val(nza) = b%val(j)
end if
enddo
else
do j=1,nzb
if ((ma + b%ia(j)) <= nr) then
nza = nza + 1
a%ia(nza) = b%ia(j)
a%ja(nza) = b%ja(j)
a%val(nza) = b%val(j)
end if
enddo
endif
call a%set_nzeros(nza)
type is (psb_z_csr_sparse_mat)
do i=1, min(nr-ma,mb)
do
jb = b%irp(i)
if (jb >= b%irp(i+1)) exit
nza = nza + 1
a%val(nza) = b%val(jb)
a%ia(nza) = ma + i
a%ja(nza) = b%ja(jb)
jb = jb + 1
end do
end do
call a%set_nzeros(nza)
class default
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
end select
call a%set_ncols(max(na,nb))
endif
call a%set_nrows(nr)
class default
info = 135
ch_err=a%get_fmt()
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_zbase_rwextd
base/serial/psb_zsymbmm.f90
+177 -95
View File
@@ -40,62 +40,107 @@
!
subroutine psb_zsymbmm(a,b,c,info)
use psb_spmat_type
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_zsymbmm
implicit none
type(psb_zspmat_type) :: a,b,c
type(psb_z_sparse_mat), intent(in) :: a,b
type(psb_z_sparse_mat), intent(out) :: c
integer, intent(out) :: info
type(psb_z_csr_sparse_mat), allocatable :: ccsr
integer :: err_act
character(len=*), parameter :: name='psb_symbmm'
call psb_erractionsave(err_act)
info = 0
if ((a%is_null()) .or.(b%is_null())) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
allocate(ccsr, stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call psb_symbmm(a%a,b%a,ccsr,info)
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call move_alloc(ccsr,c%a)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_zsymbmm
subroutine psb_zbase_symbmm(a,b,c,info)
use psb_mat_mod
use psb_serial_mod, psb_protect_name => psb_zbase_symbmm
implicit none
class(psb_z_base_sparse_mat), intent(in) :: a,b
type(psb_z_csr_sparse_mat), intent(out) :: c
integer, intent(out) :: info
integer, allocatable :: itemp(:)
integer :: nze,info
interface
subroutine symbmm (n, m, l, ia, ja, diaga, &
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
character(len=20) :: name
integer :: err_act
logical :: csra, csrb
integer :: nze, ma,na,mb,nb
character(len=20) :: name
integer :: err_act
name='psb_symbmm'
call psb_erractionsave(err_act)
info = 0
csra = (psb_toupper(a%fida(1:3))=='CSR')
csrb = (psb_toupper(b%fida(1:3))=='CSR')
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
if (b%m /= a%k) then
write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k
if ( mb /= na ) then
write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
endif
allocate(itemp(max(a%m,a%k,b%m,b%k)),stat=info)
allocate(itemp(max(ma,na,mb,nb)),stat=info)
if (info /= 0) then
return
info = 4000
call psb_Errpush(info,name)
goto 9999
endif
nze = max(a%m+1,2*a%m)
call psb_sp_reall(c,nze,info)
!
! Note: we need to test whether there is a performance impact
! in not using the original Douglas & Bank code.
!
if (csra.and.csrb) then
call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,&
& b%ia2,b%ia1,0,&
& c%ia2,c%ia1,0,itemp)
else
call inner_symbmm(a,b,c,itemp,info)
endif
call psb_realloc(size(c%ia1),c%aspk,info)
select type(a)
type is (psb_z_csr_sparse_mat)
select type(b)
type is (psb_z_csr_sparse_mat)
call csr_symbmm(a,b,c,itemp,info)
class default
call gen_symbmm(a,b,c,itemp,info)
end select
class default
call gen_symbmm(a,b,c,itemp,info)
end select
c%pl(1) = 0
c%pr(1) = 0
c%m=a%m
c%k=b%k
c%fida='CSR'
c%descra='GUN'
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call psb_realloc(size(c%ja),c%val,info)
deallocate(itemp)
call psb_erractionrestore(err_act)
return
@@ -108,82 +153,119 @@ subroutine psb_zsymbmm(a,b,c,info)
return
contains
subroutine inner_symbmm(a,b,c,index,info)
type(psb_zspmat_type) :: a,b,c
subroutine csr_symbmm(a,b,c,itemp,info)
type(psb_z_csr_sparse_mat), intent(in) :: a,b
type(psb_z_csr_sparse_mat), intent(out) :: c
integer :: itemp(:)
integer, intent(out) :: info
interface
subroutine symbmm (n, m, l, ia, ja, diaga, &
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
integer :: nze, ma,na,mb,nb
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
nze = max(ma+1,2*ma)
call c%allocate(ma,nb,nze)
call symbmm(ma,na,nb,a%irp,a%ja,0,&
& b%irp,b%ja,0,&
& c%irp,c%ja,0,itemp)
end subroutine csr_symbmm
subroutine gen_symbmm(a,b,c,index,info)
class(psb_z_base_sparse_mat), intent(in) :: a,b
type(psb_z_csr_sparse_mat), intent(out) :: c
integer :: index(:),info
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
complex(psb_dpk_), allocatable :: aval(:),bval(:)
integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
integer :: nze, ma,na,mb,nb
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
n = a%m
m = a%k
l = b%k
nze = max(ma+1,2*ma)
call c%allocate(ma,nb,nze)
n = ma
m = na
l = nb
maxlmn = max(l,m,n)
allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
& aval(maxlmn),bval(maxlmn), stat=info)
& stat=info)
if (info /= 0) then
info = 4000
return
endif
if (size(c%ia2) < n+1) then
call psb_realloc(n+1,c%ia2,info)
endif
do i=1,maxlmn
index(i)=0
end do
c%ia2(1)=1
minlm = min(l,m)
minmn = min(m,n)
c%irp(1)=1
minlm = min(l,m)
minmn = min(m,n)
main: do i=1,n
istart=-1
length=0
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
do jj=1, nazr
j=iacl(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
else
if(index(ibcl(k)) == 0) then
index(ibcl(k))=istart
istart=ibcl(k)
length=length+1
endif
endif
end do
end do
main: do i=1,n
istart=-1
length=0
call a%csget(i,i,nazr,iarw,iacl,info)
do jj=1, nazr
c%ia2(i+1)=c%ia2(i)+length
if (c%ia2(i+1) > size(c%ia1)) then
if (n > (2*i)) then
nze = max(c%ia2(i+1), c%ia2(i)*((n+i-1)/i))
j=iacl(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
info = 1
return
endif
call b%csget(j,j,nbzr,ibrw,ibcl,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
info=2
return
else
nze = max(c%ia2(i+1), nint((dble(c%ia2(i))*(dble(n)/i))) )
endif
call psb_realloc(nze,c%ia1,info)
end if
do j= c%ia2(i),c%ia2(i+1)-1
c%ia1(j)=istart
istart=index(istart)
index(c%ia1(j))=0
if(index(ibcl(k)) == 0) then
index(ibcl(k))=istart
istart=ibcl(k)
length=length+1
endif
endif
end do
call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1))
index(i) = 0
end do main
end do
end subroutine inner_symbmm
c%irp(i+1)=c%irp(i)+length
end subroutine psb_zsymbmm
if (c%irp(i+1) > size(c%ja)) then
if (n > (2*i)) then
nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i))
else
nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) )
endif
call psb_realloc(nze,c%ja,info)
end if
do j= c%irp(i),c%irp(i+1)-1
c%ja(j)=istart
istart=index(istart)
index(c%ja(j))=0
end do
call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1))
index(i) = 0
end do main
end subroutine gen_symbmm
end subroutine psb_zbase_symbmm