msolve seems not working with scalar complex

[guillepo]

Dear Development team, I found some not expected results using msolve. Please, check the following examples using PDL::LinearAlgebra:

1) error to large in the complex part

my $A=exp(i()*0.01*sequence(4,4,3)) ;
my $B=exp(i()*0.1*sequence(1,4,3));
my ($E,$er)=msolve($A,$B);
say $A x $E -$B;

2) Did not work at all in trivial case saying that found singular matrix after use getrf factorization

my $A=pdl([i(),0.0],[0,i()])
my $B=pdl([1.0],[1.0])->r2C

Regards, Guillermo

[mohawk2]

@guillepo Thanks for the report! Digging into the various pure-Perl functions in PDL::LinearAlgebra has revealed that few of them will handle native complex correctly, although if you were to use PDL::Complex objects, that would probably work. However, since working with PDL::Complex objects is so complicated and difficult, the right thing to do is to make them all handle native complex, which I am now doing.

[guillepo]

Thanks @mohawk2, OK, I might try with PDL::Complex to manage the interaction with PDL::LinearAlgebra. But, as you said native complex calculations are easier. Please, let me know about your progress or if in something I can try to help you.

Regards

[mohawk2]

Thanks for the offer! Right now I'm just fighting a problem with the mrcond function, which works with real ndarray, but crashes with a realloc error when I call it with similarly-valued PDL::Complex ndarray. This suggests an error in the P:LA:Complex functions it's calling. If you have any insights (or even can look at the current master test for mrcond in t/1.t and see if it's sensible values for complex data, that would be great!

[guillepo]

Just I am doing some tests with PDL::Complex and PDL::LinearAlgebra, and seems work if cast cplx is avoided. For example: $A=random(2,2,5)+i()*random(2,2,5); $B=(random 1,2,5)->r2C; ($E,$er)=msolve($A,$B); say $A x $E -$B; gives enough good I guess, and also $A=random(2,2,5)->r2C; .... But, if I try with $A= cplx random(2, 2,5); or $A=random(2,2,5)->complex; The info say that it is complex ndarray, but it seems did not report extra an dimension for complex as if were and scalar ! then my example crash due to incompatibilities with $B I am not sure if above lines can help you, but line 25 in t/1.t just are using cplx and also you are using PDL::Complex. Maybe I am wrong, but in my short experience using PDL::Complex is not planned for scalar (native) complex, at least in the beginning.

Regards

[mohawk2]

You are completely right about PDL::Complex! Until we brought in native complex, there was only PDL::Complex, and P:LA worked to a sufficient degree with that, at least on certain routines (notably, the ones used by Photonic). My strategy in fixing the problem you have correctly identified is to start by establishing a baseline of things working with PDL::Complex as well as with real, by adding actual tests (there were hardly any at all) for that. The bump in the road here is that is revealing at least potential problems with the complex routines.

Once that's settled, I can then add support for native-complex (with tests), and have assurance that I haven't broken anything. So, any extra test cases you can provide that aren't already covered will be of great value, and even more so if you identify mistakes I have made so far :-)

[guillepo]

OK, thank you very much for your work! By the way, call my attention that in t/1.t that deal complex there are not use of cgetrf, the version for complex number of getrf subrutine as you use in line 70 of t/1.t. Regards

[guillepo]

OK, just found it in line 87, but with a call to PDL::LA::C::cgetrf I suppose that cgetrf is accesible anyway if it is needed using PDL::LA only Regards

[mohawk2]

I believe (I'm looking at other stuff right now) that the problem lies in the shape of things being passed by P:LA's PDL::Complex::mrcond, to cgetrf (and clearly there is a problem in that it should not be possible to cause memory errors with Perl code, at all). If you have opinions on that, let me know, or I'll get to it in due course :-)

[mohawk2]

Update: I recall it may be cgecon which is blowing up, which is reassuring since cgetrf is used heavily in Photonic!

[guillepo]

Yes, you are right respect to shape. Just I can manage PDL::Complex with PDL::LA Te main idea is to try all variables as complex from the beginning to guarantee the correct threading.

[guillepo]

Not so fine, I get this message, Subroutine sumover redefined at /home/gortiz/perl5/perlbrew/perls/perl-5.26.1/lib/site_perl/5.26.1/x86_64-linux/PDL/Complex.pm line 1424. But msolve pass with error less than 1e-11 with 4x4(x50x50) complex matrix under threading.

[guillepo]

But, not stable at all. After some a few test, I found this message: Stringizing problem: Math::Complex::make: Cannot take real part of 'NaN'. at /home/gortiz/perl5/perlbrew/perls/perl-5.26.1/lib/site_perl/5.26.1/x86_64-linux/PDL/Core.pm line 2897. PDL::string(PDL=SCALAR(0x55903b79ab00)) called at /home/gortiz/perl5/perlbrew/perls/perl-5.26.1/lib/site_perl/5.26.1/x86_64-linux/PDL/Complex.pm line 1414

[guillepo]

powerfull solution to mange PDL::Complex in Photonic was to move first dimension to the last, peroform al opertaions included threading and then back last dimension to the first position. Maybe such a trick could work here also. Regards

[mohawk2]

I will note that (I did a lot of work on Photonic and the technique you describe may be necessary). Thank you for the further research! That redefinition isn't a problem.

[mohawk2]

Could you share the code that triggered the NaN problem?

[guillepo]

Yes, I know about your further contribution to Photonic, it is great ! The code below run without complain, But if you uncomment L25 (then comment L26) instead of show matrix B it return the message that I said. Of course, it is not msolve issue, but maybe a related problem with it ?

#!/usr/bin/env perl
use warnings;
use strict;

use feature qw(say);
use PDL;
use PDL::NiceSlice;
use PDL::Constants qw(PI);
use PDL::LinearAlgebra;
use PDL::Complex;
$PDL::BIGPDL=1;
$PDL::BIGPDL=1;

my $c=1;

my $nw=2;
my $nK=2;
my $w=zeroes($nw)->xlinvals(0.1,4)->r2C->dummy(1,$nK); # K,w
my $K=zeroes($nK)->xlinvals(0.1,PI)->r2C->dummy(2,$nw)+0.0001*i(); # K,w

my $epsA=1.0+0*i();
my $muA=1.0+0*i();
my $eps1=1.0+0.0*i();
my $mu1=1.0+0.0*i();
#my $eps2=drude($w,0.01);
my $eps2=12.0+0.0*i();
my $mu2=1.0+0.0*i();

my @mu=($muA,$mu1,$mu2);
my @eps=($epsA,$eps1,$eps2);
my @epsmu=($epsA*$muA,$eps1*$mu1,$eps2*$mu2);

my $d=1.0;
my @d=(0.6*$d,0.4*$d);

my @k=map{($w/$c)*(mysqrt($_)->complex)} @epsmu;

my $A=pdl(
    [exp(i()*$k[1]*$d[0]),-exp(i()*$k[2]*$d[0]),exp(-i()*$k[1]*$d[0]),-exp(-i()*$k[2]*$d[0])],
    [$k[1]*exp(i()*$k[1]*$d[0]),-$k[2]*exp(i()*$k[2]*$d[0]),-$k[1]*exp(-i()*$k[1]*$d[0]),$k[2]*exp(-i()*$k[2]*$d[0])],
    [r2C(1.0)->dummy(1,$nK)->dummy(2,$nw),-exp(i()*($k[2]-$K)*$d),r2C(1.0)->dummy(1,$nK)->dummy(2,$nw),-exp(-i()*($k[2]+$K)*$d)],
    [$k[1],-$k[2]*exp(i()*($k[2]-$K)*$d),-$k[1],$k[2]*exp(-i()*($k[2]+$K)*$d)]
    )->complex->mv(-1,1)->mv(-1,1);
say "A:\t", $A->info;
my $J0=1.0;
my $DE=4*PI*($w/$c)*$J0*(1.0/(i()*$c*($k[2]**2-$K**2))-1.0/(i()*$c*($k[1]**2-$K**2)));

my $B=($DE*(pdl([exp(i()*$K*$d[0]),$K*exp(i()*$K*$d[0]),r2C(1.0)->dummy(1,$nK)->dummy(2,$nw),$K])->complex))->mv(-1,1)->dummy(1);

say "B:\t", $B;
my ($E,$er)=msolve($A,$B);

say "OK\n" if all approx $A x $E , $B, 1e-10;

BEGIN{
    thread_define  '_mysqrt(a();[o]b())' => over {
    my ($a,$b)=@_;
    my $r=sqrt($a);
    $r=-$r if $r->im<0;
    $b.=$r;
    };

    sub mysqrt{
    my $a=shift;
    my $b=null;
    _mysqrt($a,$b);
    return $b;
    }

    sub drude{
    my $omega=shift;
    my $gamma=shift;
    my $eps=1-1/16/($omega*($omega+i()*$gamma));
    }

}

[mohawk2]

Thanks for this! I edited it to make it show as code; a single ` won't quote multiple lines; you need to start and end the whole block with ``` on lines by themselves.

[guillepo]

OK, thanks

[mohawk2]

Now the performance problem with PDL::NiceSlice is resolved (see https://github.com/PDLPorters/pdl/pull/74), we will be able to resume work on this. PDL 2.069 has just been released, so I'd suggest upgrading to that so that PDL::LinearAlgebra loads quicker.

[guillepo]

OK, thanks, I just upgrade to PDL and PDL::LA last version. Back to shape of ndarray and PDL::Complex. It seems that the problem is with thread_define in the code above. I am trying to update to PDL::Complex, using

  thread_define 'mysqrtaux(x(2);[o]s(2))', over {
    my $x=shift;
    my $s=shift;
    my $t;
    $t=sqrt($x->cplx); #complex
    $t=-$t if $t->im <0 || ($t->im==0 && $t->re < 0);
    $s.=$t;
};

    sub mysqrt {
    my $x=shift;
    my $r=null;
    mysqrtaux($x, $r);
    return $r->complex;
    }

But, I obtain the following error

threadover: need one creating flag per pdl! at /home/gortiz/.perlbrew/libs/perl-5.26.1@testeo/lib/perl5/x86_64-linux/PDL/Core.pm line 1658. Thank you for your comments Regards

[mohawk2]

Thank you, I'll investigate that too! There is apparently also a bug in main PDL which the mighty @zmughal has uncovered.

[mohawk2]

Update for those following at home: the mrcond problem was caused by a bug in LAPACK's complex gecon routines needing a "real"-sized work area, while the real ones wanted an "integer"-sized one, but P:LA's wrapper for gecon was doing the same for complex as for real. That's now fixed. There are a few other problems also fixed, and now all the P:LA:m* methods that have a real and complex version are tested, which will then make it safe to start adding support for native-complex.

[guillepo]

OK, thank you for this message. Only for the chance, I just test with P::LA and intrinsic complex but msolve is not working yet. the imaginary part of the results still is not fine. Regards

[mohawk2]

Sorry to be unclear! The above only means the repo master has some fixes. Until I've finished adding tests for all the functions in P:LA, it's not safe for me to try updating them to native-complex. It's going to take another couple of days, sorry.

If you're feeling adventurous and want to participate, you could fork and then clone your copy of this repo locally, and run cover -test then view the generated HTML coverage info in cover_db/. The latest Coveralls data shows that I've brilliantly missed these functions/methods (because they didn't have separate PDL::Complex methods, instead having internal switch on how many dims):

• mcond
• morth
• mnull
• mpinv
• mlu
• mhessen
• mschurx
• msolvex
• msymsolvex
• mpossolvex
• mlls
• mllsy
• mllss
• mlse
• meigenx
• mgeigenx
• msymeigenx
• msymgeigen
• msymgeigenx
• msvd

as well as the various select functions modes of the Schur functions I have tested. Data to trigger some of the non-converging errors would be valuable too. Basically, if you feel like PR-ing updates to t/1.t that exercise some of these (hopefully it's very obvious how, and easy) that would be great! If you simply provide on this issue, data that's either valid for these, or deliberately invalid, that's also great.

[guillepo]

OK, thanks for the invitations, but for now, maybe the bug reported can solve easily by noting that in line 77 of file t/1.t you may invoke getrf function directly, but in line 92 of file t/1.t the complex cgetrf is invoked through PDL::LinearAlgebra::Complex module. It sound like some wrap is needed for complex? Nevertheless, I understand that it is a tremendous work that you are doing according to your roadmap and times to all is needed.

Regards

[mohawk2]

Update: there is now what I feel is sufficient test coverage of PDL::LinearAlgebra functions and all their various (complicated) usage modes, so I can confidently start rearranging the code to handle complex versions more naturally. Once that is in place, handling native complex should be easy.

[mohawk2]

I haven't forgotten this! Just pushed an update to the PP code which allowed deleting a bunch of code, since it allows one to give nulls for output variables. Next, to finish unifying the three classes of tricky functions: Schur, linear systems, eigenvalues, then it's finished.

[mohawk2]

That's the Schur functions unified.

[mohawk2]

And the linear systems.

[guillepo]

Dear developers, Only for the case that "linear systems" also refers to msolve() routine. This example should be working yet? use PDL::LinearAlgebra; $B=(pdl[[1.0],[1.0]])->r2C; $A=pdl([[i(),0.0],[0.0,i()]]); $X=msolve($A,$B);

Regards

El sáb, 2 abr 2022 a la(s) 08:16, mohawk2 @.***) escribió:

And the linear systems.

— Reply to this email directly, view it on GitHub https://github.com/PDLPorters/pdl-linearalgebra/issues/10#issuecomment-1086618180, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABRRDHZGGINZRVXKFPRMOHDVDAUA3ANCNFSM5MPYPPQQ . You are receiving this because you were mentioned.Message ID: @.***>

-- Saludos, Guillermo

[mohawk2]

"Linear systems" here was just the functions whose names are like mlls. This is only available so far on the master branch. Native complex is getting close to working but a little more tidying up of code is required to make it easy and relatively painless for me to implement.

[mohawk2]

To give an insight into what has been required, when I started the PDL::LinearAlgebra module was about 6800 lines. It is now about 3500 lines. Now we're down to solving the problem of the native-complex PP-wrapped LAPACK functions passing back ndarrays that should be e.g. 1+0i 2+0i 3+0i and are instead 1+2i 3+0i 0+0i. It looks like a real ndarray gets passed back instead of a complex one. Still debugging!

[mohawk2]

Solved that real/complex confusion. Other (increasingly-small) problems remain!

[mohawk2]

With the above-linked commit, basic (non-broadcasting) functions all appear to work with native complex. However, a quick try of your original posted code, with 4x4x3 ndarrays and msolve, give odd results: stripping it back to just 4x4x2, multiplying $A by $E gets the first 4x1 right, but the second one is wrong. Something appears to be wrong with its broadcasting behaviour. Still investigating.

One observation is that the Perl wrapper functions in P:LA had some mistakes that meant they wouldn't have broadcast correctly. But part of the restructuring should have eliminated most of those problems.

[mohawk2]

I now suspect the above point is just due to numerical instability in the chosen example (the $E had, in the middle matrix, values up to 1e+33).

I have therefore released to CPAN version 0.27, with the various fixes. Please try using it with native complex, and report any problems! If it works Ok for you, please then close this issue :-)

[guillepo]

Now 0.27 version of PDL::LinearAlgebra is working as expected in my code using msolve(). Thanks to Ed for your support

[mohawk2]

Very glad to hear it! Thank you for the feedback.