Does Matlab Take Full Advantage of the BLAS Routines?

It's my understanding that Matlab uses BLAS routines under the hood for many matrix algebra operations. Is that correct?

If so, how does Matlab handle an expression like

C = alpha*A'*B' + beta*C

where A, B, and C are double, real matrices of compatible dimensions and alpha and beta are real scalars?

This assignment can be handled with a single call to LAPACK: dgemm. Is that what Matlab does? Or is the right hand side evaluated with a sequence of calls to @doc:mtimes and @doc:plus?

Based on this discussion I don't expect either transpose to be formed explicitly. Now I'm wondering about just how far the parser can take such optimizations to minimimize temporary storage and the number of operations, and the impact of such optimization on the documented order of operations.

[dpb]

"...and the impact of such optimization on the documented order of operations."

Mathworks doesn't document the operation of the JIT compiler as it can (and does) continually evolve; that's a purposeful decision so that folks don't try to write code specific for any given optimization.

Answer 1

[James Tursa]

"Does Matlab Take Full Advantage of the BLAS Routines?"

No. At least it didn't used to and I am assuming that is still the case. The parser is responsible for translating MATLAB source code into BLAS / LAPACK calls. Some syntax will be recognized and some will not. The rules for this are not published and subject to change from release to release.

"It's my understanding that Matlab uses BLAS routines under the hood for many matrix algebra operations. Is that correct?"

Yes. But as stated above, it relies on the parser recognizing certain syntax to take advantage of this, and I don't think it recognizes every possibility.

"If so, how does Matlab handle an expression like

C = alpha*A'*B' + beta*C

... I don't expect either transpose to be formed explicitly."

I don't know about that. The MATLAB * operator is left-to-right, so the alpha*A' would theoretically be done first. That means the transpose will be done explicitly. And probably the B' will be done explicitly as well. To do otherwise you would probably have to use parentheses around the (A'*B'). But ... this all depends on the parser. I don't think the parser is smart enough to recognize this pattern, but I haven't explicitly checked this either. Would have to code this up in a mex routine and look at numerical and timing results to verify all of this.

That being said, there is more going on here than you realize. From the documentation for DGEMM here:

https://netlib.org/blas/dgemm.f

* SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)

:

*> DGEMM performs one of the matrix-matrix operations

*>

*> C := alpha*op( A )*op( B ) + beta*C,

You can see that the DGEMM routine performs the operation in-place for the C matrix. MATLAB would in all likelihood not do this for the syntax you proposed because of the variable-sharing stuff going on in the background. To do so, it would have to evaluate C in real-time to see if it was shared or not and if an in-place operation was safe. If it happened to be unshared at the time you executed that line of code, then it could in theory make the DGEMM call directly. If it happened to be shared at the time you executed that line of code, then it would have to unshare it first before making the DGEMM call. Could be done, but I doubt the parser is smart enough to translate that line of code into the required code branching stuff needed to make it work properly. Because of all this, I highly doubt that MATLAB currently calls DGEMM directly for this particular line of code. A mex routine would be needed to verify what it actually does.

E.g., Windows R2021a

>> alpha = 1e300;

>> A = 1e300 * ones(2,2);

>> B = 1e-300 * ones(2,2);

>> alpha * A' * B'

ans =

Inf Inf

Inf Inf

>> alpha * (A' * B')

ans =

1.0e+300 *

2.0000 2.0000

2.0000 2.0000

It's clear that the leftmost * operation is done first. And online MATLAB:

alpha = 1e300;
A = 1e300 * ones(2,2);
B = 1e-300 * ones(2,2);
alpha * A' * B'

ans = 2×2

Inf Inf Inf Inf

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

We would expect inf for the answer if the leftmost * is done first, and that is what we got.

