title: Which dgemm call is the fastest?
body: <p>I need to do two matrix matrix multiplications to evaluate some intermediates:</p>
  
  <p><a href="https://i.stack.imgur.com/Oll4X.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Oll4X.gif" alt="enter image description here"></a></p>
  
  <p>I can do this using multiple variants of dgemms. The set I had until now uses one <code>'N','N'</code> multiplication and one <code>'T','N'</code> (N meaning normal, T meaning transformed): </p>
  
  <pre><code>       call dgemm('N','N',Pdim,Kdim,pqdim,           &amp;
                   &amp; 1.0d0,B,Pdim,D(1,1,iamthr),pqdim, &amp;
                   &amp; 0.0d0,Etilde(1,1,iamthr),Pdim     )
  
          call dgemm('T','N',pqdim,Kdim,Pdim,              &amp;
                   &amp; 1.0d0,B,Pdim,Etilde(1,1,iamthr),Pdim, &amp;
                   &amp; 0.0d0,E(1,1,iamthr),pqdim             )
  </code></pre>
  
  <p>Where <code>Pdim</code> is the dimension of P, <code>Kdim</code> is the dimension of K and <code>pqdim</code> is the dimension of pq and ij. `Kdim' is the smallest dimension, 'Pdim' is in from a couple hundred to maybe 2000 and pqdim can be everything from 1000 to 100000.</p>
  
  <p>Now, I tried 3 versions of dgemm calls:</p>
  
  <p>a) 'N','N' + 'N','N'</p>
  
  <pre><code>    call dgemm('N','N',Pdim,Kdim,pqdim,          &amp;
               &amp; 1.0d0,B,Pdim,D(1,1,iamthr),pqdim, &amp;
               &amp; 0.0d0,Etilde(1,1,iamthr),Pdim     )
  
      call dgemm('N','N',pqdim,Kdim,Pdim,                 &amp;
               &amp; 0.50d0,Bt,pqdim,Etilde(1,1,iamthr),Pdim, &amp;
               &amp; 0.0d0,E(1,1,iamthr),pqdim                )
  </code></pre>
  
  <p>b) 'T','N' + 'T','N'           </p>
  
  <pre><code>    call dgemm('T','N',Pdim,Kdim,pqdim,            &amp;
               &amp; 1.0d0,Bt,pqdim,D(1,1,iamthr),pqdim, &amp;
               &amp; 0.0d0,Etilde(1,1,iamthr),Pdim       )
  
      call dgemm('T','N',Kdim,pqdim,Kdim,              &amp;
               &amp; 1.0d0,Etilde(1,1,iamthr),Pdim,B,Pdim, &amp;
               &amp; 0.0d0,E(1,1,iamthr),Kdim              )
  </code></pre>
  
  <p>c) 'N','N' + 'T','N' (see above)</p>
  
  <p>I don't know why, but combination c) is the fastest. This does not make sense to me. In combination a) and b) I use pre-transformed matrices B (Bt) and the leading dimension of E is different of course.</p>
  
  <p>It does not seem logical to me why c) should be fastest, because either <code>'N','N'</code> is faster than <code>'T','N'</code> or the other way around. Either way, a) or b) would have to be fastest.</p>
  
  <p>So I'm left with a couple possibilities:</p>
  
  <p>Either the compiler (ifort 19) notices two dgemms after one another and somehow magically concatenates them or since the dimensions are so vastly different that does make a huge difference. In the latter case I would still guess combination b) to be the fastest because there pqdim (the biggest dimension) is the leading dimension for both matrices...</p>
  
  <p>Or maybe I just missed something essential?</p>