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/modules/psb_c_base_mat_mod.f90
 base/modules/psb_d_base_mat_mod.f90
 base/modules/psb_s_base_mat_mod.f90
 base/modules/psb_z_base_mat_mod.f90
 base/serial/impl/psb_c_base_mat_impl.F90
 base/serial/impl/psb_c_coo_impl.f90
 base/serial/impl/psb_c_csc_impl.f90
 base/serial/impl/psb_c_csr_impl.f90
 base/serial/impl/psb_d_base_mat_impl.F90
 base/serial/impl/psb_d_coo_impl.f90
 base/serial/impl/psb_d_csc_impl.f90
 base/serial/impl/psb_d_csr_impl.f90
 base/serial/impl/psb_s_base_mat_impl.F90
 base/serial/impl/psb_s_coo_impl.f90
 base/serial/impl/psb_s_csc_impl.f90
 base/serial/impl/psb_s_csr_impl.f90
 base/serial/impl/psb_z_base_mat_impl.F90
 base/serial/impl/psb_z_coo_impl.f90
 base/serial/impl/psb_z_csc_impl.f90
 base/serial/impl/psb_z_csr_impl.f90

Fix SCAL and other IS_UNIT() usages.
psblas3-final
Salvatore Filippone 12 years ago
parent a922c70f3d
commit 5510d3242b
20 changed files with 728 additions and 352 deletions
Show Stats Download Patch File Download Diff File

18
base/modules/psb_c_base_mat_mod.f90
Unescape Escape View File

@ -73,6 +73,7 @@ module psb_c_base_mat_mod
procedure, pass(a) :: mv_from_fmt => psb_c_base_mv_from_fmt procedure, pass(a) :: mv_from_fmt => psb_c_base_mv_from_fmt
procedure, pass(a) :: mold => psb_c_base_mold procedure, pass(a) :: mold => psb_c_base_mold
procedure, pass(a) :: clone => psb_c_base_clone procedure, pass(a) :: clone => psb_c_base_clone
procedure, pass(a) :: add_unit_diag => psb_c_base_add_unit_diag
! !
! Transpose methods: defined here but not implemented. ! Transpose methods: defined here but not implemented.
@ -431,6 +432,23 @@ module psb_c_base_mat_mod
end interface end interface
!
!
!> Function add_unit_diag:
!! \memberof psb_c_base_add_unit_diag
!! \brief Given a matrix for which is_unit() is true, explicitly
!! store the unit diagonal and set is_unit() to false.
!! This is needed e.g. when scaling
!
interface
subroutine psb_c_base_add_unit_diag(a)
import :: psb_c_base_sparse_mat
implicit none
class(psb_c_base_sparse_mat), intent(inout) :: a
end subroutine psb_c_base_add_unit_diag
end interface
! !
!> Function cp_to_coo: !> Function cp_to_coo:
!! \memberof psb_c_base_sparse_mat !! \memberof psb_c_base_sparse_mat

18
base/modules/psb_d_base_mat_mod.f90
Unescape Escape View File

@ -73,6 +73,7 @@ module psb_d_base_mat_mod
procedure, pass(a) :: mv_from_fmt => psb_d_base_mv_from_fmt procedure, pass(a) :: mv_from_fmt => psb_d_base_mv_from_fmt
procedure, pass(a) :: mold => psb_d_base_mold procedure, pass(a) :: mold => psb_d_base_mold
procedure, pass(a) :: clone => psb_d_base_clone procedure, pass(a) :: clone => psb_d_base_clone
procedure, pass(a) :: add_unit_diag => psb_d_base_add_unit_diag
! !
! Transpose methods: defined here but not implemented. ! Transpose methods: defined here but not implemented.
@ -431,6 +432,23 @@ module psb_d_base_mat_mod
end interface end interface
!
!
!> Function add_unit_diag:
!! \memberof psb_d_base_add_unit_diag
!! \brief Given a matrix for which is_unit() is true, explicitly
!! store the unit diagonal and set is_unit() to false.
!! This is needed e.g. when scaling
!
interface
subroutine psb_d_base_add_unit_diag(a)
import :: psb_d_base_sparse_mat
implicit none
class(psb_d_base_sparse_mat), intent(inout) :: a
end subroutine psb_d_base_add_unit_diag
end interface
! !
!> Function cp_to_coo: !> Function cp_to_coo:
!! \memberof psb_d_base_sparse_mat !! \memberof psb_d_base_sparse_mat

18
base/modules/psb_s_base_mat_mod.f90
Unescape Escape View File

@ -73,6 +73,7 @@ module psb_s_base_mat_mod
procedure, pass(a) :: mv_from_fmt => psb_s_base_mv_from_fmt procedure, pass(a) :: mv_from_fmt => psb_s_base_mv_from_fmt
procedure, pass(a) :: mold => psb_s_base_mold procedure, pass(a) :: mold => psb_s_base_mold
procedure, pass(a) :: clone => psb_s_base_clone procedure, pass(a) :: clone => psb_s_base_clone
procedure, pass(a) :: add_unit_diag => psb_s_base_add_unit_diag
! !
! Transpose methods: defined here but not implemented. ! Transpose methods: defined here but not implemented.
@ -431,6 +432,23 @@ module psb_s_base_mat_mod
end interface end interface
!
!
!> Function add_unit_diag:
!! \memberof psb_s_base_add_unit_diag
!! \brief Given a matrix for which is_unit() is true, explicitly
!! store the unit diagonal and set is_unit() to false.
!! This is needed e.g. when scaling
!
interface
subroutine psb_s_base_add_unit_diag(a)
import :: psb_s_base_sparse_mat
implicit none
class(psb_s_base_sparse_mat), intent(inout) :: a
end subroutine psb_s_base_add_unit_diag
end interface
! !
!> Function cp_to_coo: !> Function cp_to_coo:
!! \memberof psb_s_base_sparse_mat !! \memberof psb_s_base_sparse_mat

18
base/modules/psb_z_base_mat_mod.f90
Unescape Escape View File

@ -73,6 +73,7 @@ module psb_z_base_mat_mod
procedure, pass(a) :: mv_from_fmt => psb_z_base_mv_from_fmt procedure, pass(a) :: mv_from_fmt => psb_z_base_mv_from_fmt
procedure, pass(a) :: mold => psb_z_base_mold procedure, pass(a) :: mold => psb_z_base_mold
procedure, pass(a) :: clone => psb_z_base_clone procedure, pass(a) :: clone => psb_z_base_clone
procedure, pass(a) :: add_unit_diag => psb_z_base_add_unit_diag
! !
! Transpose methods: defined here but not implemented. ! Transpose methods: defined here but not implemented.
@ -431,6 +432,23 @@ module psb_z_base_mat_mod
end interface end interface
!
!
!> Function add_unit_diag:
!! \memberof psb_z_base_add_unit_diag
!! \brief Given a matrix for which is_unit() is true, explicitly
!! store the unit diagonal and set is_unit() to false.
!! This is needed e.g. when scaling
!
interface
subroutine psb_z_base_add_unit_diag(a)
import :: psb_z_base_sparse_mat
implicit none
class(psb_z_base_sparse_mat), intent(inout) :: a
end subroutine psb_z_base_add_unit_diag
end interface
! !
!> Function cp_to_coo: !> Function cp_to_coo:
!! \memberof psb_z_base_sparse_mat !! \memberof psb_z_base_sparse_mat

