Inexplicable different execution speeds when filling a matrix with complex entries

Hello everyone,

In my code, I got something disturbing regarding the execution speed. Here is the toy model where I isolated my issue:

n = 100000;
m = 10;
Me = ones(2,1); 
M = zeros(m,n);  % create a table M of length n that we partially fill with Me. The issue is the same for all m>1
count = 1; 

% I want to fill the to first rows of each column of M with coef*Me:
for i = 1 : n
    % Counter to see progress
    if i == floor(count * n/100) + 1   
        disp(['Done ', num2str(count), '/100']);
        count = count + 1;
    end
    % Filling with coef*M where coef can be either real or complex (arbitrary criterion n/2 just to highlight the issue)
    if i < n/2
        coef = 1;
    else
        coef = 1i;
    end
    M(1:2,i) = coef * Me;  % Fill in
end

The issue is the following: when , we see that the execution speed is considerably decreased, due to the fact that becomes complex.

However, if we take the opposite criterion

if i > n/2

then, the speed is constant, without issue. So why, when filling with complex numbers, do we encounter such a phenomenon direction?

I observed that when we fill with a complex, all the entries of M take the complex format "_ + _i". Hence, if we start filling with complex coefficients, all the entries of M are rewritten into the new format, whatever if we keep filling with complex or reals. But if we start with reals, then the first complex filling will retroactively convert the format of all previous entries.

If is a scalar, there is no issue.

If is complex at the beginning, then, as we would expect, it is the converse: the criterion

if i < n/2

works fine and

if i > n/2

leads to the issue.

I also tried to have when odd i and when even i (or the converse, just to switch frequently) and then, no issue. Or I think that the issue is somehow "dissolved" through the processing and we don't clearly observe the slowdown.

The fact is that my code deals with a million entries table M. It works well with the second criterion, but with the first I am finished. It seems to be very basic so I feel a bit puzzled. Any idea of what is the issue?

[Walter Roberson]

Interesting!!

Further testing shows that the speed depends upon the proportion of individual column that are being written to has a complex value with non-zero term -- even when the rows being written into are not the ones that have the complex values!

If you test the below code, and put in r = 1 instead of the random value, you can get times for the first half that exceed 1 second per iteration, but then switch to about 0.4 seconds once the complex values start getting written.

As you put in smaller r (so filling the bottom of more and more columns) the code gets faster per iteration -- even when not dealing with the columns that have had a complex value written to them. So for example for r = .9 then it will take about .9 seconds per block instead of over 1 second per iteration.

If you put in r about 1/1000 the code zips along.

n = 100000;
m = 10;
Me = ones(2,1);
Me1 = Me;
Me2 = 1i*Me;
M = zeros(m,n);  % create a table M of length n that we partially fill with Me. The issue is the same for all m>1
r = rand / 5;
fprintf('last row UNfilled fraction: %g\n', r);

M(end,floor(r*end):end) = 1i;

count = 1; 

start_time = tic;
last_loop_start_time = start_time;

% I want to fill the to first rows of each column of M with coef*Me:
for i = 1 : n
    % Counter to see progress
    if i == floor(count * n/100) + 1
        since_beginning = toc(start_time);
        since_last = toc(last_loop_start_time);
        last_loop_start_time = tic;
        fprintf('Done %d/100, this section %g seconds, total since beginning %g seconds\n', count, since_last, since_beginning);
        count = count + 1;
    end
    % Filling with coef*M where coef can be either real or complex (arbitrary criterion n/2 just to highlight the issue)
    if i < n/2
        coef = Me1;
    else
        coef = Me2;
    end
    M(1:2,i) = coef;  % Fill in
end

[david gasperini]

Hi Walter, thanks. Yes it is right.

I also tried your solution with complex(), and it always leads, when dealing with 0 entries, to 0 preallocating with "real" format, which leaves the issue unsolved.

I tried the variant complex(zeros(m,n),eps) to force the complex format while having almost zero imag. It works, but it's more a trick than something else. An also, I want something sparse after that, so before sparsifying I would have to vanish all small imag contributions from eps...

As I deal with numerous manipulation of the big matrix in my original code (with many different complex coefs to insert there or there: it is a FEM matrix, where the coefs represent different material properties)), I would like to find something straight forward for this issue

