base/modules/psb_error_mod.F90
 base/modules/psb_psblas_mod.f90
 base/newserial/psbn_base_mat_mod.f03
 base/newserial/psbn_d_base_mat_mod.f03
 base/newserial/psbn_d_coo_impl.f03
 base/newserial/psbn_d_csr_impl.f03
 base/newserial/psbn_d_csr_mat_mod.f03
 base/newserial/psbn_mat_mod.f03
 base/psblas/psb_dnrmi.f90
 base/psblas/psb_dspmm.f90
 base/psblas/psb_dspsm.f90
 base/tools/psb_dspalloc.f90
 prec/psb_dbjac_aply.f90
 prec/psb_dbjac_bld.f90
 prec/psb_dilu_fct.f90
 prec/psb_dprecbld.f90
 prec/psb_dprecinit.f90
 prec/psb_prec_mod.f90
 prec/psb_prec_type.f90
 test/pargen/ppde.f90
 test/pargen/runs/ppde.inp

Now both BJAC_BLD and CSSV work. Really! 
And initial performance is not too bad. 
Lots and lots of details to be fixed yet...........
psblas3-type-indexed
Salvatore Filippone 16 years ago
parent ffe5ab739d
commit 4ecc1b632d

@ -386,9 +386,15 @@ contains
case(30)
write (0,'("input argument n. ",i0," has an invalid value")')i_e_d(1)
write (0,'("current value is ",i0)')i_e_d(2)
case(31)
write (0,'("input argument n. ",i0," has an invalid value")')i_e_d(1)
write (0,'("current value is ",a)')a_e_d
case(35)
write (0,'("Size of input array argument n. ",i0," is invalid.")')i_e_d(1)
write (0,'("Current value is ",i0)')i_e_d(2)
case(36)
write (0,'("Size of input array argument n. ",i0," must be ")')i_e_d(1)
write (0,'("at least ",i0)')i_e_d(2)
case(40)
write (0,'("input argument n. ",i0," has an invalid value")')i_e_d(1)
write (0,'("current value is ",a)')a_e_d(2:2)

@ -771,7 +771,8 @@ module psb_psblas_mod
& diag, n, jx, jy, work)
use psb_serial_mod
use psb_descriptor_type
type(psb_dspmat_type), intent(in) :: t
use psbn_d_mat_mod
type(psbn_d_sparse_mat), intent(in) :: t
real(psb_dpk_), intent(in) :: x(:,:)
real(psb_dpk_), intent(inout) :: y(:,:)
real(psb_dpk_), intent(in) :: alpha, beta
@ -787,7 +788,8 @@ module psb_psblas_mod
& diag, work)
use psb_serial_mod
use psb_descriptor_type
type(psb_dspmat_type), intent(in) :: t
use psbn_d_mat_mod
type(psbn_d_sparse_mat), intent(in) :: t
real(psb_dpk_), intent(in) :: x(:)
real(psb_dpk_), intent(inout) :: y(:)
real(psb_dpk_), intent(in) :: alpha, beta

