%% test_lu
clear all %#ok<*CLSCR>
clc
mem = memory %#ok<*NASGU,*NOPTS>

tic, [l,u,p] = lu(9), toc %#ok<*ASGLU>
% now do this...
tic, [l,u,p] = lu(mp(9)), toc
... and look at differences ...
mem = memory
% ... wait for while ( for me about 10..40 sec, randomly !)

% ...
% and then I get out of memory error window

... here is my output ...

mem =

MaxPossibleArrayBytes: 3.3001e+09
MemAvailableAllArrays: 3.3001e+09
MemUsedMATLAB: 735469568


l =

1


u =

9


p =

1

Elapsed time is 0.024654 seconds.

l =

{3x1 cell}
{3x1 cell}
[ 0]


u =

{3x1 cell}
{3x1 cell}
[ 0]


p =

0

Elapsed time is 2.856406 seconds.

mem =

MaxPossibleArrayBytes: 3.2757e+09
MemAvailableAllArrays: 3.2757e+09
MemUsedMATLAB: 760942592

>> mp.Info
-------------------------------------------------------------------------------------------------------------
Multiprecision Computing Toolbox, (c) 2008-2015 Advanpix LLC.
Version : 3.8.4 Build 8915
Platform: 64-bit (Win64)
Release : 2015-07-22

P.S.
I reduced the virtual machine's memory to 1G and get much faster now out of memory

P.S. #2
L = lu(mp(9)) % Matlab crash

matlab code:
>> chol(36)
ans =
6
>> chol(mp(36))
ans =
3.6000e+01

Hi Pavel,
we have just released our Rational Krylov Toolbox for MATLAB version 2.0, and it now supports the use of your Multiple Precision toolbox: http://guettel.com/rktoolbox/

If you're interested to see how it's useful for us, I created one example in the collection which demonstrate its use:

http://guettel.com/rktoolbox/examples/html/example_rkfun.html

In a nutshell, our toolbox represents rational functions in matrix pencil form (H,K) associated with a discrete orthogonal basis. This basis has condition number 1 and hence we work with this representation whenever possible. For some applications, however, it is required to convert a rational function to partial fraction, quotient, or continued fraction form, which is nothing but a change of basis. As the new basis may be badly conditioned, we need multiple precision arithmetic to avoid loss of accuracy. We found your toolbox particularly useful as it's much faster than MATLAB's VPA, and also it supports EIG for matrix pencils, which VPA doesn't.

Thanks again for your great work!

Stefan

*** matlab code:
a = rand(3)
mp.ExtendConstAccuracy(0);
a0 = mp(a)
mp.ExtendConstAccuracy(1);
a1 = mp(a)
mp.Info

*** output

a =

0.2760 0.1626 0.9597
0.6797 0.1190 0.3404
0.6551 0.4984 0.5853


a0 =

0.2760 0.1626 0.9597
0.6797 0.1190 0.3404
0.6551 0.4984 0.5853


a1 =

0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000

-------------------------------------------------------------------------------------------------------------
Multiprecision Computing Toolbox, (c) 2008-2015 Advanpix LLC.
Version : 3.8.4 Build 8907
Platform: 64-bit (Win64)
Release : 2015-07-21

Trial version. Valid until 2015-08-04

Open Source Libraries Acknowledgements:
MPIR 2.7.0 C Library for Multiple Precision Integers and Rationals.
MPFR 3.2.0-dev C Library for Multiple-Precision Floating-point computations with correct Rounding.
MPC 1.1dev C Library for the Arithmetic of Complex numbers with arbitrary high precision.
LAPACK 3.5.0 State of the art software library for numerical linear algebra.
Eigen 3.2.90 C++ template library for linear algebra.
Nika 1.1.0 C++ template library for numerical mathematics.
Boost 1.49.0 Free peer-reviewed portable C++ source libraries.
KISS FFT 1.2.9 Simple and efficient mixed-radix Fast Fourier Transform library.
-------------------------------------------------------------------------------------------------------------
>>

Hello.

I tried to calculate the Hankel Function with mp objects. But as I saw it is not implemented yet. Therefore I used the definition to calculate it via Bessel functions. Following example:

First with double precision:
besselh(0,1,71.62771527592053131293674232738074i) =
0.000000000000000e+00 - 7.347258660573146e-33i

since besselh is not implemented yet, I used besselh(0,1,x)=besselj(0,x)+1i*bessely(0,x) to calculate it.

mp.Digits(34) %Quadruple precision

besselj(0,mp('71.62771527592053131293674232738074i'))+1i*bessely(0,mp('71.62771527592053131293674232738074i'))=
-6.103515625e-05 - 0.0003108494982263580776736011003369421i

which is quite different. Now with higher precision

mp.Digits(200) %Higher Precision
besselj(0,mp('71.62771527592053131293674232738074i'))+1i*bessely(0,mp('71.62771527592053131293674232738074i'))=
0 - 7.34725866057318428016340868623722...e-33i

seems better.

mp.Digits(500) %Higher Precisionbesselj(0,mp('71.62771527592053131293674232738074i'))+1i*bessely(0,mp('71.62771527592053131293674232738074i'))=
0 -7.34725866057318428016340868623722892521743368377951878300459469616931...e-33i

Interesting is the following: (using double precision)
besselj(0,71.62771527592053131293674232738074i)+1i*bessely(0,71.62771527592053131293674232738074i)=
.000000000000000e+00 - 7.347258660573146e-33i

I assume this is a bug is it?

Thank you for your answer.

I'm getting errors when I try to convert a single-precision complex variable to mp type and back. For example:
>> x = single(1 + i)
x =
1.0000 + 1.0000i
>> x = mp(x)
x =
6.94449606095542138733663817572487e-310 + 5.263544247120890312873936353273975e-315i

Due to the nature of the application I'm working on, I can't get around this by using a string, e.g. mp('1+i'); I have to convert from a single-precision variable. Presently I'm getting around this by casting the variable to double before passing it to mp, e.g. x = mp(double(x)), but I'd like to not have to do that.

In addition, when I try to cast from mp to single, the variable gets cast to a double instead. For example:
>> x = mp('1+i')
x =
1 + 1i
>> x = single(x)
x =
1.0000 + 1.0000i
>> whos
Name Size Bytes Class Attributes
x 1x1 16 double complex

Presently I'm getting around this issue by calling single() twice, e.g. x = single(single(x)), but again, I'd like to not have to do that.

I'm using MATLAB version 8.4.0.150421 (R2014b) and mp version 3.8.5 Build 9059.

Hi,
Is there any method that provides a speed up (over eig) for symmetric or Hermitian matrices? What is your current time estimate for diagonalizing a 1000x1000 real, symmetric matrix for quad precision?

Finally, how does the time to diagonalize scale with matrix size D at a fixed precision? Is it still D^3 as for double precision?