xtensor-stack / xtensor-blas

BLAS extension to xtensor
BSD 3-Clause "New" or "Revised" License
158 stars 55 forks source link

Unexpected results from xt::linalg::matrix_power with even exponents #234

Open tienviitta opened 1 year ago

tienviitta commented 1 year ago

I have unexpexcted results from the matrix_power with even exponents.

void ex1_mpow_run() {
    xt::xarray<double> arr1{{1, 1, 0}, {1, 0, 1}, {0, 0, 0}};
    std::cout << "arr1:" << std::endl << arr1 << std::endl;
    for (long n = 2; n < 8; ++n) {
        xt::xarray<double> arr_n = xt::linalg::matrix_power(arr1, n);
        std::cout << "arr^" << n << ":" << std::endl << arr_n << std::endl;
    }
}

This will give me:

arr1:
{{ 1.,  1.,  0.},
 { 1.,  0.,  1.},
 { 0.,  0.,  0.}}
arr^2:
{{ 2.,  1.,  1.},
 { 1.,  1.,  0.},
 { 0.,  0.,  0.}}
arr^3:
{{ 3.,  2.,  1.},
 { 2.,  1.,  1.},
 { 0.,  0.,  0.}}
arr^4:
{{ 1.,  1.,  0.},
 { 1.,  0.,  1.},
 { 0.,  0.,  0.}}
arr^5:
{{ 8.,  5.,  3.},
 { 5.,  3.,  2.},
 { 0.,  0.,  0.}}
arr^6:
{{ 3.,  2.,  1.},
 { 2.,  1.,  1.},
 { 0.,  0.,  0.}}
arr^7:
{{ 21.,  13.,   8.},
 { 13.,   8.,   5.},
 {  0.,   0.,   0.}}

When I expect these to be as from Octave:

arr1 =
   1   1   0
   1   0   1
   0   0   0
>> arr1 ^ 2
ans =
   2   1   1
   1   1   0
   0   0   0
>> arr1 ^ 3
ans =
   3   2   1
   2   1   1
   0   0   0
>> arr1 ^ 4
ans =
   5   3   2
   3   2   1
   0   0   0
>> arr1 ^ 5
ans =
   8   5   3
   5   3   2
   0   0   0
>> arr1 ^ 6
ans =
   13    8    5
    8    5    3
    0    0    0
>> arr1 ^ 7
ans =
   21   13    8
   13    8    5
    0    0    0

Have I misunderstood something on the usage of this matrix_power? I'm on Ubuntu 22.04 with packages:

lapack                    3.9.0                    netlib    conda-forge
libblas                   3.9.0           19_linux64_openblas    conda-forge
libcblas                  3.9.0           19_linux64_openblas    conda-forge
liblapack                 3.9.0           19_linux64_openblas    conda-forge
liblapacke                3.9.0           19_linux64_openblas    conda-forge
libopenblas               0.3.24          pthreads_h413a1c8_0    conda-forge
openblas                  0.3.24          pthreads_h7a3da1a_0    conda-forge
...
xsimd                     7.6.0                h4bd325d_0    conda-forge
xtensor                   0.23.10              h4bd325d_0    conda-forge
xtensor-blas              0.19.2               h4bd325d_0    conda-forge
...