我正在尝试制作一个有趣的项目，该项目使用 BOOST Ublas 制作矩阵幂函数。它与矩阵的这个Numpy 库幂函数非常相似，它使用矩阵求幂来计算矩阵在对数时间内的 n 次幂。计算矩阵的 n 次方有 3 种情况：
- 如果 power > 0直接使用矩阵求幂
- 如果 power = 0 检查矩阵是否有逆矩阵（检查 Using lu_factorize），如果是，则返回单位矩阵
- 如果 Power < 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，i Get 有时会得到不正确的答案。（我使用矩阵和逆矩阵的乘积是单位矩阵的想法检查了这一点）
这可能与 lu_factorize 和 lu_substitute 有关，它们进行了某些检查以确保变量类型正确。
由于它们没有可用于 lu_factorize 的文档，因此我不确定如何使用它。（我刚刚阅读了一个使用 Ublas lu_factorize 和 lu_substitute 计算矩阵逆的示例。我什至阅读了源代码，但由于缺乏经验，代码掌握得并不多。
我现在对我遇到的问题有几个问题。由于这是我第一次尝试提升，如果我问一些愚蠢的问题，请不要对我苛刻。所以这些是我的问题 -
- 不正确的答案可能是由于数据类型之间的不正确转换，或类似的东西。我该如何解决这个问题？我如何确保在每一步都使用正确的数据类型？
- 当用户输入不正确的类型时，我如何给出错误，我知道我可以使用 boost assert 但这会产生大量难以理解的编译错误。确保输入类型有效的最简单方法是什么。例如，如果用户提供输入字符串，我想给出一个错误。你能为此提供一个例子吗？
我尝试了各种方法来解决编译错误，其中一种是使用#define BOOST_UBLAS_TYPE_CHECK 0 这有助于绕过数据类型的编译错误，但后来我得到了不正确的答案。你能解释一下至少在这种情况下它是如何工作的吗？
由于这是我第一次尝试从提升中做出任何事情，我可以理解这当然可以做得更好。此头文件中必须存在哪些其他内容以确保诸如错误处理、支持多个编译器等库标准？