@ -88,7 +88,8 @@ module psbn_base_mat_mod
generic, public :: allocate => allocate_mnnz
generic, public :: reallocate => reallocate_nz
procedure, pass(a) :: print => sparse_print
procedure, pass(a) :: sizeof
end type psbn_base_sparse_mat
private :: set_nrows, set_ncols, set_dupl, set_state, &
@ -97,10 +98,21 @@ module psbn_base_mat_mod
& get_nzeros, get_size, get_state, get_dupl, is_null, is_bld, &
& is_upd, is_asb, is_sorted, is_upper, is_lower, is_triangle, &
& is_unit, get_neigh, allocate_mn, allocate_mnnz, reallocate_nz, &
& free, sparse_print, get_fmt, trim
& free, sparse_print, get_fmt, trim, sizeof
contains
function sizeof(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
integer(psb_long_int_k_) :: res
res = 8
end function sizeof
function get_fmt(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a

@ -9,13 +9,17 @@ module psbn_d_base_mat_mod
generic, public :: csmm => d_base_csmm, d_base_csmv
procedure, pass(a) :: d_base_cssv
procedure, pass(a) :: d_base_cssm
generic, public :: cssm => d_base_cssm, d_base_cssv
generic, public :: base_cssm => d_base_cssm, d_base_cssv
procedure, pass(a) :: d_cssv
procedure, pass(a) :: d_cssm
generic, public :: cssm => d_cssm, d_cssv
procedure, pass(a) :: d_scals
procedure, pass(a) :: d_scal
generic, public :: scal => d_scals, d_scal
procedure, pass(a) :: get_diag
procedure, pass(a) :: csnmi
procedure, pass(a) :: get_diag
procedure, pass(a) :: csput
procedure, pass(a) :: d_csgetrow
procedure, pass(a) :: d_csgetblk
generic, public :: csget => d_csgetrow, d_csgetblk
@ -34,7 +38,7 @@ module psbn_d_base_mat_mod
& d_scals, d_scal, csnmi, csput, d_csgetrow, d_csgetblk, &
& cp_to_coo, cp_from_coo, cp_to_fmt, cp_from_fmt, &
& mv_to_coo, mv_from_coo, mv_to_fmt, mv_from_fmt, &
& get_diag, csclip
& get_diag, csclip, d_cssv, d_cssm
type, extends(psbn_d_base_sparse_mat) :: psbn_d_coo_sparse_mat
@ -73,6 +77,7 @@ module psbn_d_base_mat_mod
procedure, pass(a) :: d_csgetrow => d_coo_csgetrow
procedure, pass(a) :: print => d_coo_print
procedure, pass(a) :: get_fmt => d_coo_get_fmt
procedure, pass(a) :: sizeof => d_coo_sizeof
end type psbn_d_coo_sparse_mat
@ -82,7 +87,7 @@ module psbn_d_base_mat_mod
& d_fix_coo, d_coo_free, d_coo_print, d_coo_get_fmt, &
& d_cp_coo_to_coo, d_cp_coo_from_coo, &
& d_cp_coo_to_fmt, d_cp_coo_from_fmt, &
& d_coo_scals, d_coo_scal, d_coo_csgetrow
& d_coo_scals, d_coo_scal, d_coo_csgetrow, d_coo_sizeof
interface
@ -838,6 +843,265 @@ contains
end subroutine d_base_cssv
subroutine d_cssm(alpha,a,x,beta,y,info,trans,side,d)
use psb_error_mod
use psb_string_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:,:)
real(psb_dpk_), intent(inout) :: y(:,:)
integer, intent(out) :: info
character, optional, intent(in) :: trans, side
real(psb_dpk_), intent(in), optional :: d(:)
real(psb_dpk_), allocatable :: tmp(:,:)
Integer :: err_act, nar,nac,nc, i
character(len=1) :: side_
character(len=20) :: name='d_cssm'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
if (.not.a%is_asb()) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
nar = a%get_nrows()
nac = a%get_ncols()
nc = min(size(x,2), size(y,2))
if (size(x,1) < nac) then
info = 36
call psb_errpush(info,name,i_err=(/3,nac,0,0,0/))
goto 9999
end if
if (size(y,1) < nar) then
info = 36
call psb_errpush(info,name,i_err=(/3,nar,0,0,0/))
goto 9999
end if
if (.not. (a%is_triangle())) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
if (present(d)) then
if (present(side)) then
side_ = side
else
side_ = 'L'
end if
if (psb_toupper(side_) == 'R') then
if (size(d,1) < nac) then
info = 36
call psb_errpush(info,name,i_err=(/9,nac,0,0,0/))
goto 9999
end if
allocate(tmp(nac,nc),stat=info)
if (info /= 0) info = 4000
if (info == 0) then
do i=1, nac
tmp(i,1:nc) = d(i)*x(i,1:nc)
end do
end if
if (info == 0)&
& call a%base_cssm(alpha,tmp,beta,y,info,trans)
if (info == 0) then
deallocate(tmp,stat=info)
if (info /= 0) info = 4000
end if
else if (psb_toupper(side_) == 'L') then
if (size(d,1) < nar) then
info = 36
call psb_errpush(info,name,i_err=(/9,nar,0,0,0/))
goto 9999
end if
allocate(tmp(nar,nc),stat=info)
if (info /= 0) info = 4000
if (info == 0)&
& call a%base_cssm(done,x,dzero,tmp,info,trans)
if (info == 0)then
do i=1, nar
tmp(i,1:nc) = d(i)*tmp(i,1:nc)
end do
end if
if (info == 0)&
& call daxpby(nar,nc,alpha,tmp,size(tmp,1),beta,y,size(y,1),info)
if (info == 0) then
deallocate(tmp,stat=info)
if (info /= 0) info = 4000
end if
else
info = 31
call psb_errpush(info,name,i_err=(/8,0,0,0,0/),a_err=side_)
goto 9999
end if
else
! Side is ignored in this case
call a%base_cssm(alpha,x,beta,y,info,trans)
end if
if (info /= 0) then
info = 4010
call psb_errpush(info,name, a_err='base_cssm')
goto 9999
end if
return
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 d_cssm
subroutine d_cssv(alpha,a,x,beta,y,info,trans,side,d)
use psb_error_mod
use psb_string_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:)
real(psb_dpk_), intent(inout) :: y(:)
integer, intent(out) :: info
character, optional, intent(in) :: trans, side
real(psb_dpk_), intent(in), optional :: d(:)
real(psb_dpk_), allocatable :: tmp(:)
Integer :: err_act, nar,nac,nc, i
character(len=1) :: side_
character(len=20) :: name='d_cssm'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
if (.not.a%is_asb()) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
nar = a%get_nrows()
nac = a%get_ncols()
nc = 1
if (size(x,1) < nac) then
info = 36
call psb_errpush(info,name,i_err=(/3,nac,0,0,0/))
goto 9999
end if
if (size(y,1) < nar) then
info = 36
call psb_errpush(info,name,i_err=(/3,nar,0,0,0/))
goto 9999
end if
if (.not. (a%is_triangle())) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
if (present(d)) then
if (present(side)) then
side_ = side
else
side_ = 'L'
end if
if (psb_toupper(side_) == 'R') then
if (size(d,1) < nac) then
info = 36
call psb_errpush(info,name,i_err=(/9,nac,0,0,0/))
goto 9999
end if
allocate(tmp(nac),stat=info)
if (info /= 0) info = 4000
if (info == 0) tmp(1:nac) = d(1:nac)*x(1:nac)
if (info == 0)&
& call a%base_cssm(alpha,tmp,beta,y,info,trans)
if (info == 0) then
deallocate(tmp,stat=info)
if (info /= 0) info = 4000
end if
else if (psb_toupper(side_) == 'L') then
if (size(d,1) < nar) then
info = 36
call psb_errpush(info,name,i_err=(/9,nar,0,0,0/))
goto 9999
end if
allocate(tmp(nar),stat=info)
if (info /= 0) info = 4000
if (info == 0)&
& call a%base_cssm(done,x,dzero,tmp,info,trans)
if (info == 0) tmp(1:nar) = d(1:nar)*tmp(1:nar)
if (info == 0)&
& call daxpby(nar,nc,alpha,tmp,size(tmp,1),beta,y,size(y,1),info)
if (info == 0) then
deallocate(tmp,stat=info)
if (info /= 0) info = 4000
end if
else
info = 31
call psb_errpush(info,name,i_err=(/8,0,0,0,0/),a_err=side_)
goto 9999
end if
else
! Side is ignored in this case
call a%base_cssm(alpha,x,beta,y,info,trans)
end if
if (info /= 0) then
info = 4010
call psb_errpush(info,name, a_err='base_cssm')
goto 9999
end if
return
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 d_cssv
subroutine d_scals(d,a,info)
use psb_error_mod
@ -908,7 +1172,7 @@ contains
! so we throw an error.
info = 700
call psb_errpush(info,name,a_err=a%get_fmt())
write(0,*) 'Got into error path',err_act,psb_act_ret_
if (err_act /= psb_act_ret_) then
call psb_error()
end if
@ -961,6 +1225,18 @@ contains
!====================================
function d_coo_sizeof(a) result(res)
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
integer(psb_long_int_k_) :: res
res = 8 + 1
res = res + psb_sizeof_dp * size(a%val)
res = res + psb_sizeof_int * size(a%ia)
res = res + psb_sizeof_int * size(a%ja)
end function d_coo_sizeof
function d_coo_get_fmt(a) result(res)
implicit none
@ -1597,7 +1873,6 @@ contains
call psb_errpush(info,name,i_err=(/3,0,0,0,0/))
goto 9999
endif
if (info == 0) call psb_realloc(nz_,a%ia,info)
if (info == 0) call psb_realloc(nz_,a%ja,info)
if (info == 0) call psb_realloc(nz_,a%val,info)
@ -1721,7 +1996,7 @@ contains
character, optional, intent(in) :: trans
character :: trans_
integer :: i,j,k,m,n, nnz, ir, jc
integer :: i,j,k,m,n, nnz, ir, jc, nac, nar
real(psb_dpk_) :: acc
logical :: tra
Integer :: err_act
@ -1735,7 +2010,19 @@ contains
call psb_errpush(info,name)
goto 9999
endif
nar = a%get_nrows()
nac = a%get_ncols()
if (size(x) < nac) then
info = 36
call psb_errpush(info,name,i_err=(/3,nac,0,0,0/))
goto 9999
end if
if (size(y) < nar) then
info = 36
call psb_errpush(info,name,i_err=(/3,nar,0,0,0/))
goto 9999
end if
call d_coo_csmm_impl(alpha,a,x,beta,y,info,trans)
@ -1765,7 +2052,7 @@ contains
character, optional, intent(in) :: trans
character :: trans_
integer :: i,j,k,m,n, nnz, ir, jc, nc
integer :: i,j,k,m,n, nnz, ir, jc, nc, nar, nac
real(psb_dpk_), allocatable :: acc(:)
logical :: tra
Integer :: err_act
@ -1775,7 +2062,24 @@ contains
call psb_erractionsave(err_act)
if (.not.a%is_asb()) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
nar = a%get_nrows()
nac = a%get_ncols()
if (size(x,1) < nac) then
info = 36
call psb_errpush(info,name,i_err=(/3,nac,0,0,0/))
goto 9999
end if
if (size(y,1) < nar) then
info = 36
call psb_errpush(info,name,i_err=(/3,nar,0,0,0/))
goto 9999
end if
call d_coo_csmm_impl(alpha,a,x,beta,y,info,trans)
if (info /= 0) goto 9999
@ -1805,7 +2109,7 @@ contains
character, optional, intent(in) :: trans
character :: trans_
integer :: i,j,k,m,n, nnz, ir, jc
integer :: i,j,k,m,n, nnz, ir, jc, nar, nac
real(psb_dpk_) :: acc
real(psb_dpk_), allocatable :: tmp(:)
logical :: tra
@ -1821,9 +2125,21 @@ contains
goto 9999
endif
nar = a%get_nrows()
nac = a%get_ncols()
if (size(x,1) < nac) then
info = 36
call psb_errpush(info,name,i_err=(/3,nac,0,0,0/))
goto 9999
end if
if (size(y,1) < nar) then
info = 36
call psb_errpush(info,name,i_err=(/3,nar,0,0,0/))
goto 9999
end if
if (.not. (a%is_triangle())) then
write(0,*) 'Called SM on a non-triangular mat!'
info = 1121
call psb_errpush(info,name)
goto 9999
@ -1859,7 +2175,7 @@ contains
character, optional, intent(in) :: trans
character :: trans_
integer :: i,j,k,m,n, nnz, ir, jc, nc
integer :: i,j,k,m,n, nnz, ir, jc, nc, nar, nac
real(psb_dpk_) :: acc
real(psb_dpk_), allocatable :: tmp(:,:)
logical :: tra
@ -1875,9 +2191,21 @@ contains
goto 9999
endif
nar = a%get_nrows()
nac = a%get_ncols()
if (size(x,1) < nac) then
info = 36
call psb_errpush(info,name,i_err=(/3,nac,0,0,0/))
goto 9999
end if
if (size(y,1) < nar) then
info = 36
call psb_errpush(info,name,i_err=(/3,nar,0,0,0/))
goto 9999
end if
if (.not. (a%is_triangle())) then
write(0,*) 'Called SM on a non-triangular mat!'
info = 1121
call psb_errpush(info,name)
goto 9999

