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

psblas3:

Browse Source 
base/newserial/Makefile
 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_impl.f03
 base/newserial/psbn_mat_mod.f03
 base/serial/csr/dcsrsm.f
 test/pargen/runs/ppde.inp
 test/serial
 test/serial/Makefile
 test/serial/d_coo_matgen.f03

Added serial test directory.
First tests of conversion COO-CSR
psblas3-type-indexed
Salvatore Filippone 15 years ago
parent 9bee709ca3
commit 95aeca09c3
12 changed files with 2539 additions and 896 deletions
Show Stats Download Patch File Download Diff File

2
base/newserial/Makefile
Unescape Escape View File

@ -19,7 +19,9 @@ lib: $(MODULES) $(OBJS) $(LIBMOD)
psbn_mat_mod.o: psbn_base_mat_mod.o
psbn_coo_mat.o psbn_csr_mat.o: psbn_d_base_mat_mod.o
psbn_d_coo_impl.o: psbn_d_base_mat_mod.o
psbn_d_csr_impl.o: psbn_d_csr_mat_mod.o
psbn_mat_impl.o: psbn_mat_mod.o
clean:

68
base/newserial/psbn_base_mat_mod.f03
Unescape Escape View File

@ -70,18 +70,22 @@ module psbn_base_mat_mod
generic, public :: allocate => allocate_mn, allocate_mnnz
generic, public :: reallocate => reallocate_nz
procedure, pass(a) :: print => sparse_print
end type psbn_base_sparse_mat
private :: 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, 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, allocate_mn, allocate_mnnz, reallocate_nz, &
& free
& free, sparse_print
contains
function get_dupl(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
integer :: res
res = a%duplicate
@ -89,18 +93,21 @@ contains
function get_state(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
integer :: res
res = a%state
end function get_state
function get_nrows(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
integer :: res
res = a%m
end function get_nrows
function get_ncols(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
integer :: res
res = a%n
@ -108,12 +115,14 @@ contains
subroutine set_nrows(m,a)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
integer, intent(in) :: m
a%m = m
end subroutine set_nrows
subroutine set_ncols(n,a)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
integer, intent(in) :: n
a%n = n
@ -121,6 +130,7 @@ contains
subroutine set_state(n,a)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
integer, intent(in) :: n
a%state = n
@ -128,36 +138,42 @@ contains
subroutine set_dupl(n,a)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
integer, intent(in) :: n
a%duplicate = n
end subroutine set_dupl
subroutine set_null(a)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
a%state = psbn_spmat_null_
end subroutine set_null
subroutine set_bld(a)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
a%state = psbn_spmat_bld_
end subroutine set_bld
subroutine set_upd(a)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
a%state = psbn_spmat_upd_
end subroutine set_upd
subroutine set_asb(a)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
a%state = psbn_spmat_asb_
end subroutine set_asb
subroutine set_sorted(a,val)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
logical, intent(in), optional :: val
@ -169,6 +185,7 @@ contains
end subroutine set_sorted
subroutine set_triangle(a,val)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
logical, intent(in), optional :: val
@ -180,6 +197,7 @@ contains
end subroutine set_triangle
subroutine set_unit(a,val)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
logical, intent(in), optional :: val
@ -191,6 +209,7 @@ contains
end subroutine set_unit
subroutine set_lower(a,val)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
logical, intent(in), optional :: val
@ -202,6 +221,7 @@ contains
end subroutine set_lower
subroutine set_upper(a,val)
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
logical, intent(in), optional :: val
@ -213,54 +233,63 @@ contains
end subroutine set_upper
function is_triangle(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
logical :: res
res = a%triangle
end function is_triangle
function is_unit(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
logical :: res
res = a%unitd
end function is_unit
function is_upper(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
logical :: res
res = a%upper
end function is_upper
function is_lower(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
logical :: res
res = .not.a%upper
end function is_lower
function is_null(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
logical :: res
res = (a%state == psbn_spmat_null_)
end function is_null
function is_bld(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
logical :: res
res = (a%state == psbn_spmat_bld_)
end function is_bld
function is_upd(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
logical :: res
res = (a%state == psbn_spmat_upd_)
end function is_upd
function is_asb(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
logical :: res
res = (a%state == psbn_spmat_asb_)
end function is_asb
function is_sorted(a) result(res)
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
logical :: res
res = a%sorted
@ -269,6 +298,7 @@ contains
function get_nzeros(a) result(res)
use psb_error_mod
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
integer :: res
@ -292,6 +322,7 @@ contains
function get_size(a) result(res)
use psb_error_mod
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
integer :: res
@ -313,9 +344,40 @@ contains
end function get_size
subroutine sparse_print(iout,a,iv,eirs,eics,head,ivr,ivc)
use psb_spmat_type
use psb_string_mod
implicit none
integer, intent(in) :: iout
class(psbn_base_sparse_mat), intent(in) :: a
integer, intent(in), optional :: iv(:)
integer, intent(in), optional :: eirs,eics
character(len=*), optional :: head
integer, intent(in), optional :: ivr(:), ivc(:)
Integer :: err_act, info
character(len=20) :: name='sparse_print'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 700
! 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)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine sparse_print
subroutine get_neigh(a,idx,neigh,n,info,lev)
use psb_error_mod
implicit none
class(psbn_base_sparse_mat), intent(in) :: a
integer, intent(in) :: idx
integer, intent(out) :: n
@ -343,6 +405,7 @@ contains
subroutine allocate_mn(m,n,a)
use psb_error_mod
implicit none
integer, intent(in) :: m,n
class(psbn_base_sparse_mat), intent(inout) :: a
@ -365,6 +428,7 @@ contains
subroutine allocate_mnnz(m,n,nz,a)
use psb_error_mod
implicit none
integer, intent(in) :: m,n,nz
class(psbn_base_sparse_mat), intent(inout) :: a
Integer :: err_act
@ -386,6 +450,7 @@ contains
subroutine reallocate_nz(nz,a)
use psb_error_mod
implicit none
integer, intent(in) :: nz
class(psbn_base_sparse_mat), intent(inout) :: a
Integer :: err_act
@ -407,6 +472,7 @@ contains
subroutine free(a)
use psb_error_mod
implicit none
class(psbn_base_sparse_mat), intent(inout) :: a
Integer :: err_act
character(len=20) :: name='free'

509
base/newserial/psbn_d_base_mat_mod.f03
Unescape Escape View File

@ -11,13 +11,18 @@ module psbn_d_base_mat_mod
procedure, pass(a) :: d_base_cssm
generic, public :: psbn_cssm => d_base_cssm, d_base_cssv
procedure, pass(a) :: csins
procedure, pass(a) :: to_coo
procedure, pass(a) :: from_coo
procedure, pass(a) :: to_fmt
procedure, pass(a) :: from_fmt
procedure, pass(a) :: cp_to_coo
procedure, pass(a) :: cp_from_coo
procedure, pass(a) :: cp_to_fmt
procedure, pass(a) :: cp_from_fmt
procedure, pass(a) :: mv_to_coo
procedure, pass(a) :: mv_from_coo
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,&
& csins, to_coo, from_coo, to_fmt, from_fmt
& csins, cp_to_coo, cp_from_coo, cp_to_fmt, cp_from_fmt, &
& mv_to_coo, mv_from_coo, mv_to_fmt, mv_from_fmt
type, extends(psbn_d_base_sparse_mat) :: psbn_d_coo_sparse_mat
@ -28,6 +33,7 @@ module psbn_d_base_mat_mod
contains
procedure, pass(a) :: get_size => d_coo_get_size
procedure, pass(a) :: get_nzeros => d_coo_get_nzeros
procedure, pass(a) :: set_nzeros => d_coo_set_nzeros
procedure, pass(a) :: d_base_csmm => d_coo_csmm
@ -38,21 +44,34 @@ module psbn_d_base_mat_mod
procedure, pass(a) :: reallocate_nz => d_coo_reallocate_nz
procedure, pass(a) :: allocate_mnnz => d_coo_allocate_mnnz
procedure, pass(a) :: allocate_mn => d_coo_allocate_mn
procedure, pass(a) :: to_coo => d_coo_to_coo
procedure, pass(a) :: from_coo => d_coo_from_coo
procedure, pass(a) :: to_fmt => d_coo_to_fmt
procedure, pass(a) :: from_fmt => d_coo_from_fmt
procedure, pass(a) :: cp_to_coo => d_cp_coo_to_coo
procedure, pass(a) :: cp_from_coo => d_cp_coo_from_coo
procedure, pass(a) :: cp_to_fmt => d_cp_coo_to_fmt
procedure, pass(a) :: cp_from_fmt => d_cp_coo_from_fmt
procedure, pass(a) :: fix => d_fix_coo
procedure, pass(a) :: free => d_coo_free
procedure, pass(a) :: print => d_coo_print
end type psbn_d_coo_sparse_mat
private :: d_coo_get_nzeros, d_coo_set_nzeros, &
& d_coo_csmm, d_coo_csmv, d_coo_cssm, d_coo_cssv, &
& d_coo_csins, d_coo_reallocate_nz, d_coo_allocate_mnnz, &
& d_coo_allocate_mn, d_coo_to_coo, d_coo_from_coo, &
& d_fix_coo, d_coo_free
& d_coo_allocate_mn, d_fix_coo, d_coo_free, &
& d_cp_coo_to_coo, d_cp_coo_from_coo, &
& d_cp_coo_to_fmt, d_cp_coo_from_fmt
interface
subroutine d_fix_coo_inner(nzin,dupl,ia,ja,val,nzout,info,idir)
use psb_const_mod
integer, intent(in) :: nzin,dupl
integer, intent(inout) :: ia(:), ja(:)
real(psb_dpk_), intent(inout) :: val(:)
integer, intent(out) :: nzout, info
integer, intent(in), optional :: idir
end subroutine d_fix_coo_inner
end interface
interface
subroutine d_fix_coo_impl(a,info,idir)
use psb_const_mod
@ -64,45 +83,86 @@ module psbn_d_base_mat_mod
end interface
interface
subroutine d_coo_to_coo_impl(a,b,info)
subroutine d_cp_coo_to_coo_impl(a,b,info)
use psb_const_mod
import psbn_d_coo_sparse_mat
class(psbn_d_coo_sparse_mat), intent(in) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer, intent(out) :: info
end subroutine d_coo_to_coo_impl
end subroutine d_cp_coo_to_coo_impl
end interface
interface
subroutine d_coo_from_coo_impl(a,b,info)
subroutine d_cp_coo_from_coo_impl(a,b,info)
use psb_const_mod
import psbn_d_coo_sparse_mat
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(in) :: b
integer, intent(out) :: info
end subroutine d_coo_from_coo_impl
end subroutine d_cp_coo_from_coo_impl
end interface
interface
subroutine d_coo_to_fmt_impl(a,b,info)
subroutine d_cp_coo_to_fmt_impl(a,b,info)
use psb_const_mod
import psbn_d_coo_sparse_mat, psbn_d_base_sparse_mat
class(psbn_d_coo_sparse_mat), intent(in) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
end subroutine d_coo_to_fmt_impl
end subroutine d_cp_coo_to_fmt_impl
end interface
interface
subroutine d_coo_from_fmt_impl(a,b,info)
subroutine d_cp_coo_from_fmt_impl(a,b,info)
use psb_const_mod
import psbn_d_coo_sparse_mat, psbn_d_base_sparse_mat
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(in) :: b
integer, intent(out) :: info
end subroutine d_coo_from_fmt_impl
end subroutine d_cp_coo_from_fmt_impl
end interface
interface
subroutine d_mv_coo_to_coo_impl(a,b,info)
use psb_const_mod
import psbn_d_coo_sparse_mat
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer, intent(out) :: info
end subroutine d_mv_coo_to_coo_impl
end interface
interface
subroutine d_mv_coo_from_coo_impl(a,b,info)
use psb_const_mod
import psbn_d_coo_sparse_mat
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
end subroutine d_mv_coo_from_coo_impl
end interface
interface
subroutine d_mv_coo_to_fmt_impl(a,b,info)
use psb_const_mod
import psbn_d_coo_sparse_mat, psbn_d_base_sparse_mat
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
end subroutine d_mv_coo_to_fmt_impl
end interface
interface
subroutine d_mv_coo_from_fmt_impl(a,b,info)
use psb_const_mod
import psbn_d_coo_sparse_mat, psbn_d_base_sparse_mat
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
end subroutine d_mv_coo_from_fmt_impl
end interface
interface
subroutine d_coo_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
use psb_const_mod
@ -160,10 +220,10 @@ module psbn_d_base_mat_mod
contains
subroutine to_coo(a,b,info)
subroutine cp_to_coo(a,b,info)
use psb_error_mod
use psb_realloc_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
@ -184,11 +244,12 @@ contains
end if
return
end subroutine to_coo
end subroutine cp_to_coo
subroutine from_coo(a,b,info)
subroutine cp_from_coo(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_base_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(in) :: b
integer, intent(out) :: info
@ -209,12 +270,13 @@ contains
end if
return
end subroutine from_coo
end subroutine cp_from_coo
subroutine to_fmt(a,b,info)
subroutine cp_to_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
@ -235,11 +297,12 @@ contains
end if
return
end subroutine to_fmt
end subroutine cp_to_fmt
subroutine from_fmt(a,b,info)
subroutine cp_from_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_base_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(in) :: b
integer, intent(out) :: info
@ -260,13 +323,120 @@ contains
end if
return
end subroutine from_fmt
end subroutine cp_from_fmt
subroutine mv_to_coo(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_base_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='to_coo'
logical, parameter :: debug=.false.
call psb_erractionsave(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)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine mv_to_coo
subroutine mv_from_coo(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_base_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='from_coo'
logical, parameter :: debug=.false.
call psb_erractionsave(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)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine mv_from_coo
subroutine mv_to_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_base_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='to_fmt'
logical, parameter :: debug=.false.
call psb_erractionsave(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)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine mv_to_fmt
subroutine mv_from_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_base_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='from_fmt'
logical, parameter :: debug=.false.
call psb_erractionsave(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)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine mv_from_fmt
subroutine d_fix_coo(a,info,idir)
use psb_error_mod
use psb_const_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
integer, intent(out) :: info
integer, intent(in), optional :: idir
@ -300,6 +470,7 @@ contains
subroutine csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_base_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: val(:)
integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax
@ -326,6 +497,7 @@ contains
subroutine d_base_csmm(alpha,a,x,beta,y,info,trans)
use psb_error_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(:,:)
@ -352,6 +524,7 @@ contains
subroutine d_base_csmv(alpha,a,x,beta,y,info,trans)
use psb_error_mod
implicit none
class(psbn_d_base_sparse_mat), intent(in) :: a
real(kind(1.d0)), intent(in) :: alpha, beta, x(:)
real(kind(1.d0)), intent(inout) :: y(:)
@ -379,6 +552,7 @@ contains
subroutine d_base_cssm(alpha,a,x,beta,y,info,trans)
use psb_error_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(:,:)
@ -405,6 +579,7 @@ contains
subroutine d_base_cssv(alpha,a,x,beta,y,info,trans)
use psb_error_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(:)
@ -431,9 +606,10 @@ contains
end subroutine d_base_cssv
subroutine d_coo_to_coo(a,b,info)
subroutine d_cp_coo_to_coo(a,b,info)
use psb_error_mod
use psb_realloc_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
@ -444,7 +620,7 @@ contains
call psb_erractionsave(err_act)
info = 0
call d_coo_to_coo_impl(a,b,info)
call d_cp_coo_to_coo_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
@ -460,11 +636,12 @@ contains
end if
return
end subroutine d_coo_to_coo
end subroutine d_cp_coo_to_coo
subroutine d_coo_from_coo(a,b,info)
subroutine d_cp_coo_from_coo(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(in) :: b
integer, intent(out) :: info
@ -475,7 +652,7 @@ contains
call psb_erractionsave(err_act)
info = 0
call d_coo_from_coo_impl(a,b,info)
call d_cp_coo_from_coo_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
@ -491,11 +668,12 @@ contains
end if
return
end subroutine d_coo_from_coo
end subroutine d_cp_coo_from_coo
subroutine d_coo_to_fmt(a,b,info)
subroutine d_cp_coo_to_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
@ -506,7 +684,7 @@ contains
call psb_erractionsave(err_act)
info = 0
call d_coo_to_fmt_impl(a,b,info)
call d_cp_coo_to_fmt_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
@ -522,11 +700,12 @@ contains
end if
return
end subroutine d_coo_to_fmt
end subroutine d_cp_coo_to_fmt
subroutine d_coo_from_fmt(a,b,info)
subroutine d_cp_coo_from_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(in) :: b
integer, intent(out) :: info
@ -537,7 +716,73 @@ contains
call psb_erractionsave(err_act)
info = 0
call d_coo_from_fmt_impl(a,b,info)
call d_cp_coo_from_fmt_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_cp_coo_from_fmt
subroutine d_mv_coo_to_coo(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='to_coo'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
call d_mv_coo_to_coo_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_mv_coo_to_coo
subroutine d_mv_coo_from_coo(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='from_coo'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
call d_mv_coo_from_coo_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
@ -553,20 +798,87 @@ contains
end if
return
end subroutine d_coo_from_fmt
end subroutine d_mv_coo_from_coo
subroutine d_mv_coo_to_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='to_coo'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
call d_mv_coo_to_fmt_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_mv_coo_to_fmt
subroutine d_mv_coo_from_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='from_coo'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
call d_mv_coo_from_fmt_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_mv_coo_from_fmt
subroutine d_coo_reallocate_nz(nz,a)
use psb_error_mod
use psb_realloc_mod
implicit none
integer, intent(in) :: nz
class(psbn_d_coo_sparse_mat), intent(inout) :: a
Integer :: err_act
Integer :: err_act, info
character(len=20) :: name='d_coo_reallocate_nz'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
call psb_realloc(nx,a%ia,a%ja,a%val,info)
call psb_realloc(nz,a%ia,a%ja,a%val,info)
if (info /= 0) then
call psb_errpush(4000,name)
@ -588,7 +900,32 @@ contains
end subroutine d_coo_reallocate_nz
function d_coo_get_size(a) result(res)
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
integer :: res
res = -1
if (allocated(a%ia)) res = size(a%ia)
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_coo_get_size
function d_coo_get_nzeros(a) result(res)
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
integer :: res
res = a%nnz
@ -596,6 +933,7 @@ contains
subroutine d_coo_set_nzeros(nz,a)
implicit none
integer, intent(in) :: nz
class(psbn_d_coo_sparse_mat), intent(inout) :: a
@ -607,6 +945,7 @@ contains
subroutine d_coo_csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
use psb_error_mod
use psb_realloc_mod
implicit none
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
@ -649,7 +988,8 @@ contains
end if
if (nz == 0) return
nza = a%get_nzeros()
!!$ write(0,*) 'On entry to csins: ',nza
call d_coo_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
if (info /= 0) goto 9999
@ -670,6 +1010,7 @@ contains
subroutine d_coo_csmv(alpha,a,x,beta,y,info,trans)
use psb_error_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:)
real(psb_dpk_), intent(inout) :: y(:)
@ -713,6 +1054,7 @@ contains
subroutine d_coo_csmm(alpha,a,x,beta,y,info,trans)
use psb_error_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:,:)
real(psb_dpk_), intent(inout) :: y(:,:)
@ -752,6 +1094,7 @@ contains
subroutine d_coo_cssv(alpha,a,x,beta,y,info,trans)
use psb_error_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:)
real(psb_dpk_), intent(inout) :: y(:)
@ -805,6 +1148,7 @@ contains
subroutine d_coo_cssm(alpha,a,x,beta,y,info,trans)
use psb_error_mod
implicit none
class(psbn_d_coo_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:,:)
real(psb_dpk_), intent(inout) :: y(:,:)
@ -854,6 +1198,7 @@ contains
subroutine d_coo_free(a)
implicit none
class(psbn_d_coo_sparse_mat), intent(inout) :: a
@ -871,6 +1216,7 @@ contains
subroutine d_coo_allocate_mnnz(m,n,nz,a)
use psb_error_mod
use psb_realloc_mod
implicit none
integer, intent(in) :: m,n,nz
class(psbn_d_coo_sparse_mat), intent(inout) :: a
Integer :: err_act, info
@ -904,8 +1250,9 @@ contains
call a%set_nzeros(0)
call a%set_bld()
call a%set_triangle(.false.)
call a%set_dupl(psbn_dupl_def_)
end if
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
@ -924,6 +1271,7 @@ contains
subroutine d_coo_allocate_mn(m,n,a)
use psb_error_mod
use psb_realloc_mod
implicit none
integer, intent(in) :: m,n
class(psbn_d_coo_sparse_mat), intent(inout) :: a
Integer :: err_act, info, nz
@ -961,6 +1309,79 @@ contains
end subroutine d_coo_allocate_mn
subroutine d_coo_print(iout,a,iv,eirs,eics,head,ivr,ivc)
use psb_spmat_type
use psb_string_mod
implicit none
integer, intent(in) :: iout
class(psbn_d_coo_sparse_mat), intent(in) :: a
integer, intent(in), optional :: iv(:)
integer, intent(in), optional :: eirs,eics
character(len=*), optional :: head
integer, intent(in), optional :: ivr(:), ivc(:)
Integer :: err_act
character(len=20) :: name='d_coo_print'
logical, parameter :: debug=.false.
character(len=80) :: frmtv
integer :: irs,ics,i,j, nmx, ni, nr, nc, nz
if (present(eirs)) then
irs = eirs
else
irs = 0
endif
if (present(eics)) then
ics = eics
else
ics = 0
endif
if (present(head)) then
write(iout,'(a)') '%%MatrixMarket matrix coordinate real general'
write(iout,'(a,a)') '% ',head
write(iout,'(a)') '%'
write(iout,'(a,a)') '% COO'
endif
nr = a%get_nrows()
nc = a%get_ncols()
nz = a%get_nzeros()
nmx = max(nr,nc,1)
ni = floor(log10(1.0*nmx)) + 1
write(frmtv,'(a,i3.3,a,i3.3,a)') '(2(i',ni,',1x),es26.18,1x,2(i',ni,',1x))'
write(iout,*) nr, nc, nz
if(present(iv)) then
do j=1,a%get_nzeros()
write(iout,frmtv) iv(a%ia(j)),iv(a%ja(j)),a%val(j)
enddo
else
if (present(ivr).and..not.present(ivc)) then
do j=1,a%get_nzeros()
write(iout,frmtv) ivr(a%ia(j)),a%ja(j),a%val(j)
enddo
else if (present(ivr).and.present(ivc)) then
do j=1,a%get_nzeros()
write(iout,frmtv) ivr(a%ia(j)),ivc(a%ja(j)),a%val(j)
enddo
else if (.not.present(ivr).and.present(ivc)) then
do j=1,a%get_nzeros()
write(iout,frmtv) a%ia(j),ivc(a%ja(j)),a%val(j)
enddo
else if (.not.present(ivr).and..not.present(ivc)) then
do j=1,a%get_nzeros()
write(iout,frmtv) a%ia(j),a%ja(j),a%val(j)
enddo
endif
endif
end subroutine d_coo_print
end module psbn_d_base_mat_mod

1028
base/newserial/psbn_d_coo_impl.f03
View File

File diff suppressed because it is too large Load Diff

597
base/newserial/psbn_d_csr_impl.f03
Unescape Escape View File

@ -2,6 +2,7 @@
subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
use psb_error_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_csr_csmv_impl
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:)
real(psb_dpk_), intent(inout) :: y(:)
@ -63,7 +64,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
if (alpha == done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -73,7 +74,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
else if (alpha == -done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -83,7 +84,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
else
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -97,7 +98,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
if (alpha == done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -107,7 +108,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
else if (alpha == -done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -117,7 +118,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
else
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -130,7 +131,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
if (alpha == done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -140,7 +141,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
else if (alpha == -done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -150,7 +151,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
else
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -163,7 +164,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
if (alpha == done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -173,7 +174,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
else if (alpha == -done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -183,7 +184,7 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans)
else
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j))
enddo
@ -265,6 +266,7 @@ end subroutine d_csr_csmv_impl
subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
use psb_error_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_csr_csmm_impl
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:,:)
real(psb_dpk_), intent(inout) :: y(:,:)
@ -279,6 +281,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
character(len=20) :: name='d_csr_csmm'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
if (present(trans)) then
@ -330,7 +333,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
if (alpha == done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -340,7 +343,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
else if (alpha == -done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -350,7 +353,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
else
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -364,7 +367,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
if (alpha == done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -374,7 +377,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
else if (alpha == -done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -384,7 +387,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
else
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -397,7 +400,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
if (alpha == done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -407,7 +410,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
else if (alpha == -done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -417,7 +420,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
else
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -430,7 +433,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
if (alpha == done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -440,7 +443,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
else if (alpha == -done) then
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -450,7 +453,7 @@ subroutine d_csr_csmm_impl(alpha,a,x,beta,y,info,trans)
else
do i=1,m
acc = zero
acc = dzero
do j=a%irp(i), a%irp(i+1)-1
acc = acc + a%val(j) * x(a%ja(j),:)
enddo
@ -533,6 +536,7 @@ end subroutine d_csr_csmm_impl
subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
use psb_error_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_csr_cssv_impl
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:)
real(psb_dpk_), intent(inout) :: y(:)
@ -548,6 +552,7 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
character(len=20) :: name='d_csr_cssv'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
if (present(trans)) then
@ -617,6 +622,7 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
contains
subroutine inner_csrsv(tra,a,x,y)
implicit none
logical, intent(in) :: tra
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: x(:)
@ -724,6 +730,7 @@ end subroutine d_csr_cssv_impl
subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
use psb_error_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_csr_cssm_impl
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:,:)
real(psb_dpk_), intent(inout) :: y(:,:)
@ -739,6 +746,7 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
character(len=20) :: name='d_base_csmm'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
if (present(trans)) then
@ -821,6 +829,7 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
contains
subroutine inner_csrsm(tra,a,x,y,info)
implicit none
logical, intent(in) :: tra
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: x(:,:)
@ -829,6 +838,7 @@ contains
integer :: i,j,k,m, ir, jc
real(psb_dpk_), allocatable :: acc(:)
info = 0
allocate(acc(size(x,2)), stat=info)
if(info /= 0) then
info=4010
@ -936,6 +946,7 @@ subroutine d_csr_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
use psb_error_mod
use psb_realloc_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_csr_csins_impl
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: val(:)
@ -949,7 +960,7 @@ subroutine d_csr_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
logical, parameter :: debug=.false.
integer :: nza, i,j,k, nzl, isza, int_err(5)
info = 0
nza = a%get_nzeros()
if (a%is_bld()) then
@ -957,20 +968,11 @@ subroutine d_csr_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
info = 1121
else if (a%is_upd()) then
if (a%is_sorted()) then
call d_csr_srch_upd(nz,ia,ja,val,a,&
& imin,imax,jmin,jmax,info,gtl)
if (info /= 0) then
!!$#ifdef FIXED_NAG_SEGV
!!$ call d_csr_srch_upd(nz,ia,ja,val,a,&
!!$ & imin,imax,jmin,jmax,info,gtl)
!!$#else
call d_csr_srch_upd(nz,ia,ja,val,&
& a%irp,a%ja,a%val,&
& a%get_dupl(),a%get_nzeros(),a%get_nrows(),&
& info,gtl)
!!$#endif
else
info = 1121
end if
@ -998,255 +1000,40 @@ subroutine d_csr_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
contains
subroutine d_csr_srch_upd(nz,ia,ja,val,a,&
& imin,imax,jmin,jmax,info,gtl)
!!$#ifdef FIXED_NAG_SEGV
!!$ subroutine d_csr_srch_upd(nz,ia,ja,val,a,&
!!$ & imin,imax,jmin,jmax,info,gtl)
!!$
!!$ use psb_const_mod
!!$ use psb_realloc_mod
!!$ use psb_string_mod
!!$ use psb_serial_mod
!!$ implicit none
!!$
!!$ class(psbn_d_csr_sparse_mat), intent(inout) :: a
!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax
!!$ integer, intent(in) :: ia(:),ja(:)
!!$ real(psb_dpk_), intent(in) :: val(:)
!!$ integer, intent(out) :: info
!!$ integer, intent(in), optional :: gtl(:)
!!$ integer :: i,ir,ic, ilr, ilc, ip, &
!!$ & i1,i2,nc,nnz,dupl,ng
!!$ integer :: debug_level, debug_unit
!!$ character(len=20) :: name='d_csr_srch_upd'
!!$
!!$ info = 0
!!$ debug_unit = psb_get_debug_unit()
!!$ debug_level = psb_get_debug_level()
!!$
!!$ dupl = a%get_dupl()
!!$
!!$ if (.not.a%is_sorted()) then
!!$ info = -4
!!$ return
!!$ end if
!!$
!!$ ilr = -1
!!$ ilc = -1
!!$ nnz = a%get_nzeros()
!!$
!!$
!!$ if (present(gtl)) then
!!$ ng = size(gtl)
!!$
!!$ select case(dupl)
!!$ case(psbn_dupl_ovwrt_,psbn_dupl_err_)
!!$ ! Overwrite.
!!$ ! Cannot test for error, should have been caught earlier.
!!$ do i=1, nz
!!$ ir = ia(i)
!!$ ic = ja(i)
!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
!!$ ir = gtl(ir)
!!$ if ((ir > 0).and.(ir <= a%m)) then
!!$ ic = gtl(ic)
!!$ if (ir /= ilr) then
!!$ i1 = psb_ibsrch(ir,nnz,a%ia)
!!$ i2 = i1
!!$ do
!!$ if (i2+1 > nnz) exit
!!$ if (a%ia(i2+1) /= a%ia(i2)) exit
!!$ i2 = i2 + 1
!!$ end do
!!$ do
!!$ if (i1-1 < 1) exit
!!$ if (a%ia(i1-1) /= a%ia(i1)) exit
!!$ i1 = i1 - 1
!!$ end do
!!$ ilr = ir
!!$ else
!!$ i1 = 1
!!$ i2 = 1
!!$ end if
!!$ nc = i2-i1+1
!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2))
!!$ if (ip>0) then
!!$ a%val(i1+ip-1) = val(i)
!!$ else
!!$ info = i
!!$ return
!!$ end if
!!$ else
!!$ if (debug_level >= psb_debug_serial_) &
!!$ & write(debug_unit,*) trim(name),&
!!$ & ': Discarding row that does not belong to us.'
!!$ endif
!!$ end if
!!$ end do
!!$ case(psbn_dupl_add_)
!!$ ! Add
!!$ do i=1, nz
!!$ ir = ia(i)
!!$ ic = ja(i)
!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then
!!$ ir = gtl(ir)
!!$ ic = gtl(ic)
!!$ if ((ir > 0).and.(ir <= a%m)) then
!!$
!!$ if (ir /= ilr) then
!!$ i1 = psb_ibsrch(ir,nnz,a%ia)
!!$ i2 = i1
!!$ do
!!$ if (i2+1 > nnz) exit
!!$ if (a%ia(i2+1) /= a%ia(i2)) exit
!!$ i2 = i2 + 1
!!$ end do
!!$ do
!!$ if (i1-1 < 1) exit
!!$ if (a%ia(i1-1) /= a%ia(i1)) exit
!!$ i1 = i1 - 1
!!$ end do
!!$ ilr = ir
!!$ else
!!$ i1 = 1
!!$ i2 = 1
!!$ end if
!!$ nc = i2-i1+1
!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2))
!!$ if (ip>0) then
!!$ a%val(i1+ip-1) = a%val(i1+ip-1) + val(i)
!!$ else
!!$ info = i
!!$ return
!!$ end if
!!$ else
!!$ if (debug_level >= psb_debug_serial_) &
!!$ & write(debug_unit,*) trim(name),&
!!$ & ': Discarding row that does not belong to us.'
!!$ end if
!!$ end if
!!$ end do
!!$
!!$ case default
!!$ info = -3
!!$ if (debug_level >= psb_debug_serial_) &
!!$ & write(debug_unit,*) trim(name),&
!!$ & ': Duplicate handling: ',dupl
!!$ end select
!!$
!!$ else
!!$
!!$ select case(dupl)
!!$ case(psbn_dupl_ovwrt_,psbn_dupl_err_)
!!$ ! Overwrite.
!!$ ! Cannot test for error, should have been caught earlier.
!!$ do i=1, nz
!!$ ir = ia(i)
!!$ ic = ja(i)
!!$ if ((ir > 0).and.(ir <= a%m)) then
!!$
!!$ if (ir /= ilr) then
!!$ i1 = psb_ibsrch(ir,nnz,a%ia)
!!$ i2 = i1
!!$ do
!!$ if (i2+1 > nnz) exit
!!$ if (a%ia(i2+1) /= a%ia(i2)) exit
!!$ i2 = i2 + 1
!!$ end do
!!$ do
!!$ if (i1-1 < 1) exit
!!$ if (a%ia(i1-1) /= a%ia(i1)) exit
!!$ i1 = i1 - 1
!!$ end do
!!$ ilr = ir
!!$ else
!!$ i1 = 1
!!$ i2 = 1
!!$ end if
!!$ nc = i2-i1+1
!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2))
!!$ if (ip>0) then
!!$ a%val(i1+ip-1) = val(i)
!!$ else
!!$ info = i
!!$ return
!!$ end if
!!$ end if
!!$ end do
!!$
!!$ case(psbn_dupl_add_)
!!$ ! Add
!!$ do i=1, nz
!!$ ir = ia(i)
!!$ ic = ja(i)
!!$ if ((ir > 0).and.(ir <= a%m)) then
!!$
!!$ if (ir /= ilr) then
!!$ i1 = psb_ibsrch(ir,nnz,a%ia)
!!$ i2 = i1
!!$ do
!!$ if (i2+1 > nnz) exit
!!$ if (a%ia(i2+1) /= a%ia(i2)) exit
!!$ i2 = i2 + 1
!!$ end do
!!$ do
!!$ if (i1-1 < 1) exit
!!$ if (a%ia(i1-1) /= a%ia(i1)) exit
!!$ i1 = i1 - 1
!!$ end do
!!$ ilr = ir
!!$ else
!!$ i1 = 1
!!$ i2 = 1
!!$ end if
!!$ nc = i2-i1+1
!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2))
!!$ if (ip>0) then
!!$ a%val(i1+ip-1) = a%val(i1+ip-1) + val(i)
!!$ else
!!$ info = i
!!$ return
!!$ end if
!!$ end if
!!$ end do
!!$
!!$ case default
!!$ info = -3
!!$ if (debug_level >= psb_debug_serial_) &
!!$ & write(debug_unit,*) trim(name),&
!!$ & ': Duplicate handling: ',dupl
!!$ end select
!!$
!!$ end if
!!$
!!$ end subroutine d_csr_srch_upd
!!$
!!$#else
subroutine d_csr_srch_upd(nz,ia,ja,val,&
& airp,aja,aval,dupl,nza,nra,&
& info,gtl)
use psb_error_mod
use psb_sort_mod
use psb_const_mod
use psb_realloc_mod
use psb_string_mod
use psb_serial_mod
implicit none
integer, intent(inout) :: airp(:),aja(:)
real(psb_dpk_), intent(inout) :: aval(:)
integer, intent(in) :: nz, dupl,nza, nra
class(psbn_d_csr_sparse_mat), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax
integer, intent(in) :: ia(:),ja(:)
real(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: info
integer, intent(in), optional :: gtl(:)
integer :: i,ir,ic, ilr, ilc, ip, &
& i1,i2,nc,nnz,dupl,ng
integer :: debug_level, debug_unit
character(len=20) :: name='d_csr_srch_upd'
integer :: i,ir,ic, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m, ng
info = 0
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
dupl = a%get_dupl()
if (.not.a%is_sorted()) then
info = -4
return
end if
ilr = -1
ilc = -1
nnz = a%get_nzeros()
if (present(gtl)) then
ng = size(gtl)
@ -1265,18 +1052,18 @@ contains
ir = gtl(ir)
ic = gtl(ic)
if ((ir > 0).and.(ir <= a%m)) then
i1 = airp(ir)
i2 = airp(ir+1)
i1 = a%irp(ir)
i2 = a%irp(ir+1)
nc=i2-i1
ip = psb_ibsrch(ic,nc,aja(i1:i2-1))
ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1))
if (ip>0) then
aval(i1+ip-1) = val(i)
a%val(i1+ip-1) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',aja(i1:i2-1)
& ' : ',a%ja(i1:i2-1)
info = i
return
end if
@ -1301,17 +1088,17 @@ contains
ir = gtl(ir)
ic = gtl(ic)
if ((ir > 0).and.(ir <= a%m)) then
i1 = airp(ir)
i2 = airp(ir+1)
i1 = a%irp(ir)
i2 = a%irp(ir+1)
nc = i2-i1
ip = psb_ibsrch(ic,nc,aja(i1:i2-1))
ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1))
if (ip>0) then
aval(i1+ip-1) = aval(i1+ip-1) + val(i)
a%val(i1+ip-1) = a%val(i1+ip-1) + val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',aja(i1:i2-1)
& ' : ',a%ja(i1:i2-1)
info = i
return
end if
@ -1346,18 +1133,18 @@ contains
if ((ir > 0).and.(ir <= a%m)) then
i1 = airp(ir)
i2 = airp(ir+1)
i1 = a%irp(ir)
i2 = a%irp(ir+1)
nc=i2-i1
ip = psb_ibsrch(ic,nc,aja(i1:i2-1))
ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1))
if (ip>0) then
aval(i1+ip-1) = val(i)
a%val(i1+ip-1) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',aja(i1:i2-1)
& ' : ',a%ja(i1:i2-1)
info = i
return
end if
@ -1378,12 +1165,12 @@ contains
ir = ia(i)
ic = ja(i)
if ((ir > 0).and.(ir <= a%m)) then
i1 = airp(ir)
i2 = airp(ir+1)
i1 = a%irp(ir)
i2 = a%irp(ir+1)
nc = i2-i1
ip = psb_ibsrch(ic,nc,aja(i1:i2-1))
ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1))
if (ip>0) then
aval(i1+ip-1) = aval(i1+ip-1) + val(i)
a%val(i1+ip-1) = a%val(i1+ip-1) + val(i)
else
info = i
return
@ -1403,6 +1190,232 @@ contains
end select
end if
end subroutine d_csr_srch_upd
!!$#endif
end subroutine d_csr_csins_impl
subroutine d_cp_csr_from_coo_impl(a,b,info)
use psb_const_mod
use psb_realloc_mod
use psbn_d_base_mat_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_cp_csr_from_coo_impl
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(in) :: b
integer, intent(out) :: info
type(psbn_d_coo_sparse_mat) :: tmp
integer, allocatable :: itemp(:)
!locals
logical :: rwshr_
Integer :: nza, nr, i,j,irw, idl,err_act, nc
Integer, Parameter :: maxtry=8
integer :: debug_level, debug_unit
character(len=20) :: name
info = 0
! This is to have fix_coo called behind the scenes
call tmp%cp_from_coo(b,info)
if (info ==0) call a%mv_from_coo(tmp,info)
end subroutine d_cp_csr_from_coo_impl
subroutine d_cp_csr_to_coo_impl(a,b,info)
use psb_const_mod
use psbn_d_base_mat_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_cp_csr_to_coo_impl
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, allocatable :: itemp(:)
!locals
logical :: rwshr_
Integer :: nza, nr, nc,i,j,irw, idl,err_act
Integer, Parameter :: maxtry=8
integer :: debug_level, debug_unit
character(len=20) :: name
info = 0
nr = a%get_nrows()
nc = a%get_ncols()
nza = a%get_nzeros()
call b%allocate(nr,nr,nza)
do i=1, nr
do j=a%irp(i),a%irp(i+1)-1
b%ia(j) = i
b%ja(j) = a%ja(j)
b%val(j) = a%val(j)
end do
end do
call b%set_nzeros(a%get_nzeros())
call b%set_nrows(a%get_nrows())
call b%set_ncols(a%get_ncols())
call b%set_dupl(a%get_dupl())
call b%set_state(a%get_state())
call b%set_triangle(a%is_triangle())
call b%set_upper(a%is_upper())
call b%fix(info)
end subroutine d_cp_csr_to_coo_impl
subroutine d_mv_csr_from_coo_impl(a,b,info)
use psb_const_mod
use psb_realloc_mod
use psbn_d_base_mat_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_mv_csr_from_coo_impl
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
integer, allocatable :: itemp(:)
!locals
logical :: rwshr_
Integer :: nza, nr, i,j,irw, idl,err_act, nc
Integer, Parameter :: maxtry=8
integer :: debug_level, debug_unit
character(len=20) :: name
info = 0
call b%fix(info)
if (info /= 0) return
nr = b%get_nrows()
nc = b%get_ncols()
nza = b%get_nzeros()
call a%set_nrows(b%get_nrows())
call a%set_ncols(b%get_ncols())
call a%set_dupl(b%get_dupl())
call a%set_state(b%get_state())
call a%set_triangle(b%is_triangle())
call a%set_upper(b%is_upper())
! Dirty trick: call move_alloc to have the new data allocated just once.
call move_alloc(b%ia,itemp)
call move_alloc(b%ja,a%ja)
call move_alloc(b%val,a%val)
call psb_realloc(max(nr+1,nc+1),a%irp,info)
call b%free()
if (nza <= 0) then
a%irp(:) = 1
else
a%irp(1) = 1
if (nr < itemp(nza)) then
write(debug_unit,*) trim(name),': RWSHR=.false. : ',&
&nr,itemp(nza),' Expect trouble!'
info = 12
end if
j = 1
i = 1
irw = itemp(j)
outer: do
inner: do
if (i >= irw) exit inner
if (i>nr) then
write(debug_unit,*) trim(name),&
& 'Strange situation: i>nr ',i,nr,j,nza,irw,idl
exit outer
end if
a%irp(i+1) = a%irp(i)
i = i + 1
end do inner
j = j + 1
if (j > nza) exit
if (itemp(j) /= irw) then
a%irp(i+1) = j
irw = itemp(j)
i = i + 1
endif
if (i>nr) exit
enddo outer
!
! Cleanup empty rows at the end
!
if (j /= (nza+1)) then
write(debug_unit,*) trim(name),': Problem from loop :',j,nza
info = 13
endif
do
if (i>nr) exit
a%irp(i+1) = j
i = i + 1
end do
endif
end subroutine d_mv_csr_from_coo_impl
subroutine d_mv_csr_to_coo_impl(a,b,info)
use psb_const_mod
use psb_realloc_mod
use psbn_d_base_mat_mod
use psbn_d_csr_mat_mod, psb_protect_name => d_mv_csr_to_coo_impl
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer, intent(out) :: info
integer, allocatable :: itemp(:)
!locals
logical :: rwshr_
Integer :: nza, nr, i,j,irw, idl,err_act, nc
Integer, Parameter :: maxtry=8
integer :: debug_level, debug_unit
character(len=20) :: name
info = 0
nr = a%get_nrows()
nc = a%get_ncols()
nza = a%get_nzeros()
call b%set_nrows(a%get_nrows())
call b%set_ncols(a%get_ncols())
call b%set_dupl(a%get_dupl())
call b%set_state(a%get_state())
call b%set_triangle(a%is_triangle())
call b%set_upper(a%is_upper())
! Dirty trick: call move_alloc to have the new data allocated just once.
call move_alloc(a%irp,itemp)
call move_alloc(a%ja,b%ja)
call move_alloc(a%val,b%val)
call psb_realloc(nza,b%ia,info)
if (info /= 0) return
if (allocated(itemp)) then
do i=1, nr
do j=itemp(i),itemp(i+1)-1
b%ia(j) = i
end do
end do
end if
call b%fix(info)
end subroutine d_mv_csr_to_coo_impl

399
base/newserial/psbn_d_csr_mat_mod.f03
Unescape Escape View File

@ -17,36 +17,112 @@ module psbn_d_csr_mat_mod
procedure, pass(a) :: csins => d_csr_csins
procedure, pass(a) :: allocate_mnnz => d_csr_allocate_mnnz
procedure, pass(a) :: allocate_mn => d_csr_allocate_mn
procedure, pass(a) :: to_coo => d_csr_to_coo
procedure, pass(a) :: from_coo => d_csr_from_coo
procedure, pass(a) :: cp_to_coo => d_cp_csr_to_coo
procedure, pass(a) :: cp_from_coo => d_cp_csr_from_coo
procedure, pass(a) :: cp_to_fmt => d_cp_csr_to_fmt
procedure, pass(a) :: cp_from_fmt => d_cp_csr_from_fmt
procedure, pass(a) :: mv_to_coo => d_mv_csr_to_coo
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) :: free => d_csr_free
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_csins, d_csr_reallocate_nz, d_csr_allocate_mnnz, &
& d_csr_allocate_mn, d_csr_to_coo, d_csr_from_coo, &
& d_csr_free
& d_csr_allocate_mn, d_csr_free, d_csr_print, &
& d_cp_csr_to_coo, d_cp_csr_from_coo, &
& d_mv_csr_to_coo, d_mv_csr_from_coo
interface
subroutine d_csr_to_coo_impl(a,b,info)
subroutine d_cp_csr_to_fmt_impl(a,b,info)
use psb_const_mod
use psbn_d_base_mat_mod
import psbn_d_csr_sparse_mat
class(psbn_d_csr_sparse_mat), intent(in) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
end subroutine d_cp_csr_to_fmt_impl
end interface
interface
subroutine d_cp_csr_from_fmt_impl(a,b,info)
use psb_const_mod
use psbn_d_base_mat_mod
import psbn_d_csr_sparse_mat
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(in) :: b
integer, intent(out) :: info
end subroutine d_cp_csr_from_fmt_impl
end interface
interface
subroutine d_cp_csr_to_coo_impl(a,b,info)
use psb_const_mod
use psbn_d_base_mat_mod
import psbn_d_csr_sparse_mat
class(psbn_d_csr_sparse_mat), intent(in) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer, intent(out) :: info
end subroutine d_csr_to_coo_impl
end subroutine d_cp_csr_to_coo_impl
end interface
interface
subroutine d_csr_from_coo_impl(a,b,info)
subroutine d_cp_csr_from_coo_impl(a,b,info)
use psb_const_mod
use psbn_d_base_mat_mod
import psbn_d_csr_sparse_mat
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(in) :: b
integer, intent(out) :: info
end subroutine d_csr_from_coo_impl
end subroutine d_cp_csr_from_coo_impl
end interface
interface
subroutine d_mv_csr_to_fmt_impl(a,b,info)
use psb_const_mod
use psbn_d_base_mat_mod
import psbn_d_csr_sparse_mat
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
end subroutine d_mv_csr_to_fmt_impl
end interface
interface
subroutine d_mv_csr_from_fmt_impl(a,b,info)
use psb_const_mod
use psbn_d_base_mat_mod
import psbn_d_csr_sparse_mat
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
end subroutine d_mv_csr_from_fmt_impl
end interface
interface
subroutine d_mv_csr_to_coo_impl(a,b,info)
use psb_const_mod
use psbn_d_base_mat_mod
import psbn_d_csr_sparse_mat
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer, intent(out) :: info
end subroutine d_mv_csr_to_coo_impl
end interface
interface
subroutine d_mv_csr_from_coo_impl(a,b,info)
use psb_const_mod
use psbn_d_base_mat_mod
import psbn_d_csr_sparse_mat
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
end subroutine d_mv_csr_from_coo_impl
end interface
interface
@ -110,6 +186,7 @@ contains
subroutine d_csr_reallocate_nz(nz,a)
use psb_error_mod
use psb_realloc_mod
implicit none
integer, intent(in) :: nz
class(psbn_d_csr_sparse_mat), intent(inout) :: a
Integer :: err_act, info
@ -118,9 +195,9 @@ contains
call psb_erractionsave(err_act)
call psb_realloc(nx,a%ja,info)
if (info == 0) call psb_realloc(nx,a%val,info)
if (info == 0) call psb_realloc(max(nx,a%m+1,a%n+1),a%irp,info)
call psb_realloc(nz,a%ja,info)
if (info == 0) call psb_realloc(nz,a%val,info)
if (info == 0) call psb_realloc(max(nz,a%m+1,a%n+1),a%irp,info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
@ -141,6 +218,7 @@ contains
end subroutine d_csr_reallocate_nz
function d_csr_get_nzeros(a) result(res)
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
integer :: res
res = a%irp(a%m+1)-1
@ -150,6 +228,7 @@ contains
subroutine d_csr_csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl)
use psb_const_mod
use psb_error_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: val(:)
integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax
@ -212,6 +291,7 @@ contains
subroutine d_csr_csmv(alpha,a,x,beta,y,info,trans)
use psb_error_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:)
real(psb_dpk_), intent(inout) :: y(:)
@ -255,6 +335,7 @@ contains
subroutine d_csr_csmm(alpha,a,x,beta,y,info,trans)
use psb_error_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:,:)
real(psb_dpk_), intent(inout) :: y(:,:)
@ -294,6 +375,7 @@ contains
subroutine d_csr_cssv(alpha,a,x,beta,y,info,trans)
use psb_error_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:)
real(psb_dpk_), intent(inout) :: y(:)
@ -347,6 +429,7 @@ contains
subroutine d_csr_cssm(alpha,a,x,beta,y,info,trans)
use psb_error_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:,:)
real(psb_dpk_), intent(inout) :: y(:,:)
@ -396,6 +479,7 @@ contains
subroutine d_csr_free(a)
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
@ -411,9 +495,10 @@ contains
end subroutine d_csr_free
subroutine d_csr_to_coo(a,b,info)
subroutine d_cp_csr_to_coo(a,b,info)
use psb_error_mod
use psb_realloc_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
@ -424,7 +509,7 @@ contains
call psb_erractionsave(err_act)
info = 0
call d_csr_to_coo_impl(a,b,info)
call d_cp_csr_to_coo_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
@ -440,11 +525,12 @@ contains
end if
return
end subroutine d_csr_to_coo
end subroutine d_cp_csr_to_coo
subroutine d_csr_from_coo(a,b,info)
subroutine d_cp_csr_from_coo(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(in) :: b
integer, intent(out) :: info
@ -455,7 +541,72 @@ contains
call psb_erractionsave(err_act)
info = 0
call d_csr_from_coo_impl(a,b,info)
call d_cp_csr_from_coo_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_cp_csr_from_coo
subroutine d_cp_csr_to_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(in) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='to_fmt'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
!!$ call d_cp_csr_to_fmt_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_cp_csr_to_fmt
subroutine d_cp_csr_from_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(in) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='from_fmt'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
!!$ call d_cp_csr_from_fmt_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
@ -471,12 +622,143 @@ contains
end if
return
end subroutine d_csr_from_coo
end subroutine d_cp_csr_from_fmt
subroutine d_mv_csr_to_coo(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(out) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='to_coo'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
call d_mv_csr_to_coo_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_mv_csr_to_coo
subroutine d_mv_csr_from_coo(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_coo_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='from_coo'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
call d_mv_csr_from_coo_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_mv_csr_from_coo
subroutine d_mv_csr_to_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='to_fmt'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
!!$ call d_mv_csr_to_fmt_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_mv_csr_to_fmt
subroutine d_mv_csr_from_fmt(a,b,info)
use psb_error_mod
use psb_realloc_mod
implicit none
class(psbn_d_csr_sparse_mat), intent(inout) :: a
class(psbn_d_base_sparse_mat), intent(inout) :: b
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='from_fmt'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = 0
!!$ call d_mv_csr_from_fmt_impl(a,b,info)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
call psb_errpush(info,name)
if (err_act /= psb_act_ret_) then
call psb_error()
end if
return
end subroutine d_mv_csr_from_fmt
subroutine d_csr_allocate_mnnz(m,n,nz,a)
use psb_error_mod
use psb_realloc_mod
implicit none
integer, intent(in) :: m,n,nz
class(psbn_d_csr_sparse_mat), intent(inout) :: a
Integer :: err_act, info
@ -530,6 +812,7 @@ contains
subroutine d_csr_allocate_mn(m,n,a)
use psb_error_mod
use psb_realloc_mod
implicit none
integer, intent(in) :: m,n
class(psbn_d_csr_sparse_mat), intent(inout) :: a
Integer :: err_act, info, nz
@ -566,4 +849,86 @@ contains
end subroutine d_csr_allocate_mn
subroutine d_csr_print(iout,a,iv,eirs,eics,head,ivr,ivc)
use psb_spmat_type
use psb_string_mod
implicit none
integer, intent(in) :: iout
class(psbn_d_csr_sparse_mat), intent(in) :: a
integer, intent(in), optional :: iv(:)
integer, intent(in), optional :: eirs,eics
character(len=*), optional :: head
integer, intent(in), optional :: ivr(:), ivc(:)
Integer :: err_act
character(len=20) :: name='d_csr_print'
logical, parameter :: debug=.false.
character(len=80) :: frmtv
integer :: irs,ics,i,j, nmx, ni, nr, nc, nz
if (present(eirs)) then
irs = eirs
else
irs = 0
endif
if (present(eics)) then
ics = eics
else
ics = 0
endif
if (present(head)) then
write(iout,'(a)') '%%MatrixMarket matrix coordinate real general'
write(iout,'(a,a)') '% ',head
write(iout,'(a)') '%'
write(iout,'(a,a)') '% COO'
endif
nr = a%get_nrows()
nc = a%get_ncols()
nz = a%get_nzeros()
nmx = max(nr,nc,1)
ni = floor(log10(1.0*nmx)) + 1
write(frmtv,'(a,i3.3,a,i3.3,a)') '(2(i',ni,',1x),es26.18,1x,2(i',ni,',1x))'
write(iout,*) nr, nc, nz
if(present(iv)) then
do i=1, nr
do j=a%irp(i),a%irp(i+1)-1
write(iout,frmtv) iv(i),iv(a%ja(j)),a%val(j)
end do
enddo
else
if (present(ivr).and..not.present(ivc)) then
do i=1, nr
do j=a%irp(i),a%irp(i+1)-1
write(iout,frmtv) ivr(i),(a%ja(j)),a%val(j)
end do
enddo
else if (present(ivr).and.present(ivc)) then
do i=1, nr
do j=a%irp(i),a%irp(i+1)-1
write(iout,frmtv) ivr(i),ivc(a%ja(j)),a%val(j)
end do
enddo
else if (.not.present(ivr).and.present(ivc)) then
do i=1, nr
do j=a%irp(i),a%irp(i+1)-1
write(iout,frmtv) (i),ivc(a%ja(j)),a%val(j)
end do
enddo
else if (.not.present(ivr).and..not.present(ivc)) then
do i=1, nr
do j=a%irp(i),a%irp(i+1)-1
write(iout,frmtv) (i),(a%ja(j)),a%val(j)
end do
enddo
endif
endif
end subroutine d_csr_print
end module psbn_d_csr_mat_mod

2
base/newserial/psbn_mat_impl.f03
Unescape Escape View File

@ -2,6 +2,7 @@ subroutine psbn_d_spcnv(a,b,info,type,mold,upd,dupl)
use psbn_d_mat_mod, psb_protect_name => psbn_d_spcnv
use psb_realloc_mod
use psb_sort_mod
implicit none
type(psbn_d_sparse_mat), intent(in) :: a
type(psbn_d_sparse_mat), intent(out) :: b
integer, intent(out) :: info
@ -15,6 +16,7 @@ subroutine psbn_d_spcnv_ip(a,info,type,mold,dupl)
use psbn_d_mat_mod, psb_protect_name => psbn_d_spcnv_ip
use psb_realloc_mod
use psb_sort_mod
implicit none
type(psbn_d_sparse_mat), intent(inout) :: a
integer, intent(out) :: info

37
base/newserial/psbn_mat_mod.f03
Unescape Escape View File

@ -77,6 +77,7 @@ contains
function get_dupl(a) result(res)
use psb_error_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
integer :: res
@ -89,6 +90,7 @@ contains
function get_state(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
integer :: res
@ -100,6 +102,7 @@ contains
end function get_state
function get_nrows(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
integer :: res
@ -112,6 +115,7 @@ contains
end function get_nrows
function get_ncols(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
integer :: res
@ -124,6 +128,7 @@ contains
end function get_ncols
function is_triangle(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
logical :: res
@ -136,6 +141,7 @@ contains
end function is_triangle
function is_unit(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
logical :: res
@ -148,6 +154,7 @@ contains
end function is_unit
function is_upper(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
logical :: res
@ -160,6 +167,7 @@ contains
end function is_upper
function is_lower(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
logical :: res
@ -172,6 +180,7 @@ contains
end function is_lower
function is_null(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
logical :: res
@ -184,6 +193,7 @@ contains
end function is_null
function is_bld(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
logical :: res
@ -196,6 +206,7 @@ contains
end function is_bld
function is_upd(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
logical :: res
@ -208,6 +219,7 @@ contains
end function is_upd
function is_asb(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
logical :: res
@ -220,6 +232,7 @@ contains
end function is_asb
function is_sorted(a) result(res)
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
logical :: res
@ -234,10 +247,11 @@ contains
function get_nzeros(a) result(res)
use psb_error_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
integer :: res
Integer :: err_act
Integer :: err_act, info
character(len=20) :: name='get_nzeros'
logical, parameter :: debug=.false.
@ -265,10 +279,11 @@ contains
function get_size(a) result(res)
use psb_error_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
integer :: res
Integer :: err_act
Integer :: err_act, info
character(len=20) :: name='get_size'
logical, parameter :: debug=.false.
@ -298,6 +313,7 @@ contains
subroutine get_neigh(a,idx,neigh,n,info,lev)
use psb_error_mod
implicit none
class(psbn_d_sparse_mat), intent(in) :: a
integer, intent(in) :: idx
integer, intent(out) :: n
@ -309,6 +325,7 @@ contains
character(len=20) :: name='get_neigh'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
if (.not.allocated(a%a)) then
info = 1121
@ -337,6 +354,7 @@ contains
subroutine allocate_mn(m,n,a,type,mold)
use psb_error_mod
use psb_string_mod
implicit none
integer, intent(in) :: m,n
class(psbn_d_sparse_mat), intent(inout) :: a
character(len=*), intent(in), optional :: type
@ -405,6 +423,7 @@ contains
subroutine allocate_mnnz(m,n,nz,a,type,mold)
use psb_error_mod
use psb_string_mod
implicit none
integer, intent(in) :: m,n,nz
class(psbn_d_sparse_mat), intent(inout) :: a
character(len=*), intent(in), optional :: type
@ -472,9 +491,10 @@ contains
subroutine reallocate_nz(nz,a)
use psb_error_mod
implicit none
integer, intent(in) :: nz
class(psbn_d_sparse_mat), intent(inout) :: a
Integer :: err_act
Integer :: err_act, info
character(len=20) :: name='reallocate_nz'
logical, parameter :: debug=.false.
@ -504,8 +524,9 @@ contains
subroutine free(a)
use psb_error_mod
implicit none
class(psbn_d_sparse_mat), intent(inout) :: a
Integer :: err_act
Integer :: err_act, info
character(len=20) :: name='free'
logical, parameter :: debug=.false.
@ -535,6 +556,7 @@ contains
subroutine d_csmm(alpha,a,x,beta,y,info,trans)
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(:,:)
@ -544,6 +566,7 @@ contains
character(len=20) :: name='psbn_csmm'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
if (.not.allocated(a%a)) then
info = 1121
@ -569,6 +592,7 @@ contains
subroutine d_csmv(alpha,a,x,beta,y,info,trans)
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(:)
@ -578,6 +602,7 @@ contains
character(len=20) :: name='psbn_csmv'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
if (.not.allocated(a%a)) then
info = 1121
@ -603,6 +628,7 @@ contains
subroutine d_cssm(alpha,a,x,beta,y,info,trans)
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(:,:)
@ -612,6 +638,7 @@ contains
character(len=20) :: name='psbn_cssm'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
if (.not.allocated(a%a)) then
info = 1121
@ -637,6 +664,7 @@ contains
subroutine d_cssv(alpha,a,x,beta,y,info,trans)
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(:)
@ -646,6 +674,7 @@ contains
character(len=20) :: name='psbn_cssv'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
if (.not.allocated(a%a)) then
info = 1121

2
base/serial/csr/dcsrsm.f
Unescape Escape View File

@ -60,7 +60,7 @@ C .. Local Arrays ..
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
END IF
IF (psb_toupper(DESCRA(3:3)).EQ.'N') DIAG = 'N'
DIAG = 'N'
IF (psb_toupper(DESCRA(3:3)).EQ.'U') DIAG = 'U'
IF(UNITD.EQ.'B') THEN
IERROR=5

2
test/pargen/runs/ppde.inp
Unescape Escape View File

@ -2,7 +2,7 @@
BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL RGMRES
BJAC Preconditioner NONE DIAG BJAC
CSR Storage format for matrix A: CSR COO JAD
080 Domain size (acutal system is this**3)
060 Domain size (acutal system is this**3)
1 Stopping criterion
200 MAXIT
-1 ITRACE

44
test/serial/Makefile
Unescape Escape View File

@ -0,0 +1,44 @@
include ../../Make.inc
#
# Libraries used
#
LIBDIR=../../lib/
PSBLAS_LIB= -L$(LIBDIR) -lpsb_util -lpsb_krylov -lpsb_prec -lpsb_base
LDLIBS=$(PSBLDLIBS)
#
# Compilers and such
#
CCOPT= -g
FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG).
EXEDIR=./runs
all: d_coo_matgen
d_coo_matgen: d_coo_matgen.o
$(F90LINK) d_coo_matgen.o -o d_coo_matgen $(PSBLAS_LIB) $(LDLIBS)
/bin/mv d_coo_matgen $(EXEDIR)
#ppde spde
ppde: ppde.o
$(F90LINK) ppde.o -o ppde $(PSBLAS_LIB) $(LDLIBS)
/bin/mv ppde $(EXEDIR)
spde: spde.o
$(F90LINK) spde.o -o spde $(PSBLAS_LIB) $(LDLIBS)
/bin/mv spde $(EXEDIR)
.f90.o:
$(MPF90) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $<
clean:
/bin/rm -f d_coo_matgen.o tpg.o ppde.o spde.o $(EXEDIR)/ppde
verycleanlib:
(cd ../..; make veryclean)
lib:
(cd ../../; make library)

475
test/serial/d_coo_matgen.f03
Unescape Escape View File

@ -0,0 +1,475 @@
!
program d_coo_matgen
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod
use psbn_d_base_mat_mod
use psbn_d_csr_mat_mod
implicit none
! input parameters
character(len=20) :: kmethd, ptype
character(len=5) :: afmt
integer :: idim
! miscellaneous
real(psb_dpk_), parameter :: one = 1.d0
real(psb_dpk_) :: t1, t2, tprec
! sparse matrix and preconditioner
type(psb_dspmat_type) :: a
type(psb_dprec_type) :: prec
! descriptor
type(psb_desc_type) :: desc_a
! dense matrices
real(psb_dpk_), allocatable :: b(:), x(:)
! blacs parameters
integer :: ictxt, iam, np
! solver parameters
integer :: iter, itmax,itrace, istopc, irst
integer(psb_long_int_k_) :: amatsize, precsize, descsize
real(psb_dpk_) :: err, eps
! other variables
integer :: info, err_act
character(len=20) :: name,ch_err
info=0
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ictxt)
stop
endif
if(psb_get_errstatus() /= 0) goto 9999
call psb_set_errverbosity(2)
!
! get parameters
!
call get_parms(ictxt,idim)
!
! allocate and fill in the coefficient matrix, rhs and initial guess
!
call psb_barrier(ictxt)
t1 = psb_wtime()
call create_matrix(idim,a,b,x,desc_a,ictxt,afmt,info)
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
if(info /= 0) then
call psb_error(ictxt)
end if
call psb_exit(ictxt)
stop
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
end if
stop
contains
!
! get iteration parameters from standard input
!
subroutine get_parms(ictxt,idim)
integer :: ictxt
integer :: idim
integer :: np, iam
integer :: intbuf(10), ip
call psb_info(ictxt, iam, np)
read(*,*) idim
return
end subroutine get_parms
!
! print an error message
!
subroutine pr_usage(iout)
integer :: iout
write(iout,*)'incorrect parameter(s) found'
write(iout,*)' usage: pde90 methd prec dim &
&[istop itmax itrace]'
write(iout,*)' where:'
write(iout,*)' methd: cgstab cgs rgmres bicgstabl'
write(iout,*)' prec : bjac diag none'
write(iout,*)' dim number of points along each axis'
write(iout,*)' the size of the resulting linear '
write(iout,*)' system is dim**3'
write(iout,*)' istop stopping criterion 1, 2 '
write(iout,*)' itmax maximum number of iterations [500] '
write(iout,*)' itrace <=0 (no tracing, default) or '
write(iout,*)' >= 1 do tracing every itrace'
write(iout,*)' iterations '
end subroutine pr_usage
!
! subroutine to allocate and fill in the coefficient matrix and
! the rhs.
!
subroutine create_matrix(idim,a,b,xv,desc_a,ictxt,afmt,info)
!
! discretize the partial diferential equation
!
! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u)
! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u
! dxdx dydy dzdz dx dy dz
!
! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1.
!
! Boundary conditions are set in a very simple way, by adding
! equations of the form
!
! u(x,y) = exp(-x^2-y^2-z^2)
!
! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation.
!
use psb_base_mod
implicit none
integer :: idim
integer, parameter :: nb=20
real(psb_dpk_), allocatable :: b(:),xv(:)
type(psb_desc_type) :: desc_a
integer :: ictxt, info
character :: afmt*5
type(psb_dspmat_type) :: a
real(psb_dpk_) :: zt(nb),glob_x,glob_y,glob_z
integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k
integer :: x,y,z,ia,indx_owner
integer :: np, iam, nr, nt,nz,isz
integer :: element
integer, allocatable :: irow(:),icol(:),myidx(:)
real(psb_dpk_), allocatable :: val(:)
type(psbn_d_coo_sparse_mat) :: acoo
type(psbn_d_csr_sparse_mat) :: acsr
! deltah dimension of each grid cell
! deltat discretization time
real(psb_dpk_) :: deltah
real(psb_dpk_),parameter :: rhs=0.d0,one=1.d0,zero=0.d0
real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcpy, tmov
real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3
external :: a1, a2, a3, a4, b1, b2, b3
integer :: err_act
character(len=20) :: name, ch_err
info = 0
name = 'create_matrix'
call psb_erractionsave(err_act)
call psb_info(ictxt, iam, np)
deltah = 1.d0/(idim-1)
! initialize array descriptor and sparse matrix storage. provide an
! estimate of the number of non zeroes
m = idim*idim*idim
n = m
nnz = ((n*9)/(np))
if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0,")...")')n
!
! Using a simple BLOCK distribution.
!
nt = (m+np-1)/np
nr = max(0,min(nt,m-(iam*nt)))
nt = nr
call psb_sum(ictxt,nt)
if (nt /= m) write(0,*) iam, 'Initialization error ',nr,nt,m
write(0,*) iam, 'Initialization ',nr,nt,m
nlr = nt
call psb_barrier(ictxt)
t0 = psb_wtime()
call acoo%allocate(nr,nr)
talc = psb_wtime()-t0
!!$ write(*,*) 'Test get size:',d_coo_get_size(acoo)
!!$ write(*,*) 'Test 2 get size:',acoo%get_size(),acoo%get_nzeros()
if (info /= 0) then
info=4010
ch_err='allocation rout.'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
! we build an auxiliary matrix consisting of one row at a
! time; just a small matrix. might be extended to generate
! a bunch of rows per call.
!
allocate(val(20*nb),irow(20*nb),&
&icol(20*nb),myidx(nlr),stat=info)
if (info /= 0 ) then
info=4000
call psb_errpush(info,name)
goto 9999
endif
! loop over rows belonging to current process in a block
! distribution.
call psb_barrier(ictxt)
t1 = psb_wtime()
do ii=1, nlr,nb
ib = min(nb,nlr-ii+1)
!!$ write(0,*) 'Row ',ii,ib
element = 1
do k=1,ib
i=ii+k-1
! local matrix pointer
glob_row=i
! compute gridpoint coordinates
if (mod(glob_row,(idim*idim)) == 0) then
x = glob_row/(idim*idim)
else
x = glob_row/(idim*idim)+1
endif
if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
y = (glob_row-(x-1)*idim*idim)/idim
else
y = (glob_row-(x-1)*idim*idim)/idim+1
endif
z = glob_row-(x-1)*idim*idim-(y-1)*idim
! glob_x, glob_y, glob_x coordinates
glob_x=x*deltah
glob_y=y*deltah
glob_z=z*deltah
! check on boundary points
zt(k) = 0.d0
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
if (x==1) then
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(k) = exp(-glob_y**2-glob_z**2)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)&
& -a1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element) = (x-2)*idim*idim+(y-1)*idim+(z)
irow(element) = glob_row
element = element+1
endif
! term depending on (x,y-1,z)
if (y==1) then
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(k) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
else
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element) = (x-1)*idim*idim+(y-2)*idim+(z)
irow(element) = glob_row
element = element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(k) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
else
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element) = (x-1)*idim*idim+(y-1)*idim+(z-1)
irow(element) = glob_row
element = element+1
endif
! term depending on (x,y,z)
val(element)=2*b1(glob_x,glob_y,glob_z)&
& +2*b2(glob_x,glob_y,glob_z)&
& +2*b3(glob_x,glob_y,glob_z)&
& +a1(glob_x,glob_y,glob_z)&
& +a2(glob_x,glob_y,glob_z)&
& +a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element) = (x-1)*idim*idim+(y-1)*idim+(z)
irow(element) = glob_row
element = element+1
! term depending on (x,y,z+1)
if (z==idim) then
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(k) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
else
val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element) = (x-1)*idim*idim+(y-1)*idim+(z+1)
irow(element) = glob_row
element = element+1
endif
! term depending on (x,y+1,z)
if (y==idim) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(k) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
else
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element) = (x-1)*idim*idim+(y)*idim+(z)
irow(element) = glob_row
element = element+1
endif
! term depending on (x+1,y,z)
if (x<idim) then
val(element)=-b3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element) = (x)*idim*idim+(y-1)*idim+(z)
irow(element) = glob_row
element = element+1
endif
end do
call acoo%csins(element-1,val,irow,icol,1,nr,1,nr,info)
end do
tgen = psb_wtime()-t1
if(info /= 0) then
info=4010
ch_err='insert rout.'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
!!$ call acoo%print(19)
t1 = psb_wtime()
!!$ write(0,*) 'out of loop ',acoo%get_nzeros()
call acoo%fix(info)
!!$ write(0,*) '2 out of loop ',acoo%get_nzeros()
if(info /= 0) then
info=4010
ch_err='asb rout.'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
tasb = psb_wtime()-t1
!!$ call acoo%print(20)
t1 = psb_wtime()
call acsr%cp_from_coo(acoo,info)
if(info /= 0) then
info=4010
ch_err='cp rout.'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
tcpy = psb_wtime()-t1
!!$ call acsr%print(21)
t1 = psb_wtime()
call acsr%mv_from_coo(acoo,info)
if(info /= 0) then
info=4010
ch_err='mv rout.'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
tmov = psb_wtime()-t1
!!$ call acsr%print(22)
if(iam == psb_root_) then
write(*,'("The matrix has been generated and assembled in ",a3," format.")')&
& a%fida(1:3)
write(*,'("-allocation time : ",es12.5)') talc
write(*,'("-coeff. gen. time : ",es12.5)') tgen
write(*,'("-assembly time : ",es12.5)') tasb
write(*,'("-copy time : ",es12.5)') tcpy
write(*,'("-move time : ",es12.5)') tmov
!!$ write(*,'("-total time : ",es12.5)') ttot
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
return
end subroutine create_matrix
end program d_coo_matgen
!
! functions parametrizing the differential equation
!
function a1(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a1
real(psb_dpk_) :: x,y,z
a1=1.d0
end function a1
function a2(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a2
real(psb_dpk_) :: x,y,z
a2=2.d1*y
end function a2
function a3(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a3
real(psb_dpk_) :: x,y,z
a3=1.d0
end function a3
function a4(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a4
real(psb_dpk_) :: x,y,z
a4=1.d0
end function a4
function b1(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b1
real(psb_dpk_) :: x,y,z
b1=1.d0
end function b1
function b2(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b2
real(psb_dpk_) :: x,y,z
b2=1.d0
end function b2
function b3(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b3
real(psb_dpk_) :: x,y,z
b3=1.d0
end function b3
Write Preview
Loading…
Cancel
Save