fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-16 19:02:09 +00:00
Code Issues Packages Projects Releases Wiki Activity

psblas3:

Browse Source 
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
 test/serial/d_matgen.f03

Added SCAL. 
Defined CSGET/CSCLIP: they work on the same inner kernel
implementation!
This commit is contained in:
Salvatore Filippone
2009-09-11 14:20:06 +00:00
parent 607a0aa949
commit fb28c925dc
7 changed files with 1776 additions and 79 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/newserial/psbn_base_mat_mod.f03
+30 -8
View File
@@ -84,6 +84,7 @@ module psbn_base_mat_mod
procedure, pass(a) :: allocate_mnnz
procedure, pass(a) :: reallocate_nz
procedure, pass(a) :: free
procedure, pass(a) :: trim
generic, public :: allocate => allocate_mnnz
generic, public :: reallocate => reallocate_nz
procedure, pass(a) :: print => sparse_print
@@ -96,7 +97,7 @@ 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
& free, sparse_print, get_fmt, trim
contains
@@ -329,7 +330,7 @@ contains
character(len=20) :: name='base_get_nzeros'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
call psb_get_erraction(err_act)
res = -1
! This is the base version. If we get here
! it means the derived class is incomplete,
@@ -353,7 +354,7 @@ contains
character(len=20) :: name='get_size'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
call psb_get_erraction(err_act)
res = -1
! This is the base version. If we get here
! it means the derived class is incomplete,
@@ -383,7 +384,7 @@ contains
character(len=20) :: name='sparse_print'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
call psb_get_erraction(err_act)
info = 700
! This is the base version. If we get here
! it means the derived class is incomplete,
@@ -412,7 +413,7 @@ contains
character(len=20) :: name='get_neigh'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
call psb_get_erraction(err_act)
info = 700
! This is the base version. If we get here
! it means the derived class is incomplete,
@@ -436,7 +437,7 @@ contains
character(len=20) :: name='allocate_mnz'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
@@ -458,7 +459,7 @@ contains
character(len=20) :: name='reallocate_nz'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
@@ -479,7 +480,7 @@ contains
character(len=20) :: name='free'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
@@ -492,5 +493,26 @@ contains
end subroutine free
subroutine trim(a)
use psb_error_mod
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
Integer :: err_act
character(len=20) :: name='trim'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
call psb_errpush(700,name,a_err=a%get_fmt())
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine trim
end module psbn_base_mat_mod
base/newserial/psbn_d_base_mat_mod.f03
+614 -6
View File
@@ -10,8 +10,16 @@ module psbn_d_base_mat_mod
procedure, pass(a) :: d_base_cssv
procedure, pass(a) :: d_base_cssm
generic, public :: cssm => d_base_cssm, d_base_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) :: csput
procedure, pass(a) :: d_csgetrow
procedure, pass(a) :: d_csgetblk
generic, public :: csget => d_csgetrow, d_csgetblk
procedure, pass(a) :: csclip
procedure, pass(a) :: cp_to_coo
procedure, pass(a) :: cp_from_coo
procedure, pass(a) :: cp_to_fmt
@@ -21,9 +29,12 @@ module psbn_d_base_mat_mod
procedure, pass(a) :: mv_to_fmt
procedure, pass(a) :: mv_from_fmt
end type psbn_d_base_sparse_mat
private :: d_base_csmv, d_base_csmm, d_base_cssv, d_base_cssm,&
& csnmi, csput, cp_to_coo, cp_from_coo, cp_to_fmt, cp_from_fmt, &
& mv_to_coo, mv_from_coo, mv_to_fmt, mv_from_fmt
& 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
type, extends(psbn_d_base_sparse_mat) :: psbn_d_coo_sparse_mat
@@ -41,8 +52,11 @@ module psbn_d_base_mat_mod
procedure, pass(a) :: d_base_csmv => d_coo_csmv
procedure, pass(a) :: d_base_cssm => d_coo_cssm
procedure, pass(a) :: d_base_cssv => d_coo_cssv
procedure, pass(a) :: d_scals => d_coo_scals
procedure, pass(a) :: d_scal => d_coo_scal
procedure, pass(a) :: csnmi => d_coo_csnmi
procedure, pass(a) :: csput => d_coo_csput
procedure, pass(a) :: get_diag => d_coo_get_diag
procedure, pass(a) :: reallocate_nz => d_coo_reallocate_nz
procedure, pass(a) :: allocate_mnnz => d_coo_allocate_mnnz
procedure, pass(a) :: cp_to_coo => d_cp_coo_to_coo
@@ -55,16 +69,23 @@ module psbn_d_base_mat_mod
procedure, pass(a) :: mv_from_fmt => d_mv_coo_from_fmt
procedure, pass(a) :: fix => d_fix_coo
procedure, pass(a) :: free => d_coo_free
procedure, pass(a) :: trim => d_coo_trim
procedure, pass(a) :: d_csgetrow => d_coo_csgetrow
procedure, pass(a) :: d_csgetblk => d_coo_csgetblk
procedure, pass(a) :: csclip => d_coo_csclip
procedure, pass(a) :: print => d_coo_print
procedure, pass(a) :: get_fmt => d_coo_get_fmt
end type psbn_d_coo_sparse_mat
private :: d_coo_get_nzeros, d_coo_set_nzeros, &
private :: d_coo_get_nzeros, d_coo_set_nzeros, d_coo_get_diag, &
& d_coo_csmm, d_coo_csmv, d_coo_cssm, d_coo_cssv, d_coo_csnmi, &
& d_coo_csput, d_coo_reallocate_nz, d_coo_allocate_mnnz, &
& 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_cp_coo_to_fmt, d_cp_coo_from_fmt, &
& d_coo_scals, d_coo_scal, d_coo_csgetrow, d_coo_csgetblk, &
& d_coo_csclip
interface
@@ -175,12 +196,33 @@ module psbn_d_base_mat_mod
import psbn_d_coo_sparse_mat
class(psbn_d_coo_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: val(:)
integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax
integer, intent(in) :: nz,ia(:), ja(:),&
& imin,imax,jmin,jmax
integer, intent(out) :: info
integer, intent(in), optional :: gtl(:)
end subroutine d_coo_csput_impl
end interface
interface
subroutine d_coo_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
use psb_const_mod
import psbn_d_coo_sparse_mat
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
integer, intent(in) :: imin,imax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale
end subroutine d_coo_csgetrow_impl
end interface
interface d_coo_cssm_impl
subroutine d_coo_cssv_impl(alpha,a,x,beta,y,info,trans)
use psb_const_mod
@@ -490,6 +532,108 @@ contains
end subroutine csput
subroutine d_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
integer, intent(in) :: imin,imax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale
Integer :: err_act
character(len=20) :: name='csget'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
info = 700
call psb_errpush(info,name,a_err=a%get_fmt())
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_csgetrow
subroutine d_csgetblk(imin,imax,a,b,info,&
& jmin,jmax,iren,append,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
class(psbn_d_coo_sparse_mat), intent(inout) :: b
integer, intent(in) :: imin,imax
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax
logical, intent(in), optional :: rscale,cscale
Integer :: err_act
character(len=20) :: name='csget'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
info = 700
call psb_errpush(info,name,a_err=a%get_fmt())
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_csgetblk
subroutine csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer,intent(out) :: info
integer, intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
Integer :: err_act
character(len=20) :: name='csclip'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
info = 700
call psb_errpush(info,name,a_err=a%get_fmt())
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine csclip
!====================================
!
!
@@ -610,9 +754,61 @@ contains
end if
return
end subroutine d_base_cssv
subroutine d_scals(d,a,info)
use psb_error_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: d
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='d_scals'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
info = 700
call psb_errpush(info,name,a_err=a%get_fmt())
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_scals
subroutine d_scal(d,a,info)
use psb_error_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: d(:)
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='d_scal'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
info = 700
call psb_errpush(info,name,a_err=a%get_fmt())
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_scal
function csnmi(a) result(res)
use psb_error_mod
use psb_const_mod
@@ -640,6 +836,33 @@ contains
end function csnmi
subroutine get_diag(a,d,info)
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(out) :: d(:)
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='get_diag'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
! This is the base version. If we get here
! it means the derived class is incomplete,
! so we throw an error.
info = 700
call psb_errpush(info,name,a_err=a%get_fmt())
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine get_diag
@@ -1125,6 +1348,242 @@ contains
end subroutine d_coo_csput
subroutine d_coo_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
integer, intent(in) :: imin,imax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale
Integer :: err_act
character(len=20) :: name='csget'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
call d_coo_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
if (info /= 0) goto 9999
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_coo_csgetrow
subroutine d_coo_csgetblk(imin,imax,a,b,info,&
& jmin,jmax,iren,append,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
class(psbn_d_coo_sparse_mat), intent(inout) :: b
integer, intent(in) :: imin,imax
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax
logical, intent(in), optional :: rscale,cscale
Integer :: err_act, nzin, nzout
character(len=20) :: name='csget'
logical :: append_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
if (present(append)) then
append_ = append
else
append_ = .false.
endif
if (append_) then
nzin = a%get_nzeros()
else
nzin = 0
endif
call a%csget(imin,imax,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin, jmax=jmax, iren=iren, append=append_, &
& nzin=nzin, rscale=rscale, cscale=cscale)
if (info /= 0) goto 9999
call b%set_nzeros(nzin+nzout)
call b%fix(info)
if (info /= 0) goto 9999
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_coo_csgetblk
subroutine d_coo_csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer,intent(out) :: info
integer, intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
Integer :: err_act, nzin, nzout, imin_, imax_, jmin_, jmax_, mb,nb
character(len=20) :: name='csget'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
nzin = 0
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = a%get_nrows() ! Should this be imax_ ??
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = a%get_ncols() ! Should this be jmax_ ??
endif
call b%allocate(mb,nb)
call a%csget(imin_,imax_,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=jmax_, append=.false., &
& nzin=nzin, rscale=rscale_, cscale=cscale_)
if (info /= 0) goto 9999
call b%set_nzeros(nzin+nzout)
call b%fix(info)
if (info /= 0) goto 9999
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_coo_csclip
!!$
!!$ subroutine d_coo_csget(irw,a,nz,ia,ja,val,info,iren,lrw,append,nzin)
!!$ ! Output is always in COO format
!!$ use psb_error_mod
!!$ use psb_const_mod
!!$ implicit none
!!$
!!$ class(psbn_d_coo_sparse_mat), intent(inout) :: a
!!$ integer, intent(in) :: irw
!!$ integer, intent(out) :: nz
!!$ integer, allocatable, intent(inout) :: ia(:), ja(:)
!!$ real(psb_dpk_), allocatable, intent(inout) :: val(:)
!!$ integer,intent(out) :: info
!!$ logical, intent(in), optional :: append
!!$ integer, intent(in), optional :: iren(:)
!!$ integer, intent(in), optional :: lrw, nzin
!!$ Integer :: err_act
!!$ character(len=20) :: name='csget'
!!$ logical, parameter :: debug=.false.
!!$
!!$ call psb_erractionsave(err_act)
!!$ info = 0
!!$
!!$ call d_coo_csget_impl(irw,a,nz,ia,ja,val,info,iren,lrw,append,nzin)
!!$ if (info /= 0) goto 9999
!!$
!!$ 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_coo_csget
!!$
subroutine d_coo_free(a)
implicit none
@@ -1141,6 +1600,38 @@ contains
end subroutine d_coo_free
subroutine d_coo_trim(a)
use psb_realloc_mod
use psb_error_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
Integer :: err_act, info, nz
character(len=20) :: name='trim'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
nz = a%get_nzeros()
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)
if (info /= 0) goto 9999
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_coo_trim
subroutine d_coo_allocate_mnnz(m,n,a,nz)
use psb_error_mod
use psb_realloc_mod
@@ -1275,6 +1766,7 @@ contains
!====================================
!
!
@@ -1494,6 +1986,122 @@ contains
end function d_coo_csnmi
subroutine d_coo_get_diag(a,d,info)
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(out) :: d(:)
integer, intent(out) :: info
Integer :: err_act,mnm, i, j
character(len=20) :: name='get_diag'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
mnm = min(a%get_nrows(),a%get_ncols())
if (size(d) < mnm) then
info=35
call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/))
goto 9999
end if
d(:) = dzero
do i=1,a%get_nzeros()
j=a%ia(i)
if ((j==a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then
d(j) = a%val(i)
endif
enddo
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_coo_get_diag
subroutine d_coo_scal(d,a,info)
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: d(:)
integer, intent(out) :: info
Integer :: err_act,mnm, i, j, m
character(len=20) :: name='scal'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
m = a%get_nrows()
if (size(d) < m) then
info=35
call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/))
goto 9999
end if
do i=1,a%get_nzeros()
j = a%ia(i)
a%val(i) = a%val(i) * d(j)
enddo
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_coo_scal
subroutine d_coo_scals(d,a,info)
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: d
integer, intent(out) :: info
Integer :: err_act,mnm, i, j, m
character(len=20) :: name='scal'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d
enddo
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_coo_scals
end module psbn_d_base_mat_mod
base/newserial/psbn_d_coo_impl.f03
+295 -11
View File
@@ -881,26 +881,309 @@ function d_coo_csnmi_impl(a) result(res)
acc = acc + abs(a%val(k))
end do
res = max(res,acc)
i = j
end do
end function d_coo_csnmi_impl
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!
! Data management
!
!
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!====================================
!
!
!
! Data management
!
!
!
!
!
!====================================
subroutine d_coo_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_error_mod
use psbn_d_base_mat_mod, psb_protect_name => d_coo_csgetrow_impl
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
integer, intent(in) :: imin,imax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale
logical :: append_, rscale_, cscale_
integer :: nzin_, jmin_, jmax_, err_act, i
character(len=20) :: name='csget'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
endif
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
endif
if ((imax<imin).or.(jmax_<jmin_)) return
if (present(append)) then
append_=append
else
append_=.false.
endif
if ((append_).and.(present(nzin))) then
nzin_ = nzin
else
nzin_ = 0
endif
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .false.
endif
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .false.
endif
if ((rscale_.or.cscale_).and.(present(iren))) then
info = 583
call psb_errpush(info,name,a_err='iren (rscale.or.cscale)')
goto 9999
end if
call coo_getrow(imin,imax,jmin_,jmax_,a,nz,ia,ja,val,nzin_,append_,info,&
& iren)
if (rscale_) then
do i=nzin_+1, nzin_+nz
ia(i) = ia(i) - imin + 1
end do
end if
if (cscale_) then
do i=nzin_+1, nzin_+nz
ja(i) = ja(i) - jmin_ + 1
end do
end if
if (info /= 0) goto 9999
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 coo_getrow(imin,imax,jmin,jmax,a,nz,ia,ja,val,nzin,append,info,&
& iren)
use psb_const_mod
use psb_error_mod
use psb_realloc_mod
use psb_sort_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
integer :: imin,imax,jmin,jmax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer, intent(in) :: nzin
logical, intent(in) :: append
integer :: info
integer, optional :: iren(:)
integer :: nzin_, nza, idx,ip,jp,i,k, nzt, irw, lrw
integer :: debug_level, debug_unit
character(len=20) :: name='coo_getrow'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nza = a%get_nzeros()
irw = imin
lrw = imax
if (irw<0) then
write(debug_unit,*) ' spgtrow Error : idx no good ',irw
info = 2
return
end if
if (append) then
nzin_ = nzin
else
nzin_ = 0
endif
if (a%is_sorted()) then
! In this case we can do a binary search.
if (debug_level >= psb_debug_serial_)&
& write(debug_unit,*) trim(name), ': srtdcoo '
do
ip = psb_ibsrch(irw,nza,a%ia)
if (ip /= -1) exit
irw = irw + 1
if (irw > imax) then
write(debug_unit,*) trim(name),&
& 'Warning : did not find any rows. Is this an error? ',&
& irw,lrw,imin
exit
end if
end do
if (ip /= -1) then
! expand [ip,jp] to contain all row entries.
do
if (ip < 2) exit
if (a%ia(ip-1) == irw) then
ip = ip -1
else
exit
end if
end do
end if
do
jp = psb_ibsrch(lrw,nza,a%ia)
if (jp /= -1) exit
lrw = lrw - 1
if (irw > lrw) then
write(debug_unit,*) trim(name),&
& 'Warning : did not find any rows. Is this an error?'
exit
end if
end do
if (jp /= -1) then
! expand [ip,jp] to contain all row entries.
do
if (jp == nza) exit
if (a%ia(jp+1) == lrw) then
jp = jp + 1
else
exit
end if
end do
end if
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),': ip jp',ip,jp,nza
if ((ip /= -1) .and.(jp /= -1)) then
! Now do the copy.
nzt = jp - ip +1
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
do i=ip,jp
if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then
nzin_ = nzin_ + 1
nz = nz + 1
val(nzin_) = a%val(i)
ia(nzin_) = iren(a%ia(i))
ja(nzin_) = iren(a%ja(i))
end if
enddo
else
do i=ip,jp
if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then
nzin_ = nzin_ + 1
nz = nz + 1
val(nzin_) = a%val(i)
ia(nzin_) = a%ia(i)
ja(nzin_) = a%ja(i)
end if
enddo
end if
else
nz = 0
end if
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),': unsorted '
nzt = (nza*(lrw-irw+1))/max(a%get_nrows(),1)
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
k = 0
do i=1, a%get_nzeros()
if ((a%ia(i)>=irw).and.(a%ia(i)<=lrw).and.&
& (jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then
k = k + 1
if (k > nzt) then
nzt = k
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
end if
val(nzin_+k) = a%val(i)
ia(nzin_+k) = iren(a%ia(i))
ja(nzin_+k) = iren(a%ja(i))
endif
enddo
else
k = 0
do i=1,a%get_nzeros()
if ((a%ia(i)>=irw).and.(a%ia(i)<=lrw).and.&
& (jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then
k = k + 1
if (k > nzt) then
nzt = k
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
end if
val(nzin_+k) = a%val(i)
ia(nzin_+k) = (a%ia(i))
ja(nzin_+k) = (a%ja(i))
endif
enddo
nzin_=nzin_+k
end if
nz = k
end if
end subroutine coo_getrow
end subroutine d_coo_csgetrow_impl
subroutine d_coo_csput_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
use psb_error_mod
use psb_realloc_mod
@@ -965,6 +1248,7 @@ subroutine d_coo_csput_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
call psb_inner_ins(nz,ia,ja,val,nza,a%ia,a%ja,a%val,isza,&
& imin,imax,jmin,jmax,info,gtl)
call a%set_nzeros(nza)
call a%set_sorted(.false.)
else if (a%is_upd()) then
base/newserial/psbn_d_csr_impl.f03
+201 -23
View File
@@ -1,16 +1,16 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!
! Computational routines
!
!
!
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!=====================================
!
!
!
! Computational routines
!
!
!
!
!
!
!=====================================
subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
use psb_error_mod
@@ -980,17 +980,195 @@ function d_csr_csnmi_impl(a) result(res)
end function d_csr_csnmi_impl
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!
! Data management
!
!
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!=====================================
!
!
!
! Data management
!
!
!
!
!
!=====================================
subroutine d_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_error_mod
use psbn_d_base_mat_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_csr_csgetrow_impl
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
integer, intent(in) :: imin,imax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale
logical :: append_, rscale_, cscale_
integer :: nzin_, jmin_, jmax_, err_act, i
character(len=20) :: name='csget'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
endif
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
endif
if ((imax<imin).or.(jmax_<jmin_)) return
if (present(append)) then
append_=append
else
append_=.false.
endif
if ((append_).and.(present(nzin))) then
nzin_ = nzin
else
nzin_ = 0
endif
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .false.
endif
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .false.
endif
if ((rscale_.or.cscale_).and.(present(iren))) then
info = 583
call psb_errpush(info,name,a_err='iren (rscale.or.cscale)')
goto 9999
end if
call csr_getrow(imin,imax,jmin_,jmax_,a,nz,ia,ja,val,nzin_,append_,info,&
& iren)
if (rscale_) then
do i=nzin_+1, nzin_+nz
ia(i) = ia(i) - imin + 1
end do
end if
if (cscale_) then
do i=nzin_+1, nzin_+nz
ja(i) = ja(i) - jmin_ + 1
end do
end if
if (info /= 0) goto 9999
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 csr_getrow(imin,imax,jmin,jmax,a,nz,ia,ja,val,nzin,append,info,&
& iren)
use psb_const_mod
use psb_error_mod
use psb_realloc_mod
use psb_sort_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
integer :: imin,imax,jmin,jmax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer, intent(in) :: nzin
logical, intent(in) :: append
integer :: info
integer, optional :: iren(:)
integer :: nzin_, nza, idx,i,j,k, nzt, irw, lrw
integer :: debug_level, debug_unit
character(len=20) :: name='coo_getrow'
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nza = a%get_nzeros()
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
if (append) then
nzin_ = nzin
else
nzin_ = 0
endif
nzt = a%irp(lrw)-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
do i=irw, lrw
do j=a%irp(i), a%irp(i+1) - 1
if ((jmin <= a%ja(j)).and.(a%ja(j)<=jmax)) then
nzin_ = nzin_ + 1
nz = nz + 1
val(nzin_) = a%val(j)
ia(nzin_) = iren(i)
ja(nzin_) = iren(a%ja(j))
end if
enddo
end do
else
do i=irw, lrw
do j=a%irp(i), a%irp(i+1) - 1
if ((jmin <= a%ja(j)).and.(a%ja(j)<=jmax)) then
nzin_ = nzin_ + 1
nz = nz + 1
val(nzin_) = a%val(j)
ia(nzin_) = (i)
ja(nzin_) = (a%ja(j))
end if
enddo
end do
end if
end subroutine csr_getrow
end subroutine d_csr_csgetrow_impl
base/newserial/psbn_d_csr_mat_mod.f03
+405 -5
View File
@@ -10,10 +10,13 @@ module psbn_d_csr_mat_mod
contains
procedure, pass(a) :: get_nzeros => d_csr_get_nzeros
procedure, pass(a) :: get_fmt => d_csr_get_fmt
procedure, pass(a) :: get_diag => d_csr_get_diag
procedure, pass(a) :: d_base_csmm => d_csr_csmm
procedure, pass(a) :: d_base_csmv => d_csr_csmv
procedure, pass(a) :: d_base_cssm => d_csr_cssm
procedure, pass(a) :: d_base_cssv => d_csr_cssv
procedure, pass(a) :: d_scals => d_csr_scals
procedure, pass(a) :: d_scal => d_csr_scal
procedure, pass(a) :: csnmi => d_csr_csnmi
procedure, pass(a) :: reallocate_nz => d_csr_reallocate_nz
procedure, pass(a) :: csput => d_csr_csput
@@ -26,16 +29,23 @@ module psbn_d_csr_mat_mod
procedure, pass(a) :: mv_from_coo => d_mv_csr_from_coo
procedure, pass(a) :: mv_to_fmt => d_mv_csr_to_fmt
procedure, pass(a) :: mv_from_fmt => d_mv_csr_from_fmt
procedure, pass(a) :: d_csgetrow => d_csr_csgetrow
procedure, pass(a) :: d_csgetblk => d_csr_csgetblk
procedure, pass(a) :: csclip => d_csr_csclip
procedure, pass(a) :: get_size => d_csr_get_size
procedure, pass(a) :: free => d_csr_free
procedure, pass(a) :: trim => d_csr_trim
procedure, pass(a) :: print => d_csr_print
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, &
& d_csr_free, d_csr_print, d_csr_get_fmt, d_csr_csnmi, &
& d_csr_free, d_csr_print, d_csr_get_fmt, d_csr_csnmi, get_diag, &
& d_cp_csr_to_coo, d_cp_csr_from_coo, &
& 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_mv_csr_to_fmt, d_mv_csr_from_fmt, &
& d_csr_scals, d_csr_scal, d_csr_trim, d_csr_csgetrow, d_csr_csgetblk, &
& d_csr_csclip, d_csr_get_size
interface
@@ -140,6 +150,26 @@ module psbn_d_csr_mat_mod
end subroutine d_csr_csput_impl
end interface
interface
subroutine d_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
use psb_const_mod
import psbn_d_csr_sparse_mat
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
integer, intent(in) :: imin,imax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale
end subroutine d_csr_csgetrow_impl
end interface
interface d_csr_cssm_impl
subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
use psb_const_mod
@@ -221,8 +251,31 @@ contains
res = a%irp(a%m+1)-1
end function d_csr_get_nzeros
function d_csr_get_size(a) result(res)
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
integer :: res
!=====================================
res = -1
if (allocated(a%ja)) then
if (res >= 0) then
res = min(res,size(a%ja))
else
res = size(a%ja)
end if
end if
if (allocated(a%val)) then
if (res >= 0) then
res = min(res,size(a%val))
else
res = size(a%val)
end if
end if
end function d_csr_get_size
!=====================================
!
!
!
@@ -232,7 +285,7 @@ contains
!
!
!
!=====================================
!=====================================
subroutine d_csr_reallocate_nz(nz,a)
@@ -332,6 +385,198 @@ contains
return
end subroutine d_csr_csput
subroutine d_csr_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
integer, intent(in) :: imin,imax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale
Integer :: err_act
character(len=20) :: name='csget'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
call d_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
if (info /= 0) goto 9999
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_csr_csgetrow
subroutine d_csr_csgetblk(imin,imax,a,b,info,&
& jmin,jmax,iren,append,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
class(psbn_d_coo_sparse_mat), intent(inout) :: b
integer, intent(in) :: imin,imax
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax
logical, intent(in), optional :: rscale,cscale
Integer :: err_act, nzin, nzout
character(len=20) :: name='csget'
logical :: append_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
if (present(append)) then
append_ = append
else
append_ = .false.
endif
if (append_) then
nzin = a%get_nzeros()
else
nzin = 0
endif
call a%csget(imin,imax,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin, jmax=jmax, iren=iren, append=append_, &
& nzin=nzin, rscale=rscale, cscale=cscale)
if (info /= 0) goto 9999
call b%set_nzeros(nzin+nzout)
call b%fix(info)
if (info /= 0) goto 9999
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_csr_csgetblk
subroutine d_csr_csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer,intent(out) :: info
integer, intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
Integer :: err_act, nzin, nzout, imin_, imax_, jmin_, jmax_, mb,nb
character(len=20) :: name='csget'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
nzin = 0
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = a%get_nrows() ! Should this be imax_ ??
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = a%get_ncols() ! Should this be jmax_ ??
endif
call b%allocate(mb,nb)
call a%csget(imin_,imax_,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=jmax_, append=.false., &
& nzin=nzin, rscale=rscale_, cscale=cscale_)
if (info /= 0) goto 9999
call b%set_nzeros(nzin+nzout)
call b%fix(info)
if (info /= 0) goto 9999
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_csr_csclip
subroutine d_csr_free(a)
implicit none
@@ -349,6 +594,38 @@ contains
end subroutine d_csr_free
subroutine d_csr_trim(a)
use psb_realloc_mod
use psb_error_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
Integer :: err_act, info, nz, m
character(len=20) :: name='trim'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
m = a%get_nrows()
nz = a%get_nzeros()
if (info == 0) call psb_realloc(m+1,a%irp,info)
if (info == 0) call psb_realloc(nz,a%ja,info)
if (info == 0) call psb_realloc(nz,a%val,info)
if (info /= 0) goto 9999
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_csr_trim
subroutine d_cp_csr_to_coo(a,b,info)
use psb_error_mod
@@ -965,13 +1242,136 @@ contains
character(len=20) :: name='csnmi'
logical, parameter :: debug=.false.
res = d_csr_csnmi_impl(a)
return
end function d_csr_csnmi
subroutine d_csr_get_diag(a,d,info)
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(out) :: d(:)
integer, intent(out) :: info
Integer :: err_act, mnm, i, j, k
character(len=20) :: name='get_diag'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
mnm = min(a%get_nrows(),a%get_ncols())
if (size(d) < mnm) then
info=35
call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/))
goto 9999
end if
do i=1, mnm
do k=a%irp(i),a%irp(i+1)-1
j=a%ja(k)
if ((j==i) .and.(j <= mnm )) then
d(i) = a%val(k)
endif
enddo
end do
do i=mnm+1,size(d)
d(i) = dzero
end do
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_csr_get_diag
subroutine d_csr_scal(d,a,info)
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: d(:)
integer, intent(out) :: info
Integer :: err_act,mnm, i, j, m
character(len=20) :: name='scal'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
m = a%get_nrows()
if (size(d) < m) then
info=35
call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/))
goto 9999
end if
do i=1, m
do j = a%irp(i), a%irp(i+1) -1
a%val(j) = a%val(j) * d(i)
end do
enddo
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_csr_scal
subroutine d_csr_scals(d,a,info)
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: d
integer, intent(out) :: info
Integer :: err_act,mnm, i, j, m
character(len=20) :: name='scal'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d
enddo
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_csr_scals
end module psbn_d_csr_mat_mod
base/newserial/psbn_mat_mod.f03
+207 -23
View File
@@ -44,8 +44,12 @@ module psbn_d_mat_mod
! Memory/data management
procedure, pass(a) :: csall
procedure, pass(a) :: free
procedure, pass(a) :: trim
procedure, pass(a) :: csput
procedure, pass(a) :: csget
procedure, pass(a) :: d_csgetrow
procedure, pass(a) :: d_csgetblk
generic, public :: csget => d_csgetrow, d_csgetblk
procedure, pass(a) :: csclip
procedure, pass(a) :: reall => reallocate_nz
procedure, pass(a) :: get_neigh
procedure, pass(a) :: d_cscnv
@@ -54,6 +58,7 @@ module psbn_d_mat_mod
procedure, pass(a) :: print => sparse_print
! Computational routines
procedure, pass(a) :: get_diag
procedure, pass(a) :: csnmi
procedure, pass(a) :: d_csmv
procedure, pass(a) :: d_csmm
@@ -68,11 +73,12 @@ module psbn_d_mat_mod
private :: get_nrows, get_ncols, 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, csall, csput, csget, d_cscnv, d_cscnv_ip, &
& reallocate_nz, free, d_csmv, d_csmm, d_cssv, d_cssm, sparse_print, &
& 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, &
& 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
& set_unit, csnmi, get_diag
contains
@@ -904,6 +910,7 @@ contains
character(len=20) :: name='reallocate_nz'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
if (.not.allocated(a%a)) then
info = 1121
call psb_errpush(info,name)
@@ -933,6 +940,7 @@ contains
character(len=20) :: name='free'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
if (.not.allocated(a%a)) then
info = 1121
call psb_errpush(info,name)
@@ -944,7 +952,6 @@ contains
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
@@ -954,6 +961,36 @@ contains
end subroutine free
subroutine trim(a)
use psb_error_mod
implicit none
class(psbn_d_sparse_mat), intent(inout) :: a
Integer :: err_act, info
character(len=20) :: name='trim'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
if (.not.allocated(a%a)) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
call a%a%trim()
return
9999 continue
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine trim
subroutine csput(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
use psbn_d_base_mat_mod
use psb_error_mod
@@ -993,19 +1030,27 @@ contains
end subroutine csput
subroutine csget(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
use psbn_d_base_mat_mod
subroutine d_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
! Output is always in COO format
use psb_error_mod
implicit none
class(psbn_d_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: nz, ia(:), ja(:)
integer, intent(in) :: imin,imax,jmin,jmax
integer, intent(out) :: info
integer, intent(in), optional :: gtl(:)
use psb_const_mod
use psbn_d_base_mat_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
integer, intent(in) :: imin,imax
integer, intent(out) :: nz
integer, allocatable, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), allocatable, intent(inout) :: val(:)
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale
Integer :: err_act
character(len=20) :: name='csput'
character(len=20) :: name='csget'
logical, parameter :: debug=.false.
info = 0
@@ -1016,14 +1061,10 @@ contains
goto 9999
endif
info = 700
call psb_errpush(info,name,a_err='CSGET')
goto 9999
!!$
!!$ call a%a%csget(nz,val,ia,ja,imin,imax,jmin,jmax,info,gtl)
!!$ if (info /= 0) goto 9999
call a%a%csget(imin,imax,nz,ia,ja,val,info,&
& jmin,jmax,iren,append,nzin,rscale,cscale)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
@@ -1036,7 +1077,109 @@ contains
return
end if
end subroutine csget
end subroutine d_csgetrow
subroutine d_csgetblk(imin,imax,a,b,info,&
& jmin,jmax,iren,append,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psbn_d_base_mat_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
class(psbn_d_sparse_mat), intent(out) :: b
integer, intent(in) :: imin,imax
integer,intent(out) :: info
logical, intent(in), optional :: append
integer, intent(in), optional :: iren(:)
integer, intent(in), optional :: jmin,jmax
logical, intent(in), optional :: rscale,cscale
Integer :: err_act
character(len=20) :: name='csget'
logical, parameter :: debug=.false.
type(psbn_d_coo_sparse_mat), allocatable :: acoo
info = 0
call psb_erractionsave(err_act)
if (a%is_null()) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == 0) call a%a%csget(imin,imax,acoo,info,&
& jmin,jmax,iren,append,rscale,cscale)
if (info == 0) call move_alloc(acoo,b%a)
if (info /= 0) goto 9999
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
end subroutine d_csgetblk
subroutine csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psbn_d_base_mat_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
class(psbn_d_sparse_mat), intent(out) :: b
integer,intent(out) :: info
integer, intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
Integer :: err_act
character(len=20) :: name='csclip'
logical, parameter :: debug=.false.
type(psbn_d_coo_sparse_mat), allocatable :: acoo
info = 0
call psb_erractionsave(err_act)
if (a%is_null()) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == 0) call a%a%csclip(acoo,info,&
& imin,imax,jmin,jmax,rscale,cscale)
if (info == 0) call move_alloc(acoo,b%a)
if (info /= 0) goto 9999
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
end subroutine csclip
subroutine d_cscnv(a,b,info,type,mold,upd,dupl)
@@ -1399,6 +1542,7 @@ contains
character(len=20) :: name='csnmi'
logical, parameter :: debug=.false.
call psb_get_erraction(err_act)
if (.not.allocated(a%a)) then
info = 1121
call psb_errpush(info,name)
@@ -1420,5 +1564,45 @@ contains
end function csnmi
subroutine get_diag(a,d,info)
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(out) :: d(:)
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='csnmi'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
if (.not.allocated(a%a)) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
call a%a%get_diag(d,info)
if (info /= 0) goto 9999
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 get_diag
end module psbn_d_mat_mod
test/serial/d_matgen.f03
+24 -3
View File
@@ -1,5 +1,5 @@
!
program d_coo_matgen
program d_matgen
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod
@@ -154,7 +154,7 @@ contains
integer :: np, iam, nr, nt,nz,isz
integer :: element
integer, allocatable :: irow(:),icol(:),myidx(:)
real(psb_dpk_), allocatable :: val(:)
real(psb_dpk_), allocatable :: val(:), diag(:)
type(psbn_d_sparse_mat) :: a_n
type(psbn_d_coo_sparse_mat) :: acoo
type(psbn_d_csr_sparse_mat) :: acsr
@@ -382,6 +382,26 @@ contains
call a_n%print(20)
anorm = a_n%csnmi()
write(0,*) 'Nrm infinity ',anorm
call a_n%csget(2,3,element,irow,icol,val,info)
write(0,*) 'From csget ',element,info
if (info == 0) then
do i=1,element
write(0,*) irow(i),icol(i),val(i)
end do
end if
isz = a_n%get_size()
write(0,*) 'Size 1: ',isz
call a_n%trim()
isz = a_n%get_size()
write(0,*) 'Size 2: ',isz
allocate(diag(nlr),stat=info)
if (info == 0) then
call a_n%get_diag(diag,info)
end if
!!$
t1 = psb_wtime()
call a_n%cscnv(info,mold=acxx)
@@ -396,6 +416,7 @@ contains
call a_n%print(21)
anorm = a_n%csnmi()
write(0,*) 'Nrm infinity ',anorm
!!$
if(iam == psb_root_) then
@@ -419,7 +440,7 @@ contains
end if
return
end subroutine create_matrix
end program d_coo_matgen
end program d_matgen
!
! functions parametrizing the differential equation
!