Napisałem sobie program do wyliczania rozwiązań układu równań liniowych. Napisałem sobie najpierw pierwszy program, który działał dobrze w konsoli. Chciałem sobie zrobić program w qt (aby się lepiej obsługiwało). Napisałem sobie kodzik, wszystko ok. Tylko jeden szkopuł. Wyniki algorytmu nie są poprawne (chociaż je bezpośrednio zaimportowałem z 1. programu (plik nagłówkowy z szablonami funkcji)). Hmm... gdzie robię błąd, bo nie podejrzewam szczerze takiej specyfiki biblioteki qt.
Kod:
ortho.h
#ifndef _ORTHO_
#define _ORTHO_
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
namespace ublas = boost::numeric::ublas;
using namespace ublas;
namespace myNam
{
template <class S>
int QR_Ortho_Decomposition (ublas::matrix<S> A, ublas::matrix<S> &Q, ublas::matrix<S> &R)
{
const int n = A.size1();
const int m = A.size2();
ublas::vector<S> *u = new ublas::vector<S>[n];
ublas::vector<S> *e = new ublas::vector<S>[n];
//std::cout << "Macierz A: " << A << std::endl;
for (int k = 0; k < m; k++)
{
u[k] = column (A, k);
ublas::vector<S> temp (n);
temp.clear();
for (int j = 0; j < k ; j++)
{
// std::cout << "ALA" << std::endl;
vector<S> tmp = ortho_projection(e[j], u[k]);
// std::cout << "Projection: " << tmp << std::endl;
temp += tmp;
}
// std::cout << "Temp(" << k << ")" << temp << std::endl;
u[k] -= temp;
e[k] = u[k]/norm_2(u[k]);
}
/*
for (int i = 0; i < n; i++)
{
std::cout << "E(" << i << "): " << e[i] << std::endl;
std::cout << "U(" << i << "): " << u[i] << std::endl;
}
*/
for (int i = 0; i < n; i++)
for (int j = i; j < m; j++)
R(i, j) = inner_prod(e[i], column(A, j));
for (int i = 0; i < n; i++)
for (int j = 0; j < m; j++)
Q(i, j) = e[j](i);
delete[] e;
delete[] u;
return 0;
}
template <class S>
ublas::vector<S> ortho_projection (ublas::vector<S> e, ublas::vector<S> a)
{
ublas::vector<S> ret = inner_prod(e,a)/inner_prod(e,e)*e;
//std::cout << "Projection: " << ret << std::endl;
return ret;
}
}
template <class S>
ublas::vector<S> ortho_solve (ublas::matrix<S> A, ublas::vector<S> b)
{
ublas::vector<S> x(b.size());
if (A.size1() != A.size2()) throw (-1);
int n = A.size1();
matrix<S> Q(n,n);
matrix<S> R(n,n);
int result = myNam::QR_Ortho_Decomposition(A, Q, R);
matrix<S> D = prod(Q, trans(Q));
// odwrócenie D
for (int i = 0; i < n; i++)
D(i,i) = 1/D(i,i);
//std::cout << "Odwrócone D - D^-1" << D << std::endl;
b = prod(prod(D, trans(Q)), b); // nowe b jest źle obliczane
std::cout << "DATA" << std::endl;
std::cout << "R: " << R << std::endl;
std::cout << "B new: " << b << std::endl;
std::cout << "ENDDATA" << std::endl;
//std::cout << prod (Q, trans(Q)) << std::endl;
x = matrix_upper_solve (R, b); // rozwiązywanie z macierzami diagonalnymi
//std::cout << x << std::endl;
return x;
}
template <class S>
ublas::vector<S> matrix_upper_solve (matrix<S> U, vector<S> b)
{
// dobrze napisane
const int n = b.size();
vector<S> x(n);
x.clear();
for (int i = n-1; i>=0; i--)
{
S suma = 0;
for (int j = i+1; j <= n-1 ; j++)
//std::cout << U(i,j)*x(j) << std::endl;
suma = suma + U(i,j)*x(j);
x(i) = (b(i) - suma)/U(i,i);
}
return x;
}
#endif _ORTHO_
mainwindow.cpp
#include "mainwindow.h"
#include "ui_mainwindow.h"
#include "ortho.h"
#include <iostream>
MainWindow::MainWindow(QWidget *parent) :
QMainWindow(parent),
ui(new Ui::MainWindow)
{
ui->setupUi(this);
}
MainWindow::~MainWindow()
{
delete ui;
}
void MainWindow::changeEvent(QEvent *e)
{
QMainWindow::changeEvent(e);
switch (e->type()) {
case QEvent::LanguageChange:
ui->retranslateUi(this);
break;
default:
break;
}
}
void MainWindow::on_pushButton_clicked()
{
int size = this->ui->tableWidget->rowCount();
ublas::matrix<long double> A(size,size);
ublas::vector<long double> b(size);
for (int i = 0; i < size; i++)
for (int j = 0; j < size; j++)
A(i,j) = this->ui->tableWidget->item(i,j)->text().toLongLong(NULL, 10);
std::cout << A << std::endl;
for (int i = 0; i < size; i++)
b(i) = this->ui->tableWidget_3->item(i,0)->text().toLongLong(NULL, 10);
std::cout << b << std::endl;
ublas::vector<long double> x = ortho_solve (A, b);
std::cout << x << std::endl;
for (int i = 0; i < size; i++)
{
QTableWidgetItem item;
QString s = QString::number(x(i), 10);
std::cout << s.toStdString() << std::endl;
item.setTextAlignment(Qt::AlignCenter);
item.setText(s);
item.setTextColor(QColor(0, 255, 0));
this->ui->tableWidget_4->setItem(i, 0, &item);
}
}
void MainWindow::on_lineEdit_textChanged(QString s)
{
int size = s.toInt();
if (size < 1) return;
ui->tableWidget->setRowCount(size);
ui->tableWidget->setColumnCount(size);
ui->tableWidget_3->setRowCount(size);
ui->tableWidget_4->setRowCount(size);
}
Reszta kodu jest mało istotna. (main.cpp to same app.exec() itd.)