[Bruno Luong]

If the speed is important, then don't create a full temporary matrix. This is a waste of time. Create directly your sparse FEM matrix.

[david gasperini]

The fact is that I build with the full matrix all the entries of my elementary fem matrices (hence this matrix is almost full, but some coefs, maybe one half, will remain zero). Then I spread it to my sparse global matrix using sparse. It is the best way to write it, since otherwise I would have to reallocate n a loop each entry of the global sparse matrix, which is the worst thing to do regarding time

Answer 1

[Bruno Luong]

I also dig deeper into this "issue" of slowness.

The strange issue I discover is that when filling with real values into a complex array, the runtime seems to depend on the density of complex numbers in the preallocation array, as showed as this code and plot below.

The code initialize M with zero(m,n) and distribute p complex numbers. The begining and end are forced to be complexes. So the array should remain complex in the entire for-loop regardless the direction of scan/populate.

I bench for various values of the density (~1/p in the plot) and provide the result on my 2-year old laptop, R2020b update 3.

[ptab,t]= fillM(10,100000);
loglog(ptab,t)
xlabel('p number of complexes when array is initialized');
ylabel('execution time [s]')
grid on
title('real filling')

%%
function [ptab,t] = fillM(m,n)
%
dummynan = nan*1i;
ptab = 2.^(0:floor(log2(m*n)));
ptab = union(ptab,m*n); % the last point, who array is complex
t = zeros(size(ptab));
for k=1:length(ptab) % loop on density
    p = ptab(k);
    M = zeros(m,n);
    Me = ones(2,1);
    count = 1;
    
    % distribute p complex numbers in the allocate array
    inan = round(linspace(1,m*n,p));
    M(inan) = dummynan;
    
    tic
    % I want to fill the to first rows of each column of M with coef*Me:
    for i = 1 : n
        % Counter to see progress
        if i == floor(count * n/100) + 1
            disp(['Done ', num2str(count), '/100']);
            count = count + 1;
        end
        % Filling with coef*M where coef is real
        M(1:2,i) = 1*Me;  % Fill in, change 1 to 1i to switch to complex filling
    end
    % replace the dummmy nan by 0
    b = isnan(M(inan));
    M(inan(b)) = 0;
    t(k) = toc;
end % for

end

If we fill with complex

M(1:2,i) = 1i*Me; % instead of M(1:2,i) = 1*Me

the CPU time is fast and does not depend on the density of complex numbers distributed in the preallocation array.

It looks to me that when filling with REAL on COMPLEX array, the data deepcoying is (bi)multi-level. It proceed by block of around 100-1000 elements, then progressivey gather together to the entire array,

Not of my trick of filling NaN complex then undo later does not have significant speed penalty.

So there we go. It seems a candidate for workaround.

[Walter Roberson]

Bruno, is it possible that the block size is 128 ? That would correspond to 1024 bytes of real-only, 2048 bytes for inline-complex

[Bruno Luong]

Possibly Walter. 100 is really a gross estimation, 128 makes more sense.

Answer 2

[Bruno Luong]

You should preallocate M as complex

M = 1i+zeros(m,n);

if you intend to populate it with complex number later in the for-loop.

Everytime a real array becomes complex it will take time. The opposite is also true.

I think it's better keep the array complex from start to finish.

[Walter Roberson]

My tests show that it does not depend upon whether the entire array is complex or not. If it did, then

M = complex(zeros(m,n), 0);

would be enough.

[david gasperini]

Hi Bruno Luong, thanks for this answer

However, proceeding like this, we obtain all entries as 0+1i, whereas I expect the non filled entries to be kept at 0 or 0+0i.

Hence this solution would imply to rewrite all vanishing entries, which leads to significant additionnal costs

[Bruno Luong]

Walter:

"

M = complex(zeros(m,n), 0);

would be enough."

No, as soon as a real number is set to single element M becomes real. Useless.

To David : you can populate M with 1i*NaN then later replace them with 0s.

EDIT: I see James already propose this solution.

Another suggestion is populate the array with 1i*realmin and leave it as it is. The appreciation whereas 1i*realmin can be considered as practically 0s is leave to you to decide because you know how they will be used downstream.