31
base/serial/impl/psb_c_base_mat_impl.F90
Unescape Escape View File

@ -580,6 +580,37 @@ subroutine psb_c_base_clone(a,b,info)
end subroutine psb_c_base_clone end subroutine psb_c_base_clone
subroutine psb_c_base_add_unit_diag(a)
use psb_c_base_mat_mod, psb_protect_name => psb_c_base_add_unit_diag
use psb_error_mod
implicit none
class(psb_c_base_sparse_mat), intent(inout) :: a
type(psb_c_coo_sparse_mat) :: tmp
integer(psb_ipk_) :: i, j, m, n, nz, mnm, info
if (a%is_unit()) then
call a%mv_to_coo(tmp,info)
if (info /= 0) return
m = tmp%get_nrows()
n = tmp%get_ncols()
mnm = min(m,n)
nz = tmp%get_nzeros()
call tmp%reallocate(nz+mnm)
do i=1, mnm
tmp%val(nz+i) = cone
tmp%ia(nz+i) = i
tmp%ja(nz+i) = i
end do
call tmp%set_nzeros(nz+mnm)
call tmp%set_unit(.false.)
call tmp%fix(info)
if (info /= 0) &
& call a%mv_from_coo(tmp,info)
end if
end subroutine psb_c_base_add_unit_diag
subroutine psb_c_base_mold(a,b,info) subroutine psb_c_base_mold(a,b,info)
use psb_c_base_mat_mod, psb_protect_name => psb_c_base_mold use psb_c_base_mat_mod, psb_protect_name => psb_c_base_mold
use psb_error_mod use psb_error_mod

77
base/serial/impl/psb_c_coo_impl.f90
Unescape Escape View File

@ -54,11 +54,11 @@ subroutine psb_c_coo_get_diag(a,d,info)
call psb_errpush(info,name,i_err=ierr) call psb_errpush(info,name,i_err=ierr)
goto 9999 goto 9999
end if end if
d(:) = czero
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = cone d(1:mnm) = cone
else else
d(1:mnm) = czero
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
j=a%ia(i) j=a%ia(i)
if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then
@ -101,6 +101,10 @@ subroutine psb_c_coo_scal(d,a,info,side)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
side_ = 'L' side_ = 'L'
if (present(side)) then if (present(side)) then
side_ = psb_toupper(side) side_ = psb_toupper(side)
@ -167,6 +171,9 @@ subroutine psb_c_coo_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d
@ -1313,7 +1320,7 @@ subroutine psb_c_coo_csmv(alpha,a,x,beta,y,info,trans)
endif endif
return return
else else
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
if (beta == czero) then if (beta == czero) then
do i = 1, min(m,n) do i = 1, min(m,n)
y(i) = alpha*x(i) y(i) = alpha*x(i)
@ -1524,7 +1531,7 @@ subroutine psb_c_coo_csmm(alpha,a,x,beta,y,info,trans)
endif endif
return return
else else
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
if (beta == czero) then if (beta == czero) then
do i = 1, min(m,n) do i = 1, min(m,n)
y(i,1:nc) = alpha*x(i,1:nc) y(i,1:nc) = alpha*x(i,1:nc)
@ -1661,12 +1668,17 @@ function psb_c_coo_maxval(a) result(res)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
if (a%is_unit()) then
res = sone
else
res = szero res = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
if (allocated(a%val)) then if (allocated(a%val)) then
nnz = min(nnz,size(a%val)) nnz = min(nnz,size(a%val))
res = maxval(abs(a%val(1:nnz))) res = maxval(abs(a%val(1:nnz)))
end if end if
end function psb_c_coo_maxval end function psb_c_coo_maxval
function psb_c_coo_csnmi(a) result(res) function psb_c_coo_csnmi(a) result(res)
@ -1696,7 +1708,11 @@ function psb_c_coo_csnmi(a) result(res)
do while ((a%ia(j) == a%ia(i)).and. (j <= nnz)) do while ((a%ia(j) == a%ia(i)).and. (j <= nnz))
j = j+1 j = j+1
enddo enddo
if (a%is_unit()) then
acc = sone
else
acc = szero acc = szero
end if
do k=i, j-1 do k=i, j-1
acc = acc + abs(a%val(k)) acc = acc + abs(a%val(k))
end do end do
@ -1707,7 +1723,11 @@ function psb_c_coo_csnmi(a) result(res)
m = a%get_nrows() m = a%get_nrows()
allocate(vt(m),stat=info) allocate(vt(m),stat=info)
if (info /= 0) return if (info /= 0) return
vt(:) = szero if (a%is_unit()) then
vt = sone
else
vt = szero
end if
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
vt(i) = vt(i) + abs(a%val(j)) vt(i) = vt(i) + abs(a%val(j))
@ -1743,7 +1763,11 @@ function psb_c_coo_csnm1(a) result(res)
n = a%get_ncols() n = a%get_ncols()
allocate(vt(n),stat=info) allocate(vt(n),stat=info)
if (info /= 0) return if (info /= 0) return
vt(:) = szero if (a%is_unit()) then
vt = sone
else
vt = szero
end if
do j=1, nnz do j=1, nnz
i = a%ja(j) i = a%ja(j)
vt(i) = vt(i) + abs(a%val(j)) vt(i) = vt(i) + abs(a%val(j))
@ -1781,18 +1805,17 @@ subroutine psb_c_coo_rowsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = cone
else
d = czero d = czero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
d(i) = d(i) + a%val(j) d(i) = d(i) + a%val(j)
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, m
d(i) = d(i) + cone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -1835,19 +1858,17 @@ subroutine psb_c_coo_arwsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = sone
else
d = szero d = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
d(i) = d(i) + abs(a%val(j)) d(i) = d(i) + abs(a%val(j))
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, m
d(i) = d(i) + sone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1889,19 +1910,18 @@ subroutine psb_c_coo_colsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = cone
else
d = czero d = czero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
k = a%ja(j) k = a%ja(j)
d(k) = d(k) + a%val(j) d(k) = d(k) + a%val(j)
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, n
d(i) = d(i) + cone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1943,19 +1963,18 @@ subroutine psb_c_coo_aclsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = sone
else
d = szero d = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
k = a%ja(j) k = a%ja(j)
d(k) = d(k) + abs(a%val(j)) d(k) = d(k) + abs(a%val(j))
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, n
d(i) = d(i) + sone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

71
base/serial/impl/psb_c_csc_impl.f90
Unescape Escape View File

@ -304,7 +304,7 @@ subroutine psb_c_csc_csmv(alpha,a,x,beta,y,info,trans)
endif endif
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, min(m,n) do i=1, min(m,n)
y(i) = y(i) + alpha*x(i) y(i) = y(i) + alpha*x(i)
end do end do
@ -590,7 +590,7 @@ subroutine psb_c_csc_csmm(alpha,a,x,beta,y,info,trans)
endif endif
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, min(m,n) do i=1, min(m,n)
y(i,:) = y(i,:) + alpha*x(i,:) y(i,:) = y(i,:) + alpha*x(i,:)
end do end do
@ -1081,7 +1081,12 @@ function psb_c_csc_maxval(a) result(res)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
if (a%is_unit()) then
res = sone
else
res = szero res = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
if (allocated(a%val)) then if (allocated(a%val)) then
nnz = min(nnz,size(a%val)) nnz = min(nnz,size(a%val))
@ -1112,7 +1117,11 @@ function psb_c_csc_csnmi(a) result(res)
if (info /= psb_success_) then if (info /= psb_success_) then
return return
end if end if
acc(:) = szero if (a%is_unit()) then
acc = sone
else
acc = szero
end if
do i=1, nc do i=1, nc
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j)) acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j))
@ -1149,7 +1158,11 @@ function psb_c_csc_csnm1(a) result(res)
m = a%get_nrows() m = a%get_nrows()
n = a%get_ncols() n = a%get_ncols()
do j=1, n do j=1, n
if (a%is_unit()) then
acc = sone
else
acc = szero acc = szero
end if
do k=a%icp(j),a%icp(j+1)-1 do k=a%icp(j),a%icp(j+1)-1
acc = acc + abs(a%val(k)) acc = acc + abs(a%val(k))
end do end do
@ -1187,19 +1200,17 @@ subroutine psb_c_csc_colsum(d,a)
end if end if
do i = 1, a%get_ncols() do i = 1, a%get_ncols()
if (a%is_unit()) then
d(i) = cone
else
d(i) = czero d(i) = czero
end if
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
d(i) = d(i) + (a%val(j)) d(i) = d(i) + (a%val(j))
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_ncols()
d(i) = d(i) + cone
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1242,13 +1253,18 @@ subroutine psb_c_csc_aclsum(d,a)
do i = 1, a%get_ncols() do i = 1, a%get_ncols()
if (a%is_unit()) then
d(i) = sone
else
d(i) = szero d(i) = szero
end if
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
d(i) = d(i) + abs(a%val(j)) d(i) = d(i) + abs(a%val(j))
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, a%get_ncols() do i=1, a%get_ncols()
d(i) = d(i) + sone d(i) = d(i) + sone
end do end do
@ -1295,7 +1311,11 @@ subroutine psb_c_csc_rowsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = cone
else
d = czero d = czero
end if
do i=1, m do i=1, m
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
@ -1304,13 +1324,6 @@ subroutine psb_c_csc_rowsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_nrows()
d(i) = d(i) + cone
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1352,7 +1365,12 @@ subroutine psb_c_csc_arwsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = sone
else
d = szero d = szero
end if
do i=1, m do i=1, m
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
k = a%ia(j) k = a%ia(j)
@ -1360,13 +1378,6 @@ subroutine psb_c_csc_arwsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_nrows()
d(i) = d(i) + sone
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1408,7 +1419,7 @@ subroutine psb_c_csc_get_diag(a,d,info)
end if end if
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = cone d(1:mnm) = cone
else else
do i=1, mnm do i=1, mnm
@ -1450,6 +1461,7 @@ subroutine psb_c_csc_scal(d,a,info,side)
character, intent(in), optional :: side character, intent(in), optional :: side
integer(psb_ipk_) :: err_act,mnm, i, j, n integer(psb_ipk_) :: err_act,mnm, i, j, n
type(psb_c_coo_sparse_mat) :: tmp
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='scal' character(len=20) :: name='scal'
character :: side_ character :: side_
@ -1464,6 +1476,10 @@ subroutine psb_c_csc_scal(d,a,info,side)
side_ = psb_toupper(side) side_ = psb_toupper(side)
end if end if
if (a%is_unit()) then
call a%add_unit_diag()
end if
left = (side_ == 'L') left = (side_ == 'L')
if (left) then if (left) then
@ -1524,6 +1540,9 @@ subroutine psb_c_csc_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d

21
base/serial/impl/psb_c_csr_impl.f90
Unescape Escape View File

@ -380,7 +380,7 @@ contains
endif endif
if (is_triangle.and.is_unit) then if (is_unit) then
do i=1, min(m,n) do i=1, min(m,n)
y(i) = y(i) + alpha*x(i) y(i) = y(i) + alpha*x(i)
end do end do
@ -733,7 +733,7 @@ contains
endif endif
if (is_triangle.and.is_unit) then if (is_unit) then
do i=1, min(m,n) do i=1, min(m,n)
y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc) y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc)
end do end do
@ -1398,7 +1398,7 @@ subroutine psb_c_csr_rowsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, m do i=1, m
d(i) = d(i) + cone d(i) = d(i) + cone
end do end do
@ -1453,7 +1453,7 @@ subroutine psb_c_csr_arwsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, m do i=1, m
d(i) = d(i) + sone d(i) = d(i) + sone
end do end do
@ -1509,7 +1509,7 @@ subroutine psb_c_csr_colsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, n do i=1, n
d(i) = d(i) + cone d(i) = d(i) + cone
end do end do
@ -1566,7 +1566,7 @@ subroutine psb_c_csr_aclsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, n do i=1, n
d(i) = d(i) + sone d(i) = d(i) + sone
end do end do
@ -1613,7 +1613,7 @@ subroutine psb_c_csr_get_diag(a,d,info)
end if end if
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = cone d(1:mnm) = cone
else else
do i=1, mnm do i=1, mnm
@ -1665,6 +1665,10 @@ subroutine psb_c_csr_scal(d,a,info,side)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
side_ = 'L' side_ = 'L'
if (present(side)) then if (present(side)) then
side_ = psb_toupper(side) side_ = psb_toupper(side)
@ -1734,6 +1738,9 @@ subroutine psb_c_csr_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d

31
base/serial/impl/psb_d_base_mat_impl.F90
Unescape Escape View File

@ -580,6 +580,37 @@ subroutine psb_d_base_clone(a,b,info)
end subroutine psb_d_base_clone end subroutine psb_d_base_clone
subroutine psb_d_base_add_unit_diag(a)
use psb_d_base_mat_mod, psb_protect_name => psb_d_base_add_unit_diag
use psb_error_mod
implicit none
class(psb_d_base_sparse_mat), intent(inout) :: a
type(psb_d_coo_sparse_mat) :: tmp
integer(psb_ipk_) :: i, j, m, n, nz, mnm, info
if (a%is_unit()) then
call a%mv_to_coo(tmp,info)
if (info /= 0) return
m = tmp%get_nrows()
n = tmp%get_ncols()
mnm = min(m,n)
nz = tmp%get_nzeros()
call tmp%reallocate(nz+mnm)
do i=1, mnm
tmp%val(nz+i) = done
tmp%ia(nz+i) = i
tmp%ja(nz+i) = i
end do
call tmp%set_nzeros(nz+mnm)
call tmp%set_unit(.false.)
call tmp%fix(info)
if (info /= 0) &
& call a%mv_from_coo(tmp,info)
end if
end subroutine psb_d_base_add_unit_diag
subroutine psb_d_base_mold(a,b,info) subroutine psb_d_base_mold(a,b,info)
use psb_d_base_mat_mod, psb_protect_name => psb_d_base_mold use psb_d_base_mat_mod, psb_protect_name => psb_d_base_mold
use psb_error_mod use psb_error_mod

77
base/serial/impl/psb_d_coo_impl.f90
Unescape Escape View File

