Ошибка в использовании Boost Ublas lu_Factorize

Я пытаюсь сделать забавный проект, который делает функцию мощности матрицы
используя BOOST Ublas. Это очень похоже на это Степенная функция библиотеки Numpy для матрицы Он использует матричное возведение в степень для вычисления n-й степени матрицы в логарифмическом времени.
Он имеет 3 случая для вычисления n-й степени матрицы:

1. Если мощность> 0, используйте матрицу возведения в степень непосредственно
2. Если power = 0, проверьте, есть ли у матрицы обратное (проверьте используя lu_factorize), и если да, верните единичная матрица
3. Если сила < 0 найти обратный (если он существует)
   затем использовать матрицу возведения в степень на нем

Я хорошо с алгоритмами и реализацией, но это
в первый раз я использую любую библиотеку с открытым исходным кодом, я хочу
чтобы узнать это, чтобы я мог в конечном итоге внести свой вклад в повышение.

Это мой заголовочный файл

//  Distributed under the Boost Software License, Version 1.0. (See
//  accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)
//
#ifndef BOOST_UBLAS_POW_HPP
#define BOOST_UBLAS_POW_HPP
#include <iostream>
#include <vector>
#include <iomanip>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/mpl/set.hpp>
#include <boost/mpl/assert.hpp>
#include <boost/multiprecision/cpp_int.hpp>
using namespace boost::numeric::ublas;
namespace boost { namespace numeric { namespace ublas {

typedef permutation_matrix<std::size_t> pmatrix;
template< typename T,typename T2 >
matrix<T> matrix_power(const matrix<T> input, T2 exponent)
{
matrix<T> resultant=input;
BOOST_ASSERT_MSG(input.size1()==input.size2(),"Not a square matrix\n");
if(exponent>0)
resultant=matrix_exponent(input,exponent);// this where you could directly use matrix exponentiation

else if(exponent==0)
{
pmatrix pm(input.size1());
matrix<T> A(input);
BOOST_ASSERT_MSG(lu_factorize(A, pm)==0,"Attempted to compute a  power of 0 in a matrix without inverse\n");
resultant.assign(identity_matrix<T> (input.size1()));// if the matrix is invertible output identity matrix for the 0t hpower
}
else {
matrix<T> A(input);
pmatrix pm(A.size1());
BOOST_ASSERT_MSG(lu_factorize(A, pm)==0,"Attempted to compute inverse in a singular matrix\n");
resultant.assign(identity_matrix<T> (A.size1()));
lu_substitute(A, pm, resultant);
resultant=matrix_exponent(resultant,std::abs(exponent));
}// in last case first we compute the inverse of the matrix and then we do matrix exponentiation
return resultant;
}template< typename T,typename T2 >
matrix<T> matrix_exponent(matrix<T> base,T2 exponent)
{
matrix<T> resultant=base;
exponent-=2;
while(exponent>0)
{
if(exponent%2==1)resultant=prod(resultant,base);
exponent=exponent >> 1;
base=prod(base,base);
}
return resultant;
}
}
}
}
#endif

Я тестирую этот заголовочный файл, используя этот

#include "power.hpp"using namespace boost::numeric::ublas;
using namespace std;
typedef permutation_matrix<std::size_t> pmatrix;
int main()
{
matrix<double> m (3, 3);
for(int i=0;i<3;i++)for(int j=0;j<3;j++)m(i,j)=3*i+j+8;
m(0,0)=11;
matrix<double> c(3, 3);
int h=matrix_power(m,c,-1);// c stores -1 power of m
if(h)// h tells whether the power exists or not
std:: cout << c << "\n\n\n";}

Функции отлично работают при мощности> 0, потому что используют матрицу
Возведение в степень непосредственно. Работает намного быстрее, чем повторяющееся умножение, я видел, что на протяжении примерно 1000 итераций время его выполнения и использование цикла различаются в 100-1000 раз. Вы можете наблюдать это
Но для власти <= 0, иногда получаю неправильный ответ. (Я проверил это, используя идею о том, что произведение матрицы и инверсии является единичной матрицей)

Это, вероятно, как-то связано с
lu_factorize и lu_substitute который имеет определенные проверки для обеспечения правильности типов переменных.

Поскольку их нет в наличии для lu_factorize, я не уверен, как его использовать. (Я только что прочитал пример, доступный здесь используя Ublas lu_factorize и lu_substitute для вычисления обратной матрицы Я даже читаю исходный код, но код не очень-то разбирается из-за недостатка опыта.

Теперь у меня есть несколько вопросов относительно проблемы, с которой я сталкиваюсь.
Так как это моя первая попытка повышения, если я спрошу что-нибудь глупое, пожалуйста, не будь настойчивым со мной.
Так что это мои вопросы —

1. Возможно, неправильный ответ связан с неправильным преобразованием типа данных или чем-то подобным. Как я могу решить это? Как обеспечить использование правильных типов данных на каждом этапе?
2. Как выдавать ошибки, когда пользователь вводит неверные типы, я знаю, что могу использовать boost assert, но это
   дает множество непонятных ошибок компиляции. каковы
   простейшие способы убедиться, что тип ввода является действительным. например я хочу
   выдавать ошибку, если пользователь вводит строку для ввода.
   Можете ли вы привести пример для этого?

3. Я пробовал различные методы для устранения ошибок компиляции
   один из которых должен был использовать #define BOOST_UBLAS_TYPE_CHECK 0
   это помогло обойти ошибки компиляции для типов данных, но затем я получил неправильные ответы. Можете ли вы объяснить, как это работает по крайней мере в этом случае?

4. Поскольку это моя первая попытка сделать что-либо из ускорения, я могу понять, что это, безусловно, можно сделать намного лучше. Какие другие вещи должны присутствовать в этом заголовочном файле для обеспечения таких стандартов библиотеки, как
   Обработка ошибок, поддержка нескольких компиляторов и т. Д.?