Another option is to allocate M wih one more column with some dummy complex value, but never use it.

There is no really a good solution.

>> M=1i*realmin+zeros(5)

M =

  1.0e-307 *

   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i
   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i
   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i
   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i
   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i   0.0000 + 0.2225i

>> M(1)=1

M =

   1.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

>> 

[david gasperini]

I would say that would be ok but as I said, when sparsing it will probably consider e-307 as a non zero entry, which then would need to be replaced by zero, the same applies for the NaN solution I think: need to replace all, which is time consuming

[Bruno Luong]

That's why my comment about SPARSE. Why don't you build directly the sparse matrix?

Believe me I have build many FEM matrices unsing MATLAB and it's always a good idea to tackle directly the sparse matrix with the construction

M = sparse(i, j, s, ...)

Don't do any full temporary matrix, assign elements blocks in the loop. Otherwise you'll run into inefficient code.

[david gasperini]

Interesting

I am not sure to understand fully your approach, can I ask you to develop a bit?

Basically with a loop you compute the full elementary ones for each element, so how do you spread it over the sparse global one according to the nodes, without rewritting it all the time? (due to sparse format)

[Bruno Luong]

Building FEM matrix;

1. Preallocate long vectors i, j, s depending on number of elements (and the FEM order)
2. you loop on elements (or phi_k the test functions, which is, numbered by nodes k) and populate i, j, s; where i, j are node numbers and s is integral A(phi_i,phi_j) ds; A is your FEM operator.
3. You call SPARSE once with the syntax

S = sparse(i,j,s,...)

[david gasperini]

Ah yes exactly, i and j gives the position where to spread (so local to global for the nodes) and s my table M of all elementary matrices , right?

If I well understood, it is already how I proceed:

M = zeros(9,msh.nb_TRI); % Initialize tables of elementary matrices

for i = 1 : msh.nb_TRI
    Me = getMassElem(i); % Get the elementary for element i
    M(1 : 9,i) = Me;    % Load it in the table
end
I_glob = msh.TRIANGLES(:,[1 2 3 1 2 3 1 2 3]).';  % Connectivity list of triangles, rows
J_glob = msh.TRIANGLES(:,[1 1 1 2 2 2 3 3 3]).';  % Connectivity list of triangles, columns

M_glob = sparse(I_glob,J_glob,M,msh.nbNod,msh.nbNod); % Build global

Is it what you mean? I precise that I cannot vectorize everything since I have to consider some complex coefs that have to be computed for each element, in a way that I cannot avoid the loop. And the 's, which are sometimes complex, will be multiply, for i, to the corresponding .

The fact is that in reality, I have for example up to N "other" sparse mass matrices where the elementary ones are differently coefficiented (due to a coupling in my PDE). Then instead of creating N tables M, I create a big one of size (9N,n) and use sparsify on a part of M, like for instance for the 2nd table which gives the second mass :

M_glob2 = sparse(I_glob,J_glob,M(1+9*2 : 9*3,:))

N could be 3 or either 20, so in this last case it is too large and Walter's technique does not work

[Bruno Luong]

... M(1+9*2 : 9*3,:) ...

Accessing row wise is inefficient. Better working with the transpose your matrix M.

I just can't believe you work with such inefficient accessing and yet complain about extra cpu time to filter out NaN entry.

And again as long as you have as many as half elements that is 0 in your M(:) data, the matrix assembly method seems wasteful. You should be able not storing unnessary any 0 from the topology of your mesh. It seems all that is links to the fact that you put the mass matrix elements in some big intermediate matrix rather than work them as separate vectors.

[david gasperini]

I see your point. Indeed, it would be better to create N tables for each of the N global matrices, and fit them to avoid any 0. And then we don't use row wise access, right?

Then with such smaller tables, each of the techniques work.

But N is an input parameter, and isn't fixed. So should I create a cell array like

M = cell(N)

and then treat each of the N cells as my mass tables of elementary matrices, with adapted size to avoid the 0's?

[Bruno Luong]

Yes I would use CELL.

But only you know the characteristics of your mesh/equation.

What is N? The number of subdomains?

[david gasperini]