@ -29,7 +29,6 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
if (.not. (a%is_triangle())) then
write(0,*) 'Called SM on a non-triangular mat!'
info = 1121
call psb_errpush(info,name)
goto 9999
@ -48,20 +47,20 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
if (alpha == dzero) then
if (beta == dzero) then
do i = 1, m
y(i,:) = dzero
y(i,1:nc) = dzero
enddo
else
do i = 1, m
y(i,:) = beta*y(i,:)
y(i,1:nc) = beta*y(i,1:nc)
end do
endif
return
end if
if (beta == dzero) then
call inner_coosm(tra,a,x,y,info)
call inner_coosm(tra,a,x(:,1:nc),y(:,1:nc),info)
do i = 1, m
y(i,:) = alpha*y(i,:)
y(i,1:nc) = alpha*y(i,1:nc)
end do
else
allocate(tmp(m,nc), stat=info)
@ -71,10 +70,10 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
goto 9999
end if
tmp(1:m,:) = x(1:m,:)
call inner_coosm(tra,a,tmp,y,info)
tmp(1:m,1:nc) = x(1:m,1:nc)
call inner_coosm(tra,a,tmp(:,1:nc),y(:,1:nc),info)
do i = 1, m
y(i,:) = alpha*tmp(i,:) + beta*y(i,:)
y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc)
end do
end if
@ -135,7 +134,7 @@ contains
do
if (j > nnz) exit
if (a%ia(j) > i) exit
acc = acc + a%val(j)*x(a%ja(j),:)
acc = acc + a%val(j)*y(a%ja(j),:)
j = j + 1
end do
y(i,:) = x(i,:) - acc
@ -152,7 +151,7 @@ contains
j = j + 1
exit
end if
acc = acc + a%val(j)*x(a%ja(j),:)
acc = acc + a%val(j)*y(a%ja(j),:)
j = j + 1
end do
end do
@ -185,7 +184,7 @@ contains
j = j - 1
exit
end if
acc = acc + a%val(j)*x(a%ja(j),:)
acc = acc + a%val(j)*y(a%ja(j),:)
j = j - 1
end do
end do
@ -396,7 +395,7 @@ contains
do
if (j > nnz) exit
if (a%ia(j) > i) exit
acc = acc + a%val(j)*x(a%ja(j))
acc = acc + a%val(j)*y(a%ja(j))
j = j + 1
end do
y(i) = x(i) - acc
@ -413,7 +412,7 @@ contains
j = j + 1
exit
end if
acc = acc + a%val(j)*x(a%ja(j))
acc = acc + a%val(j)*y(a%ja(j))
j = j + 1
end do
end do
@ -427,7 +426,7 @@ contains
do
if (j < 1) exit
if (a%ia(j) < i) exit
acc = acc + a%val(j)*x(a%ja(j))
acc = acc + a%val(j)*y(a%ja(j))
j = j - 1
end do
y(i) = x(i) - acc
@ -446,7 +445,7 @@ contains
j = j - 1
exit
end if
acc = acc + a%val(j)*x(a%ja(j))
acc = acc + a%val(j)*y(a%ja(j))
j = j - 1
end do
end do
@ -747,11 +746,11 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans)
if (alpha == dzero) then
if (beta == dzero) then
do i = 1, m
y(i,:) = dzero
y(i,1:nc) = dzero
enddo
else
do i = 1, m
y(i,:) = beta*y(i,:)
y(i,1:nc) = beta*y(i,1:nc)
end do
endif
return
@ -759,28 +758,28 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans)
if (.not.a%is_unit()) then
if (beta == dzero) then
do i = 1, m
y(i,:) = dzero
y(i,1:nc) = dzero
enddo
else
do i = 1, m
y(i,:) = beta*y(i,:)
y(i,1:nc) = beta*y(i,1:nc)
end do
endif
else if (a%is_unit()) then
if (beta == dzero) then
do i = 1, min(m,n)
y(i,:) = alpha*x(i,:)
y(i,1:nc) = alpha*x(i,1:nc)
enddo
do i = min(m,n)+1, m
y(i,:) = dzero
y(i,1:nc) = dzero
enddo
else
do i = 1, min(m,n)
y(i,:) = beta*y(i,:) + alpha*x(i,:)
y(i,1:nc) = beta*y(i,1:nc) + alpha*x(i,1:nc)
end do
do i = min(m,n)+1, m
y(i,:) = beta*y(i,:)
y(i,1:nc) = beta*y(i,1:nc)
enddo
endif
@ -796,15 +795,15 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans)
acc = dzero
do
if (i>nnz) then
y(ir,:) = y(ir,:) + alpha * acc
y(ir,1:nc) = y(ir,1:nc) + alpha * acc
exit
endif
if (a%ia(i) /= ir) then
y(ir,:) = y(ir,:) + alpha * acc
y(ir,1:nc) = y(ir,1:nc) + alpha * acc
ir = a%ia(i)
acc = dzero
endif
acc = acc + a%val(i) * x(a%ja(i),:)
acc = acc + a%val(i) * x(a%ja(i),1:nc)
i = i + 1
enddo
end if
@ -815,7 +814,7 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans)
do i=1,nnz
ir = a%ja(i)
jc = a%ia(i)
y(ir,:) = y(ir,:) + a%val(i)*x(jc,:)
y(ir,1:nc) = y(ir,1:nc) + a%val(i)*x(jc,1:nc)
enddo
else if (alpha.eq.-done) then
@ -823,7 +822,7 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans)
do i=1,nnz
ir = a%ja(i)
jc = a%ia(i)
y(ir,:) = y(ir,:) - a%val(i)*x(jc,:)
y(ir,1:nc) = y(ir,1:nc) - a%val(i)*x(jc,1:nc)
enddo
else
@ -831,7 +830,7 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans)
do i=1,nnz
ir = a%ja(i)
jc = a%ia(i)
y(ir,:) = y(ir,:) + alpha*a%val(i)*x(jc,:)
y(ir,1:nc) = y(ir,1:nc) + alpha*a%val(i)*x(jc,1:nc)
enddo
end if !.....end testing on alpha
@ -861,7 +860,6 @@ function d_coo_csnmi_impl(a) result(res)
integer :: i,j,k,m,n, nnz, ir, jc, nc
real(psb_dpk_) :: acc
real(psb_dpk_), allocatable :: tmp(:,:)
logical :: tra
Integer :: err_act
character(len=20) :: name='d_base_csnmi'
@ -1028,7 +1026,6 @@ contains
irw = imin
lrw = imax
if (irw<0) then
write(debug_unit,*) ' spgtrow Error : idx no good ',irw
info = 2
return
end if
@ -1237,13 +1234,11 @@ subroutine d_coo_csput_impl(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)
nza = a%get_nzeros()
isza = a%get_size()
!!$ write(0,*) 'On entry to csput_impl: ',nza
if (a%is_bld()) then
! Build phase. Must handle reallocations in a sensible way.
if (isza < (nza+nz)) then
call a%reallocate(max(nza+nz,int(1.5*isza)))
isza = a%get_size()
!!$ write(0,*) 'isza: ',isza,nza+nz
endif
call psb_inner_ins(nz,ia,ja,val,nza,a%ia,a%ja,a%val,isza,&
@ -1327,7 +1322,6 @@ contains
if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then
nza = nza + 1
if (nza > maxsz) then
write(0,*) 'err -92 ',nza,maxsz
info = -92
return
endif

@ -40,7 +40,6 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
end if
if (.not.a%is_asb()) then
write(0,*) 'Error: csmv called on an unassembled mat'
info = 1121
call psb_errpush(info,name)
goto 9999
@ -567,7 +566,6 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
info = 0
call psb_erractionsave(err_act)
if (present(trans)) then
trans_ = trans
else
@ -588,7 +586,7 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
goto 9999
end if
if (alpha == dzero) then
if (beta == dzero) then
do i = 1, m
@ -603,18 +601,29 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
end if
if (beta == dzero) then
call inner_csrsv(tra,a,x,y)
do i = 1, m
y(i) = alpha*y(i)
end do
!!$ call inner_csrsv(tra,a,x,y)
call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),&
& a%irp,a%ja,a%val,x,y)
if (alpha == done) then
! do nothing
else if (alpha == -done) then
do i = 1, m
y(i) = -y(i)
end do
else
do i = 1, m
y(i) = alpha*y(i)
end do
end if
else
allocate(tmp(m), stat=info)
if (info /= 0) then
write(0,*) 'Memory allocation error in CSRSV '
return
end if
tmp(1:m) = x(1:m)
call inner_csrsv(tra,a,tmp,y)
call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),&
& a%irp,a%ja,a%val,tmp,y)
!!$ call inner_csrsv(tra,a,tmp,y)
do i = 1, m
y(i) = alpha*tmp(i) + beta*y(i)
end do
@ -634,53 +643,56 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
contains
subroutine inner_csrsv(tra,a,x,y)
subroutine inner_csrsv(tra,lower,unit,n,irp,ja,val,x,y)
implicit none
logical, intent(in) :: tra
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: x(:)
real(psb_dpk_), intent(out) :: y(:)
logical, intent(in) :: tra,lower,unit
!!$ class(psbn_d_csr_sparse_mat), intent(in) :: a
integer, intent(in) :: irp(*), ja(*),n
real(psb_dpk_), intent(in) :: val(*)
real(psb_dpk_), intent(in) :: x(*)
real(psb_dpk_), intent(out) :: y(*)
integer :: i,j,k,m, ir, jc
real(psb_dpk_) :: acc
if (.not.tra) then
if (a%is_lower()) then
if (a%is_unit()) then
do i=1, a%get_nrows()
if (lower) then
if (unit) then
do i=1, n
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j)*x(a%ja(j))
do j=irp(i), irp(i+1)-1
acc = acc + val(j)*y(ja(j))
end do
y(i) = x(i) - acc
end do
else if (.not.a%is_unit()) then
do i=1, a%get_nrows()
else if (.not.unit) then
do i=1, n
acc = dzero
do j=a%irp(i), a%irp(i+1)-2
acc = acc + a%val(j)*x(a%ja(j))
do j=irp(i), irp(i+1)-2
acc = acc + val(j)*y(ja(j))
end do
y(i) = (x(i) - acc)/a%val(a%irp(i+1)-1)
y(i) = (x(i) - acc)/val(irp(i+1)-1)
end do
end if
else if (a%is_upper()) then
else if (.not.lower) then
if (a%is_unit()) then
do i=a%get_nrows(), 1, -1
if (unit) then
do i=n, 1, -1
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j)*x(a%ja(j))
do j=irp(i), irp(i+1)-1
acc = acc + val(j)*y(ja(j))
end do
y(i) = x(i) - acc
end do
else if (.not.a%is_unit()) then
do i=a%get_nrows(), 1, -1
else if (.not.unit) then
do i=n, 1, -1
acc = dzero
do j=a%irp(i)+1, a%irp(i+1)-1
acc = acc + a%val(j)*x(a%ja(j))
do j=irp(i)+1, irp(i+1)-1
acc = acc + val(j)*y(ja(j))
end do
y(i) = (x(i) - acc)/a%val(a%irp(i))
y(i) = (x(i) - acc)/val(irp(i))
end do
end if
@ -688,46 +700,46 @@ contains
else if (tra) then
do i=1, a%get_nrows()
do i=1, n
y(i) = x(i)
end do
if (a%is_lower()) then
if (a%is_unit()) then
do i=a%get_nrows(), 1, -1
if (lower) then
if (unit) then
do i=n, 1, -1
acc = y(i)
do j=a%irp(i), a%irp(i+1)-1
jc = a%ja(j)
y(jc) = y(jc) - a%val(j)*acc
do j=irp(i), irp(i+1)-1
jc = ja(j)
y(jc) = y(jc) - val(j)*acc
end do
end do
else if (.not.a%is_unit()) then
do i=a%get_nrows(), 1, -1
y(i) = y(i)/a%val(a%irp(i+1)-1)
else if (.not.unit) then
do i=n, 1, -1
y(i) = y(i)/val(irp(i+1)-1)
acc = y(i)
do j=a%irp(i), a%irp(i+1)-2
jc = a%ja(j)
y(jc) = y(jc) - a%val(j)*acc
do j=irp(i), irp(i+1)-2
jc = ja(j)
y(jc) = y(jc) - val(j)*acc
end do
end do
end if
else if (a%is_upper()) then
else if (.not.lower) then
if (a%is_unit()) then
do i=1, a%get_nrows()
if (unit) then
do i=1, n
acc = y(i)
do j=a%irp(i), a%irp(i+1)-1
jc = a%ja(j)
y(jc) = y(jc) - a%val(j)*acc
do j=irp(i), irp(i+1)-1
jc = ja(j)
y(jc) = y(jc) - val(j)*acc
end do
end do
else if (.not.a%is_unit()) then
do i=1, a%get_nrows()
y(i) = y(i)/a%val(a%irp(i))
else if (.not.unit) then
do i=1, n
y(i) = y(i)/val(irp(i))
acc = y(i)
do j=a%irp(i)+1, a%irp(i+1)-1
jc = a%ja(j)
y(jc) = y(jc) - a%val(j)*acc
do j=irp(i)+1, irp(i+1)-1
jc = ja(j)
y(jc) = y(jc) - val(j)*acc
end do
end do
end if
@ -779,7 +791,6 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
nc = min(size(x,2) , size(y,2))
if (.not. (a%is_triangle())) then
write(0,*) 'Called SM on a non-triangular mat!'
info = 1121
call psb_errpush(info,name)
goto 9999
@ -800,9 +811,11 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
end if
if (beta == dzero) then
call inner_csrsm(tra,a,x,y,info)
!!$ call inner_csrsm(tra,a,x,y,info)
call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,&
& a%irp,a%ja,a%val,x,size(x,1),y,size(y,1),info)
do i = 1, m
y(i,:) = alpha*y(i,:)
y(i,1:nc) = alpha*y(i,1:nc)
end do
else
allocate(tmp(m,nc), stat=info)
@ -812,10 +825,12 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
goto 9999
end if
tmp(1:m,:) = x(1:m,:)
call inner_csrsm(tra,a,tmp,y,info)
tmp(1:m,:) = x(1:m,1:nc)
call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,&
& a%irp,a%ja,a%val,tmp,size(tmp,1),y,size(y,1),info)
!!$ call inner_csrsm(tra,a,tmp,y,info)
do i = 1, m
y(i,:) = alpha*tmp(i,:) + beta*y(i,:)
y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc)
end do
end if
@ -841,18 +856,19 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
contains
subroutine inner_csrsm(tra,a,x,y,info)
subroutine inner_csrsm(tra,lower,unit,nr,nc,&
& irp,ja,val,x,ldx,y,ldy,info)
implicit none
logical, intent(in) :: tra
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: x(:,:)
real(psb_dpk_), intent(out) :: y(:,:)
logical, intent(in) :: tra,lower,unit
integer, intent(in) :: nr,nc,ldx,ldy,irp(*),ja(*)
real(psb_dpk_), intent(in) :: val(*), x(ldx,*)
real(psb_dpk_), intent(out) :: y(ldy,*)
integer, intent(out) :: info
integer :: i,j,k,m, ir, jc
real(psb_dpk_), allocatable :: acc(:)
info = 0
allocate(acc(size(x,2)), stat=info)
allocate(acc(nc), stat=info)
if(info /= 0) then
info=4010
return
@ -861,41 +877,41 @@ contains
if (.not.tra) then
if (a%is_lower()) then
if (a%is_unit()) then
do i=1, a%get_nrows()
if (lower) then
if (unit) then
do i=1, nr
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j)*x(a%ja(j),:)
acc = acc + a%val(j)*y(a%ja(j),1:nc)
end do
y(i,:) = x(i,:) - acc
y(i,1:nc) = x(i,1:nc) - acc
end do
else if (.not.a%is_unit()) then
do i=1, a%get_nrows()
else if (.not.unit) then
do i=1, nr
acc = dzero
do j=a%irp(i), a%irp(i+1)-2
acc = acc + a%val(j)*x(a%ja(j),:)
acc = acc + a%val(j)*y(a%ja(j),1:nc)
end do
y(i,:) = (x(i,:) - acc)/a%val(a%irp(i+1)-1)
y(i,1:nc) = (x(i,1:nc) - acc)/a%val(a%irp(i+1)-1)
end do
end if
else if (a%is_upper()) then
else if (.not.lower) then
if (a%is_unit()) then
do i=a%get_nrows(), 1, -1
if (unit) then
do i=nr, 1, -1
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j)*x(a%ja(j),:)
acc = acc + a%val(j)*y(a%ja(j),1:nc)
end do
y(i,:) = x(i,:) - acc
y(i,1:nc) = x(i,1:nc) - acc
end do
else if (.not.a%is_unit()) then
do i=a%get_nrows(), 1, -1
else if (.not.unit) then
do i=nr, 1, -1
acc = dzero
do j=a%irp(i)+1, a%irp(i+1)-1
acc = acc + a%val(j)*x(a%ja(j),:)
acc = acc + a%val(j)*y(a%ja(j),1:nc)
end do
y(i,:) = (x(i,:) - acc)/a%val(a%irp(i))
y(i,1:nc) = (x(i,1:nc) - acc)/a%val(a%irp(i))
end do
end if
@ -903,46 +919,46 @@ contains
else if (tra) then
do i=1, a%get_nrows()
y(i,:) = x(i,:)
do i=1, nr
y(i,1:nc) = x(i,1:nc)
end do
if (a%is_lower()) then
if (a%is_unit()) then
do i=a%get_nrows(), 1, -1
acc = y(i,:)
if (lower) then
if (unit) then
do i=nr, 1, -1
acc = y(i,1:nc)
do j=a%irp(i), a%irp(i+1)-1
jc = a%ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc
y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc
end do
end do
else if (.not.a%is_unit()) then
do i=a%get_nrows(), 1, -1
y(i,:) = y(i,:)/a%val(a%irp(i+1)-1)
acc = y(i,:)
else if (.not.unit) then
do i=nr, 1, -1
y(i,1:nc) = y(i,1:nc)/a%val(a%irp(i+1)-1)
acc = y(i,1:nc)
do j=a%irp(i), a%irp(i+1)-2
jc = a%ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc
y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc
end do
end do
end if
else if (a%is_upper()) then
else if (.not.lower) then
if (a%is_unit()) then
do i=1, a%get_nrows()
acc = y(i,:)
if (unit) then
do i=1, nr
acc = y(i,1:nc)
do j=a%irp(i), a%irp(i+1)-1
jc = a%ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc
y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc
end do
end do
else if (.not.a%is_unit()) then
do i=1, a%get_nrows()
y(i,:) = y(i,:)/a%val(a%irp(i))
acc = y(i,:)
else if (.not.unit) then
do i=1, nr
y(i,1:nc) = y(i,1:nc)/a%val(a%irp(i))
acc = y(i,1:nc)
do j=a%irp(i)+1, a%irp(i+1)-1
jc = a%ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc
y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc
end do
end do
end if
@ -1120,7 +1136,6 @@ contains
irw = imin
lrw = min(imax,a%get_nrows())
if (irw<0) then
write(debug_unit,*) ' spgtrow Error : idx no good ',irw
info = 2
return
end if
@ -1131,13 +1146,14 @@ contains
nzin_ = 0
endif
nzt = a%irp(lrw)-a%irp(irw)
nzt = a%irp(lrw+1)-a%irp(irw)
nz = 0
call psb_ensure_size(nzin_+nzt,ia,info)
if (info==0) call psb_ensure_size(nzin_+nzt,ja,info)
if (info==0) call psb_ensure_size(nzin_+nzt,val,info)
if (info /= 0) return
if (present(iren)) then
@ -1671,6 +1687,9 @@ subroutine d_mv_csr_to_fmt_impl(a,b,info)
select type (b)
class is (psbn_d_coo_sparse_mat)
call a%mv_to_coo(b,info)
! Need to fix trivial copies!
!!$ class is (psbn_d_csr_sparse_mat)
!!$ call a%mv_to_coo(b,info)
class default
call tmp%mv_from_fmt(a,info)
if (info == 0) call b%mv_from_coo(tmp,info)

