fortran-lapack

fortran-lapack copied to clipboard

Some random thoughts:

eye is just a special case of diag with a scalar, and has not a simpler interface (because most of time the mold 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 with merge are likely not well optimized by most of compilers.

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