[John D'Errico]

Beyond the good answer here, it is also true that any answer we could give might change next year. The parser has undergone many changes with time.

[Paul]

@James Tursa,

Thanks for the answer.

"I don't know about that. ... "

As we've seen in the thread that I linked, Matlab doesn't always obey the documented order of operations. So assuming that multiplication is strictly left-to-right and relying on parentheses is at least a little bit dubious.

For example

rng(100);
alpha = rand;
n = 5e3;
A = rand(n);
B = rand(n);
B1 = B(:,1);

I would expect these expressions to be equivalent, but they aren't.

isequal(transpose(A)*B1,(A.')*B1)  % (1)

ans = logical
   0

In the second term, the parentheses are ignored, as shown by

isequal(A.'*B1,(A.')*B1)  % (2)

ans = logical
   1

Also as stated in that linked thread, the "operator" '* is treated as a "fused ... operation" presumably taking advantage of the transpose flag input to DGEMV. That might be a smart thing to do, but the documentation on Matlab's operator precedence is pretty clear that parentheses is 1, transpose is 2, and multiplication is 5, so line (1) should have evaulated to true based on that documentation.

But, if we do matrix-matrix instead of matrix-vector a different result obtains

isequal(transpose(A)*B,A.'*B)

ans = logical
   1

Is this just coincidence for these specific inputs? Or in this case is Matlab actually forming the explicit transpose first and then NOT using the transpose input to DGEMM? I have no idea. If the latter, I wonder why it's smart with DGEMV but not so much with DGEMM.

So I don't know what to rely on as far as operator precedence.

"there is more going on here than you realize."

No, I realize that dgemm works on C in-place. That was actually part of my question, i.e., would Matlab under any cirumstance take advantage of the in-place operation to minimize storage.

"... for the syntax you proposed ..."

To be clear, I'm not proposing anything. Just curious about how Matlab works because we can't always rely on the documentation for operator precedence (amazingly enough).

"I highly doubt that MATLAB currently calls DGEMM directly for this particular line of code."

So do I. I highly doubt it even calls it to evaluate the RHS into a temp variable before doing the assignment.

In this line of code (assuming A and B are matrices):

% C = alpha*A'*B' + beta*C

how many times do you think DGEMM is called and how many times do you think an explict tranpose is temporarily stored? I wonder if Matlab works differently on this line if B is a vector or if B is not transposed.

[dpb]

I Think but haven't tested that you could find out what MATLAB does with these expressions by writing a function containing such and compiling it. It will generate C calls instead of Fortran, but would translate to the library calls I believe.

I don't have the Complier SDK so don't have access to codegen. But, I thought might as well see what the AI 'bot thinks...here's it's response

"If you put C = alpha*A'*B' + beta*C; into a MATLAB function and run it through codegen, the engine collapses the entire line into a single optimized function:"

/* Portion of generated C code */
void optimizedMatrixOps(double alpha, const double A[10000], const double B[10000], 
                        double beta, double C[10000])
{
  /* 
     MATLAB passes 'CblasTrans' for both arrays.
     BLAS performs the transpose on-the-fly via cache stride manipulation, 
     saving massive memory allocation overhead.
  */
  cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, 100, 100, 100, 
              alpha, &A[0], 100, &B[0], 100, beta, &C[0], 100);

So, it claims MATLAB can and will recognize the general form. It went on at length about ways can rewrite the expression to break the computation order by requiring precompution of explicit transposes, etc., ...

One must, of course, take the above with a healthy dose of skepticism and verify it is correct before drawing a firm conclusion.

[Paul]

DISCLAIMER: This is my first foray into Matlab's code generation tools.

Sounds like you have an excellent idea. Unfortunately, I couldn't show results consistent with the AI claim (see disclaimer).

When I tried the Matlab Compiler I couldn't figuire out how to view the source code (see disclaimer).

I was able to generate see the source code using codegen in Matlab Coder. That code implemented the the matrix multiplies with a double for loop, pretty much as shown in the example in the middle of the page at Input Type Specification for Code Generation - MATLAB & Simulink (see disclaimer).

Maybe someone else with more expertise than I with the code generation tools can get further than I did (see disclaimer).

I can post the results if of interest.

[dpb]

How large were the arrays you defined? For small arrays (for an unknown definition of "small") MATLAB will unroll the loops.

I believe you also have to set some codegen directives to ensure it will generate BLAS code; I might try to download a trial copy later on.

[Paul]

The arrays were small. But the loops were not unrolled. There were two for loops.

I check the codegen doc page didn't find anything to direct the generator to generage calls to BLAS (see previous disclaimer). Would be happy to do it if someone can tell/show how or point to relevant doc page.

[dpb]

I should have added "or write explicit loops if deems function overhead exceeds gains", sorry.

See if this link helps any https://www.mathworks.com/help/coder/ug/blas-calls-in-generated-code.html

As always, I couldn't recreate the exact sequence in whch I asked AI before and didn't save anything except the code, but that version showed a way to set code generation for BLAS without the class using a set syntax but I'm too unfamiliar with codegen to remember just exactly how it was written.

It does seem as though one has to go to some report format to actually see the code; it doesn't appear there's a "-c" switch to just compie/no link used to with C/Fortran compilers.

[Paul]

Thanks for the link. I'm having trouble getting it all to work. I might try again later.

[dpb]

UPDATED 9:30 AM, 6/6/26 --dpb

It appears thad my presumption that codegen would be set up to generate code consistent with the shipped JIT code generator is not true. codegen will implement matrix operations via (nested) for loops (or inline if small enough) unless one does set up the environment to use BLAS/LAPACK and links directly. Hence, using it to find out how the JIT code generator maps a given expression isn't as straightforward as I had hoped. One can probably assume that the code generated using BLAS/LSPACK will indicate which routines are called for a given MATLAB expression, but that wouldn't necessarily be definite.

However, I did find the following recent Mathworks Staff generated article on Answers that is germane to the posed question although it doesn't provide guidance on which common matrix operations are recognized and how they are mapped to library calls. <use the BLAS and LAPACK implementations included in AMD Optimizing CPU Libraries>.

[Paul]

I should have some results to report from codegen in a bit.

In the meantime, based on @James Tursa's answer it occurred to me that we might have some other oddities:

rng(100);
alpha = rand;
n = 5e3;
A = rand(n);
B = rand(n);
B1 = B(:,1);
isequal(1*A.'*B1,A.'*B1)

ans = logical
   0

So multiplication by 1 isn't really multiplication by 1.

The second term uses the "fused operation" and doesn't explicitly compute the transpose, whereas for the first term the parser decides to take the transpose explicitly before the first multiplication.

[Paul]

codegen will call BLAS routines when used to generate a .mex file and the source code can be saved. I found someone nice enough to do just that for the following four functions:

function out = lapacktest1(alpha,A,B,beta,C) %#codegen
arguments
    alpha(1,1)     double
    A    (500,500) double
    B    (500,500) double
    beta (1,1)     double
    C    (500,500) double
end

% 1. C = alpha*A*B + beta*C;
C = alpha*A*B + beta*C;

out = C;

end

function out = lapacktest2(alpha,A,B,beta,C) %#codegen
arguments
    alpha(1,1)     double
    A    (500,500) double
    B    (500,500) double
    beta (1,1)     double
    C    (500,500) double
end

% 2. C = alpha*A.'*B.' + beta*C;
C = alpha*A.'*B.' + beta*C;

out = C;

end

function out = lapacktest3(alpha,A,beta,x,y) %#codegen
arguments
    alpha(1,1) double
    A    (500,500) double
    beta (1,1)     double
    x    (500,1)   double
    y    (500,1)   double
end

% 3. y = alpha*A*x + beta*y;
y = alpha*A*x + beta*y;

out = y;

end

function out = lapacktest4(alpha,A,beta,x,y) %#codegen
arguments
    alpha(1,1)     double
    A    (500,500) double
    beta (1,1)     double
    x    (500,1)   double
    y    (500,1)   double
end

% 4. y = alpha*A.'*x + beta*y;
y = alpha*A.'*x + beta*y;

out = y;

end

I observed the following (R2025B):

1. For all four cases the m-function and the mex function yielded identical results using isequal for one set of input data.

2. Generated code for cases 3 and 4 makes calls to DGEMM, not DGEMV. Of course DGEMM can be used for matrix-vector, but I was surprised to not see DGEMV called. Does DGEMV offer any advantage over DGEMM for matrix-vector?

3. In all four cases, the generated code call to DGEMM had alpha = 1 and beta = 0. The actual multiplications by the input arguments alpha and beta is done outside the call to DGEMM.

[dpb]

Just out of curiosity, can you attach the generated code?

I'm not terribly surprised about (2), but (3) seems like a peculiar implementation decision, indeed, particullarly the beta=0; that seems as though that would produce the wrong result unless the sum is also done internally separately which means must be.

[Paul]

I'm not comfortable posting the generated code absent some guidance from Mathworks saying that it would be o.k. to do so.

For case 1, the basic flow is:

temp = alpha*A

call dgemm to compute: out = temp*B

compute out = out + beta*C

[dpb]

OK, hadn't really thought about a potential IP issue...

That code structure seems bizarre, indeed, and if that really does duplicate what the JIT code generator does, the answer to the Q? as posed would be "No".and perhaps "No, not even close"

It would then be interesting to see what a mex routine to call native LAPACK would compare performance-wise. I don't write C and presently don't have Fortran reinstalled after rebuilding this machine after a drive failure...

[Paul]

Why does that code structure seem bizarre?

[dpb]

Going to the trouble to set up to call DGEMM and then externally precompute the alpha scaling in order to pass it artificial constants as arguments and then have to go compute externally what you kept it from doing internally really seems bizarre to me.

If you can (and do) pass 1 and 0 for alpha and beta, why can't/doesn't it pass the actual application constant values? Then the psuedo code would be

call dgemm to compute: out = alpha*A*B + beta*C

as your initial Q? notes it could be done. Would have to parse and pass the transpose flag status for A,B or does (2) also do the transpose outside first so the (op) flags are always false for transpose?

I suppose it may have something to do with MATLAB not being typed so that the "constants" aren't but 1x1 arrays--but certainly it's possible to pass by value; just the code generator apparently doesn't know how to do that.

[Paul]

We don't know that the code generator actually recognizes the pattern that could be implemented as a single call to DGEMM but chooses not to do so.

I should have have said before that for case2 the transb op flag is set to 'T'.

So case 2 is basically:

temp = alpha*transpose(A)
out = dgemm on temp and B with TransA = 'N' and TransB = 'T'
out = out + beta*C

So it looks like the code generator is only using DGEMM to implement some, but not all, cases of mtimes, *, and we can infer that the op flags will be set for either or both arguments as needed.

AFAICT, the implementations of all four cases respect operator precedence, which is a good thing. Maybe that was also a motivating factor at one point, though we've seen at present that respect for operator precedence is not 100%.

[dpb]

No, we don't have any way to know what JIT does understand, but to me it seems highly unlikely it recognized this pattern and then spat out the code it did. I'd guess it's a better bet it just uses Level 3 BLAS/LAPACK for mtimes and that's it everywhere(*) and at the present more complex expressions such as these simply aren't recognized but the hooks are there for future extensions because Level 3 matrix-matrix operations were written to be generic and so have this prototype, not just a simple matrix multiply routine alone.

If I find some time I may go poking at the mex sample code and see if there's a BLAS/LAPACK example to start from and see if could compare performance. I never was much good in C and it's now been well over 20 years since did any at all.

Maybe the precedence thingie is the reason for step 1 but I can't see any logic at all in the arbitrary setting of beta to zero unless the above supposition is true and the extended capabilites are never tapped.

(*) With the known but undocumented internal heuristics that decree when it will either inline or use explicit loops instead of the BLAS overhead call, of course.

[James Tursa]

I'm just getting caught back up on this thread today. Will look over all of this to see if I have any additional comments.

I will point out that in the beta = 0 case, I have never been able to uncover any official rules for this. Is the library allowed to test for beta == 0 and skip the multiply in this case and just do a straight assignment? Or is it forced to do the 0 multiply and add? E.g., you would get a different answer if there are any NaN or inf values in C. Some time ago I looked into the original BLAS interface definitions on one of the academic journals (I don't recall which one), and they didn't comment on this one way or the other. Frankly, a library should state in their user's guide how they deal with this case so that the user knows, but I have never seen that either.

[Paul]

My reading of the DGEMM source code at netlib.org is that it takes shortcuts if beta = 0.

If alpha = beta = 0, then assign C(i,j) = 0 and then return to caller

If alpha ~= 0 and beta = 0, then assign C(i,j) = 0 and proceed with the rest of the algorithm.

I don't see any case where beta*C is performed if beta = 0.

If I'm reading that correctly, then that's a good reason for Matlab to not use a single DGEMM call (if Matlab is usin the netlib implementation).

[dpb]

I don't see that makes any real difference; by the pseudo code it's already calculated alpha*A so the first case is optimal but immaterial because it doesn't get the chance.

But, setting beta=0 artificially, while it short-circuits calculating 0, the algorithm pseudo code then went ahead and did beta*C.

I suppose it boils down to what internal specifications for behavior that aren't published rule -- something like @James Tursa's note above about carrying through non-finite results that is typical MATLAB behavior rather than simply using the LAPACK implementation.

ADDENDUM

But, if this is the way LAPACK is going to be used, it appears to me would be more advantageous to just implement an optimized matrix multiply function and forget all the extra overhead...unless there is thought of eventually extending use.

ADDENDUM SECOND

By using LAPACK as shown by the pseudo code structure, alpha always =1 and so unless there's also a branch inside that skips the "multiply by one" case as well as the 0, there's an extra multiply operation introduced -- and that may be the source of the mismatch in the other thread(?)(+) --. but which, if not bypassed, would add to instead of reduce the calculation speed (relative to just a straight matrix multiply routine).

(+) Would (1*A')*1 to simulate the internal 1 after externally having calculated alpha*A make the isequal test return true, maybe?

[James Tursa]

The Fortran source code listed at netlib.org does not have a shortcut for the alpha = 1 case, surprisingly. Since that case is unambiguous (you would never get a different result by doing or not doing the multiply other than maybe a NaN payload changing), it seems they missed an optimization here. It would be easy to modify that source code for this, but the real question is what happens in actual optimized libraries such as the one MATLAB uses or Intel publishes. I would assume they take every optimization they can think of, but really there is no way to tell via argument calls since the case is unambiguous. You would have to ask the writers or disassemble the microcode.

[dpb]

I was wondering if it were possible since alpha is a variable and not a constant if it might not be possible for a rounding bit get flipped by doing another multiplication operation on load/storing from the FPP back to memory.(presuming the library doesn't have the specific optimization, either).

[James Tursa]

My understanding of IEEE floating point rules is that the primitive operations +, -, *, and / must do the operations as if there is infinite precision and then rounded according to rounding mode in effect. For a multiply by 1.0, there is no chance you will get a bit flip and a different result if this rule is followed. For library functions (e.g., sin, cos, exp, etc.) that is not the rule and one should not necessarily expect that you will always get the closest floating point result to the actual value. The library functions are typically very good and may in fact satisfy the strict rule, but they aren't held to it by the standard.

[dpb]

@Paul's example is even more puzzling, then...

""There are more things in Heaven and Earth, Horatio, than are dreamt of in your philosophy."

[James Tursa]

Which example specifically?

[dpb]

I wasn't wanting to drag that in, but it was a case from thread that evolved into a discussion of MATLAB optimization with the transpose operator <answered by Christine Tobler> as related to how MATLAB does an internal optimization to not explicitly compute the transpose in certain cases. Paul introduced an example with what seems to be unusual result of

rng(1001);
alpha = rand;
n = 5e3;
A = rand(n);
B = rand(n);
B1 = B(:,1);
isequal(1*A.'*B1,A.'*B1)

ans = logical
   0

But, I failed to see the complete follow up that it is all owing to operator precedence in conjunction with the "fused optimization" of the matrix multiply with a transpose that doesn't compute the transpose.

isequal(1*A.'*B1,(1*A.')*B1)

ans = logical
   1

isequal(1*(A.'*B1),A.'*B1)

ans = logical
   1

It's not actually the "1*" that's the culprit so it was a red herring....my bad for not having caught that earlier.