Added implementation of A = alpha A + beta B for sparse matrices

This commit is contained in:
Cirdans-Home
2020-04-08 09:19:03 +02:00
parent dc55713541
commit 3b36c2196b
20 changed files with 972 additions and 8 deletions
+52
View File
@@ -1812,6 +1812,58 @@ subroutine psb_c_base_get_diag(a,d,info)
end subroutine psb_c_base_get_diag
subroutine psb_c_base_spaxpby(alpha,a,beta,b,info)
use psb_error_mod
use psb_const_mod
use psb_c_base_mat_mod, psb_protect_name => psb_c_base_spaxpby
complex(psb_spk_), intent(in) :: alpha
class(psb_c_base_sparse_mat), intent(inout) :: a
complex(psb_spk_), intent(in) :: beta
class(psb_c_base_sparse_mat), intent(inout) :: b
integer(psb_ipk_), intent(out) :: info
! Auxiliary
integer(psb_ipk_) :: err_act
character(len=20) :: name='spaxpby'
logical, parameter :: debug=.false.
type(psb_c_coo_sparse_mat) :: acoo
call psb_erractionsave(err_act)
if((a%get_ncols() /= b%get_ncols()).or.(a%get_nrows() /= b%get_nrows())) then
info = psb_err_from_subroutine_
call psb_errpush(info,name)
goto 9999
end if
call a%mv_to_coo(acoo,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_coo')
goto 9999
end if
call acoo%spaxpby(alpha,beta,b,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='spaxby')
goto 9999
end if
call acoo%mv_to_fmt(a,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_fmt')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_c_base_spaxpby
! == ==================================
!
+57
View File
@@ -219,6 +219,63 @@ subroutine psb_c_coo_scalplusidentity(d,a,info)
end subroutine psb_c_coo_scalplusidentity
subroutine psb_c_coo_spaxpby(alpha,a,beta,b,info)
use psb_c_base_mat_mod, psb_protect_name => psb_c_coo_spaxpby
use psb_error_mod
use psb_const_mod
implicit none
class(psb_c_coo_sparse_mat), intent(inout) :: a
class(psb_c_base_sparse_mat), intent(inout) :: b
complex(psb_spk_), intent(in) :: alpha
complex(psb_spk_), intent(in) :: beta
integer(psb_ipk_), intent(out) :: info
!Local
integer(psb_ipk_) :: err_act
character(len=20) :: name='c_coo_spaxpby'
type(psb_c_coo_sparse_mat) :: tcoo,bcoo
integer(psb_ipk_) :: nza, nzb, M, N
call psb_erractionsave(err_act)
! Copy (whatever) b format to coo
call b%cp_to_coo(bcoo,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='cp_to_coo')
goto 9999
end if
! Get information on the matrix
M = a%get_nrows()
N = a%get_ncols()
nza = a%get_nzeros()
nzb = b%get_nzeros()
! Allocate (temporary) space for the solution
call tcoo%allocate(M,N,(nza+nzb))
! Compute the sum
tcoo%ia(1:nza) = a%ia(1:nza)
tcoo%ja(1:nza) = a%ja(1:nza)
tcoo%val(1:nza) = alpha*a%val(1:nza)
tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
tcoo%val(nza+1:nza+nzb) = beta*bcoo%val(1:nzb)
! Fix the indexes
call tcoo%fix(info)
! Move to correct output format
call tcoo%mv_to_coo(a,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_coo')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_c_coo_spaxpby
subroutine psb_c_coo_reallocate_nz(nz,a)
use psb_c_base_mat_mod, psb_protect_name => psb_c_coo_reallocate_nz
+52
View File
@@ -1812,6 +1812,58 @@ subroutine psb_d_base_get_diag(a,d,info)
end subroutine psb_d_base_get_diag
subroutine psb_d_base_spaxpby(alpha,a,beta,b,info)
use psb_error_mod
use psb_const_mod
use psb_d_base_mat_mod, psb_protect_name => psb_d_base_spaxpby
real(psb_dpk_), intent(in) :: alpha
class(psb_d_base_sparse_mat), intent(inout) :: a
real(psb_dpk_), intent(in) :: beta
class(psb_d_base_sparse_mat), intent(inout) :: b
integer(psb_ipk_), intent(out) :: info
! Auxiliary
integer(psb_ipk_) :: err_act
character(len=20) :: name='spaxpby'
logical, parameter :: debug=.false.
type(psb_d_coo_sparse_mat) :: acoo
call psb_erractionsave(err_act)
if((a%get_ncols() /= b%get_ncols()).or.(a%get_nrows() /= b%get_nrows())) then
info = psb_err_from_subroutine_
call psb_errpush(info,name)
goto 9999
end if
call a%mv_to_coo(acoo,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_coo')
goto 9999
end if
call acoo%spaxpby(alpha,beta,b,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='spaxby')
goto 9999
end if
call acoo%mv_to_fmt(a,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_fmt')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_d_base_spaxpby
! == ==================================
!
+57
View File
@@ -219,6 +219,63 @@ subroutine psb_d_coo_scalplusidentity(d,a,info)
end subroutine psb_d_coo_scalplusidentity
subroutine psb_d_coo_spaxpby(alpha,a,beta,b,info)
use psb_d_base_mat_mod, psb_protect_name => psb_d_coo_spaxpby
use psb_error_mod
use psb_const_mod
implicit none
class(psb_d_coo_sparse_mat), intent(inout) :: a
class(psb_d_base_sparse_mat), intent(inout) :: b
real(psb_dpk_), intent(in) :: alpha
real(psb_dpk_), intent(in) :: beta
integer(psb_ipk_), intent(out) :: info
!Local
integer(psb_ipk_) :: err_act
character(len=20) :: name='d_coo_spaxpby'
type(psb_d_coo_sparse_mat) :: tcoo,bcoo
integer(psb_ipk_) :: nza, nzb, M, N
call psb_erractionsave(err_act)
! Copy (whatever) b format to coo
call b%cp_to_coo(bcoo,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='cp_to_coo')
goto 9999
end if
! Get information on the matrix
M = a%get_nrows()
N = a%get_ncols()
nza = a%get_nzeros()
nzb = b%get_nzeros()
! Allocate (temporary) space for the solution
call tcoo%allocate(M,N,(nza+nzb))
! Compute the sum
tcoo%ia(1:nza) = a%ia(1:nza)
tcoo%ja(1:nza) = a%ja(1:nza)
tcoo%val(1:nza) = alpha*a%val(1:nza)
tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
tcoo%val(nza+1:nza+nzb) = beta*bcoo%val(1:nzb)
! Fix the indexes
call tcoo%fix(info)
! Move to correct output format
call tcoo%mv_to_coo(a,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_coo')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_d_coo_spaxpby
subroutine psb_d_coo_reallocate_nz(nz,a)
use psb_d_base_mat_mod, psb_protect_name => psb_d_coo_reallocate_nz
+52
View File
@@ -1812,6 +1812,58 @@ subroutine psb_s_base_get_diag(a,d,info)
end subroutine psb_s_base_get_diag
subroutine psb_s_base_spaxpby(alpha,a,beta,b,info)
use psb_error_mod
use psb_const_mod
use psb_s_base_mat_mod, psb_protect_name => psb_s_base_spaxpby
real(psb_spk_), intent(in) :: alpha
class(psb_s_base_sparse_mat), intent(inout) :: a
real(psb_spk_), intent(in) :: beta
class(psb_s_base_sparse_mat), intent(inout) :: b
integer(psb_ipk_), intent(out) :: info
! Auxiliary
integer(psb_ipk_) :: err_act
character(len=20) :: name='spaxpby'
logical, parameter :: debug=.false.
type(psb_s_coo_sparse_mat) :: acoo
call psb_erractionsave(err_act)
if((a%get_ncols() /= b%get_ncols()).or.(a%get_nrows() /= b%get_nrows())) then
info = psb_err_from_subroutine_
call psb_errpush(info,name)
goto 9999
end if
call a%mv_to_coo(acoo,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_coo')
goto 9999
end if
call acoo%spaxpby(alpha,beta,b,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='spaxby')
goto 9999
end if
call acoo%mv_to_fmt(a,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_fmt')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_s_base_spaxpby
! == ==================================
!
+57
View File
@@ -219,6 +219,63 @@ subroutine psb_s_coo_scalplusidentity(d,a,info)
end subroutine psb_s_coo_scalplusidentity
subroutine psb_s_coo_spaxpby(alpha,a,beta,b,info)
use psb_s_base_mat_mod, psb_protect_name => psb_s_coo_spaxpby
use psb_error_mod
use psb_const_mod
implicit none
class(psb_s_coo_sparse_mat), intent(inout) :: a
class(psb_s_base_sparse_mat), intent(inout) :: b
real(psb_spk_), intent(in) :: alpha
real(psb_spk_), intent(in) :: beta
integer(psb_ipk_), intent(out) :: info
!Local
integer(psb_ipk_) :: err_act
character(len=20) :: name='s_coo_spaxpby'
type(psb_s_coo_sparse_mat) :: tcoo,bcoo
integer(psb_ipk_) :: nza, nzb, M, N
call psb_erractionsave(err_act)
! Copy (whatever) b format to coo
call b%cp_to_coo(bcoo,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='cp_to_coo')
goto 9999
end if
! Get information on the matrix
M = a%get_nrows()
N = a%get_ncols()
nza = a%get_nzeros()
nzb = b%get_nzeros()
! Allocate (temporary) space for the solution
call tcoo%allocate(M,N,(nza+nzb))
! Compute the sum
tcoo%ia(1:nza) = a%ia(1:nza)
tcoo%ja(1:nza) = a%ja(1:nza)
tcoo%val(1:nza) = alpha*a%val(1:nza)
tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
tcoo%val(nza+1:nza+nzb) = beta*bcoo%val(1:nzb)
! Fix the indexes
call tcoo%fix(info)
! Move to correct output format
call tcoo%mv_to_coo(a,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_coo')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_s_coo_spaxpby
subroutine psb_s_coo_reallocate_nz(nz,a)
use psb_s_base_mat_mod, psb_protect_name => psb_s_coo_reallocate_nz
+52
View File
@@ -1812,6 +1812,58 @@ subroutine psb_z_base_get_diag(a,d,info)
end subroutine psb_z_base_get_diag
subroutine psb_z_base_spaxpby(alpha,a,beta,b,info)
use psb_error_mod
use psb_const_mod
use psb_z_base_mat_mod, psb_protect_name => psb_z_base_spaxpby
complex(psb_dpk_), intent(in) :: alpha
class(psb_z_base_sparse_mat), intent(inout) :: a
complex(psb_dpk_), intent(in) :: beta
class(psb_z_base_sparse_mat), intent(inout) :: b
integer(psb_ipk_), intent(out) :: info
! Auxiliary
integer(psb_ipk_) :: err_act
character(len=20) :: name='spaxpby'
logical, parameter :: debug=.false.
type(psb_z_coo_sparse_mat) :: acoo
call psb_erractionsave(err_act)
if((a%get_ncols() /= b%get_ncols()).or.(a%get_nrows() /= b%get_nrows())) then
info = psb_err_from_subroutine_
call psb_errpush(info,name)
goto 9999
end if
call a%mv_to_coo(acoo,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_coo')
goto 9999
end if
call acoo%spaxpby(alpha,beta,b,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='spaxby')
goto 9999
end if
call acoo%mv_to_fmt(a,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_fmt')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_z_base_spaxpby
! == ==================================
!
+57
View File
@@ -219,6 +219,63 @@ subroutine psb_z_coo_scalplusidentity(d,a,info)
end subroutine psb_z_coo_scalplusidentity
subroutine psb_z_coo_spaxpby(alpha,a,beta,b,info)
use psb_z_base_mat_mod, psb_protect_name => psb_z_coo_spaxpby
use psb_error_mod
use psb_const_mod
implicit none
class(psb_z_coo_sparse_mat), intent(inout) :: a
class(psb_z_base_sparse_mat), intent(inout) :: b
complex(psb_dpk_), intent(in) :: alpha
complex(psb_dpk_), intent(in) :: beta
integer(psb_ipk_), intent(out) :: info
!Local
integer(psb_ipk_) :: err_act
character(len=20) :: name='z_coo_spaxpby'
type(psb_z_coo_sparse_mat) :: tcoo,bcoo
integer(psb_ipk_) :: nza, nzb, M, N
call psb_erractionsave(err_act)
! Copy (whatever) b format to coo
call b%cp_to_coo(bcoo,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='cp_to_coo')
goto 9999
end if
! Get information on the matrix
M = a%get_nrows()
N = a%get_ncols()
nza = a%get_nzeros()
nzb = b%get_nzeros()
! Allocate (temporary) space for the solution
call tcoo%allocate(M,N,(nza+nzb))
! Compute the sum
tcoo%ia(1:nza) = a%ia(1:nza)
tcoo%ja(1:nza) = a%ja(1:nza)
tcoo%val(1:nza) = alpha*a%val(1:nza)
tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
tcoo%val(nza+1:nza+nzb) = beta*bcoo%val(1:nzb)
! Fix the indexes
call tcoo%fix(info)
! Move to correct output format
call tcoo%mv_to_coo(a,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info,name, a_err='mv_to_coo')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_z_coo_spaxpby
subroutine psb_z_coo_reallocate_nz(nz,a)
use psb_z_base_mat_mod, psb_protect_name => psb_z_coo_reallocate_nz