@ -34,6 +34,7 @@ module psbn_d_csr_mat_mod
procedure, pass(a) :: free => d_csr_free
procedure, pass(a) :: trim => d_csr_trim
procedure, pass(a) :: print => d_csr_print
procedure, pass(a) :: sizeof => d_csr_sizeof
end type psbn_d_csr_sparse_mat
private :: d_csr_get_nzeros, d_csr_csmm, d_csr_csmv, d_csr_cssm, d_csr_cssv, &
& d_csr_csput, d_csr_reallocate_nz, d_csr_allocate_mnnz, &
@ -42,7 +43,8 @@ module psbn_d_csr_mat_mod
& d_mv_csr_to_coo, d_mv_csr_from_coo, &
& d_cp_csr_to_fmt, d_cp_csr_from_fmt, &
& d_mv_csr_to_fmt, d_mv_csr_from_fmt, &
& d_csr_scals, d_csr_scal, d_csr_trim, d_csr_csgetrow, d_csr_get_size
& d_csr_scals, d_csr_scal, d_csr_trim, d_csr_csgetrow, d_csr_get_size, &
& d_csr_sizeof
interface
@ -234,6 +236,18 @@ contains
!
!=====================================
function d_csr_sizeof(a) result(res)
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
integer(psb_long_int_k_) :: res
res = 8
res = res + psb_sizeof_dp * size(a%val)
res = res + psb_sizeof_int * size(a%irp)
res = res + psb_sizeof_int * size(a%ja)
end function d_csr_sizeof
function d_csr_get_fmt(a) result(res)
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
@ -1151,7 +1165,6 @@ contains
if (.not. (a%is_triangle())) then
write(0,*) 'Called SM on a non-triangular mat!'
info = 1121
call psb_errpush(info,name)
goto 9999
@ -1205,7 +1218,6 @@ contains
if (.not. (a%is_triangle())) then
write(0,*) 'Called SM on a non-triangular mat!'
info = 1121
call psb_errpush(info,name)
goto 9999

@ -40,6 +40,8 @@ module psbn_d_mat_mod
procedure, pass(a) :: is_triangle
procedure, pass(a) :: is_unit
procedure, pass(a) :: get_fmt => sparse_get_fmt
procedure, pass(a) :: sizeof => d_sizeof
! Memory/data management
procedure, pass(a) :: csall
@ -78,12 +80,32 @@ module psbn_d_mat_mod
& is_unit, get_neigh, csall, csput, d_csgetrow,&
& d_csgetblk, csclip, d_cscnv, d_cscnv_ip, &
& reallocate_nz, free, trim, &
& d_csmv, d_csmm, d_cssv, d_cssm, sparse_print, &
& sparse_print, &
& set_nrows, set_ncols, set_dupl, &
& set_state, set_null, set_bld, &
& set_upd, set_asb, set_sorted, &
& set_upper, set_lower, set_triangle, &
& set_unit, csnmi, get_diag, d_scals, d_scal
& set_unit, get_diag
interface psb_sizeof
module procedure d_sizeof
end interface
interface psbn_csmm
module procedure d_csmm, d_csmv
end interface
interface psbn_cssm
module procedure d_cssm, d_cssv
end interface
interface psbn_csnmi
module procedure csnmi
end interface
interface psbn_scal
module procedure d_scals, d_scal
end interface
contains
@ -100,6 +122,20 @@ contains
!
!=====================================
function d_sizeof(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
integer(psb_long_int_k_) :: res
res = 0
if (allocated(a%a)) then
res = a%a%sizeof()
end if
end function d_sizeof
function sparse_get_fmt(a) result(res)
implicit none
@ -1265,7 +1301,7 @@ contains
call move_alloc(altmp,b%a)
call b%set_asb()
call b%trim()
call psb_erractionrestore(err_act)
return
@ -1357,7 +1393,7 @@ contains
call move_alloc(altmp,a%a)
call a%set_asb()
call a%trim()
call psb_erractionrestore(err_act)
return
@ -1460,14 +1496,15 @@ contains
end subroutine d_csmv
subroutine d_cssm(alpha,a,x,beta,y,info,trans)
subroutine d_cssm(alpha,a,x,beta,y,info,trans,side,d)
use psb_error_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
real(kind(1.d0)), intent(in) :: alpha, beta, x(:,:)
real(kind(1.d0)), intent(inout) :: y(:,:)
integer, intent(out) :: info
character, optional, intent(in) :: trans
character, optional, intent(in) :: trans, side
real(psb_dpk_), intent(in), optional :: d(:)
Integer :: err_act
character(len=20) :: name='psbn_cssm'
logical, parameter :: debug=.false.
@ -1480,7 +1517,7 @@ contains
goto 9999
endif
call a%a%cssm(alpha,x,beta,y,info,trans)
call a%a%cssm(alpha,x,beta,y,info,trans,side,d)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
@ -1497,14 +1534,15 @@ contains
end subroutine d_cssm
subroutine d_cssv(alpha,a,x,beta,y,info,trans)
subroutine d_cssv(alpha,a,x,beta,y,info,trans,side,d)
use psb_error_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
real(kind(1.d0)), intent(in) :: alpha, beta, x(:)
real(kind(1.d0)), intent(inout) :: y(:)
integer, intent(out) :: info
character, optional, intent(in) :: trans
character, optional, intent(in) :: trans, side
real(psb_dpk_), intent(in), optional :: d(:)
Integer :: err_act
character(len=20) :: name='psbn_cssv'
logical, parameter :: debug=.false.
@ -1516,8 +1554,8 @@ contains
call psb_errpush(info,name)
goto 9999
endif
call a%a%cssm(alpha,x,beta,y,info,trans)
call a%a%cssm(alpha,x,beta,y,info,trans,side,d)
if (info /= 0) goto 9999

@ -95,7 +95,7 @@ function psb_dnrmi(a,desc_a,info)
end if
if ((m /= 0).and.(n /= 0)) then
nrmi = a%csnmi()
nrmi = psbn_csnmi(a)
if(info /= 0) then
info=4010
ch_err='psb_csnmi'

@ -251,7 +251,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
if(info /= 0) exit blk
! local Matrix-vector product
call a%csmm(alpha,x(:,jjx+i-1:jjx+i-1+ib-1),&
call psbn_csmm(alpha,a,x(:,jjx+i-1:jjx+i-1+ib-1),&
& beta,y(:,jjy+i-1:jjy+i-1+ib-1),info,trans=trans_)
if(info /= 0) exit blk
@ -266,7 +266,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
if (doswap_)&
& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& ib1,dzero,x(:,1:ik),desc_a,iwork,info)
if (info == 0) call a%csmm(alpha,x(:,1:ik),beta,y(:,1:ik),info)
if (info == 0) call psbn_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info)
end if
if(info /= 0) then
info = 4011
@ -313,7 +313,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
if (info == 0) call psi_ovrl_upd(x,desc_a,psb_avg_,info)
y(nrow+1:ncol,1:ik) = dzero
if (info == 0) call a%csmm(alpha,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_)
if (info == 0) call psbn_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_)
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' csmm ', info
if (info /= 0) then
@ -584,7 +584,7 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
& dzero,x,desc_a,iwork,info,data=psb_comm_halo_)
end if
! Just for fun
call a%csmm(alpha,x,beta,y,info)
call psbn_csmm(alpha,a,x,beta,y,info)
if(info /= 0) then
info = 4011
@ -633,7 +633,7 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
yp(nrow+1:ncol) = dzero
! local Matrix-vector product
if (info == 0) call a%csmm(alpha,x,beta,y,info,trans=trans_)
if (info == 0) call psbn_csmm(alpha,a,x,beta,y,info,trans=trans_)
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' csmm ', info