@ -54,11 +54,11 @@ subroutine psb_d_coo_get_diag(a,d,info)
call psb_errpush(info,name,i_err=ierr) call psb_errpush(info,name,i_err=ierr)
goto 9999 goto 9999
end if end if
d(:) = dzero
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = done d(1:mnm) = done
else else
d(1:mnm) = dzero
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
j=a%ia(i) j=a%ia(i)
if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then
@ -101,6 +101,10 @@ subroutine psb_d_coo_scal(d,a,info,side)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
side_ = 'L' side_ = 'L'
if (present(side)) then if (present(side)) then
side_ = psb_toupper(side) side_ = psb_toupper(side)
@ -167,6 +171,9 @@ subroutine psb_d_coo_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d
@ -1313,7 +1320,7 @@ subroutine psb_d_coo_csmv(alpha,a,x,beta,y,info,trans)
endif endif
return return
else else
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
if (beta == dzero) then if (beta == dzero) then
do i = 1, min(m,n) do i = 1, min(m,n)
y(i) = alpha*x(i) y(i) = alpha*x(i)
@ -1524,7 +1531,7 @@ subroutine psb_d_coo_csmm(alpha,a,x,beta,y,info,trans)
endif endif
return return
else else
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
if (beta == dzero) then if (beta == dzero) then
do i = 1, min(m,n) do i = 1, min(m,n)
y(i,1:nc) = alpha*x(i,1:nc) y(i,1:nc) = alpha*x(i,1:nc)
@ -1661,12 +1668,17 @@ function psb_d_coo_maxval(a) result(res)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
if (a%is_unit()) then
res = done
else
res = dzero res = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
if (allocated(a%val)) then if (allocated(a%val)) then
nnz = min(nnz,size(a%val)) nnz = min(nnz,size(a%val))
res = maxval(abs(a%val(1:nnz))) res = maxval(abs(a%val(1:nnz)))
end if end if
end function psb_d_coo_maxval end function psb_d_coo_maxval
function psb_d_coo_csnmi(a) result(res) function psb_d_coo_csnmi(a) result(res)
@ -1696,7 +1708,11 @@ function psb_d_coo_csnmi(a) result(res)
do while ((a%ia(j) == a%ia(i)).and. (j <= nnz)) do while ((a%ia(j) == a%ia(i)).and. (j <= nnz))
j = j+1 j = j+1
enddo enddo
if (a%is_unit()) then
acc = done
else
acc = dzero acc = dzero
end if
do k=i, j-1 do k=i, j-1
acc = acc + abs(a%val(k)) acc = acc + abs(a%val(k))
end do end do
@ -1707,7 +1723,11 @@ function psb_d_coo_csnmi(a) result(res)
m = a%get_nrows() m = a%get_nrows()
allocate(vt(m),stat=info) allocate(vt(m),stat=info)
if (info /= 0) return if (info /= 0) return
vt(:) = dzero if (a%is_unit()) then
vt = done
else
vt = dzero
end if
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
vt(i) = vt(i) + abs(a%val(j)) vt(i) = vt(i) + abs(a%val(j))
@ -1743,7 +1763,11 @@ function psb_d_coo_csnm1(a) result(res)
n = a%get_ncols() n = a%get_ncols()
allocate(vt(n),stat=info) allocate(vt(n),stat=info)
if (info /= 0) return if (info /= 0) return
vt(:) = dzero if (a%is_unit()) then
vt = done
else
vt = dzero
end if
do j=1, nnz do j=1, nnz
i = a%ja(j) i = a%ja(j)
vt(i) = vt(i) + abs(a%val(j)) vt(i) = vt(i) + abs(a%val(j))
@ -1781,18 +1805,17 @@ subroutine psb_d_coo_rowsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = done
else
d = dzero d = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
d(i) = d(i) + a%val(j) d(i) = d(i) + a%val(j)
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, m
d(i) = d(i) + done
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -1835,19 +1858,17 @@ subroutine psb_d_coo_arwsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = done
else
d = dzero d = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
d(i) = d(i) + abs(a%val(j)) d(i) = d(i) + abs(a%val(j))
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, m
d(i) = d(i) + done
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1889,19 +1910,18 @@ subroutine psb_d_coo_colsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = done
else
d = dzero d = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
k = a%ja(j) k = a%ja(j)
d(k) = d(k) + a%val(j) d(k) = d(k) + a%val(j)
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, n
d(i) = d(i) + done
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1943,19 +1963,18 @@ subroutine psb_d_coo_aclsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = done
else
d = dzero d = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
k = a%ja(j) k = a%ja(j)
d(k) = d(k) + abs(a%val(j)) d(k) = d(k) + abs(a%val(j))
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, n
d(i) = d(i) + done
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

71
base/serial/impl/psb_d_csc_impl.f90
Unescape Escape View File

@ -304,7 +304,7 @@ subroutine psb_d_csc_csmv(alpha,a,x,beta,y,info,trans)
endif endif
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, min(m,n) do i=1, min(m,n)
y(i) = y(i) + alpha*x(i) y(i) = y(i) + alpha*x(i)
end do end do
@ -590,7 +590,7 @@ subroutine psb_d_csc_csmm(alpha,a,x,beta,y,info,trans)
endif endif
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, min(m,n) do i=1, min(m,n)
y(i,:) = y(i,:) + alpha*x(i,:) y(i,:) = y(i,:) + alpha*x(i,:)
end do end do
@ -1081,7 +1081,12 @@ function psb_d_csc_maxval(a) result(res)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
if (a%is_unit()) then
res = done
else
res = dzero res = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
if (allocated(a%val)) then if (allocated(a%val)) then
nnz = min(nnz,size(a%val)) nnz = min(nnz,size(a%val))
@ -1112,7 +1117,11 @@ function psb_d_csc_csnmi(a) result(res)
if (info /= psb_success_) then if (info /= psb_success_) then
return return
end if end if
acc(:) = dzero if (a%is_unit()) then
acc = done
else
acc = dzero
end if
do i=1, nc do i=1, nc
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j)) acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j))
@ -1149,7 +1158,11 @@ function psb_d_csc_csnm1(a) result(res)
m = a%get_nrows() m = a%get_nrows()
n = a%get_ncols() n = a%get_ncols()
do j=1, n do j=1, n
if (a%is_unit()) then
acc = done
else
acc = dzero acc = dzero
end if
do k=a%icp(j),a%icp(j+1)-1 do k=a%icp(j),a%icp(j+1)-1
acc = acc + abs(a%val(k)) acc = acc + abs(a%val(k))
end do end do
@ -1187,19 +1200,17 @@ subroutine psb_d_csc_colsum(d,a)
end if end if
do i = 1, a%get_ncols() do i = 1, a%get_ncols()
if (a%is_unit()) then
d(i) = done
else
d(i) = dzero d(i) = dzero
end if
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
d(i) = d(i) + (a%val(j)) d(i) = d(i) + (a%val(j))
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_ncols()
d(i) = d(i) + done
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1242,13 +1253,18 @@ subroutine psb_d_csc_aclsum(d,a)
do i = 1, a%get_ncols() do i = 1, a%get_ncols()
if (a%is_unit()) then
d(i) = done
else
d(i) = dzero d(i) = dzero
end if
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
d(i) = d(i) + abs(a%val(j)) d(i) = d(i) + abs(a%val(j))
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, a%get_ncols() do i=1, a%get_ncols()
d(i) = d(i) + done d(i) = d(i) + done
end do end do
@ -1295,7 +1311,11 @@ subroutine psb_d_csc_rowsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = done
else
d = dzero d = dzero
end if
do i=1, m do i=1, m
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
@ -1304,13 +1324,6 @@ subroutine psb_d_csc_rowsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_nrows()
d(i) = d(i) + done
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1352,7 +1365,12 @@ subroutine psb_d_csc_arwsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = done
else
d = dzero d = dzero
end if
do i=1, m do i=1, m
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
k = a%ia(j) k = a%ia(j)
@ -1360,13 +1378,6 @@ subroutine psb_d_csc_arwsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_nrows()
d(i) = d(i) + done
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1408,7 +1419,7 @@ subroutine psb_d_csc_get_diag(a,d,info)
end if end if
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = done d(1:mnm) = done
else else
do i=1, mnm do i=1, mnm
@ -1450,6 +1461,7 @@ subroutine psb_d_csc_scal(d,a,info,side)
character, intent(in), optional :: side character, intent(in), optional :: side
integer(psb_ipk_) :: err_act,mnm, i, j, n integer(psb_ipk_) :: err_act,mnm, i, j, n
type(psb_d_coo_sparse_mat) :: tmp
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='scal' character(len=20) :: name='scal'
character :: side_ character :: side_
@ -1464,6 +1476,10 @@ subroutine psb_d_csc_scal(d,a,info,side)
side_ = psb_toupper(side) side_ = psb_toupper(side)
end if end if
if (a%is_unit()) then
call a%add_unit_diag()
end if
left = (side_ == 'L') left = (side_ == 'L')
if (left) then if (left) then
@ -1524,6 +1540,9 @@ subroutine psb_d_csc_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d

