! 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 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