Roughly it is a coupled system of N PDE with N corresponding mass matrices. So for each element I compute the elementary mass, apply the good coef and load it in the corresponding rows 1+9i : 9(i+1) (associated to the mass of equation i) of the corresponding column of the table M

[Bruno Luong]

So if I understand correctly the 0s in M are due to PDEs that are not directly coupled?

[david gasperini]

Yes, exactly

Answer 3

[James Tursa]

When an array switches from real to complex in R2018a or later, a new memory allocation and deep data copy must occur to change the data from packed real-only to interleaved complex. So, extra memory gets allocated and all of the current values must get copied over as well. If you start the variable as complex to begin with, as Bruno suggests, this extra allocation and deep copy can be avoided. I would need to research what the fastest way to pre-allocate a complex variable is. Basically need to find a syntax that the MATLAB parser understands to avoid unnecessary copies etc.

I don't know the behind-the-scenes rules for when MATLAB checks to see if a variable should switch from one format to the other, or even the method/order it uses for the checking. Maybe by forcing a complex number to be in the first and last spots until the very end of your algorithm you might be able to short-circuit this checking midstream and avoid midstream changes. E.g.,

>> tic;C = 1i+zeros(m);toc
Elapsed time is 0.070991 seconds.
>> tic;D = zeros(m,m,'like',1i);toc
Elapsed time is 0.008061 seconds.
>> tic;clear M;M(m,m) = 1i;toc
Elapsed time is 0.304390 seconds.

C has complex 1i in every spot, but D and M do not. If we then do this:

>> tic;C(1,1)=0;toc % C remains complex
Elapsed time is 0.000113 seconds.
>> tic;D(1,1)=0;toc % D gets turned into packed real-only, so deep data copy
Elapsed time is 0.320612 seconds.

If we get rid of all of the complex values except the last one in C, watch what happens when we get rid of the last one:

>> tic;C(1:end-1)=0;toc % the last element is still complex
Elapsed time is 0.080873 seconds.
>> tic;C(end)=0;toc % the last element (and all of C) is now pure real
Elapsed time is 0.118038 seconds.

That last action on C(end) caused a deep data copy when it got converted from complex to pure real. There may or may not be an additional memory allocation when shrinking from complex to pure real ... again these rules are not published.

[david gasperini]

Yes, but as said, Bruno's code leads to non zeros entries where I expect vanishing ones

Your solution using zeros(m,m,'like',1i) works pretty well, thanks!

However I don't see how the other tests with C solve the issue since in C we don't have vanishing imaginary parts

And your right: when changing an entry of D with a real number a, all become in real format, even using "complex(a,0). Sh_t!

[david gasperini]

So it seems that we don't have a clear solution. Do you think we can call the helpdesk?

[James Tursa]

I suppose if you want 0's left over from this process, you can start them with complex NaN and then turn those spots to 0's afterwards. E.g.,

C = NaN*1i + zeros(m,m); % start off with complex NaN's
% run your algorithm
C(isnan(C)) = 0; % replace any leftover NaN spots with 0's

[Walter Roberson]

The results of my tests say that if you start with an all-zero matrix but pad it with a single row of imaginary values, then the calculation is fast. You could then sparse() all except the last row.

[david gasperini]

Yes but it only works if the size is not to large since it covers only a block around the entry, and not the whole matrix

[Bruno Luong]

Not sure what you are talking here.

MATLAB matrix is either real or complex. It is complex where there is at least one element with non-zero imaginary part (or by complex function, but change dynamically with any sucessively operation).

There is no such thing of "covers a block".

Please illustrate with a concrete small example.

[david gasperini]

My english might be unaccurate I know

I mean:

>> M=zeros(3,3)

M =

     0     0     0     
     0     0     0     
     0     0     0
     
>> M(2,2)=1i

M =

   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   
   0.0000 + 0.0000i   0.0000 + 1.0000i   0.0000 + 0.0000i   
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i      

and conversely, from complex to real. This is ok, but it does not convert into complex the whole matrix, but only a box of size around 50. And since M is smaller than 50, it is fully converted. But on this example:

   M=zeros(3,1000);
   M(2,500)=1i;

then you can see that:

• on the command window it appears to be full complex, which is nice, but
• on the variable window you see that it is in complex format only for the columns 449 to 512

