-
Notifications
You must be signed in to change notification settings - Fork 3
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Misc about basic utilities (eye, diag,...) #45
Comments
Some code (to be polished) for set/get diagonals (work for arbitrary matrices -not necessarily square- and arbitrary diagonal number) !----------------------------------------
function sget_diagonal(A,k) result(d)
integer, parameter :: wp=kind(0.0)
real(wp), intent(in) :: A(:,:)
integer, intent(in) :: k
real(wp), allocatable :: d(:)
integer :: nd, i
!----------------------------------------
nd = get_diagonal_size(size(A,1),size(A,2),k)
allocate( d(nd) )
if (k >= 0) then
do i = 1, nd
d(i) = A(i,i+k)
end do
else
do i = 1, nd
d(i) = A(i-k,i)
end do
end if
end function sget_diagonal
!----------------------------------------
subroutine sset_diagonal_0(A,k,d,alpha)
integer, parameter :: wp=kind(0.0)
real(wp), intent(inout) :: A(:,:)
integer, intent(in) :: k
real(wp), intent(in) :: d
real(wp), intent(in), optional :: alpha
integer :: n, m, nd, i
real(wp) :: alpha___
!----------------------------------------
alpha___ = 0.0 ; if (present(alpha)) alpha___ = alpha
nd = get_diagonal_size(size(A,1),size(A,2),k)
if (k >= 0) then
do i = 1, nd
A(i,i+k) = alpha___*A(i,i+k) + d
end do
else
do i = 1, nd
A(i-k,i) = alpha___*A(i-k,i) + d
end do
end if
end subroutine sset_diagonal_0
!----------------------------------------
subroutine sset_diagonal_1(A,k,d,alpha)
integer, parameter :: wp=kind(0.0)
real(wp), intent(inout) :: A(:,:)
integer, intent(in) :: k
real(wp), intent(in) :: d(:)
real(wp), intent(in), optional :: alpha
integer :: n, m, nd, i
real(wp) :: alpha___
!----------------------------------------
alpha___ = 0.0 ; if (present(alpha)) alpha___ = alpha
nd = get_diagonal_size(size(A,1),size(A,2),k)
if( size(d) /= nd ) error stop "ERROR set_diagonal, non conformant shapes"
if (k >= 0) then
do i = 1, nd
A(i,i+k) = alpha___*A(i,i+k) + d(i)
end do
else
do i = 1, nd
A(i-k,i) = alpha___*A(i-k,i) + d(i)
end do
end if
end subroutine sset_diagonal_1
!----------------------------------------
integer function get_diagonal_size(n,m,k) result(nd)
!----------------------------------------
! get the size of the k-th diagonal of a n*m matrix
!----------------------------------------
integer, intent(in) :: n, m, k
!----------------------------------------
if (n >= m) then
if (k >= 0) then; nd = m - k
else if (k >= m-n) then; nd = m
else ; nd = n + k
end if
else
if (k <= 0) then; nd = n + k
else if (k <= m-n) then; nd = n
else ; nd = m - k
end if
end if
nd = max(nd,0)
end function get_diagonal_size |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Some random thoughts:
eye
is just a special case ofdiag
with a scalar, and has not a simpler interface (because most of time themold
parameter will be coded I think): is it really worth having it?I propose adding
set_diag
/get_diag
procedures, in order to manipulate arbitrary diagonals (not only the central one) of already existing matrices.Implementations detail that is not important at this point: what are the performance expectations for these utilities? I'm asking because the use of
do concurrent
constructs with branches or withmerge
are likely not well optimized by most of compilers.The text was updated successfully, but these errors were encountered: