Efficient computation of a matrix times its transpose

Question

Let m be a huge matrix and say I want to calculate m.Transpose[m].

The result will be a symmetric matrix, and so it feels like the standard matrix multiplication might just do too much work (or does Mathematica recognize such a product and use an optimized algorithm?). The result from m.Transpose[m] would also be a standard matrix, while a SymmetrizedArray would probably be more suitable (at least memory-wise). So, therefore, the following

Question: What would be the best way to calculate a product of a matrix with its own transpose, both memory-wise and time-wise?

---

EDIT: I tried the SYRK solution and 'Benchmarked' it (see code below). The code was run on a linux machine where Mathematica had privileged access to 4 cores (via https://unix.stackexchange.com/questions/326579/how-to-ensure-exclusive-cpu-availability-for-a-running-process)

The results are the following: |Null| Integer m| Integer Packed m| Float m| Float Packed m| |-|-|-|-|-| |m.Transpose[m]| 8.62596| 8.45888| 5.18896| 5.38446| |SYRK Integer b| 0.055395| 0.057715| 0.055435| 0.055741| |SYRK Integer Packed b| 0.200872| 0.194433| 0.195106| 0.195942| |SYRK Float b| 0.055659| 0.055549| 0.055976| 0.05589| |SYRK Float Packed b| 0.196128| 0.193704| 0.192634| 0.195042|

Conclusions:

• SYRK is indeed much faster but one should take care not to pack the array to which the result is written.
• It seems SYRK doesn't care whether m is packed
• It seems SYRK doesn't care whether any of the arrays has integers or floats

These conclusions are in line with the answer by user293787.

(* For PackedArray and BLAS functionality *)
<<Developer`
<<LinearAlgebra`BLAS`

(* Matrix Dimensions *)
{ d1, d2 } = { 40000, 10 };

(* Number of random matrices for which to test *)
nLoops = 10;

(* Generate tables of timings *)
tables = 
Reap[ # ][[2,1]]& @
Do[ (* loop over Seeds and Sow tables of timings *)
   SeedRandom[s];
   (* Define the matrices for multiplication *)   
   mPacked =      RandomInteger[ { -100, 100 }, { d1, d2 } ];       
   m =            FromPackedArray @ mPacked;   
   mFloat =       FromPackedArray @ N @ m; 
   mFloatPacked = ToPackedArray @ mFloat;
   testMats =     { m, mPacked, mFloat, mFloatPacked };

   (* Define the matrices to be filled. Use Hold to prevent memory overflow *)
   integerPackedB = Hold @ ConstantArray[ 0, { d1, d1 } ]; 
   integerB =       Hold @ FromPackedArray @ ConstantArray[ 0, { d1, d1 } ];
   floatPackedB =   Hold @ ConstantArray[ 0., { d1, d1 } ];
   floatB =         Hold @ FromPackedArray[ConstantArray[ 0., { d1, d1 } ] ];
   testBMats =      Hold @ { integerB, integerPackedB, floatB, floatPackedB };


   Do[ (* timings for Mathematica matrix multiplication *)
      ClearSystemCache[];
      t[1,i] = First @ AbsoluteTiming[ testMats[[i]].Transpose[ testMats[[i]] ]; ],
      {i,4}
   ];

   Do[ (* timings for SYRK matrix multiplication *)
      Do[
         ClearSystemCache[];
         b = ReleaseHold[ testBMats[[1,i]] ];
         t[j+1,i] = First @ AbsoluteTiming[ SYRK[ "U", "N", 1,testMats[[j]], 0, b ] ];
         b=.;,
         {i,4}
      ],
      {j,4}
   ];
   Sow[ Array[ t, {5,4} ] ];
   ClearAll[t];
,{ s, nLoops } ];

TableForm[
   Plus @@ tables / nLoops,
   TableHeadings -> 
      { 
         { "m.Transpose[m]", "SYRK Integer b", "SYRK Integer Packed b", "SYRK Float b", "SYRK Float Packed b" }, 
         {"Integer m", "Integer Packed m","Float m", "Float Packed m"}       
      }
]
```

Answer 1

Edit and Warning: I noticed that my timing may not have been quite correct, or at least has to be interpreted carefully, please see the edit.

As @Roman points out, one should check if one really has to perform this calculation, but I will not question that here. OP mentions size $37000 \times 8$ in the comments. I will use size $10000 \times 10$ for illustration.

BLAS call. Motivated by this stackoverflow post let me use SYRK from the BLAS library, and see here for other BLAS commands in Mathematica:

(* sample data *)
n=10000; k=10; 
m=RandomReal[{-1,1},{n,k}];

(* container for the result, must be a symbol *)
out=ConstantArray[N[0],{n,n}];

(* calculate and store upper triangular part of m.Transpose[m] in out *)
LinearAlgebra`BLAS`SYRK["U","N",1,m,0,out];

Check. Beware that the code above only calculates the upper triangular part of m.Transpose[m]. So let us check that the upper triangular part of out was calculated correctly:

(* calculate in usual way *)
out2=m.Transpose[m];

(* this gives zero *)
Chop[Norm[Flatten[UpperTriangularize[out-out2]]]]

Timing.

RepeatedTiming[LinearAlgebra`BLAS`SYRK["U","N",1,m,0,out];]
(* about 0.126 seconds *)

RepeatedTiming[out2=m.Transpose[m];]
(* about 0.387 seconds *)

So there does seem to be a gain.

Timing (Edit). If I reset the out container to zero each time, but still a packed array, then the SYRK calls take longer:

Table[
  out=ConstantArray[N[0],{n,n}]; (* packed! *)
  First[AbsoluteTiming[LinearAlgebra`BLAS`SYRK["U","N",1,m,0,out];]]
,{10}]//Mean
(* about 0.410 seconds *)

Note that the timing is only for the SYRK call. I do not know why it takes longer. I can only guess that a BLAS call somehow puts out in a more BLAS-optimized state, better than a mere packed array. Perhaps someone can clarify. In any case, this may not be a problem if one simply does not reset out.

Memory. If m is square, which I understand it may not be in OP's case, and if one can live with destroying m, then I believe one can use m as the output container also, assuming m is square, as in LinearAlgebra`BLAS`SYRK["U","N",1,m,0,m]. Whether there is some way to directly get a structured array of some sort, I do not know.