Commit c102f059 authored by Tobias WEBER's avatar Tobias WEBER
Browse files

started with lapack functions

parent 72fa47a1
......@@ -3855,6 +3855,89 @@ requires m::is_basic_mat<t_mat> && m::is_dyn_mat<t_mat>
// ----------------------------------------------------------------------------
}
// ----------------------------------------------------------------------------
// lapack wrappers
// ----------------------------------------------------------------------------
#ifdef USE_LAPACK
extern "C"
{
#define lapack_complex_double std::complex<double>
#define lapack_complex_double_real(z) (z.real())
#define lapack_complex_double_imag(z) (z.imag())
#define lapack_complex_float std::complex<float>
#define lapack_complex_float_real(z) (z.real())
#define lapack_complex_float_imag(z) (z.imag())
#include <lapacke.h>
}
namespace m_la {
/**
* QR decomposition of a matrix mat, mat = QR
* http://www.math.utah.edu/software/lapack/lapack-d/dgeqrf.html
* return [ok, Q, R]
*/
template<class t_mat>
std::tuple<bool, t_mat, t_mat> qr(const t_mat& mat)
requires m::is_mat<t_mat>
{
using namespace m_ops;
using T = typename t_mat::value_type;
using t_vec = std::vector<T>;
const std::size_t rows = mat.size1();
const std::size_t cols = mat.size2();
const std::size_t minor = std::min(rows, cols);
const t_mat I = m::unit<t_mat>(minor);
t_mat Q = I, R = mat;
t_vec outmat(rows*cols), outvec(minor);
for(std::size_t i=0; i<rows; ++i)
for(std::size_t j=0; j<cols; ++j)
outmat[i*cols + j] = mat(i, j);
int err = -1;
if constexpr(std::is_same_v<T, float>)
err = LAPACKE_sgeqrf(LAPACK_ROW_MAJOR, rows, cols, outmat.data(), cols, outvec.data());
else if constexpr(std::is_same_v<T, double>)
err = LAPACKE_dgeqrf(LAPACK_ROW_MAJOR, rows, cols, outmat.data(), cols, outvec.data());
for(std::size_t i=0; i<rows; ++i)
for(std::size_t j=0; j<cols; ++j)
R(i, j) = (j>=i ? outmat[i*cols + j] : T{0});
t_vec v = m::zero<t_vec>(minor);
for(std::size_t k=1; k<=minor; ++k)
{
for(std::size_t i=1; i<=k-1; ++i)
v[i-1] = T{0};
v[k-1] = T{1};
for(std::size_t i=k+1; i<=minor; ++i)
v[i-1] = outmat[(i-1)*cols + (k-1)];
Q = Q * (I - outvec[k-1]*m::outer<t_mat, t_vec>(v, v));
}
return std::make_tuple(err == 0, Q, R);
}
}
#endif
// ----------------------------------------------------------------------------
#endif
/**
......
......@@ -2094,6 +2094,7 @@ void MagStructFactDlg::PlotMouseDown(bool left, bool mid, bool right)
void MagStructFactDlg::AfterGLInitialisation()
{
if(!m_plot) return;
SetGLInfos();
// reference sphere and arrow for linked objects
m_sphere = m_plot->GetImpl()->AddSphere(0.05, 0.,0.,0., 1.,1.,1.,1.);
......@@ -2106,13 +2107,6 @@ void MagStructFactDlg::AfterGLInitialisation()
// add all 3d objects
Add3DItem(-1);
// GL device info
auto [strGlVer, strGlShaderVer, strGlVendor, strGlRenderer] = m_plot->GetImpl()->GetGlDescr();
m_labelGlInfos[0]->setText(QString("GL Version: ") + strGlVer.c_str() + QString("."));
m_labelGlInfos[1]->setText(QString("GL Shader Version: ") + strGlShaderVer.c_str() + QString("."));
m_labelGlInfos[2]->setText(QString("GL Vendor: ") + strGlVendor.c_str() + QString("."));
m_labelGlInfos[3]->setText(QString("GL Device: ") + strGlRenderer.c_str() + QString("."));
}
......@@ -2122,6 +2116,7 @@ void MagStructFactDlg::AfterGLInitialisation()
void MagStructFactDlg::AfterGLInitialisationSC()
{
if(!m_plotSC) return;
SetGLInfos();
// reference sphere and arrow for linked objects
m_sphereSC = m_plotSC->GetImpl()->AddSphere(0.05, 0.,0.,0., 1.,1.,1.,1.);
......@@ -2137,6 +2132,31 @@ void MagStructFactDlg::AfterGLInitialisationSC()
}
/**
* set descriptions of the gl device in the info tab
*/
void MagStructFactDlg::SetGLInfos()
{
static bool already_set = 0;
if(already_set) return;
// try whichever gl plotter is available first
for(auto* plot : { m_plot.get(), m_plotSC.get() })
{
if(!plot) continue;
auto [strGlVer, strGlShaderVer, strGlVendor, strGlRenderer] = plot->GetImpl()->GetGlDescr();
m_labelGlInfos[0]->setText(QString("GL Version: ") + strGlVer.c_str() + QString("."));
m_labelGlInfos[1]->setText(QString("GL Shader Version: ") + strGlShaderVer.c_str() + QString("."));
m_labelGlInfos[2]->setText(QString("GL Vendor: ") + strGlVendor.c_str() + QString("."));
m_labelGlInfos[3]->setText(QString("GL Device: ") + strGlRenderer.c_str() + QString("."));
already_set = 1;
break;
}
}
void MagStructFactDlg::closeEvent(QCloseEvent *evt)
{
if(m_sett)
......
......@@ -146,6 +146,7 @@ protected:
void PickerIntersectionSC(const t_vec3_gl* pos, std::size_t objIdx, const t_vec3_gl* posSphere);
void AfterGLInitialisation();
void AfterGLInitialisationSC();
void SetGLInfos();
virtual void closeEvent(QCloseEvent *evt) override;
......
/**
* lapack test
* @author Tobias Weber <tweber@ill.fr>
* @date feb-19
* @license GPLv3, see 'LICENSE' file
*
* g++-8 -std=c++17 -fconcepts -DUSE_LAPACK -I/usr/local/opt/lapack/include -L/usr/local/opt/lapack/lib -o la la.cpp -llapacke
*/
#include <iostream>
#include <vector>
#include "../../libs/_cxx20/math_algos.h"
using namespace m_ops;
//using t_real = double;
using t_real = float;
using t_vec = std::vector<t_real>;
using t_mat = m::mat<t_real, std::vector>;
int main()
{
auto M = m::create<t_mat>({1, 2, 3, 3, 2, 6, 4, 2, 4});
std::cout << "M = " << M << std::endl;
{
auto [Q, R] = m::qr<t_mat, t_vec>(M);
std::cout << "\nQ = " << Q << std::endl;
std::cout << "R = " << R << std::endl;
std::cout << "QR = " << Q*R << std::endl;
}
{
auto [ok, Q, R] = m_la::qr<t_mat>(M);
std::cout << "\nok = " << std::boolalpha << ok << std::endl;
std::cout << "Q = " << Q << std::endl;
std::cout << "R = " << R << std::endl;
std::cout << "QR = " << Q*R << std::endl;
}
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment