1c Copyright (C) 2009-2021 The Octave Project Developers 2c 3c See the file COPYRIGHT.md in the top-level directory of this 4c distribution or <https://octave.org/copyright/>. 5c 6c This file is part of Octave. 7c 8c Octave is free software: you can redistribute it and/or modify it 9c under the terms of the GNU General Public License as published by 10c the Free Software Foundation, either version 3 of the License, or 11c (at your option) any later version. 12c 13c Octave is distributed in the hope that it will be useful, but 14c WITHOUT ANY WARRANTY; without even the implied warranty of 15c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16c GNU General Public License for more details. 17c 18c You should have received a copy of the GNU General Public License 19c along with Octave; see the file COPYING. If not, see 20c <https://www.gnu.org/licenses/>. 21c 22 subroutine ddot3(m,n,k,a,b,c) 23c purpose: a 3-dimensional dot product. 24c c = sum (a .* b, 2), where a and b are 3d arrays. 25c arguments: 26c m,n,k (in) the dimensions of a and b 27c a,b (in) double prec. input arrays of size (m,k,n) 28c c (out) double prec. output array, size (m,n) 29 integer m,n,k,i,j,l 30 double precision a(m,k,n),b(m,k,n) 31 double precision c(m,n) 32 33 double precision ddot 34 external ddot 35 36 37c quick return if possible. 38 if (m <= 0 .or. n <= 0) return 39 40 if (m == 1) then 41c the column-major case. 42 do j = 1,n 43 c(1,j) = ddot(k,a(1,1,j),1,b(1,1,j),1) 44 end do 45 else 46c We prefer performance here, because that's what we generally 47c do by default in reduction functions. Besides, the accuracy 48c of xDOT is questionable. Hence, do a cache-aligned nested loop. 49 do j = 1,n 50 do i = 1,m 51 c(i,j) = 0d0 52 end do 53 do l = 1,k 54 do i = 1,m 55 c(i,j) = c(i,j) + a(i,l,j)*b(i,l,j) 56 end do 57 end do 58 end do 59 end if 60 61 end subroutine 62