! Copyright 2019 Wojciech Kosior
  
! This is free and unencumbered software released into the public domain.

! Anyone is free to copy, modify, publish, use, compile, sell, or
! distribute this software, either in source code form or as a compiled
! binary, for any purpose, commercial or non-commercial, and by any
! means.

! In jurisdictions that recognize copyright laws, the author or authors
! of this software dedicate any and all copyright interest in the
! software to the public domain. We make this dedication for the benefit
! of the public at large and to the detriment of our heirs and
! successors. We intend this dedication to be an overt act of
! relinquishment in perpetuity of all present and future rights to this
! software under copyright law.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
! IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
! OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
! ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
! OTHER DEALINGS IN THE SOFTWARE.

! For more information, please refer to <http://unlicense.org/>

MODULE blockmat
  IMPLICIT none
  PRIVATE
  
  PUBLIC :: blockmull
  PRIVATE :: blockmull_4, blockmull_8, blockmull_16
  
  INTERFACE blockmull
     procedure blockmull_4, blockmull_8, blockmull_16
  END INTERFACE blockmull

CONTAINS

  FUNCTION blockmull_4(A, B) result(C)
    IMPLICIT none
    real(kind=4), intent(in), dimension(1:,1:) :: A, B
    real(kind=4), dimension(size(A, 1), size(B, 2)) :: C
    integer :: i, ii, ii_max, iiend, j, k, kk, kk_max, kkend
    integer, parameter :: blocksize = 400

    C = 0
    
    kk_max = (size(A, 2) / blocksize) * blocksize + 1
    ii_max = (size(A, 1) / blocksize) * blocksize + 1
    
    DO kk = 1, kk_max, blocksize
       IF (kk .lt. kk_max) kkend = kk + blocksize - 1
       IF (kk .eq. kk_max) THEN
          kkend = size(A, 2)
          IF (kk .gt. kkend) EXIT
       END IF

       DO ii = 1, ii_max, blocksize
          IF (ii .lt. ii_max) iiend = ii + blocksize - 1
          IF (ii .eq. ii_max) THEN
             iiend = size(A, 1)
             IF (ii .gt. iiend) EXIT
          END IF

          DO j = 1, size(B, 2)
             DO k = kk, kkend

                C(ii:iiend,j) = C(ii:iiend,j) + A(ii:iiend,k) * B(k,j)
             END DO
          END DO
       END DO
    END DO

  END FUNCTION blockmull_4  

  FUNCTION blockmull_8(A, B) result(C)
    IMPLICIT none
    real(kind=8), intent(in), dimension(1:,1:) :: A, B
    real(kind=8), dimension(size(A, 1), size(B, 2)) :: C
    integer :: i, ii, ii_max, iiend, j, k, kk, kk_max, kkend
    integer, parameter :: blocksize = 44

    C = 0
    
    kk_max = (size(A, 2) / blocksize) * blocksize + 1
    ii_max = (size(A, 1) / blocksize) * blocksize + 1
    
    DO kk = 1, kk_max, blocksize
       IF (kk .lt. kk_max) kkend = kk + blocksize - 1
       IF (kk .eq. kk_max) THEN
          kkend = size(A, 2)
          IF (kk .gt. kkend) EXIT
       END IF

       DO ii = 1, ii_max, blocksize
          IF (ii .lt. ii_max) iiend = ii + blocksize - 1
          IF (ii .eq. ii_max) THEN
             iiend = size(A, 1)
             IF (ii .gt. iiend) EXIT
          END IF

          DO j = 1, size(B, 2)
             DO k = kk, kkend

                C(ii:iiend,j) = C(ii:iiend,j) + A(ii:iiend,k) * B(k,j)
             END DO
          END DO
       END DO
    END DO

  END FUNCTION blockmull_8  

  FUNCTION blockmull_16(A, B) result(C)
    IMPLICIT none
    real(kind=16), intent(in), dimension(1:,1:) :: A, B
    real(kind=16), dimension(size(A, 1), size(B, 2)) :: C
    integer :: i, ii, ii_max, iiend, j, k, kk, kk_max, kkend
    integer, parameter :: blocksize = 28

    C = 0
    
    kk_max = (size(A, 2) / blocksize) * blocksize + 1
    ii_max = (size(A, 1) / blocksize) * blocksize + 1
    
    DO kk = 1, kk_max, blocksize
       IF (kk .lt. kk_max) kkend = kk + blocksize - 1
       IF (kk .eq. kk_max) THEN
          kkend = size(A, 2)
          IF (kk .gt. kkend) EXIT
       END IF
       
       DO ii = 1, ii_max, blocksize
          IF (ii .lt. ii_max) iiend = ii + blocksize - 1
          IF (ii .eq. ii_max) THEN
             iiend = size(A, 1)
             IF (ii .gt. iiend) EXIT
          END IF
          
          DO j = 1, size(B, 2)
             DO k = kk, kkend

                C(ii:iiend,j) = C(ii:iiend,j) + A(ii:iiend,k) * B(k,j)
             END DO
          END DO
       END DO
    END DO

  END FUNCTION blockmull_16  

END MODULE blockmat