@ -64,7 +64,7 @@
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! unitd - character(optional). Specify some type of operation with
! side - character(optional). Specify some type of operation with
! the diagonal matrix D.
! choice - integer(optional). The kind of update to perform on overlap elements.
! d(:) - real , optional Matrix for diagonal scaling.
@ -75,7 +75,7 @@
!
!
subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
& trans, unitd, choice, diag, k, jx, jy, work)
& trans, side, choice, diag, k, jx, jy, work)
use psb_spmat_type
use psb_serial_mod
@ -86,17 +86,18 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
use psb_error_mod
use psb_string_mod
use psb_penv_mod
use psbn_d_mat_mod
implicit none
real(psb_dpk_), intent(in) :: alpha, beta
real(psb_dpk_), intent(in), target :: x(:,:)
real(psb_dpk_), intent(inout), target :: y(:,:)
type (psb_dspmat_type), intent(in) :: a
type(psbn_d_sparse_mat), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
real(psb_dpk_), intent(in), optional, target :: diag(:)
real(psb_dpk_), optional, target :: work(:)
character, intent(in), optional :: trans, unitd
character, intent(in), optional :: trans, side
integer, intent(in), optional :: choice
integer, intent(in), optional :: k, jx, jy
@ -106,7 +107,7 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
& ix, iy, ik, ijx, ijy, i, lld,&
& m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm
character :: lunitd
character :: lside
integer, parameter :: nb=4
real(psb_dpk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
character :: itrans
@ -158,10 +159,10 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
choice_ = psb_avg_
endif
if (present(unitd)) then
lunitd = psb_toupper(unitd)
if (present(side)) then
lside = psb_toupper(side)
else
lunitd = 'U'
lside = 'U'
endif
if (present(trans)) then
@ -192,8 +193,6 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
! check for presence/size of a work area
iwork => null()
liwork= 2*ncol
if (a%pr(1) /= 0) llwork = liwork + m * ik
if (a%pl(1) /= 0) llwork = llwork + m * ik
if (present(work)) then
if (size(work) >= liwork) then
aliw =.false.
@ -259,7 +258,7 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
! Perform local triangular system solve
xp => x(iix:lldx,jjx:jjx+ik-1)
yp => y(iiy:lldy,jjy:jjy+ik-1)
call a%cssm(alpha,xp,beta,yp,info,unitd=lunitd,d=id,trans=itrans)
call psbn_cssm(alpha,a,xp,beta,yp,info,side=side,d=diag,trans=trans)
if(info /= 0) then
info = 4010
@ -357,14 +356,14 @@ end subroutine psb_dspsm
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! trans - character(optional). Whether A or A'. If not present 'N' is assumed.
! unitd - character(optional). Specify some type of operation with
! side - character(optional). Specify some type of operation with
! the diagonal matrix D.
! choice - integer(optional). The kind of update to perform on overlap elements.
! d(:) - real , optional Matrix for diagonal scaling.
! work(:) - real , optional Working area.
!
subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
& trans, unitd, choice, diag, work)
& trans, side, choice, diag, work)
use psb_spmat_type
use psb_serial_mod
use psb_descriptor_type
@ -374,17 +373,18 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
use psb_error_mod
use psb_string_mod
use psb_penv_mod
use psbn_d_mat_mod
implicit none
real(psb_dpk_), intent(in) :: alpha, beta
real(psb_dpk_), intent(in), target :: x(:)
real(psb_dpk_), intent(inout), target :: y(:)
type(psb_dspmat_type), intent(in) :: a
type(psbn_d_sparse_mat), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
real(psb_dpk_), intent(in), optional, target :: diag(:)
real(psb_dpk_), optional, target :: work(:)
character, intent(in), optional :: trans, unitd
character, intent(in), optional :: trans, side
integer, intent(in), optional :: choice
! locals
@ -393,7 +393,7 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
& ix, iy, ik, jx, jy, i, lld,&
& m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm
character :: lunitd
character :: lside
integer, parameter :: nb=4
real(psb_dpk_),pointer :: iwork(:), xp(:), yp(:), id(:)
character :: itrans
@ -429,10 +429,10 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
choice_ = psb_avg_
endif
if (present(unitd)) then
lunitd = psb_toupper(unitd)
if (present(side)) then
lside = psb_toupper(side)
else
lunitd = 'U'
lside = 'U'
endif
if (present(trans)) then
@ -529,7 +529,8 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
! Perform local triangular system solve
xp => x(iix:lldx)
yp => y(iiy:lldy)
call a%cssm(alpha,xp,beta,yp,info,unitd=lunitd,d=id,trans=itrans)
call psbn_cssm(alpha,a,xp,beta,yp,info,side=side,d=diag,trans=trans)
!!$ call psbn_cssm(alpha,a,xp,beta,yp,info,side=side,d=id,trans=itrans)
if(info /= 0) then
info = 4010

