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

Added new optional arguments to TRIL and TRIU. Document same.

Browse Source
This commit is contained in:
Salvatore Filippone
2017-08-09 10:43:38 +01:00
parent a0a0e96a96
commit 7218ae86d6
19 changed files with 808 additions and 346 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/serial/impl/psb_c_base_mat_impl.F90
+133 -42
View File
@@ -630,8 +630,8 @@ end subroutine psb_c_base_csclip
! this is just based on the getrow.
! If performance is critical it can be overridden.
!
subroutine psb_c_base_tril(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
subroutine psb_c_base_tril(a,l,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,u)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
@@ -639,12 +639,16 @@ subroutine psb_c_base_tril(a,b,info,&
implicit none
class(psb_c_base_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(out) :: b
class(psb_c_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_c_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_spk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
@@ -700,28 +704,70 @@ subroutine psb_c_base_tril(a,b,info,&
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
call b%fix(info)
nzout = b%get_nzeros()
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
if (present(u)) then
nzlin = l%get_nzeros() ! At this point it should be 0
call u%allocate(mb,nb,nz)
nzuin = u%get_nzeros() ! At this point it should be 0
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = ia(k)
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end do
end do
call l%set_nzeros(nzlin)
call u%set_nzeros(nzuin)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= -1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
nzin = l%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,l%ia,l%ja,l%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call l%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
end if
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_lower(.true.)
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
@@ -735,8 +781,8 @@ subroutine psb_c_base_tril(a,b,info,&
end subroutine psb_c_base_tril
subroutine psb_c_base_triu(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
subroutine psb_c_base_triu(a,u,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,l)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
@@ -744,12 +790,16 @@ subroutine psb_c_base_triu(a,b,info,&
implicit none
class(psb_c_base_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(out) :: b
class(psb_c_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_c_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_spk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
@@ -805,28 +855,69 @@ subroutine psb_c_base_triu(a,b,info,&
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
call b%fix(info)
nzout = b%get_nzeros()
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
if (present(l)) then
nzuin = u%get_nzeros() ! At this point it should be 0
call l%allocate(mb,nb,nz)
nzlin = l%get_nzeros() ! At this point it should be 0
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = ia(k)
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end do
end do
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros()
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,u%ia,u%ja,u%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call u%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
end if
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_upper(.true.)
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
base/serial/impl/psb_c_mat_impl.F90
+44 -28
View File
@@ -906,24 +906,24 @@ subroutine psb_c_csgetblk(imin,imax,a,b,info,&
end subroutine psb_c_csgetblk
subroutine psb_c_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
subroutine psb_c_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
use psb_error_mod
use psb_const_mod
use psb_c_base_mat_mod
use psb_c_mat_mod, psb_protect_name => psb_c_tril
implicit none
class(psb_cspmat_type), intent(in) :: a
class(psb_cspmat_type), intent(inout) :: b
class(psb_cspmat_type), intent(inout) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_cspmat_type), optional, intent(inout) :: u
integer(psb_ipk_) :: err_act
character(len=20) :: name='tril'
logical, parameter :: debug=.false.
type(psb_c_coo_sparse_mat), allocatable :: acoo
type(psb_c_coo_sparse_mat), allocatable :: lcoo, ucoo
info = psb_success_
call psb_erractionsave(err_act)
@@ -932,23 +932,30 @@ subroutine psb_c_tril(a,b,info,diag,imin,imax,&
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
call b%free()
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
allocate(lcoo,stat=info)
call l%free()
if (present(u)) then
if (info == psb_success_) allocate(ucoo,stat=info)
call u%free()
if (info == psb_success_) call a%a%tril(lcoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,ucoo)
if (info == psb_success_) call move_alloc(ucoo,u%a)
if (info == psb_success_) call u%cscnv(info,mold=a%a)
else
info = psb_err_alloc_dealloc_
if (info == psb_success_) then
call a%a%tril(lcoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info == psb_success_) call move_alloc(lcoo,l%a)
if (info == psb_success_) call l%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
@@ -956,23 +963,24 @@ subroutine psb_c_tril(a,b,info,diag,imin,imax,&
end subroutine psb_c_tril
subroutine psb_c_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
subroutine psb_c_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
use psb_error_mod
use psb_const_mod
use psb_c_base_mat_mod
use psb_c_mat_mod, psb_protect_name => psb_c_triu
implicit none
class(psb_cspmat_type), intent(in) :: a
class(psb_cspmat_type), intent(inout) :: b
class(psb_cspmat_type), intent(inout) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_cspmat_type), optional, intent(inout) :: l
integer(psb_ipk_) :: err_act
character(len=20) :: name='triu'
logical, parameter :: debug=.false.
type(psb_c_coo_sparse_mat), allocatable :: acoo
type(psb_c_coo_sparse_mat), allocatable :: lcoo, ucoo
info = psb_success_
@@ -983,23 +991,31 @@ subroutine psb_c_triu(a,b,info,diag,imin,imax,&
goto 9999
endif
allocate(acoo,stat=info)
call b%free()
allocate(ucoo,stat=info)
call u%free()
if (info == psb_success_) then
call a%a%triu(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
if (present(l)) then
if (info == psb_success_) allocate(lcoo,stat=info)
call l%free()
if (info == psb_success_) call a%a%triu(ucoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,lcoo)
if (info == psb_success_) call move_alloc(lcoo,l%a)
if (info == psb_success_) call l%cscnv(info,mold=a%a)
else
info = psb_err_alloc_dealloc_
if (info == psb_success_) then
call a%a%triu(ucoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info == psb_success_) call move_alloc(ucoo,u%a)
if (info == psb_success_) call u%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
base/serial/impl/psb_d_base_mat_impl.F90
+133 -42
View File
@@ -630,8 +630,8 @@ end subroutine psb_d_base_csclip
! this is just based on the getrow.
! If performance is critical it can be overridden.
!
subroutine psb_d_base_tril(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
subroutine psb_d_base_tril(a,l,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,u)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
@@ -639,12 +639,16 @@ subroutine psb_d_base_tril(a,b,info,&
implicit none
class(psb_d_base_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(out) :: b
class(psb_d_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_d_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_dpk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
@@ -700,28 +704,70 @@ subroutine psb_d_base_tril(a,b,info,&
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
call b%fix(info)
nzout = b%get_nzeros()
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
if (present(u)) then
nzlin = l%get_nzeros() ! At this point it should be 0
call u%allocate(mb,nb,nz)
nzuin = u%get_nzeros() ! At this point it should be 0
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = ia(k)
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end do
end do
call l%set_nzeros(nzlin)
call u%set_nzeros(nzuin)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= -1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
nzin = l%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,l%ia,l%ja,l%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call l%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
end if
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_lower(.true.)
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
@@ -735,8 +781,8 @@ subroutine psb_d_base_tril(a,b,info,&
end subroutine psb_d_base_tril
subroutine psb_d_base_triu(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
subroutine psb_d_base_triu(a,u,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,l)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
@@ -744,12 +790,16 @@ subroutine psb_d_base_triu(a,b,info,&
implicit none
class(psb_d_base_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(out) :: b
class(psb_d_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_d_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_dpk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
@@ -805,28 +855,69 @@ subroutine psb_d_base_triu(a,b,info,&
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
call b%fix(info)
nzout = b%get_nzeros()
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
if (present(l)) then
nzuin = u%get_nzeros() ! At this point it should be 0
call l%allocate(mb,nb,nz)
nzlin = l%get_nzeros() ! At this point it should be 0
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = ia(k)
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end do
end do
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros()
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,u%ia,u%ja,u%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call u%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
end if
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_upper(.true.)
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
base/serial/impl/psb_d_mat_impl.F90
+44 -28
View File
@@ -906,24 +906,24 @@ subroutine psb_d_csgetblk(imin,imax,a,b,info,&
end subroutine psb_d_csgetblk
subroutine psb_d_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
subroutine psb_d_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
use psb_error_mod
use psb_const_mod
use psb_d_base_mat_mod
use psb_d_mat_mod, psb_protect_name => psb_d_tril
implicit none
class(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(inout) :: b
class(psb_dspmat_type), intent(inout) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_dspmat_type), optional, intent(inout) :: u
integer(psb_ipk_) :: err_act
character(len=20) :: name='tril'
logical, parameter :: debug=.false.
type(psb_d_coo_sparse_mat), allocatable :: acoo
type(psb_d_coo_sparse_mat), allocatable :: lcoo, ucoo
info = psb_success_
call psb_erractionsave(err_act)
@@ -932,23 +932,30 @@ subroutine psb_d_tril(a,b,info,diag,imin,imax,&
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
call b%free()
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
allocate(lcoo,stat=info)
call l%free()
if (present(u)) then
if (info == psb_success_) allocate(ucoo,stat=info)
call u%free()
if (info == psb_success_) call a%a%tril(lcoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,ucoo)
if (info == psb_success_) call move_alloc(ucoo,u%a)
if (info == psb_success_) call u%cscnv(info,mold=a%a)
else
info = psb_err_alloc_dealloc_
if (info == psb_success_) then
call a%a%tril(lcoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info == psb_success_) call move_alloc(lcoo,l%a)
if (info == psb_success_) call l%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
@@ -956,23 +963,24 @@ subroutine psb_d_tril(a,b,info,diag,imin,imax,&
end subroutine psb_d_tril
subroutine psb_d_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
subroutine psb_d_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
use psb_error_mod
use psb_const_mod
use psb_d_base_mat_mod
use psb_d_mat_mod, psb_protect_name => psb_d_triu
implicit none
class(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(inout) :: b
class(psb_dspmat_type), intent(inout) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_dspmat_type), optional, intent(inout) :: l
integer(psb_ipk_) :: err_act
character(len=20) :: name='triu'
logical, parameter :: debug=.false.
type(psb_d_coo_sparse_mat), allocatable :: acoo
type(psb_d_coo_sparse_mat), allocatable :: lcoo, ucoo
info = psb_success_
@@ -983,23 +991,31 @@ subroutine psb_d_triu(a,b,info,diag,imin,imax,&
goto 9999
endif
allocate(acoo,stat=info)
call b%free()
allocate(ucoo,stat=info)
call u%free()
if (info == psb_success_) then
call a%a%triu(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
if (present(l)) then
if (info == psb_success_) allocate(lcoo,stat=info)
call l%free()
if (info == psb_success_) call a%a%triu(ucoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,lcoo)
if (info == psb_success_) call move_alloc(lcoo,l%a)
if (info == psb_success_) call l%cscnv(info,mold=a%a)
else
info = psb_err_alloc_dealloc_
if (info == psb_success_) then
call a%a%triu(ucoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info == psb_success_) call move_alloc(ucoo,u%a)
if (info == psb_success_) call u%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
base/serial/impl/psb_s_base_mat_impl.F90
+133 -42
View File
@@ -630,8 +630,8 @@ end subroutine psb_s_base_csclip
! this is just based on the getrow.
! If performance is critical it can be overridden.
!
subroutine psb_s_base_tril(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
subroutine psb_s_base_tril(a,l,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,u)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
@@ -639,12 +639,16 @@ subroutine psb_s_base_tril(a,b,info,&
implicit none
class(psb_s_base_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(out) :: b
class(psb_s_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_s_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_spk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
@@ -700,28 +704,70 @@ subroutine psb_s_base_tril(a,b,info,&
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
call b%fix(info)
nzout = b%get_nzeros()
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
if (present(u)) then
nzlin = l%get_nzeros() ! At this point it should be 0
call u%allocate(mb,nb,nz)
nzuin = u%get_nzeros() ! At this point it should be 0
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = ia(k)
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end do
end do
call l%set_nzeros(nzlin)
call u%set_nzeros(nzuin)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= -1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
nzin = l%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,l%ia,l%ja,l%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call l%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
end if
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_lower(.true.)
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
@@ -735,8 +781,8 @@ subroutine psb_s_base_tril(a,b,info,&
end subroutine psb_s_base_tril
subroutine psb_s_base_triu(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
subroutine psb_s_base_triu(a,u,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,l)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
@@ -744,12 +790,16 @@ subroutine psb_s_base_triu(a,b,info,&
implicit none
class(psb_s_base_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(out) :: b
class(psb_s_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_s_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_spk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
@@ -805,28 +855,69 @@ subroutine psb_s_base_triu(a,b,info,&
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
call b%fix(info)
nzout = b%get_nzeros()
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
if (present(l)) then
nzuin = u%get_nzeros() ! At this point it should be 0
call l%allocate(mb,nb,nz)
nzlin = l%get_nzeros() ! At this point it should be 0
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = ia(k)
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end do
end do
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros()
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,u%ia,u%ja,u%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call u%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
end if
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_upper(.true.)
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
base/serial/impl/psb_s_mat_impl.F90
+44 -28
View File
@@ -906,24 +906,24 @@ subroutine psb_s_csgetblk(imin,imax,a,b,info,&
end subroutine psb_s_csgetblk
subroutine psb_s_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
subroutine psb_s_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
use psb_error_mod
use psb_const_mod
use psb_s_base_mat_mod
use psb_s_mat_mod, psb_protect_name => psb_s_tril
implicit none
class(psb_sspmat_type), intent(in) :: a
class(psb_sspmat_type), intent(inout) :: b
class(psb_sspmat_type), intent(inout) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_sspmat_type), optional, intent(inout) :: u
integer(psb_ipk_) :: err_act
character(len=20) :: name='tril'
logical, parameter :: debug=.false.
type(psb_s_coo_sparse_mat), allocatable :: acoo
type(psb_s_coo_sparse_mat), allocatable :: lcoo, ucoo
info = psb_success_
call psb_erractionsave(err_act)
@@ -932,23 +932,30 @@ subroutine psb_s_tril(a,b,info,diag,imin,imax,&
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
call b%free()
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
allocate(lcoo,stat=info)
call l%free()
if (present(u)) then
if (info == psb_success_) allocate(ucoo,stat=info)
call u%free()
if (info == psb_success_) call a%a%tril(lcoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,ucoo)
if (info == psb_success_) call move_alloc(ucoo,u%a)
if (info == psb_success_) call u%cscnv(info,mold=a%a)
else
info = psb_err_alloc_dealloc_
if (info == psb_success_) then
call a%a%tril(lcoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info == psb_success_) call move_alloc(lcoo,l%a)
if (info == psb_success_) call l%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
@@ -956,23 +963,24 @@ subroutine psb_s_tril(a,b,info,diag,imin,imax,&
end subroutine psb_s_tril
subroutine psb_s_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
subroutine psb_s_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
use psb_error_mod
use psb_const_mod
use psb_s_base_mat_mod
use psb_s_mat_mod, psb_protect_name => psb_s_triu
implicit none
class(psb_sspmat_type), intent(in) :: a
class(psb_sspmat_type), intent(inout) :: b
class(psb_sspmat_type), intent(inout) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_sspmat_type), optional, intent(inout) :: l
integer(psb_ipk_) :: err_act
character(len=20) :: name='triu'
logical, parameter :: debug=.false.
type(psb_s_coo_sparse_mat), allocatable :: acoo
type(psb_s_coo_sparse_mat), allocatable :: lcoo, ucoo
info = psb_success_
@@ -983,23 +991,31 @@ subroutine psb_s_triu(a,b,info,diag,imin,imax,&
goto 9999
endif
allocate(acoo,stat=info)
call b%free()
allocate(ucoo,stat=info)
call u%free()
if (info == psb_success_) then
call a%a%triu(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
if (present(l)) then
if (info == psb_success_) allocate(lcoo,stat=info)
call l%free()
if (info == psb_success_) call a%a%triu(ucoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,lcoo)
if (info == psb_success_) call move_alloc(lcoo,l%a)
if (info == psb_success_) call l%cscnv(info,mold=a%a)
else
info = psb_err_alloc_dealloc_
if (info == psb_success_) then
call a%a%triu(ucoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info == psb_success_) call move_alloc(ucoo,u%a)
if (info == psb_success_) call u%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
base/serial/impl/psb_z_base_mat_impl.F90
+133 -42
View File
@@ -630,8 +630,8 @@ end subroutine psb_z_base_csclip
! this is just based on the getrow.
! If performance is critical it can be overridden.
!
subroutine psb_z_base_tril(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
subroutine psb_z_base_tril(a,l,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,u)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
@@ -639,12 +639,16 @@ subroutine psb_z_base_tril(a,b,info,&
implicit none
class(psb_z_base_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(out) :: b
class(psb_z_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_z_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_dpk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
@@ -700,28 +704,70 @@ subroutine psb_z_base_tril(a,b,info,&
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
call b%fix(info)
nzout = b%get_nzeros()
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
if (present(u)) then
nzlin = l%get_nzeros() ! At this point it should be 0
call u%allocate(mb,nb,nz)
nzuin = u%get_nzeros() ! At this point it should be 0
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = ia(k)
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end do
end do
call l%set_nzeros(nzlin)
call u%set_nzeros(nzuin)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= -1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
nzin = l%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,l%ia,l%ja,l%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call l%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
end if
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_lower(.true.)
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
@@ -735,8 +781,8 @@ subroutine psb_z_base_tril(a,b,info,&
end subroutine psb_z_base_tril
subroutine psb_z_base_triu(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
subroutine psb_z_base_triu(a,u,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,l)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
@@ -744,12 +790,16 @@ subroutine psb_z_base_triu(a,b,info,&
implicit none
class(psb_z_base_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(out) :: b
class(psb_z_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_z_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_dpk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
@@ -805,28 +855,69 @@ subroutine psb_z_base_triu(a,b,info,&
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
call b%fix(info)
nzout = b%get_nzeros()
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
if (present(l)) then
nzuin = u%get_nzeros() ! At this point it should be 0
call l%allocate(mb,nb,nz)
nzlin = l%get_nzeros() ! At this point it should be 0
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = ia(k)
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end do
end do
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros()
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,u%ia,u%ja,u%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call u%set_nzeros(nzin+nzout)
nzin = nzin+nzout
end do
end if
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_upper(.true.)
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
base/serial/impl/psb_z_mat_impl.F90
+44 -28
View File
@@ -906,24 +906,24 @@ subroutine psb_z_csgetblk(imin,imax,a,b,info,&
end subroutine psb_z_csgetblk
subroutine psb_z_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
subroutine psb_z_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
use psb_error_mod
use psb_const_mod
use psb_z_base_mat_mod
use psb_z_mat_mod, psb_protect_name => psb_z_tril
implicit none
class(psb_zspmat_type), intent(in) :: a
class(psb_zspmat_type), intent(inout) :: b
class(psb_zspmat_type), intent(inout) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_zspmat_type), optional, intent(inout) :: u
integer(psb_ipk_) :: err_act
character(len=20) :: name='tril'
logical, parameter :: debug=.false.
type(psb_z_coo_sparse_mat), allocatable :: acoo
type(psb_z_coo_sparse_mat), allocatable :: lcoo, ucoo
info = psb_success_
call psb_erractionsave(err_act)
@@ -932,23 +932,30 @@ subroutine psb_z_tril(a,b,info,diag,imin,imax,&
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
call b%free()
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
allocate(lcoo,stat=info)
call l%free()
if (present(u)) then
if (info == psb_success_) allocate(ucoo,stat=info)
call u%free()
if (info == psb_success_) call a%a%tril(lcoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,ucoo)
if (info == psb_success_) call move_alloc(ucoo,u%a)
if (info == psb_success_) call u%cscnv(info,mold=a%a)
else
info = psb_err_alloc_dealloc_
if (info == psb_success_) then
call a%a%tril(lcoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info == psb_success_) call move_alloc(lcoo,l%a)
if (info == psb_success_) call l%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
@@ -956,23 +963,24 @@ subroutine psb_z_tril(a,b,info,diag,imin,imax,&
end subroutine psb_z_tril
subroutine psb_z_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
subroutine psb_z_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
use psb_error_mod
use psb_const_mod
use psb_z_base_mat_mod
use psb_z_mat_mod, psb_protect_name => psb_z_triu
implicit none
class(psb_zspmat_type), intent(in) :: a
class(psb_zspmat_type), intent(inout) :: b
class(psb_zspmat_type), intent(inout) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_zspmat_type), optional, intent(inout) :: l
integer(psb_ipk_) :: err_act
character(len=20) :: name='triu'
logical, parameter :: debug=.false.
type(psb_z_coo_sparse_mat), allocatable :: acoo
type(psb_z_coo_sparse_mat), allocatable :: lcoo, ucoo
info = psb_success_
@@ -983,23 +991,31 @@ subroutine psb_z_triu(a,b,info,diag,imin,imax,&
goto 9999
endif
allocate(acoo,stat=info)
call b%free()
allocate(ucoo,stat=info)
call u%free()
if (info == psb_success_) then
call a%a%triu(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
if (present(l)) then
if (info == psb_success_) allocate(lcoo,stat=info)
call l%free()
if (info == psb_success_) call a%a%triu(ucoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,lcoo)
if (info == psb_success_) call move_alloc(lcoo,l%a)
if (info == psb_success_) call l%cscnv(info,mold=a%a)
else
info = psb_err_alloc_dealloc_
if (info == psb_success_) then
call a%a%triu(ucoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info == psb_success_) call move_alloc(ucoo,u%a)
if (info == psb_success_) call u%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return