and the execution speed shows that indeed, it is not fully complex and as you said we have deep reallocating and copying

[Bruno Luong]

Don't get fooled by what you see in the display screen or variable window. The matrix is either complex or real, the WHOLE thing, period.

I (and James) with mex programming many years and we pretty know how internally the data are stored.

You might also interested in reading doc of the command ISREAL that applies for the entire matrix and returns a single boolean result.

[david gasperini]

I trust you

But then I don't understand why, with large number of columns and small number of rows, when I preallocate a complex row with a complex, I get correct speed, whereas when I preallocate a column, it slows down

[Paul]

James,

Referring to the code in the original question, doesn't the deep copy from real to complex occur only once, i.e., the first time that i > n/2? But when I run the code I see it slow down for all iterations where i > n/2. So it seems that this issue relates more to how elements of M are accessed and assigned to, rather than whether or not M is preallocated real or complex?

There's a lot to digest in this thread, I'm just trying to understand it all ....

Answer 4

[david gasperini]

As far I understand, Bruno's and James's solution work, but would need to replace NaN entries by 0, or every small entries if the goal is to sparsify.

Walter's one seems to be a good trade off: my table is of large length and relatively small width so I create an additional last row with 1i, which seems to be enough to force the format along the whole matrix, then I populate as wanted, and sparsify everything but the last row

M = [zeros(m,n,'like',1i); 1i + zeros(n,1)];
%apply the fill in
% then sparsify M(1:end-1,:) according to the expected size

Here we hope it will be sufficient for the different widths of my matrices, since we see that it forces the complex format only for a certain block around the complex entry. If to large, then the same issue appears.

It seems that if we put a complex entry A at (a1,a2), it forces the format for a block like (a1-50 : a1+10, a2-50 : a2+10)

It is definitely weird

We would need something to force in the whole matrix, regardless to its size...

[Bruno Luong]

"We would need something to force in the whole matrix, regardless to its size..."

I second that. Currently there is no-way to control the complexity and leave it persistent for a given array.

The automatic data deep copying sometime is anoying, command such as COMPLEX(A) is useless due to dynamic MATLAB decision to allocate internal data as real or complex. The matter of efficiency gets worse with recent complex data interleave.

[James Tursa]

Really, the only guaranteed way to keep a variable complex when assigning elements is to essentially bypass MATLAB processing and modify variables inplace via a mex routine. Not official, not recommended, has potential side effects, etc. etc. etc.

There are plenty of unused bits left in the flags int of the mxArray for a new "persistently complex" bit. This is where the "is 2D", "is complex", "is numeric", etc. bits are stored. But it would take an enhancement request to TMW to accomplish this of course.

Answer 5

[Walter Roberson]

Test code included.

On the first line, set only_do_subset to false to run the larger suite of tests.

[1:2]no complex: real half min = 0.00167084, max = 0.0968766, mean = 0.00534692

[1:2]no complex: complex half min = 0.457809, max = 0.540868, mean = 0.486138

[1:2](1,:): real half min = 0.0214277, max = 0.539432, mean = 0.21429

[1:2](1,:): complex half min = 0.45741, max = 0.569145, mean = 0.489417

[1:2](:,1): real half min = 0.00175333, max = 0.12962, mean = 0.00614126

[1:2](:,1): complex half min = 0.00141391, max = 0.00303835, mean = 0.00224259

[1:2](end,:): real half min = 0.0014808, max = 0.090944, mean = 0.00523808

[1:2](end,:): complex half min = 0.00156304, max = 0.0023197, mean = 0.00202867

[1:2](:,end): real half min = 1.11345, max = 1.55058, mean = 1.17526

[1:2](:,end): complex half min = 0.458782, max = 0.709971, mean = 0.494252

[1](end,:): real half min = 0.00172286, max = 0.149381, mean = 0.00841403

[1](end,:): complex half min = 0.00109197, max = 0.00278375, mean = 0.00192156

[1,4](end,:): real half min = 0.0019306, max = 0.0885242, mean = 0.00633371

[1,4](end,:): complex half min = 0.00206219, max = 0.00495296, mean = 0.00286406

[1:4](end,:): real half min = 0.00207948, max = 0.0889896, mean = 0.00749856