21
base/serial/impl/psb_d_csr_impl.f90
Unescape Escape View File

@ -380,7 +380,7 @@ contains
endif endif
if (is_triangle.and.is_unit) then if (is_unit) then
do i=1, min(m,n) do i=1, min(m,n)
y(i) = y(i) + alpha*x(i) y(i) = y(i) + alpha*x(i)
end do end do
@ -733,7 +733,7 @@ contains
endif endif
if (is_triangle.and.is_unit) then if (is_unit) then
do i=1, min(m,n) do i=1, min(m,n)
y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc) y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc)
end do end do
@ -1398,7 +1398,7 @@ subroutine psb_d_csr_rowsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, m do i=1, m
d(i) = d(i) + done d(i) = d(i) + done
end do end do
@ -1453,7 +1453,7 @@ subroutine psb_d_csr_arwsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, m do i=1, m
d(i) = d(i) + done d(i) = d(i) + done
end do end do
@ -1509,7 +1509,7 @@ subroutine psb_d_csr_colsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, n do i=1, n
d(i) = d(i) + done d(i) = d(i) + done
end do end do
@ -1566,7 +1566,7 @@ subroutine psb_d_csr_aclsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, n do i=1, n
d(i) = d(i) + done d(i) = d(i) + done
end do end do
@ -1613,7 +1613,7 @@ subroutine psb_d_csr_get_diag(a,d,info)
end if end if
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = done d(1:mnm) = done
else else
do i=1, mnm do i=1, mnm
@ -1665,6 +1665,10 @@ subroutine psb_d_csr_scal(d,a,info,side)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
side_ = 'L' side_ = 'L'
if (present(side)) then if (present(side)) then
side_ = psb_toupper(side) side_ = psb_toupper(side)
@ -1734,6 +1738,9 @@ subroutine psb_d_csr_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d

31
base/serial/impl/psb_s_base_mat_impl.F90
Unescape Escape View File

@ -580,6 +580,37 @@ subroutine psb_s_base_clone(a,b,info)
end subroutine psb_s_base_clone end subroutine psb_s_base_clone
subroutine psb_s_base_add_unit_diag(a)
use psb_s_base_mat_mod, psb_protect_name => psb_s_base_add_unit_diag
use psb_error_mod
implicit none
class(psb_s_base_sparse_mat), intent(inout) :: a
type(psb_s_coo_sparse_mat) :: tmp
integer(psb_ipk_) :: i, j, m, n, nz, mnm, info
if (a%is_unit()) then
call a%mv_to_coo(tmp,info)
if (info /= 0) return
m = tmp%get_nrows()
n = tmp%get_ncols()
mnm = min(m,n)
nz = tmp%get_nzeros()
call tmp%reallocate(nz+mnm)
do i=1, mnm
tmp%val(nz+i) = sone
tmp%ia(nz+i) = i
tmp%ja(nz+i) = i
end do
call tmp%set_nzeros(nz+mnm)
call tmp%set_unit(.false.)
call tmp%fix(info)
if (info /= 0) &
& call a%mv_from_coo(tmp,info)
end if
end subroutine psb_s_base_add_unit_diag
subroutine psb_s_base_mold(a,b,info) subroutine psb_s_base_mold(a,b,info)
use psb_s_base_mat_mod, psb_protect_name => psb_s_base_mold use psb_s_base_mat_mod, psb_protect_name => psb_s_base_mold
use psb_error_mod use psb_error_mod

77
base/serial/impl/psb_s_coo_impl.f90
Unescape Escape View File

@ -54,11 +54,11 @@ subroutine psb_s_coo_get_diag(a,d,info)
call psb_errpush(info,name,i_err=ierr) call psb_errpush(info,name,i_err=ierr)
goto 9999 goto 9999
end if end if
d(:) = szero
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = sone d(1:mnm) = sone
else else
d(1:mnm) = szero
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
j=a%ia(i) j=a%ia(i)
if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then
@ -101,6 +101,10 @@ subroutine psb_s_coo_scal(d,a,info,side)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
side_ = 'L' side_ = 'L'
if (present(side)) then if (present(side)) then
side_ = psb_toupper(side) side_ = psb_toupper(side)
@ -167,6 +171,9 @@ subroutine psb_s_coo_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d
@ -1313,7 +1320,7 @@ subroutine psb_s_coo_csmv(alpha,a,x,beta,y,info,trans)
endif endif
return return
else else
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
if (beta == szero) then if (beta == szero) then
do i = 1, min(m,n) do i = 1, min(m,n)
y(i) = alpha*x(i) y(i) = alpha*x(i)
@ -1524,7 +1531,7 @@ subroutine psb_s_coo_csmm(alpha,a,x,beta,y,info,trans)
endif endif
return return
else else
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
if (beta == szero) then if (beta == szero) then
do i = 1, min(m,n) do i = 1, min(m,n)
y(i,1:nc) = alpha*x(i,1:nc) y(i,1:nc) = alpha*x(i,1:nc)
@ -1661,12 +1668,17 @@ function psb_s_coo_maxval(a) result(res)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
if (a%is_unit()) then
res = sone
else
res = szero res = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
if (allocated(a%val)) then if (allocated(a%val)) then
nnz = min(nnz,size(a%val)) nnz = min(nnz,size(a%val))
res = maxval(abs(a%val(1:nnz))) res = maxval(abs(a%val(1:nnz)))
end if end if
end function psb_s_coo_maxval end function psb_s_coo_maxval
function psb_s_coo_csnmi(a) result(res) function psb_s_coo_csnmi(a) result(res)
@ -1696,7 +1708,11 @@ function psb_s_coo_csnmi(a) result(res)
do while ((a%ia(j) == a%ia(i)).and. (j <= nnz)) do while ((a%ia(j) == a%ia(i)).and. (j <= nnz))
j = j+1 j = j+1
enddo enddo
if (a%is_unit()) then
acc = sone
else
acc = szero acc = szero
end if
do k=i, j-1 do k=i, j-1
acc = acc + abs(a%val(k)) acc = acc + abs(a%val(k))
end do end do
@ -1707,7 +1723,11 @@ function psb_s_coo_csnmi(a) result(res)
m = a%get_nrows() m = a%get_nrows()
allocate(vt(m),stat=info) allocate(vt(m),stat=info)
if (info /= 0) return if (info /= 0) return
vt(:) = szero if (a%is_unit()) then
vt = sone
else
vt = szero
end if
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
vt(i) = vt(i) + abs(a%val(j)) vt(i) = vt(i) + abs(a%val(j))
@ -1743,7 +1763,11 @@ function psb_s_coo_csnm1(a) result(res)
n = a%get_ncols() n = a%get_ncols()
allocate(vt(n),stat=info) allocate(vt(n),stat=info)
if (info /= 0) return if (info /= 0) return
vt(:) = szero if (a%is_unit()) then
vt = sone
else
vt = szero
end if
do j=1, nnz do j=1, nnz
i = a%ja(j) i = a%ja(j)
vt(i) = vt(i) + abs(a%val(j)) vt(i) = vt(i) + abs(a%val(j))
@ -1781,18 +1805,17 @@ subroutine psb_s_coo_rowsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = sone
else
d = szero d = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
d(i) = d(i) + a%val(j) d(i) = d(i) + a%val(j)
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, m
d(i) = d(i) + sone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -1835,19 +1858,17 @@ subroutine psb_s_coo_arwsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = sone
else
d = szero d = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
d(i) = d(i) + abs(a%val(j)) d(i) = d(i) + abs(a%val(j))
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, m
d(i) = d(i) + sone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1889,19 +1910,18 @@ subroutine psb_s_coo_colsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = sone
else
d = szero d = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
k = a%ja(j) k = a%ja(j)
d(k) = d(k) + a%val(j) d(k) = d(k) + a%val(j)
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, n
d(i) = d(i) + sone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1943,19 +1963,18 @@ subroutine psb_s_coo_aclsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = sone
else
d = szero d = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
k = a%ja(j) k = a%ja(j)
d(k) = d(k) + abs(a%val(j)) d(k) = d(k) + abs(a%val(j))
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, n
d(i) = d(i) + sone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

