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/serial/f03/psb_d_coo_impl.f03
 base/serial/f03/psb_d_csr_impl.f03
 base/serial/f03/psb_s_coo_impl.f03
 base/serial/f03/psb_s_csr_impl.f03

Fixed various implementation details for SV/SM.
psblas3-type-indexed
Salvatore Filippone 16 years ago
parent e9579cc9a4
commit 3a69bef4f6
4 changed files with 286 additions and 271 deletions
Show Stats Download Patch File Download Diff File

273
base/serial/f03/psb_d_coo_impl.f03
Unescape Escape View File

@ -44,7 +44,7 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C') tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C')
m = a%get_nrows() m = a%get_nrows()
nc = min(size(x,2) , size(y,2)) nc = min(size(x,2) , size(y,2))
nnz = a%get_nzeros()
if (alpha == dzero) then if (alpha == dzero) then
if (beta == dzero) then if (beta == dzero) then
@ -60,7 +60,9 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
end if end if
if (beta == dzero) then if (beta == dzero) then
call inner_coosm(tra,a,x(:,1:nc),y(:,1:nc),info) call inner_coosm(tra,a%is_lower(),a%is_unit(),a%is_sorted(),&
& m,nc,nnz,a%ia,a%ja,a%val,&
& x,size(x,1),y,size(y,1),info)
do i = 1, m do i = 1, m
y(i,1:nc) = alpha*y(i,1:nc) y(i,1:nc) = alpha*y(i,1:nc)
end do end do
@ -72,8 +74,9 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
goto 9999 goto 9999
end if end if
tmp(1:m,1:nc) = x(1:m,1:nc) call inner_coosm(tra,a%is_lower(),a%is_unit(),a%is_sorted(),&
call inner_coosm(tra,a,tmp(:,1:nc),y(:,1:nc),info) & m,nc,nnz,a%ia,a%ja,a%val,&
& x,size(x,1),tmp,size(tmp,1),info)
do i = 1, m do i = 1, m
y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc)
end do end do
@ -101,92 +104,94 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
contains contains
subroutine inner_coosm(tra,a,x,y,info) subroutine inner_coosm(tra,lower,unit,sorted,nr,nc,nz,&
& ia,ja,val,x,ldx,y,ldy,info)
implicit none implicit none
logical, intent(in) :: tra logical, intent(in) :: tra,lower,unit,sorted
class(psb_d_coo_sparse_mat), intent(in) :: a integer, intent(in) :: nr,nc,nz,ldx,ldy,ia(*),ja(*)
real(psb_dpk_), intent(in) :: x(:,:) real(psb_dpk_), intent(in) :: val(*), x(ldx,*)
real(psb_dpk_), intent(out) :: y(:,:) real(psb_dpk_), intent(out) :: y(ldy,*)
integer, intent(out) :: info integer, intent(out) :: info
integer :: i,j,k,m, ir, jc integer :: i,j,k,m, ir, jc
real(psb_dpk_), allocatable :: acc(:) real(psb_dpk_), allocatable :: acc(:)
info = 0 info = 0
allocate(acc(size(x,2)), stat=info) allocate(acc(nc), stat=info)
if(info /= 0) then if(info /= 0) then
info=4010 info=4010
return return
end if end if
if (.not.a%is_sorted()) then if (.not.sorted) then
info = 1121 info = 1121
return return
end if end if
nnz = a%get_nzeros() nnz = nz
if (.not.tra) then if (.not.tra) then
if (a%is_lower()) then if (lower) then
if (a%is_unit()) then if (unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = dzero acc(1:nc) = dzero
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
acc = acc + a%val(j)*y(a%ja(j),:) acc(1:nc) = acc(1:nc) + val(j)*y(ja(j),1:nc)
j = j + 1 j = j + 1
end do end do
y(i,:) = x(i,:) - acc y(i,1:nc) = x(i,1:nc) - acc(1:nc)
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = dzero acc(1:nc) = dzero
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
if (a%ja(j) == i) then if (ja(j) == i) then
y(i,:) = (x(i,:) - acc)/a%val(j) y(i,1:nc) = (x(i,1:nc) - acc(1:nc))/val(j)
j = j + 1 j = j + 1
exit exit
end if end if
acc = acc + a%val(j)*y(a%ja(j),:) acc(1:nc) = acc(1:nc) + val(j)*y(ja(j),1:nc)
j = j + 1 j = j + 1
end do end do
end do end do
end if end if
else if (a%is_upper()) then else if (.not.lower) then
if (a%is_unit()) then if (unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = dzero acc(1:nc) = dzero
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
acc = acc + a%val(j)*x(a%ja(j),:) acc(1:nc) = acc(1:nc) + val(j)*x(ja(j),1:nc)
j = j - 1 j = j - 1
end do end do
y(i,:) = x(i,:) - acc y(i,1:nc) = x(i,1:nc) - acc(1:nc)
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = dzero acc(1:nc) = dzero
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
if (a%ja(j) == i) then if (ja(j) == i) then
y(i,:) = (x(i,:) - acc)/a%val(j) y(i,1:nc) = (x(i,1:nc) - acc(1:nc))/val(j)
j = j - 1 j = j - 1
exit exit
end if end if
acc = acc + a%val(j)*y(a%ja(j),:) acc(1:nc) = acc(1:nc) + val(j)*y(ja(j),1:nc)
j = j - 1 j = j - 1
end do end do
end do end do
@ -196,66 +201,66 @@ contains
else if (tra) then else if (tra) then
do i=1, a%get_nrows() do i=1, nr
y(i,:) = x(i,:) y(i,1:nc) = x(i,1:nc)
end do end do
if (a%is_lower()) then if (lower) then
if (a%is_unit()) then if (unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = y(i,:) acc(1:nc) = y(i,1:nc)
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
jc = a%ja(j) jc = ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc y(jc,1:nc) = y(jc,1:nc) - val(j)*acc(1:nc)
j = j - 1 j = j - 1
end do end do
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
if (a%ja(j) == i) then if (ja(j) == i) then
y(i,:) = y(i,:) /a%val(j) y(i,1:nc) = y(i,1:nc) /val(j)
j = j - 1 j = j - 1
end if end if
acc = y(i,:) acc(1:nc) = y(i,1:nc)
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
jc = a%ja(j) jc = ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc y(jc,1:nc) = y(jc,1:nc) - val(j)*acc(1:nc)
j = j - 1 j = j - 1
end do end do
end do end do
else if (a%is_upper()) then else if (.not.lower) then
if (a%is_unit()) then if (unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = y(i,:) acc(1:nc) = y(i,1:nc)
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
jc = a%ja(j) jc = ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc y(jc,1:nc) = y(jc,1:nc) - val(j)*acc(1:nc)
j = j + 1 j = j + 1
end do end do
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
if (a%ja(j) == i) then if (ja(j) == i) then
y(i,:) = y(i,:) /a%val(j) y(i,1:nc) = y(i,1:nc) /val(j)
j = j + 1 j = j + 1
end if end if
acc = y(i,:) acc(1:nc) = y(i,1:nc)
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
jc = a%ja(j) jc = ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc y(jc,1:nc) = y(jc,1:nc) - val(j)*acc(1:nc)
j = j + 1 j = j + 1
end do end do
end do end do
@ -328,7 +333,9 @@ subroutine d_coo_cssv_impl(alpha,a,x,beta,y,info,trans)
end if end if
if (beta == dzero) then if (beta == dzero) then
call inner_coosv(tra,a,x,y,info) call inner_coosv(tra,a%is_lower(),a%is_unit(),a%is_sorted(),&
& a%get_nrows(),a%get_nzeros(),a%ia,a%ja,a%val,&
& x,y,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
@ -344,8 +351,9 @@ subroutine d_coo_cssv_impl(alpha,a,x,beta,y,info,trans)
goto 9999 goto 9999
end if end if
tmp(1:m) = x(1:m) call inner_coosv(tra,a%is_lower(),a%is_unit(),a%is_sorted(),&
call inner_coosv(tra,a,tmp,y,info) & a%get_nrows(),a%get_nzeros(),a%ia,a%ja,a%val,&
& x,tmp,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
@ -369,86 +377,87 @@ subroutine d_coo_cssv_impl(alpha,a,x,beta,y,info,trans)
contains contains
subroutine inner_coosv(tra,a,x,y,info) subroutine inner_coosv(tra,lower,unit,sorted,nr,nz,&
& ia,ja,val,x,y,info)
implicit none implicit none
logical, intent(in) :: tra logical, intent(in) :: tra,lower,unit,sorted
class(psb_d_coo_sparse_mat), intent(in) :: a integer, intent(in) :: nr,nz,ia(*),ja(*)
real(psb_dpk_), intent(in) :: x(:) real(psb_dpk_), intent(in) :: val(*), x(*)
real(psb_dpk_), intent(out) :: y(:) real(psb_dpk_), intent(out) :: y(*)
integer, intent(out) :: info integer, intent(out) :: info
integer :: i,j,k,m, ir, jc, nnz integer :: i,j,k,m, ir, jc, nnz
real(psb_dpk_) :: acc real(psb_dpk_) :: acc
info = 0 info = 0
if (.not.a%is_sorted()) then if (.not.sorted) then
info = 1121 info = 1121
return return
end if end if
nnz = a%get_nzeros() nnz = nz
if (.not.tra) then if (.not.tra) then
if (a%is_lower()) then if (lower) then
if (a%is_unit()) then if (unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = dzero acc = dzero
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
acc = acc + a%val(j)*y(a%ja(j)) acc = acc + val(j)*y(ja(j))
j = j + 1 j = j + 1
end do end do
y(i) = x(i) - acc y(i) = x(i) - acc
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = dzero acc = dzero
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
if (a%ja(j) == i) then if (ja(j) == i) then
y(i) = (x(i) - acc)/a%val(j) y(i) = (x(i) - acc)/val(j)
j = j + 1 j = j + 1
exit exit
end if end if
acc = acc + a%val(j)*y(a%ja(j)) acc = acc + val(j)*y(ja(j))
j = j + 1 j = j + 1
end do end do
end do end do
end if end if
else if (a%is_upper()) then else if (.not.lower) then
if (a%is_unit()) then if (unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = dzero acc = dzero
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
acc = acc + a%val(j)*y(a%ja(j)) acc = acc + val(j)*y(ja(j))
j = j - 1 j = j - 1
end do end do
y(i) = x(i) - acc y(i) = x(i) - acc
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = dzero acc = dzero
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
if (a%ja(j) == i) then if (ja(j) == i) then
y(i) = (x(i) - acc)/a%val(j) y(i) = (x(i) - acc)/val(j)
j = j - 1 j = j - 1
exit exit
end if end if
acc = acc + a%val(j)*y(a%ja(j)) acc = acc + val(j)*y(ja(j))
j = j - 1 j = j - 1
end do end do
end do end do
@ -458,66 +467,66 @@ contains
else if (tra) then else if (tra) then
do i=1, a%get_nrows() do i=1, nr
y(i) = x(i) y(i) = x(i)
end do end do
if (a%is_lower()) then if (lower) then
if (a%is_unit()) then if (unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = y(i) acc = y(i)
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
jc = a%ja(j) jc = ja(j)
y(jc) = y(jc) - a%val(j)*acc y(jc) = y(jc) - val(j)*acc
j = j - 1 j = j - 1
end do end do
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
if (a%ja(j) == i) then if (ja(j) == i) then
y(i) = y(i) /a%val(j) y(i) = y(i) /val(j)
j = j - 1 j = j - 1
end if end if
acc = y(i) acc = y(i)
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
jc = a%ja(j) jc = ja(j)
y(jc) = y(jc) - a%val(j)*acc y(jc) = y(jc) - val(j)*acc
j = j - 1 j = j - 1
end do end do
end do end do
else if (a%is_upper()) then else if (.not.lower) then
if (a%is_unit()) then if (unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = y(i) acc = y(i)
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
jc = a%ja(j) jc = ja(j)
y(jc) = y(jc) - a%val(j)*acc y(jc) = y(jc) - val(j)*acc
j = j + 1 j = j + 1
end do end do
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
if (a%ja(j) == i) then if (ja(j) == i) then
y(i) = y(i) /a%val(j) y(i) = y(i) /val(j)
j = j + 1 j = j + 1
end if end if
acc = y(i) acc = y(i)
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
jc = a%ja(j) jc = ja(j)
y(jc) = y(jc) - a%val(j)*acc y(jc) = y(jc) - val(j)*acc
j = j + 1 j = j + 1
end do end do
end do end do

7
base/serial/f03/psb_d_csr_impl.f03
Unescape Escape View File

@ -662,9 +662,9 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
if (info /= 0) then if (info /= 0) then
return return
end if end if
tmp(1:m) = x(1:m)
call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),& call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),&
& a%irp,a%ja,a%val,tmp,y) & a%irp,a%ja,a%val,x,tmp)
do i = 1, m do i = 1, m
y(i) = alpha*tmp(i) + beta*y(i) y(i) = alpha*tmp(i) + beta*y(i)
end do end do
@ -863,9 +863,8 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
goto 9999 goto 9999
end if end if
tmp(1:m,:) = x(1:m,1:nc)
call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,&
& a%irp,a%ja,a%val,tmp,size(tmp,1),y,size(y,1),info) & a%irp,a%ja,a%val,x,size(x,1),tmp,size(tmp,1),info)
do i = 1, m do i = 1, m
y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc)
end do end do

271
base/serial/f03/psb_s_coo_impl.f03
Unescape Escape View File

@ -44,6 +44,7 @@ subroutine s_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C') tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C')
m = a%get_nrows() m = a%get_nrows()
nc = min(size(x,2) , size(y,2)) nc = min(size(x,2) , size(y,2))
nnz = a%get_nzeros()
if (alpha == szero) then if (alpha == szero) then
@ -60,7 +61,9 @@ subroutine s_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
end if end if
if (beta == szero) then if (beta == szero) then
call inner_coosm(tra,a,x(:,1:nc),y(:,1:nc),info) call inner_coosm(tra,a%is_lower(),a%is_unit(),a%is_sorted(),&
& m,nc,nnz,a%ia,a%ja,a%val,&
& x,size(x,1),y,size(y,1),info)
do i = 1, m do i = 1, m
y(i,1:nc) = alpha*y(i,1:nc) y(i,1:nc) = alpha*y(i,1:nc)
end do end do
@ -72,8 +75,9 @@ subroutine s_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
goto 9999 goto 9999
end if end if
tmp(1:m,1:nc) = x(1:m,1:nc) call inner_coosm(tra,a%is_lower(),a%is_unit(),a%is_sorted(),&
call inner_coosm(tra,a,tmp(:,1:nc),y(:,1:nc),info) & m,nc,nnz,a%ia,a%ja,a%val,&
& x,size(x,1),tmp,size(tmp,1),info)
do i = 1, m do i = 1, m
y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc)
end do end do
@ -101,92 +105,93 @@ subroutine s_coo_cssm_impl(alpha,a,x,beta,y,info,trans)
contains contains
subroutine inner_coosm(tra,a,x,y,info) subroutine inner_coosm(tra,lower,unit,sorted,nr,nc,nz,&
& ia,ja,val,x,ldx,y,ldy,info)
implicit none implicit none
logical, intent(in) :: tra logical, intent(in) :: tra,lower,unit,sorted
class(psb_s_coo_sparse_mat), intent(in) :: a integer, intent(in) :: nr,nc,nz,ldx,ldy,ia(*),ja(*)
real(psb_spk_), intent(in) :: x(:,:) real(psb_spk_), intent(in) :: val(*), x(ldx,*)
real(psb_spk_), intent(out) :: y(:,:) real(psb_spk_), intent(out) :: y(ldy,*)
integer, intent(out) :: info integer, intent(out) :: info
integer :: i,j,k,m, ir, jc integer :: i,j,k,m, ir, jc
real(psb_spk_), allocatable :: acc(:) real(psb_spk_), allocatable :: acc(:)
info = 0 info = 0
allocate(acc(size(x,2)), stat=info) allocate(acc(nc), stat=info)
if(info /= 0) then if(info /= 0) then
info=4010 info=4010
return return
end if end if
if (.not.a%is_sorted()) then if (.not.sorted) then
info = 1121 info = 1121
return return
end if end if
nnz = a%get_nzeros() nnz = nz
if (.not.tra) then if (.not.tra) then
if (a%is_lower()) then if (lower) then
if (a%is_unit()) then if (unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = szero acc(1:nc) = szero
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
acc = acc + a%val(j)*y(a%ja(j),:) acc(1:nc) = acc(1:nc) + val(j)*y(ja(j),1:nc)
j = j + 1 j = j + 1
end do end do
y(i,:) = x(i,:) - acc y(i,1:nc) = x(i,1:nc) - acc(1:nc)
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = szero acc(1:nc) = szero
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
if (a%ja(j) == i) then if (ja(j) == i) then
y(i,:) = (x(i,:) - acc)/a%val(j) y(i,1:nc) = (x(i,1:nc) - acc(1:nc))/val(j)
j = j + 1 j = j + 1
exit exit
end if end if
acc = acc + a%val(j)*y(a%ja(j),:) acc(1:nc) = acc(1:nc) + val(j)*y(ja(j),1:nc)
j = j + 1 j = j + 1
end do end do
end do end do
end if end if
else if (a%is_upper()) then else if (.not.lower) then
if (a%is_unit()) then if (unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = szero acc(1:nc) = szero
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
acc = acc + a%val(j)*x(a%ja(j),:) acc(1:nc) = acc(1:nc) + val(j)*x(ja(j),1:nc)
j = j - 1 j = j - 1
end do end do
y(i,:) = x(i,:) - acc y(i,1:nc) = x(i,1:nc) - acc(1:nc)
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = szero acc(1:nc) = szero
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
if (a%ja(j) == i) then if (ja(j) == i) then
y(i,:) = (x(i,:) - acc)/a%val(j) y(i,1:nc) = (x(i,1:nc) - acc(1:nc))/val(j)
j = j - 1 j = j - 1
exit exit
end if end if
acc = acc + a%val(j)*y(a%ja(j),:) acc(1:nc) = acc(1:nc) + val(j)*y(ja(j),1:nc)
j = j - 1 j = j - 1
end do end do
end do end do
@ -196,66 +201,66 @@ contains
else if (tra) then else if (tra) then
do i=1, a%get_nrows() do i=1, nr
y(i,:) = x(i,:) y(i,1:nc) = x(i,1:nc)
end do end do
if (a%is_lower()) then if (lower) then
if (a%is_unit()) then if (unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = y(i,:) acc(1:nc) = y(i,1:nc)
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
jc = a%ja(j) jc = ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc y(jc,1:nc) = y(jc,1:nc) - val(j)*acc(1:nc)
j = j - 1 j = j - 1
end do end do
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
if (a%ja(j) == i) then if (ja(j) == i) then
y(i,:) = y(i,:) /a%val(j) y(i,1:nc) = y(i,1:nc) /val(j)
j = j - 1 j = j - 1
end if end if
acc = y(i,:) acc(1:nc) = y(i,1:nc)
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
jc = a%ja(j) jc = ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc y(jc,1:nc) = y(jc,1:nc) - val(j)*acc(1:nc)
j = j - 1 j = j - 1
end do end do
end do end do
else if (a%is_upper()) then else if (.not.lower) then
if (a%is_unit()) then if (unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = y(i,:) acc(1:nc) = y(i,1:nc)
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
jc = a%ja(j) jc = ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc y(jc,1:nc) = y(jc,1:nc) - val(j)*acc(1:nc)
j = j + 1 j = j + 1
end do end do
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
if (a%ja(j) == i) then if (ja(j) == i) then
y(i,:) = y(i,:) /a%val(j) y(i,1:nc) = y(i,1:nc) /val(j)
j = j + 1 j = j + 1
end if end if
acc = y(i,:) acc(1:nc) = y(i,1:nc)
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
jc = a%ja(j) jc = ja(j)
y(jc,:) = y(jc,:) - a%val(j)*acc y(jc,1:nc) = y(jc,1:nc) - val(j)*acc(1:nc)
j = j + 1 j = j + 1
end do end do
end do end do
@ -328,7 +333,9 @@ subroutine s_coo_cssv_impl(alpha,a,x,beta,y,info,trans)
end if end if
if (beta == szero) then if (beta == szero) then
call inner_coosv(tra,a,x,y,info) call inner_coosv(tra,a%is_lower(),a%is_unit(),a%is_sorted(),&
& a%get_nrows(),a%get_nzeros(),a%ia,a%ja,a%val,&
& x,y,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
@ -344,8 +351,9 @@ subroutine s_coo_cssv_impl(alpha,a,x,beta,y,info,trans)
goto 9999 goto 9999
end if end if
tmp(1:m) = x(1:m) call inner_coosv(tra,a%is_lower(),a%is_unit(),a%is_sorted(),&
call inner_coosv(tra,a,tmp,y,info) & a%get_nrows(),a%get_nzeros(),a%ia,a%ja,a%val,&
& x,tmp,info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
@ -369,86 +377,87 @@ subroutine s_coo_cssv_impl(alpha,a,x,beta,y,info,trans)
contains contains
subroutine inner_coosv(tra,a,x,y,info) subroutine inner_coosv(tra,lower,unit,sorted,nr,nz,&
& ia,ja,val,x,y,info)
implicit none implicit none
logical, intent(in) :: tra logical, intent(in) :: tra,lower,unit,sorted
class(psb_s_coo_sparse_mat), intent(in) :: a integer, intent(in) :: nr,nz,ia(*),ja(*)
real(psb_spk_), intent(in) :: x(:) real(psb_spk_), intent(in) :: val(*), x(*)
real(psb_spk_), intent(out) :: y(:) real(psb_spk_), intent(out) :: y(*)
integer, intent(out) :: info integer, intent(out) :: info
integer :: i,j,k,m, ir, jc, nnz integer :: i,j,k,m, ir, jc, nnz
real(psb_spk_) :: acc real(psb_spk_) :: acc
info = 0 info = 0
if (.not.a%is_sorted()) then if (.not.sorted) then
info = 1121 info = 1121
return return
end if end if
nnz = a%get_nzeros() nnz = nz
if (.not.tra) then if (.not.tra) then
if (a%is_lower()) then if (lower) then
if (a%is_unit()) then if (unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = szero acc = szero
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
acc = acc + a%val(j)*y(a%ja(j)) acc = acc + val(j)*y(ja(j))
j = j + 1 j = j + 1
end do end do
y(i) = x(i) - acc y(i) = x(i) - acc
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = szero acc = szero
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
if (a%ja(j) == i) then if (ja(j) == i) then
y(i) = (x(i) - acc)/a%val(j) y(i) = (x(i) - acc)/val(j)
j = j + 1 j = j + 1
exit exit
end if end if
acc = acc + a%val(j)*y(a%ja(j)) acc = acc + val(j)*y(ja(j))
j = j + 1 j = j + 1
end do end do
end do end do
end if end if
else if (a%is_upper()) then else if (.not.lower) then
if (a%is_unit()) then if (unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = szero acc = szero
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
acc = acc + a%val(j)*y(a%ja(j)) acc = acc + val(j)*y(ja(j))
j = j - 1 j = j - 1
end do end do
y(i) = x(i) - acc y(i) = x(i) - acc
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = szero acc = szero
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
if (a%ja(j) == i) then if (ja(j) == i) then
y(i) = (x(i) - acc)/a%val(j) y(i) = (x(i) - acc)/val(j)
j = j - 1 j = j - 1
exit exit
end if end if
acc = acc + a%val(j)*y(a%ja(j)) acc = acc + val(j)*y(ja(j))
j = j - 1 j = j - 1
end do end do
end do end do
@ -458,66 +467,66 @@ contains
else if (tra) then else if (tra) then
do i=1, a%get_nrows() do i=1, nr
y(i) = x(i) y(i) = x(i)
end do end do
if (a%is_lower()) then if (lower) then
if (a%is_unit()) then if (unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
acc = y(i) acc = y(i)
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
jc = a%ja(j) jc = ja(j)
y(jc) = y(jc) - a%val(j)*acc y(jc) = y(jc) - val(j)*acc
j = j - 1 j = j - 1
end do end do
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = nnz j = nnz
do i=a%get_nrows(), 1, -1 do i=nr, 1, -1
if (a%ja(j) == i) then if (ja(j) == i) then
y(i) = y(i) /a%val(j) y(i) = y(i) /val(j)
j = j - 1 j = j - 1
end if end if
acc = y(i) acc = y(i)
do do
if (j < 1) exit if (j < 1) exit
if (a%ia(j) < i) exit if (ia(j) < i) exit
jc = a%ja(j) jc = ja(j)
y(jc) = y(jc) - a%val(j)*acc y(jc) = y(jc) - val(j)*acc
j = j - 1 j = j - 1
end do end do
end do end do
else if (a%is_upper()) then else if (.not.lower) then
if (a%is_unit()) then if (unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
acc = y(i) acc = y(i)
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
jc = a%ja(j) jc = ja(j)
y(jc) = y(jc) - a%val(j)*acc y(jc) = y(jc) - val(j)*acc
j = j + 1 j = j + 1
end do end do
end do end do
else if (.not.a%is_unit()) then else if (.not.unit) then
j = 1 j = 1
do i=1, a%get_nrows() do i=1, nr
if (a%ja(j) == i) then if (ja(j) == i) then
y(i) = y(i) /a%val(j) y(i) = y(i) /val(j)
j = j + 1 j = j + 1
end if end if
acc = y(i) acc = y(i)
do do
if (j > nnz) exit if (j > nnz) exit
if (a%ia(j) > i) exit if (ia(j) > i) exit
jc = a%ja(j) jc = ja(j)
y(jc) = y(jc) - a%val(j)*acc y(jc) = y(jc) - val(j)*acc
j = j + 1 j = j + 1
end do end do
end do end do

6
base/serial/f03/psb_s_csr_impl.f03
Unescape Escape View File

@ -674,9 +674,8 @@ subroutine s_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
if (info /= 0) then if (info /= 0) then
return return
end if end if
tmp(1:m) = x(1:m)
call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),& call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),&
& a%irp,a%ja,a%val,tmp,y) & a%irp,a%ja,a%val,x,tmp)
do i = 1, m do i = 1, m
y(i) = alpha*tmp(i) + beta*y(i) y(i) = alpha*tmp(i) + beta*y(i)
end do end do
@ -875,9 +874,8 @@ subroutine s_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
goto 9999 goto 9999
end if end if
tmp(1:m,:) = x(1:m,1:nc)
call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,&
& a%irp,a%ja,a%val,tmp,size(tmp,1),y,size(y,1),info) & a%irp,a%ja,a%val,x,size(x,1),tmp,size(tmp,1),info)
do i = 1, m do i = 1, m
y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc)
end do end do

Write Preview
Loading…
Cancel
Save