使用Boost Ublas lu_Factorize时出错

我正在尝试制作一个有趣的项目，它使用BOOST Ublas制作矩阵幂函数 . 它非常类似于Numpy library power function for matrix它使用矩阵求幂来计算矩阵在对数时间内的n次幂 . 它有3个案例来计算矩阵的n次幂：

• 如果power> 0使用矩阵求幂 directly

• 如果power = 0，检查矩阵是否有反转（检查使用lu_factorize），如果是，则返回 identity matrix

• 如果Power <0找到 inverse(if it exists) ，则对其使用矩阵求幂

我很擅长算法和实现，但这是我第一次使用任何开源库，我想学习这个，这样我最终可以为提升做出贡献 .

这是我的头文件

//  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 and lu_substitute 有关，它有一定的检查以确保变量类型是正确的 .

由于他们没有可用于lu_factorize的文档，我不知道如何使用它 . （我刚刚阅读了这里可用的一个例子using Ublas lu_factorize and lu_substitute to compute matrix inverse） . 我甚至阅读了源代码，但由于缺乏经验，代码没有太多把握 .

我现在对我遇到的问题有几个问题 . 由于这是我第一次尝试提升，如果我问一些愚蠢的事情，请不要对我不好 . 所以这些是我的问题 -

• 错误答案可能是由于数据类型之间的不正确转换或类似情况 . 我该如何解决这个问题？如何确保在每一步使用正确的数据类型？

• 当用户输入不正确的类型时，我如何给出错误，我知道我可以使用boost断言，但这会给出一些难以理解的编译错误 . 确保输入类型有效的最简单方法是什么 . 例如，如果用户给出输入字符串，我想给出一个错误 . 你能为此提供一个例子吗？

• 我已经尝试了各种方法来解决编译错误，其中一个是使用 #define BOOST_UBLAS_TYPE_CHECK 0 这有助于绕过数据类型的编译错误，但后来我得到了错误的答案 . 在这种情况下，你能解释一下这至少是如何运作的吗？

• 由于这是我第一次尝试通过提升做出任何事情，我可以理解，这肯定可以做得更好 . 此头文件中必须存在哪些其他内容以确保库标准，如错误处理，支持多个编译器等？