71
base/serial/impl/psb_s_csc_impl.f90
Unescape Escape View File

@ -304,7 +304,7 @@ subroutine psb_s_csc_csmv(alpha,a,x,beta,y,info,trans)
endif endif
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, min(m,n) do i=1, min(m,n)
y(i) = y(i) + alpha*x(i) y(i) = y(i) + alpha*x(i)
end do end do
@ -590,7 +590,7 @@ subroutine psb_s_csc_csmm(alpha,a,x,beta,y,info,trans)
endif endif
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, min(m,n) do i=1, min(m,n)
y(i,:) = y(i,:) + alpha*x(i,:) y(i,:) = y(i,:) + alpha*x(i,:)
end do end do
@ -1081,7 +1081,12 @@ function psb_s_csc_maxval(a) result(res)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
if (a%is_unit()) then
res = sone
else
res = szero res = szero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
if (allocated(a%val)) then if (allocated(a%val)) then
nnz = min(nnz,size(a%val)) nnz = min(nnz,size(a%val))
@ -1112,7 +1117,11 @@ function psb_s_csc_csnmi(a) result(res)
if (info /= psb_success_) then if (info /= psb_success_) then
return return
end if end if
acc(:) = szero if (a%is_unit()) then
acc = sone
else
acc = szero
end if
do i=1, nc do i=1, nc
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j)) acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j))
@ -1149,7 +1158,11 @@ function psb_s_csc_csnm1(a) result(res)
m = a%get_nrows() m = a%get_nrows()
n = a%get_ncols() n = a%get_ncols()
do j=1, n do j=1, n
if (a%is_unit()) then
acc = sone
else
acc = szero acc = szero
end if
do k=a%icp(j),a%icp(j+1)-1 do k=a%icp(j),a%icp(j+1)-1
acc = acc + abs(a%val(k)) acc = acc + abs(a%val(k))
end do end do
@ -1187,19 +1200,17 @@ subroutine psb_s_csc_colsum(d,a)
end if end if
do i = 1, a%get_ncols() do i = 1, a%get_ncols()
if (a%is_unit()) then
d(i) = sone
else
d(i) = szero d(i) = szero
end if
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
d(i) = d(i) + (a%val(j)) d(i) = d(i) + (a%val(j))
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_ncols()
d(i) = d(i) + sone
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1242,13 +1253,18 @@ subroutine psb_s_csc_aclsum(d,a)
do i = 1, a%get_ncols() do i = 1, a%get_ncols()
if (a%is_unit()) then
d(i) = sone
else
d(i) = szero d(i) = szero
end if
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
d(i) = d(i) + abs(a%val(j)) d(i) = d(i) + abs(a%val(j))
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, a%get_ncols() do i=1, a%get_ncols()
d(i) = d(i) + sone d(i) = d(i) + sone
end do end do
@ -1295,7 +1311,11 @@ subroutine psb_s_csc_rowsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = sone
else
d = szero d = szero
end if
do i=1, m do i=1, m
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
@ -1304,13 +1324,6 @@ subroutine psb_s_csc_rowsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_nrows()
d(i) = d(i) + sone
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1352,7 +1365,12 @@ subroutine psb_s_csc_arwsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = sone
else
d = szero d = szero
end if
do i=1, m do i=1, m
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
k = a%ia(j) k = a%ia(j)
@ -1360,13 +1378,6 @@ subroutine psb_s_csc_arwsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_nrows()
d(i) = d(i) + sone
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1408,7 +1419,7 @@ subroutine psb_s_csc_get_diag(a,d,info)
end if end if
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = sone d(1:mnm) = sone
else else
do i=1, mnm do i=1, mnm
@ -1450,6 +1461,7 @@ subroutine psb_s_csc_scal(d,a,info,side)
character, intent(in), optional :: side character, intent(in), optional :: side
integer(psb_ipk_) :: err_act,mnm, i, j, n integer(psb_ipk_) :: err_act,mnm, i, j, n
type(psb_s_coo_sparse_mat) :: tmp
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='scal' character(len=20) :: name='scal'
character :: side_ character :: side_
@ -1464,6 +1476,10 @@ subroutine psb_s_csc_scal(d,a,info,side)
side_ = psb_toupper(side) side_ = psb_toupper(side)
end if end if
if (a%is_unit()) then
call a%add_unit_diag()
end if
left = (side_ == 'L') left = (side_ == 'L')
if (left) then if (left) then
@ -1524,6 +1540,9 @@ subroutine psb_s_csc_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d

21
base/serial/impl/psb_s_csr_impl.f90
Unescape Escape View File

@ -380,7 +380,7 @@ contains
endif endif
if (is_triangle.and.is_unit) then if (is_unit) then
do i=1, min(m,n) do i=1, min(m,n)
y(i) = y(i) + alpha*x(i) y(i) = y(i) + alpha*x(i)
end do end do
@ -733,7 +733,7 @@ contains
endif endif
if (is_triangle.and.is_unit) then if (is_unit) then
do i=1, min(m,n) do i=1, min(m,n)
y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc) y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc)
end do end do
@ -1398,7 +1398,7 @@ subroutine psb_s_csr_rowsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, m do i=1, m
d(i) = d(i) + sone d(i) = d(i) + sone
end do end do
@ -1453,7 +1453,7 @@ subroutine psb_s_csr_arwsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, m do i=1, m
d(i) = d(i) + sone d(i) = d(i) + sone
end do end do
@ -1509,7 +1509,7 @@ subroutine psb_s_csr_colsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, n do i=1, n
d(i) = d(i) + sone d(i) = d(i) + sone
end do end do
@ -1566,7 +1566,7 @@ subroutine psb_s_csr_aclsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, n do i=1, n
d(i) = d(i) + sone d(i) = d(i) + sone
end do end do
@ -1613,7 +1613,7 @@ subroutine psb_s_csr_get_diag(a,d,info)
end if end if
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = sone d(1:mnm) = sone
else else
do i=1, mnm do i=1, mnm
@ -1665,6 +1665,10 @@ subroutine psb_s_csr_scal(d,a,info,side)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
side_ = 'L' side_ = 'L'
if (present(side)) then if (present(side)) then
side_ = psb_toupper(side) side_ = psb_toupper(side)
@ -1734,6 +1738,9 @@ subroutine psb_s_csr_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d

31
base/serial/impl/psb_z_base_mat_impl.F90
Unescape Escape View File

@ -580,6 +580,37 @@ subroutine psb_z_base_clone(a,b,info)
end subroutine psb_z_base_clone end subroutine psb_z_base_clone
subroutine psb_z_base_add_unit_diag(a)
use psb_z_base_mat_mod, psb_protect_name => psb_z_base_add_unit_diag
use psb_error_mod
implicit none
class(psb_z_base_sparse_mat), intent(inout) :: a
type(psb_z_coo_sparse_mat) :: tmp
integer(psb_ipk_) :: i, j, m, n, nz, mnm, info
if (a%is_unit()) then
call a%mv_to_coo(tmp,info)
if (info /= 0) return
m = tmp%get_nrows()
n = tmp%get_ncols()
mnm = min(m,n)
nz = tmp%get_nzeros()
call tmp%reallocate(nz+mnm)
do i=1, mnm
tmp%val(nz+i) = zone
tmp%ia(nz+i) = i
tmp%ja(nz+i) = i
end do
call tmp%set_nzeros(nz+mnm)
call tmp%set_unit(.false.)
call tmp%fix(info)
if (info /= 0) &
& call a%mv_from_coo(tmp,info)
end if
end subroutine psb_z_base_add_unit_diag
subroutine psb_z_base_mold(a,b,info) subroutine psb_z_base_mold(a,b,info)
use psb_z_base_mat_mod, psb_protect_name => psb_z_base_mold use psb_z_base_mat_mod, psb_protect_name => psb_z_base_mold
use psb_error_mod use psb_error_mod

77
base/serial/impl/psb_z_coo_impl.f90
Unescape Escape View File