@ -109,7 +109,7 @@ subroutine psb_dspalloc(a, desc_a, info, nnz)
& write(debug_unit,*) me,' ',trim(name),':allocating size:',length_ia1
!....allocate aspk, ia1, ia2.....
call a%csall(loc_row,loc_col,length_ia1,info)
call a%csall(loc_row,loc_col,info,nz=length_ia1)
if(info /= 0) then
info=4010
ch_err='sp_all'

@ -38,6 +38,7 @@ subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
use psb_base_mod
use psbn_d_mat_mod
use psb_prec_mod, psb_protect_name => psb_dbjac_aply
implicit none

@ -30,6 +30,7 @@
!!$
!!$
subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
use psbn_d_mat_mod
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_dbjac_bld
implicit none
@ -37,7 +38,7 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
! .. Scalar Arguments ..
integer, intent(out) :: info
! .. array Arguments ..
type(psb_dspmat_type), intent(in), target :: a
type(psbn_d_sparse_mat), intent(in), target :: a
type(psb_dprec_type), intent(inout) :: p
type(psb_desc_type), intent(in) :: desc_a
character, intent(in) :: upd
@ -47,6 +48,7 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
integer :: int_err(5)
character :: trans, unitd
type(psb_dspmat_type) :: atmp
type(psbn_d_csr_sparse_mat), allocatable :: lf, uf
real(psb_dpk_) :: t1,t2,t3,t4,t5,t6, t7, t8
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
@ -61,7 +63,7 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
m = a%m
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
@ -89,12 +91,12 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
if (allocated(p%av)) then
if (size(p%av) < psb_bp_ilu_avsz) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
! Actually, we don't care here about this.
! Just let it go.
! return
end if
call p%av(i)%free()
!!$ if (info /= 0) then
!!$ ! Actually, we don't care here about this.
!!$ ! Just let it go.
!!$ ! return
!!$ end if
enddo
deallocate(p%av,stat=info)
endif
@ -108,17 +110,19 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
endif
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = psb_sp_get_nnzeros(a)
nztota = a%get_nzeros()
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
n_row = p%desc_data%matrix_data(psb_n_row_)
p%av(psb_l_pr_)%m = n_row
p%av(psb_l_pr_)%k = n_row
p%av(psb_u_pr_)%m = n_row
p%av(psb_u_pr_)%k = n_row
call psb_sp_all(n_row,n_row,p%av(psb_l_pr_),nztota,info)
if (info == 0) call psb_sp_all(n_row,n_row,p%av(psb_u_pr_),nztota,info)
allocate(lf,uf,stat=info)
if (info == 0) call lf%allocate(n_row,n_row,nztota)
if (info == 0) call uf%allocate(n_row,n_row,nztota)
!!$ call p%av(psb_l_pr_)%csall(n_row,n_row,info,nztota)
!!$ if (info == 0) call p%av(psb_u_pr_)%csall(n_row,n_row,info,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
@ -140,24 +144,34 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
endif
t3 = psb_wtime()
! This is where we have mo renumbering, thus no need
! This is where we have no renumbering, thus no need
! for ATMP
call psb_ilu_fct(a,p%av(psb_l_pr_),p%av(psb_u_pr_),p%d,info)
if(info/=0) then
!!$ call p%av(psb_l_pr_)%cscnv(info,type='CSR')
!!$ call p%av(psb_u_pr_)%cscnv(info,type='CSR')
call psb_ilu_fct(a,lf,uf,p%d,info)
if(info==0) then
call move_alloc(lf,p%av(psb_l_pr_)%a)
call move_alloc(uf,p%av(psb_u_pr_)%a)
call p%av(psb_l_pr_)%set_asb()
call p%av(psb_u_pr_)%set_asb()
call p%av(psb_l_pr_)%trim()
call p%av(psb_u_pr_)%trim()
else
info=4010
ch_err='psb_ilu_fct'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (psb_sp_getifld(psb_upd_,p%av(psb_u_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(psb_u_pr_),info)
endif
if (psb_sp_getifld(psb_upd_,p%av(psb_l_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(psb_l_pr_),info)
endif
!!$
!!$ if (psb_sp_getifld(psb_upd_,p%av(psb_u_pr_),info) /= psb_upd_perm_) then
!!$ call psb_sp_trim(p%av(psb_u_pr_),info)
!!$ endif
!!$
!!$ if (psb_sp_getifld(psb_upd_,p%av(psb_l_pr_),info) /= psb_upd_perm_) then
!!$ call psb_sp_trim(p%av(psb_l_pr_),info)
!!$ endif
case(psb_f_none_)

@ -37,12 +37,13 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck)
!
!
use psb_base_mod
use psbn_d_mat_mod
implicit none
! .. Scalar Arguments ..
integer, intent(out) :: info
! .. Array Arguments ..
type(psb_dspmat_type),intent(in) :: a
type(psb_dspmat_type),intent(inout) :: l,u
type(psbn_d_sparse_mat),intent(in) :: a
type(psbn_d_csr_sparse_mat),intent(inout) :: l,u
type(psb_dspmat_type),intent(in), optional, target :: blck
real(psb_dpk_), intent(inout) :: d(:)
! .. Local Scalars ..
@ -76,25 +77,26 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck)
blck_%m=0
endif
call psb_dilu_fctint(m,a%m,a,blck_%m,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
call psb_dilu_fctint(m,a%get_nrows(),a,blck_%m,blck_,&
& d,l%val,l%ja,l%irp,u%val,u%ja,u%irp,l1,l2,info)
if(info /= 0) then
info=4010
ch_err='psb_dilu_fctint'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
l%infoa(1) = l1
l%fida = 'CSR'
l%descra = 'TLU'
u%infoa(1) = l2
u%fida = 'CSR'
u%descra = 'TUU'
l%m = m
l%k = m
u%m = m
u%k = m
call l%set_triangle()
call l%set_lower()
call l%set_unit()
call u%set_triangle()
call u%set_upper()
call u%set_unit()
call l%set_nrows(m)
call l%set_ncols(m)
call u%set_nrows(m)
call u%set_ncols(m)
if (present(blck)) then
blck_ => null()
else
@ -124,17 +126,22 @@ contains
& d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
implicit none
type(psb_dspmat_type) :: a,b
type(psbn_d_sparse_mat) :: a
type(psb_dspmat_type) :: b
integer :: m,ma,mb,l1,l2,info
integer, dimension(:) :: lia1,lia2,uia1,uia2
real(psb_dpk_), dimension(:) :: laspk,uaspk,d
integer :: i,j,k,l,low1,low2,kk,jj,ll, irb, ktrw,err_act
integer :: i,j,k,l,low1,low2,kk,jj,ll, irb, ktrw,err_act, nz
real(psb_dpk_) :: dia,temp
integer, parameter :: nrb=16
type(psb_dspmat_type) :: trw
integer, allocatable :: irow(:), icol(:)
real(psb_dpk_), allocatable :: val(:)
integer :: int_err(5)
character(len=20) :: name, ch_err
name='psb_dilu_fctint'
if(psb_get_errstatus() /= 0) return
@ -162,61 +169,23 @@ contains
d(i) = dzero
!
! Here we take a fast shortcut if possible, otherwise
! use spgtblk, slower but able (in principle) to handle
! anything.
!
if (a%fida=='CSR') then
do j = a%ia2(i), a%ia2(i+1) - 1
k = a%ia1(j)
! write(0,*)'KKKKK',k
if ((k < i).and.(k >= 1)) then
l1 = l1 + 1
laspk(l1) = a%aspk(j)
lia1(l1) = k
else if (k == i) then
d(i) = a%aspk(j)
else if ((k > i).and.(k <= m)) then
l2 = l2 + 1
uaspk(l2) = a%aspk(j)
uia1(l2) = k
end if
enddo
else
if ((mod(i,nrb) == 1).or.(nrb==1)) then
irb = min(ma-i+1,nrb)
call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
if(info /= 0) then
info=4010
ch_err='psb_sp_getblk'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
ktrw=1
call a%csget(i,i,nz,irow,icol,val,info)
do j=1, nz
k = icol(j)
! write(0,*)'KKKKK',k
if ((k < i).and.(k >= 1)) then
l1 = l1 + 1
laspk(l1) = val(j)
lia1(l1) = k
else if (k == i) then
d(i) = val(j)
else if ((k > i).and.(k <= m)) then
l2 = l2 + 1
uaspk(l2) = val(j)
uia1(l2) = k
end if
do
if (ktrw > trw%infoa(psb_nnz_)) exit
if (trw%ia1(ktrw) > i) exit
k = trw%ia2(ktrw)
if ((k < i).and.(k >= 1)) then
l1 = l1 + 1
laspk(l1) = trw%aspk(ktrw)
lia1(l1) = k
else if (k == i) then
d(i) = trw%aspk(ktrw)
else if ((k > i).and.(k <= m)) then
l2 = l2 + 1
uaspk(l2) = trw%aspk(ktrw)
uia1(l2) = k
end if
ktrw = ktrw + 1
enddo
end if
end do
!!$
lia2(i+1) = l1 + 1

@ -115,7 +115,7 @@ subroutine psb_dprecbld(aa,desc_a,p,info,upd)
call psb_check_def(p%iprcparm(psb_f_type_),'fact',&
& psb_f_ilu_n_,is_legal_ml_fact)
call psb_bjac_bld(a,desc_a,p,upd_,info)
call psb_bjac_bld(aa,desc_a,p,upd_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_bjac_bld')

@ -55,12 +55,12 @@ subroutine psb_dprecinit(p,ptype,info)
p%iprcparm(psb_p_type_) = psb_diag_
p%iprcparm(psb_f_type_) = psb_f_none_
!!$ case ('BJAC')
!!$ p%iprcparm(:) = 0
!!$ p%iprcparm(psb_p_type_) = psb_bjac_
!!$ p%iprcparm(psb_f_type_) = psb_f_ilu_n_
!!$ p%iprcparm(psb_ilu_fill_in_) = 0
!!$
case ('BJAC')
p%iprcparm(:) = 0
p%iprcparm(psb_p_type_) = psb_bjac_
p%iprcparm(psb_f_type_) = psb_f_ilu_n_
p%iprcparm(psb_ilu_fill_in_) = 0
case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2

@ -330,9 +330,12 @@ module psb_prec_mod
end subroutine psb_silu_fct
subroutine psb_dilu_fct(a,l,u,d,info,blck)
use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_
use psbn_d_mat_mod
integer, intent(out) :: info
type(psb_dspmat_type),intent(in) :: a
type(psb_dspmat_type),intent(inout) :: l,u
type(psbn_d_sparse_mat),intent(in) :: a
type(psbn_d_csr_sparse_mat),intent(inout) :: l,u
!!$ type(psb_dspmat_type),intent(in) :: a
!!$ type(psb_dspmat_type),intent(inout) :: l,u
type(psb_dspmat_type),intent(in), optional, target :: blck
real(psb_dpk_), intent(inout) :: d(:)
end subroutine psb_dilu_fct
@ -365,10 +368,11 @@ module psb_prec_mod
character, intent(in) :: upd
end subroutine psb_sbjac_bld
subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
use psbn_d_mat_mod
use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_
use psb_prec_type, only : psb_dprec_type
integer, intent(out) :: info
type(psb_dspmat_type), intent(in), target :: a
type(psbn_d_sparse_mat), intent(in), target :: a
type(psb_dprec_type), intent(inout) :: p
type(psb_desc_type), intent(in) :: desc_a
character, intent(in) :: upd

@ -41,6 +41,7 @@ module psb_prec_type
& psb_dspmat_type, psb_zspmat_type, psb_dpk_, psb_spk_, psb_long_int_k_,&
& psb_desc_type, psb_sizeof, psb_sp_free, psb_cdfree,&
& psb_erractionsave, psb_erractionrestore, psb_error, psb_get_errstatus
use psbn_d_mat_mod, only : psbn_d_sparse_mat
integer, parameter :: psb_min_prec_=0, psb_noprec_=0, psb_diag_=1, &
& psb_bjac_=2, psb_max_prec_=2
@ -75,7 +76,7 @@ module psb_prec_type
end type psb_sprec_type
type psb_dprec_type
type(psb_dspmat_type), allocatable :: av(:)
type(psbn_d_sparse_mat), allocatable :: av(:)
real(psb_dpk_), allocatable :: d(:)
type(psb_desc_type) :: desc_data
integer, allocatable :: iprcparm(:)
@ -402,7 +403,8 @@ contains
if (allocated(p%av)) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
!!$ call psb_sp_free(p%av(i),info)
call p%av(i)%free()
if (info /= 0) then
! Actually, we don't care here about this.
! Just let it go.
@ -600,6 +602,7 @@ contains
function psb_dprec_sizeof(prec) result(val)
use psbn_d_mat_mod
type(psb_dprec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
integer :: i

@ -94,7 +94,7 @@ program ppde
real(psb_dpk_) :: err, eps
! other variables
integer :: info
integer :: info, i
character(len=20) :: name,ch_err
info=0
@ -131,7 +131,6 @@ program ppde
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es12.5)')t2
if (iam == psb_root_) write(*,'(" ")')
!
@ -156,7 +155,6 @@ program ppde
if (iam == psb_root_) write(*,'("Preconditioner time : ",es12.5)')tprec
if (iam == psb_root_) write(*,'(" ")')
!
! iterative method parameters
!
@ -178,8 +176,7 @@ program ppde
t2 = psb_wtime() - t1
call psb_amx(ictxt,t2)
!!$ amatsize = psb_sizeof(a)
amatsize = 0
amatsize = psb_sizeof(a)
descsize = psb_sizeof(desc_a)
precsize = psb_sizeof(prec)
call psb_sum(ictxt,amatsize)

@ -1,11 +1,11 @@
7 Number of entries below this
BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL RGMRES
DIAG Preconditioner NONE DIAG BJAC
CSR Storage format for matrix A: CSR COO JAD
BJAC Preconditioner NONE DIAG BJAC
COO Storage format for matrix A: CSR COO JAD
060 Domain size (acutal system is this**3)
1 Stopping criterion
2 Stopping criterion
0400 MAXIT
-40 ITRACE
-01 ITRACE
20 IRST restart for RGMRES and BiCGSTABL

Loading…
Cancel
Save