[1:4](end,:): complex half min = 0.00206808, max = 0.0133852, mean = 0.00356371

[1:2](diagonal): real half min = 0.00146681, max = 0.0696898, mean = 0.00450054

[1:2](diagonal): complex half min = 0.00196183, max = 0.00275398, mean = 0.00212438

[1:2](end,1): real half min = 0.00156037, max = 0.0687157, mean = 0.00459935

[1:2](end,1): complex half min = 0.0014585, max = 0.00262812, mean = 0.00203858

The notation here is that the leading [] shows the row indices of each column that are written to. The part immediately after that talks which parts of M were pre-initialized to be complex-valued -- for example (end,:) meaning M(end,:)=1i was used for initialization. (diagonal) was linear (m+1:m+1:end) which runs down the +1 diagonal but wrapping as often as needed

The "real half" is the segments up to n/2 that have m1 written into the array locations. The "complex half" is the segments beyond that have m2 written into the array locations.

The "min", "max", and "mean" refer to timings in second for each of the groups of n/100. With n = 100000, each segment corresponds to 1000 columns.

Row index [1:2] corresponds to the original test M(1:2,i) = coef * Me where the first two rows of each column were written into. The subset also shows [1] and [1:4] to explore how much the number of rows written to affects the results, and shows [1,4] to compare to [1:2] to explore how much effect there is if the locations being written to are not consecutive rows.

Discussion

The worst case by far is [1:2](:,end) which is the case where each row has a complex value in it. Writing real values into slots that are 0 inside such an array is the slowest operation -- over 1 second per 1000 columns. Any discussion that claims that the writing speed should depend only on the global properties cannot account for this case.

Agggh... I just added a test that appears to account for major aspects of the time cases (though not necessarily the fine details)

After each potential assignment, MATLAB needs to test to see whether the new matrix is functionally complex-valued or not, in order to know whether to convert to real-only or complex-only. Imagine it does that by starting at the beginning of the matrix and scanning through the complex components in linear order until it encounters a non-zero complex value, or gets to the end. If it gets to the end and no values had non-zero complex component, it needs to convert the matrix to real-valued; otherwise if it was real-valued and a complex component is encountered, it needs to convert the array to complex-valued.

Now, if this is correct... then the portion of the array it has to travel through before it finds the non-zero complex value could have a major impact on the speed of the assignments.

It just so happens that the first half of the columns, David happened to write only real values. With M initialized to all 0, and with the case involving arrays both marked as real optimized, writing to the first half of the columns is fast. Then when you start writing the complex values in, each time the assignment hypothesizes that maybe the array is no longer complex valued, and starts scanning through it... and in this particular code, has to get about half way through the array before it finds a non-zero imaginary component and says "Okay, guess it is still complex-valued".

It is clear to me that there are heuristics that could be used to optimize the process. If MATLAB kept an approximate count of non-zero imaginary components, and decremented it once for each real-only value assigned in, but did not increment it for each imaginary component assigned in, and re-scanned the array each time the approximate count reached zero, then you could get away with a lot fewer checks.

My first thought had been to increase the approximate count for each complex value written in, but I realized that if you kept assigning to the same location, you could increase the count arbitrarily; if it was the only location that was complex valued then even just a single assignment of a real value could turn the entire array real-valued. But if you know you have 10 non-zero components imaginary components, and you have only written 5 real-only components so far, then no matter whether you wrote those to 5 different locations or to the same location, you would know logically that there must be at least 5 non-zero components out there somewhere. You might have more than 5 -- A(1)=0;A(1)=1i repeated 5 times would leave either 10 or 11 locations with non-zero components, so maybe there are more efficient ways that do not involve having to track each individual location, but this heuristic should be fairly cheap.

[Paul]

Wouldn't it be quite surprising if Matlab actually scanned through all of the imaginary components either before or after the assignment to determine if the result needs to be converted from complex to double? I mean, wouldn't it be better to keep track of whether or not the imaginary part has any non-zero element as the matrix is created and then operated on?

On a related note, is there any way to make a complex matrix with zero imaginay part other than

M = complex(realpart,0)
M = complex(realpart)

In other words, if M already exists and is complex is there any way to force the imaginary part of M to zero but have M maintain its complex attribute w/o rewiting M?