@ -54,11 +54,11 @@ subroutine psb_z_coo_get_diag(a,d,info)
call psb_errpush(info,name,i_err=ierr) call psb_errpush(info,name,i_err=ierr)
goto 9999 goto 9999
end if end if
d(:) = zzero
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = zone d(1:mnm) = zone
else else
d(1:mnm) = zzero
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
j=a%ia(i) j=a%ia(i)
if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then if ((j == a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then
@ -101,6 +101,10 @@ subroutine psb_z_coo_scal(d,a,info,side)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
side_ = 'L' side_ = 'L'
if (present(side)) then if (present(side)) then
side_ = psb_toupper(side) side_ = psb_toupper(side)
@ -167,6 +171,9 @@ subroutine psb_z_coo_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d
@ -1313,7 +1320,7 @@ subroutine psb_z_coo_csmv(alpha,a,x,beta,y,info,trans)
endif endif
return return
else else
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
if (beta == zzero) then if (beta == zzero) then
do i = 1, min(m,n) do i = 1, min(m,n)
y(i) = alpha*x(i) y(i) = alpha*x(i)
@ -1524,7 +1531,7 @@ subroutine psb_z_coo_csmm(alpha,a,x,beta,y,info,trans)
endif endif
return return
else else
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
if (beta == zzero) then if (beta == zzero) then
do i = 1, min(m,n) do i = 1, min(m,n)
y(i,1:nc) = alpha*x(i,1:nc) y(i,1:nc) = alpha*x(i,1:nc)
@ -1661,12 +1668,17 @@ function psb_z_coo_maxval(a) result(res)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
if (a%is_unit()) then
res = done
else
res = dzero res = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
if (allocated(a%val)) then if (allocated(a%val)) then
nnz = min(nnz,size(a%val)) nnz = min(nnz,size(a%val))
res = maxval(abs(a%val(1:nnz))) res = maxval(abs(a%val(1:nnz)))
end if end if
end function psb_z_coo_maxval end function psb_z_coo_maxval
function psb_z_coo_csnmi(a) result(res) function psb_z_coo_csnmi(a) result(res)
@ -1696,7 +1708,11 @@ function psb_z_coo_csnmi(a) result(res)
do while ((a%ia(j) == a%ia(i)).and. (j <= nnz)) do while ((a%ia(j) == a%ia(i)).and. (j <= nnz))
j = j+1 j = j+1
enddo enddo
if (a%is_unit()) then
acc = done
else
acc = dzero acc = dzero
end if
do k=i, j-1 do k=i, j-1
acc = acc + abs(a%val(k)) acc = acc + abs(a%val(k))
end do end do
@ -1707,7 +1723,11 @@ function psb_z_coo_csnmi(a) result(res)
m = a%get_nrows() m = a%get_nrows()
allocate(vt(m),stat=info) allocate(vt(m),stat=info)
if (info /= 0) return if (info /= 0) return
vt(:) = dzero if (a%is_unit()) then
vt = done
else
vt = dzero
end if
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
vt(i) = vt(i) + abs(a%val(j)) vt(i) = vt(i) + abs(a%val(j))
@ -1743,7 +1763,11 @@ function psb_z_coo_csnm1(a) result(res)
n = a%get_ncols() n = a%get_ncols()
allocate(vt(n),stat=info) allocate(vt(n),stat=info)
if (info /= 0) return if (info /= 0) return
vt(:) = dzero if (a%is_unit()) then
vt = done
else
vt = dzero
end if
do j=1, nnz do j=1, nnz
i = a%ja(j) i = a%ja(j)
vt(i) = vt(i) + abs(a%val(j)) vt(i) = vt(i) + abs(a%val(j))
@ -1781,18 +1805,17 @@ subroutine psb_z_coo_rowsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = zone
else
d = zzero d = zzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
d(i) = d(i) + a%val(j) d(i) = d(i) + a%val(j)
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, m
d(i) = d(i) + zone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -1835,19 +1858,17 @@ subroutine psb_z_coo_arwsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = done
else
d = dzero d = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
i = a%ia(j) i = a%ia(j)
d(i) = d(i) + abs(a%val(j)) d(i) = d(i) + abs(a%val(j))
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, m
d(i) = d(i) + done
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1889,19 +1910,18 @@ subroutine psb_z_coo_colsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = zone
else
d = zzero d = zzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
k = a%ja(j) k = a%ja(j)
d(k) = d(k) + a%val(j) d(k) = d(k) + a%val(j)
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, n
d(i) = d(i) + zone
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1943,19 +1963,18 @@ subroutine psb_z_coo_aclsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = done
else
d = dzero d = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
do j=1, nnz do j=1, nnz
k = a%ja(j) k = a%ja(j)
d(k) = d(k) + abs(a%val(j)) d(k) = d(k) + abs(a%val(j))
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, n
d(i) = d(i) + done
end do
end if
return return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

71
base/serial/impl/psb_z_csc_impl.f90
Unescape Escape View File

@ -304,7 +304,7 @@ subroutine psb_z_csc_csmv(alpha,a,x,beta,y,info,trans)
endif endif
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, min(m,n) do i=1, min(m,n)
y(i) = y(i) + alpha*x(i) y(i) = y(i) + alpha*x(i)
end do end do
@ -590,7 +590,7 @@ subroutine psb_z_csc_csmm(alpha,a,x,beta,y,info,trans)
endif endif
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, min(m,n) do i=1, min(m,n)
y(i,:) = y(i,:) + alpha*x(i,:) y(i,:) = y(i,:) + alpha*x(i,:)
end do end do
@ -1081,7 +1081,12 @@ function psb_z_csc_maxval(a) result(res)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
if (a%is_unit()) then
res = done
else
res = dzero res = dzero
end if
nnz = a%get_nzeros() nnz = a%get_nzeros()
if (allocated(a%val)) then if (allocated(a%val)) then
nnz = min(nnz,size(a%val)) nnz = min(nnz,size(a%val))
@ -1112,7 +1117,11 @@ function psb_z_csc_csnmi(a) result(res)
if (info /= psb_success_) then if (info /= psb_success_) then
return return
end if end if
acc(:) = dzero if (a%is_unit()) then
acc = done
else
acc = dzero
end if
do i=1, nc do i=1, nc
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j)) acc(a%ia(j)) = acc(a%ia(j)) + abs(a%val(j))
@ -1149,7 +1158,11 @@ function psb_z_csc_csnm1(a) result(res)
m = a%get_nrows() m = a%get_nrows()
n = a%get_ncols() n = a%get_ncols()
do j=1, n do j=1, n
if (a%is_unit()) then
acc = done
else
acc = dzero acc = dzero
end if
do k=a%icp(j),a%icp(j+1)-1 do k=a%icp(j),a%icp(j+1)-1
acc = acc + abs(a%val(k)) acc = acc + abs(a%val(k))
end do end do
@ -1187,19 +1200,17 @@ subroutine psb_z_csc_colsum(d,a)
end if end if
do i = 1, a%get_ncols() do i = 1, a%get_ncols()
if (a%is_unit()) then
d(i) = zone
else
d(i) = zzero d(i) = zzero
end if
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
d(i) = d(i) + (a%val(j)) d(i) = d(i) + (a%val(j))
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_ncols()
d(i) = d(i) + zone
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1242,13 +1253,18 @@ subroutine psb_z_csc_aclsum(d,a)
do i = 1, a%get_ncols() do i = 1, a%get_ncols()
if (a%is_unit()) then
d(i) = done
else
d(i) = dzero d(i) = dzero
end if
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
d(i) = d(i) + abs(a%val(j)) d(i) = d(i) + abs(a%val(j))
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, a%get_ncols() do i=1, a%get_ncols()
d(i) = d(i) + done d(i) = d(i) + done
end do end do
@ -1295,7 +1311,11 @@ subroutine psb_z_csc_rowsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = zone
else
d = zzero d = zzero
end if
do i=1, m do i=1, m
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
@ -1304,13 +1324,6 @@ subroutine psb_z_csc_rowsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_nrows()
d(i) = d(i) + zone
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1352,7 +1365,12 @@ subroutine psb_z_csc_arwsum(d,a)
goto 9999 goto 9999
end if end if
if (a%is_unit()) then
d = done
else
d = dzero d = dzero
end if
do i=1, m do i=1, m
do j=a%icp(i),a%icp(i+1)-1 do j=a%icp(i),a%icp(i+1)-1
k = a%ia(j) k = a%ia(j)
@ -1360,13 +1378,6 @@ subroutine psb_z_csc_arwsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then
do i=1, a%get_nrows()
d(i) = d(i) + done
end do
end if
return
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -1408,7 +1419,7 @@ subroutine psb_z_csc_get_diag(a,d,info)
end if end if
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = zone d(1:mnm) = zone
else else
do i=1, mnm do i=1, mnm
@ -1450,6 +1461,7 @@ subroutine psb_z_csc_scal(d,a,info,side)
character, intent(in), optional :: side character, intent(in), optional :: side
integer(psb_ipk_) :: err_act,mnm, i, j, n integer(psb_ipk_) :: err_act,mnm, i, j, n
type(psb_z_coo_sparse_mat) :: tmp
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='scal' character(len=20) :: name='scal'
character :: side_ character :: side_
@ -1464,6 +1476,10 @@ subroutine psb_z_csc_scal(d,a,info,side)
side_ = psb_toupper(side) side_ = psb_toupper(side)
end if end if
if (a%is_unit()) then
call a%add_unit_diag()
end if
left = (side_ == 'L') left = (side_ == 'L')
if (left) then if (left) then
@ -1524,6 +1540,9 @@ subroutine psb_z_csc_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d

21
base/serial/impl/psb_z_csr_impl.f90
Unescape Escape View File

@ -380,7 +380,7 @@ contains
endif endif
if (is_triangle.and.is_unit) then if (is_unit) then
do i=1, min(m,n) do i=1, min(m,n)
y(i) = y(i) + alpha*x(i) y(i) = y(i) + alpha*x(i)
end do end do
@ -733,7 +733,7 @@ contains
endif endif
if (is_triangle.and.is_unit) then if (is_unit) then
do i=1, min(m,n) do i=1, min(m,n)
y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc) y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc)
end do end do
@ -1398,7 +1398,7 @@ subroutine psb_z_csr_rowsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, m do i=1, m
d(i) = d(i) + zone d(i) = d(i) + zone
end do end do
@ -1453,7 +1453,7 @@ subroutine psb_z_csr_arwsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, m do i=1, m
d(i) = d(i) + done d(i) = d(i) + done
end do end do
@ -1509,7 +1509,7 @@ subroutine psb_z_csr_colsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, n do i=1, n
d(i) = d(i) + zone d(i) = d(i) + zone
end do end do
@ -1566,7 +1566,7 @@ subroutine psb_z_csr_aclsum(d,a)
end do end do
end do end do
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
do i=1, n do i=1, n
d(i) = d(i) + done d(i) = d(i) + done
end do end do
@ -1613,7 +1613,7 @@ subroutine psb_z_csr_get_diag(a,d,info)
end if end if
if (a%is_triangle().and.a%is_unit()) then if (a%is_unit()) then
d(1:mnm) = zone d(1:mnm) = zone
else else
do i=1, mnm do i=1, mnm
@ -1665,6 +1665,10 @@ subroutine psb_z_csr_scal(d,a,info,side)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
side_ = 'L' side_ = 'L'
if (present(side)) then if (present(side)) then
side_ = psb_toupper(side) side_ = psb_toupper(side)
@ -1734,6 +1738,9 @@ subroutine psb_z_csr_scals(d,a,info)
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (a%is_unit()) then
call a%add_unit_diag()
end if
do i=1,a%get_nzeros() do i=1,a%get_nzeros()
a%val(i) = a%val(i) * d a%val(i) = a%val(i) * d

Write Preview
Loading…
Cancel
Save