Yes, unfortunately, this is expected behavior, at least for now. MATLAB doesn't allow overloading sparse data types, so we must use workarounds to add our own sparse matrix type. This comes with copy overhead in each operation (even in indexed access).

That is why, for maximum speed, please try to avoid element-wise access operations for sparse matrices (assignment/reading).

For example, when assembling sparse matrix, generate arrays of non-zeros and their indices separately, and then convert them into sparse matrix using command "sparse".

***

Starting from September, we have been working on new ideas on how to reduce copy-overhead further (now we are using undocumented functions in MATLAB, etc.). I hope this will allow us to reduce timings in future versions. 

Not directly, but you can convert required number of bits to decimal digits using simple formula: 

>> d = floor(bits*log10(2));  % precision in decimal digits
>> mp.Digits(d)

This is confirmed bug in MATLAB R2019b, see bug-report: https://www.mathworks.com/support/bugreports/details/2088279 

I have just updated toolbox bundle for GNU Linux to cope with this mess from TMW. Please re-download toolbox one more time.

Now it includes only proper version of cross, compatible with pre-R2019b and post-R2019b Update 2.

And it is removed from mp.Test to avoid confusion with R2019b.

Oh, I see the difference now, thank you (never check mails on the go...).

So, TWM added conjugation to cross function in R2019b. 

In Update 3 they reverted it back to pre-R2019b state (without conjugation).

Will add workaround in next version of toolbox.

Thank you for the report.

Strange, the "cross.m" was not changed in Update 3. What version of toolbox do you use on Windows? 

Please update it if it is too old.

Looks like TWM just fixed the bug they introduced in initial release of R2019b.

While I am installing it, could you please share code of cross.m distributed with Update 3 (email or here)?

Any follow up on this? 

Thank you for your response.

I have done some quick tests using ill-conditioned problems from NLEVP library. The 'polyeig' provided in our toolbox works quite well on them.

I wrote simple function 'test_polyeig' which returns absolute error of eigen-decomposition for the the problem specified (see its full code at the bottom) for 'double' and extended precision ('mp') cases:

>> test_polyeig('power_plant','double')

ans =

   1.8567e+05

>> test_polyeig('power_plant','mp')

ans = 

    4.091947679918651639316533910037398e-13


>> test_polyeig('cd_player','double')

ans =

   1.1410e+04

>> test_polyeig('cd_player','mp')

ans = 

    4.953565354331632662743256645399015e-15

>> test_polyeig('damped_beam','double')  % use 'double' precision

ans =

   34.2276

>> test_polyeig('damped_beam','mp') % use extended precision

ans = 

    7.547834125835157737664345305865886e-17

Our 'polyeig' gives much more accurate result out-of-the-box. The only possible downside compared to 'quadeig' is that 'polyeig' is (much) slower. But it is perfectly usable for research right now. 

Could you please justify why you think "...without 'quadeig' is supported, the toolbox is scarcely useful for a large number of people who work on quadratic eigenvalue problems." Quadratic (and other) eigenproblems can be handled by toolbox right now using 'polyeig'. Or am I missing something? 

The code of 'test_polyeig' below. It supports 'double' and 'mp' computations:

function abserror = test_polyeig(name,float_t,varargin)

[COEFFS,~] = nlevp(name,varargin{:});

p = length(COEFFS);
A = cell(p,1);
for k = 1:p
    A{k} = cast(COEFFS{k},float_t);
end

n = size(A{1},1);

[X,E] = polyeig(A{:});
error = zeros(n*(p-1),float_t);

for j = 1:n*(p-1)
   lambda = E(j);
   x = X(:,j);
   
   if lambda ~=Inf
      T = A{1};
      for k = 2:p
          T = T + lambda^(k-1) * A{k};
      end
      
      error(j) = norm(T*x,Inf);
   else
      error(j) = norm(A{1}*x,Inf);       
   end
end

abserror = norm(error,Inf)

P.S.

To get even higher accuracy, NLEVP code must be converted to support 'mp' computations.

Because right now, NLEVP computes matrices in 'double' precision and then we use 'mp' to compute eigenvalues.

If matrices would be computed with high precision from the onset, than overall accuracy would be higher.

Toolbox uses binary representation of floating-point numbers. 

When we extend double precision number with zeroes (in binary representation) - it doesn't refer to zeros in decimal representation (please refer to some guides on floating-point numbers & their formats for deeper understanding).

Do not use mp(1/3) because 1/3 is computed in double precision and then converted to quadruple - final value has double precision accuracy (not the intended quadruple).

The correct way to compute constants to full accuracy is:

>> mp.Digits(34);
>> mp('1/3')
ans = 
    0.3333333333333333333333333333333333   

>> mp.Digits(50);
>> mp('1/3')
ans = 

    0.33333333333333333333333333333333333333333333333333

>> mp.Digits(70);
>> mp('1/3')

ans = 

    0.3333333333333333333333333333333333333333333333333333333333333333333333

In this case, double precision constant is not extended by zeros (in binary format) but computed with full accuracy from the onset (within limits of specified precision).

***

Could you please reply to another thread you started yesterday: https://mct.userecho.com/communities/1/topics/179-it-would-be-fantastic-if-quadeig-is-supported

We prioritize only well-motivated requests to add new functionality.