Hi Guys,

Nice catch!

As for disp - it is funny and easy to fix right away.
Just replace display(x) with disp(mpimpl(2,x,inputname(1))); in disp(x) function in mp.m file (line #903)

As for compact - need to implement its support in C++. Have been thinking about this, but was waiting until somebody asks. I guess time has come :).

Both issues will be resolved in next version.

Thank you,
Pavel.

Hi again,

Sorry for delayed reply (the forum didn't send me notification on your reply).

Please send me your name & email address - so that I will be able to issue trial license and provide download link to MacOSX version of toolbox. My direct e-mail: pavel@advanpix.com

Thank you,
Pavel.

Also in case of symmetric/Hermitian matrix you can choose what method to use, by using the second argument to "eig":

>> A = mp(rand(1000)-0.5);
>> A = A+A';
>> tic; eig(A); toc;          % use default, which is MRRR
Elapsed time is 32.662471 seconds.

>> tic; eig(A,'dc'); toc;     % use divide & conquer
Elapsed time is 31.053027 seconds.

>> tic; eig(A,'mr'); toc;     % use MRRR
Elapsed time is 32.703078 seconds.

MRRR and divide & conquer are both optimized for symmetric matrices. Generally timings are comparable, but we use MRRR by default as it shows better results over wide range of matrices.

Second parameter of "eig" is toolbox-specific, don't try using it with built-in routine :).

Hi,

Thank you for your questions.

1. Of course, the "eig" routine in toolbox uses different algorithms depending on matrix properties.
In particular, it uses MRRR in symmetric/Hermitian case. It is current state of the art delivering near O(n^2) complexity.
(see references on the page: http://goo.gl/L2Nzrv and also this link: http://goo.gl/v9iSxq).

2. Time estimate depends on CPU/number of cores, etc. On my computer with Core i7 (4 HW cores) I see following results for quadruple precision:

>> A = mp(rand(1000)-0.5);
>> A = A+A';
>> tic; eig(A); toc;
Elapsed time is 32.662471 seconds.

Unsymmetric case for comparison:

>> A = mp(rand(1000)-0.5);
>> tic; eig(A); toc;
Elapsed time is 134.259668 seconds.

3. MATLAB/VPA, Maple and Mathematica need hours (days in case of VPA) to handle these problems, as all are using the standard QR iteration which is O(n^3).

4. Our trial version is freely accessible from our website - you can test the toolbox in your environment in a minutes.

Let me know if you have any other questions.
Pavel.

I think you are right. For some reason I assumed that second output from chol should also reflect if matrix Hermitian or not. But it does check only positive definiteness.

This clears things out. Thank you!
(I used documentation for 2011a, which has a less strict wording - source of my confusion).

Yes, your matrix is Hermitian, of course.

But in my example matrix was not Hermitian. Yet MATLAB/zpotrf was not able to recognize this right away and returned p=3, indicating that A(2,2) minor is HPD! Which is not true (matrix A is not even Hermitian!).

Moreover returned R gives:

>> norm(R'*R-A(1:p-1,1:p-1),1)
ans =
       2.8084

Which contradicts the documentation: When A is full, R is an upper triangular matrix of order q=p-1 such that R'*R=A(1:q,1:q).

That was my whole point. I still think correct result should be R=[], p=1, meaning A has no HPD minors :).
What do you think?

(I have traced this behavior down to zpotrf).

New version has been released, please download it from our website.

Just completed adding the feature - not building the release.

Here is some example you might find funny (this is pure MATLAB):

>> A =
          100 +          0i      -1.0151 +     1.4042i     -0.12828 +    0.94526i
      -1.0151 +     1.4042i          100 +          0i     -0.15982 +      1.622i
     -0.12828 +    0.94526i     -0.15982 +      1.622i     -0.12722 +    0.70364i

>> [R,p] = chol(A)
R =
           10 +          0i     -0.10151 +    0.14042i
            0 +          0i       9.9985 +          0i
p =
     3

>> norm(R'*R-A(1:p-1,1:p-1),1)
ans =
       2.8084

Matrix A has no positive definite minors, but real positive diagonal tricked MATLAB and it generates unexpected result.
Correct result should be R = [], p =1. Or am I missing something?

Yes, up to this moment we supported only the following syntax:

R = chol(A)
L = chol(A,'lower')
R = chol(A,'upper')

I understand that the two-outputs version is important and it will be ready shortly :).

***
I have been focusing on covering the MATLAB functionality in breadth, e.g. main functions in prime syntax.
I am always open and happy to do any refinements requested by users.

I believe this user-driven style of development naturally shapes the toolbox functionality, and get audience in different fields.
Besides, that way I am able to actually get some work done:), as some of the functions might take months/years to implement in all possible variations.

I hope for understanding and I am always welcome any requests. Please continue sending your excellent bug reports, comments and suggestions. It is highly appreciated!
***

I plan to release the updated "chol" today or tomorrow the latest....