title: How to efficiently calculate matrix inner product in Fortran?
body: <p>I am trying to calculate something similar to a weighted matrix inner product in Fortran. The current script that I am using for calculating the inner product is as follows</p>
  
  <pre><code>! --&gt; In
  real(kind=8), intent(in), dimension(ni, nj, nk, nVar) :: U1, U2
  real(kind=8), intent(in), dimension(ni, nj, nk) :: intW
  
  ! --&gt; Out
  real(kind=8), intent(out) :: innerProd
  
  ! --&gt; Local
  integer :: ni, nj, nk, nVar, iVar
  
  ! --&gt; Computing inner product
  do iVar = 1, nVar
      innerProd = innerProd + sum(U1(:,:,:,iVar)*U2(:,:,:,iVar)*intW)
  enddo
  </code></pre>
  
  <p>But I found that the above script that I am currently using is not very efficient. The same operation can be performed in Python using NumPy as follows,</p>
  
  <pre><code>import numpy as np 
  import os
  
  # --&gt; Preventing numpy from multi-threading
  os.environ['OPENBLAS_NUM_THREADS'] = '1'
  os.environ['MKL_NUM_THREADS'] = '1'   
  
  innerProd = 0
  
  # --&gt; Toy matrices
  U1 = np.random.random((ni,nj,nk,nVar))
  U2 = np.random.random((ni,nj,nk,nVar))
  intW = np.random.random((ni,nj,nk))
  
  # --&gt; Reshaping 
  U1 = np.reshape(np.ravel(U1), (ni*nj*nk, nVar))
  U2 = np.reshape(np.ravel(U1), (ni*nj*nk, nVar))
  intW = np.reshape(np.ravel(intW), (ni*nj*nk))
  
  # --&gt; Calculating inner product
  for iVar in range(nVar):
      innerProd = innerProd + np.dot(U1[:, iVar], U2[:, iVar]*intW)
  </code></pre>
  
  <p>The second method using Numpy seems to be far more faster than the method using Fortran. For a specific case of <code>ni = nj = nk = nVar = 130</code>, the time taken by the two methods are as follows</p>
  
  <pre><code> fortran_time = 25.8641 s
   numpy_time = 6.8924 s
  </code></pre>
  
  <p>I tried improving my Fortran code with <code>ddot</code> from BLAS as follows,</p>
  
  <pre><code>do iVar = 1, nVar
      do k = 1, nk
          do j = 1, nj
              innerProd = innerProd + ddot(ni, U1(:,j,k,iVar), 1, U2(:,j,k,iVar)*intW(:,j,k), 1)
          enddo
      enddo
  enddo
  </code></pre>
  
  <p>But there was no considerable improvement in time. The time taken by the above method for the case of <code>ni = nj = nk = nVar = 130</code> is <code>~24s</code>. (I forgot to mention that I compiled the Fortran code with '-O2' option for optimizing the performance).</p>
  
  <p>Unfortunately, there is no BLAS function for element-wise matrix multiplication in Fortran. And I don't want to use reshape in Fortran because unlike python reshaping in Fortran will lead to copying my array to a new array leading to more RAM usage.</p>
  
  <p>Is there any way to speed up the performance in Fortran so as to get close to the performance of Numpy? </p>