forum.alglib.net

I'm having problems with rmatrixgemm that I don't understand. Here's some code that demonstrates the error I'm seeing. I'm trying to compute Vetst = Jtst.VstdsTst.Transpose[Jtst].
I break it up as VtmpTst = Jtst.VstdsTst then Vetst = VtmpTst . Transpose[Jtst]. VtmpTst gets computed correctly, but Vetst returns a matrix with the elements mostly equal to
1.5762493960899296e+134 with a few equal to 0. When I compute the result in Mathematica I get a reasonable result. I tried different data for Jurf2r with the same result. Any suggestions?
(Note, the different data for Jurf2r was created by replacing scaling values less than 0.009 by a factor of 100.)
double Jurf2r[240] = {
1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,
0., 0., 0., 0., 0., 0., 0.000283063, 0.000248748, 0.0133716, 0.0146574, 1., 0.,
0., 0., 0., 0., 0., 0., -0.000248748, 0.000283063, -0.0146574, 0.0133716, 0., 1.,
0., 0., 0.5556, 0.364916, 0., 0., 0.00876054, -0.0113775, 0.262187, -0.272316, 0., 0.,
0., 0., -0.364916, 0.5556, 0., 0., 0.0113775, 0.00876054, 0.272316, 0.262187, 0., 0.,
-0.000150607, -0.0000509269, 0.0165818, -0.00413272, 0.262187, -0.272316, 0., 0., 0., 0., 0., 0.,
0.0000509269, -0.000150607, 0.00413272, 0.0165818, 0.272316, 0.262187, 0., 0., 0., 0., 0., 0.,
-0.0023959, -0.0120558, 0.458209, 0.475855, 0., 0., 0., 0., 0.313991, -0.206255, 0., 0.,
0.0120558, -0.0023959, -0.475855, 0.458209, 0., 0., 0., 0., 0.206255, 0.313991, 0., 0.};

double VrfInDta[144] = {
7.409346220569623e-6, -4.235164736271502e-22, 1.3909263055141595e-6, 4.935692141050101e-7, 9.32878743417026e-6, 9.765893771249042e-6, 0., 0., 0., 0., 0., 0.,
4.235164736271502e-22, 7.409346220569623e-6, -4.935692141050101e-7, 1.3909263055141595e-6, -9.765893771249043e-6, 9.32878743417026e-6, 0., 0., 0., 0., 0., 0.,
1.3909263055141595e-6, -4.935692141050101e-7, 8.14680793030778e-6, 0., 1.822267912898904e-6, -1.0270045827287743e-7, 0., 0., 0., 0., 0., 0.,
4.935692141050101e-7, 1.3909263055141595e-6, 0., 8.14680793030778e-6, 1.0270045827287752e-7, 1.8222679128989042e-6, 0., 0., 0., 0., 0., 0.,
9.328787434170258e-6, -9.765893771249043e-6, 1.822267912898904e-6, 1.0270045827287817e-7, 0.00005137656878429575, 2.646977960169689e-23, 0., 0., 0., 0., 0., 0.,
9.765893771249042e-6, 9.328787434170258e-6, -1.0270045827287808e-7, 1.8222679128989042e-6, -2.646977960169689e-23, 0.00005137656878429575, 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0.00002291046759316638, 0., -1.2254910215411243e-6, 7.71989946091886e-6, -0.000014649761118660017, 0.000019756720382045573,
0., 0., 0., 0., 0., 0., 0., 0.00002291046759316638, -7.71989946091886e-6, -1.2254910215411243e-6, -0.000019756720382045573, -0.00001464976111866002,
0., 0., 0., 0., 0., 0., -1.2254910215411248e-6, -7.719899460918858e-6, 0.000026948568219827768, 8.470329472543003e-22, 0.00001049646805839054, 6.943211629142092e-6,
0., 0., 0., 0., 0., 0., 7.719899460918858e-6, -1.2254910215411248e-6, -8.470329472543003e-22, 0.000026948568219827768, -6.943211629142091e-6, 0.00001049646805839054,
0., 0., 0., 0., 0., 0., -0.000014649761118660017, -0.000019756720382045573, 0.000010496468058390536, -6.943211629142092e-6, 0.000053672058978352506, -3.9704669402545333e-22,
0., 0., 0., 0., 0., 0., 0.000019756720382045573, -0.00001464976111866002, 6.943211629142092e-6, 0.000010496468058390536, 3.9704669402545333e-22, 0.00005367205897835251};

alglib::real_2d_array Jtst, VstdsTst, VtmpTst, VeTst;
Jtst.setcontent(20,12,Jurf2r);
VstdsTst.setcontent(12,12,VrfInDta);
VtmpTst.setlength(20,12);
VeTst.setlength(20,20);
rmatrixgemm(20, 12, 12, 1.0, Jtst, 0, 0, 0, VeRFIn, 0, 0, 0, 0.0, VtmpTst, 0, 0); // this one works
rmatrixgemm(20, 12, 20, 1.0, VtmpTst, 0, 0, 0, Jtst, 0, 0, 1, 0.0, VeTst, 0, 0); //this one doesn't

Hi!

I see that you made a mistake when specifying matrix sizes. The first pair of numbers in rmatrixgemm() is a result size, so the first GEMM produces a 20x12 matrix, and the second one should also produce a 20x12 one. Try replacing "20, 12, 20" in the second call by "20, 20, 12", does it work?

Sergey.

Yes that works. Thanks for the help.