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

1

Решение

Ваш matrix_power() Функция не возвращает значение из последнего блока else. Если я добавлю return true в конце этого блока ваш код запускается для меня:

$ clang++ --version
Apple LLVM version 7.0.0 (clang-700.1.76)
Target: x86_64-apple-darwin14.5.0
Thread model: posix

$ clang++ -std=c++11 -Wall -Wextra -g -I/opt/local/include lu.cpp -o lu
$ ./lu
[3,3]((-0.0644531,-0.00195312,0.861328),(1.09896,-0.0677083,-11.1406),(-0.9375,0.0625,9.4375))

Я также смог удалить это BOOST_UBLAS_TYPE_CHECK определить без получения ошибок во время компиляции (или во время выполнения). Я использую Boost 1.59.

1

Другие решения

Других решений пока нет …

А ты уже прошел курс программирования? Супер скидка!
Прокачать скилл $$$
×