HyperHDG
lapack.hxx
Go to the documentation of this file.
1 /*!*************************************************************************************************
2  * \file lapack.hxx
3  * \brief This file provides the function lapack_solve.
4  *
5  * This is a wrapper file to provide LAPACK based functions that have the ability to solve (dense)
6  * local systems of linear equations in an efficient way. The functions \c daxpy_ and \c dnrm2_ are
7  * needed to provide the three functions that solve local systems of equations. From these three
8  * functions, one is chosen to be used in the remainder of the code (i.e., \c dgesv_ cf. LAPACK
9  * manual for further details).
10  *
11  * \authors Guido Kanschat, Heidelberg University, 2019--2020.
12  * \authors Andreas Rupp, Heidelberg University, 2019--2020.
13  **************************************************************************************************/
14 
15 #pragma once // Ensure that file is included only once in a single compilation.
16 
17 #include <array>
18 // #include <exception>
19 
20 namespace Wrapper
21 {
22 /*!*************************************************************************************************
23  * \brief Exception to be thrown if LAPACK's solve fails.
24  **************************************************************************************************/
25 struct LAPACKexception : public std::exception
26 {
27  const char* what() const throw() { return "LAPACK's function failed!"; }
28 };
29 
30 // -------------------------------------------------------------------------------------------------
31 // -------------------------------------------------------------------------------------------------
32 //
33 // MAIN FUNCTIONS THAT IMPLEMENT MATRIX OPERATIONS:
34 // Only the functions within this section are supposed to be used outside of this file!
35 //
36 // -------------------------------------------------------------------------------------------------
37 // -------------------------------------------------------------------------------------------------
38 
39 /*!*************************************************************************************************
40  * \brief Solve local system of equations.
41  *
42  * Solve linear (dense) system of equations \f$Ax=b\f$, where \$A\$ is an \f$n \times n\f$ square
43  * matrix, which enters as a \c std::array of \c double, \f$n\f$ is provided via \c system_size, and
44  * the \c std::array of \c double \c rhs_b is both, input (i.e., \f$b\f$) and output (i.e., \f$x\f$)
45  * of the function.
46  *
47  * Independent of \c const expressions of functions using \c lapack_solve one should not use
48  * \c mat_a after calling this function. The input that has been within \c rhs_b will have been
49  * replaced by the solution of the system of equations.
50  *
51  * \tparam system_size Size of the system of equations.
52  * \tparam n_rhs_cols Number of columns of the right hand side matrix. Defaults to 1.
53  * \tparam float_t Floating type which this function should be executed with. Only \c float and
54  * \c double are supported.
55  * \param dense_mat Array comprising the matrix describing the linear system of equations.
56  * \param rhs Array comprising the right-hand side of the system.
57  * \retval rhs_b Array comprising the solution of the system of equations.
58  **************************************************************************************************/
59 template <unsigned int system_size, unsigned int n_rhs_cols = 1, typename lapack_float_t>
60 std::array<lapack_float_t, system_size * n_rhs_cols> lapack_solve(
61  std::array<lapack_float_t, system_size * system_size>& dense_mat,
62  std::array<lapack_float_t, system_size * n_rhs_cols>& rhs);
63 /*!*************************************************************************************************
64  * \brief Determinant of a rectangular system.
65  *
66  * Calculate the generalized determinant of a rectangular matrix. If the matrix is square, this is
67  * the standard determinant. The determinant is determined by doing a QR decomposition based on
68  * Householder transformations. Thus det(Q) = +1, if the number of rows if even and det(Q) = -1, if
69  * the number of rows is odd. This number is multiplied by the diagonal entries of R.
70  *
71  * The matrix is destroyed using this function. Its entries might be used to generate matrices Q and
72  * R of the QR descomposition.
73  *
74  * \tparam n_rows Number of rows of the matrix whose determinant should be calculated.
75  * \tparam n_cols Number of columns of the matrix whose determinant should be calculated.
76  * \tparam float_t Floating type which this function should be executed with. Only \c float and
77  * \c double are supported.
78  * \param dense_mat Array comprising the matrix describing the linear system of equations.
79  * \retval determinant Generalized determinant of the matrix.
80  **************************************************************************************************/
81 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
82 lapack_float_t lapack_det(std::array<lapack_float_t, n_rows * n_cols>& dense_mat);
83 /*!*************************************************************************************************
84  * \brief Matrix Q of QR decomposition.
85  *
86  * Return matrix Q of the Householder QR decomposition of the matrix. The matrix is destroyed using
87  * this function. Its entries might be used to generate matrices Q and R of the QR descomposition.
88  *
89  * \tparam n_rows Number of rows of the matrix whose determinant should be calculated.
90  * \tparam n_cols Number of columns of the matrix whose determinant should be calculated.
91  * \tparam float_t Floating type which this function should be executed with. Only \c float and
92  * \c double are supported.
93  * \param dense_mat Array comprising the matrix describing the linear system of equations.
94  * \retval mat_q Matrix Q of Householder QR decomposition.
95  **************************************************************************************************/
96 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
97 std::array<lapack_float_t, n_rows * n_rows> lapack_qr_decomp_q(
98  std::array<lapack_float_t, n_rows * n_cols>& dense_mat);
99 /*!*************************************************************************************************
100  * \brief Matrix R of QR decomposition.
101  *
102  * Return matrix R of the Householder QR decomposition of the matrix. The matrix is destroyed using
103  * this function. Its entries are the same as the ones of R of the QR descomposition.
104  *
105  * \tparam n_rows Number of rows of the matrix whose determinant should be calculated.
106  * \tparam n_cols Number of columns of the matrix whose determinant should be calculated.
107  * \tparam float_t Floating type which this function should be executed with. Only \c float and
108  * \c double are supported.
109  * \param dense_mat Array comprising the matrix describing the linear system of equations.
110  * \retval dense_mat Matrix R of Householder QR decomposition.
111  * \retval mat_r Matrix R of Householder QR decomposition.
112  **************************************************************************************************/
113 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
114 std::array<lapack_float_t, n_rows * n_cols>& lapack_qr_decomp_r(
115  std::array<lapack_float_t, n_rows * n_cols>& dense_mat);
116 /*!*************************************************************************************************
117  * \brief Matrices Q and R of QR decomposition.
118  *
119  * Return matriices Q and R of the Householder QR decomposition of the matrix.
120  *
121  * \tparam n_rows Number of rows of the matrix whose determinant should be calculated.
122  * \tparam n_cols Number of columns of the matrix whose determinant should be calculated.
123  * \tparam float_t Floating type which this function should be executed with. Only \c float and
124  * \c double are supported.
125  * \param dense_mat Array comprising the matrix describing the linear system of equations.
126  * \param mat_q Reference to empy matrix to be filled.
127  * \retval dense_mat Matrix R of Householder QR decomposition.
128  * \retval mat_q Matrix Q of Householder QR decomposition.
129  **************************************************************************************************/
130 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
131 void lapack_qr_decomp(std::array<lapack_float_t, n_rows * n_cols>& dense_mat,
132  std::array<lapack_float_t, n_rows * n_rows>& mat_q);
133 /*!*************************************************************************************************
134  * \brief Matrices Q and R of QR decomposition.
135  *
136  * Return matriices Q and R of the Householder QR decomposition of the matrix.
137  *
138  * \tparam n_rows Number of rows of the matrix whose determinant should be calculated.
139  * \tparam n_cols Number of columns of the matrix whose determinant should be calculated.
140  * \tparam float_t Floating type which this function should be executed with. Only \c float and
141  * \c double are supported.
142  * \param dense_mat Array comprising the matrix describing the linear system of equations.
143  * \param mat_q Reference to empy matrix to be filled.
144  * \param mat_r Reference to empy matrix to be filled.
145  * \retval dense_mat Matrix R of Householder QR decomposition.
146  * \retval mat_q Matrix Q of Householder QR decomposition.
147  * \retval mat_r Square system of size n_cols that contains respective part of R.
148  **************************************************************************************************/
149 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
150 void lapack_qr_decomp(std::array<lapack_float_t, n_rows * n_cols>& dense_mat,
151  std::array<lapack_float_t, n_rows * n_rows>& mat_q,
152  std::array<lapack_float_t, n_cols * n_cols>& mat_r);
153 
154 // -------------------------------------------------------------------------------------------------
156 // -------------------------------------------------------------------------------------------------
157 
158 // -------------------------------------------------------------------------------------------------
159 // -------------------------------------------------------------------------------------------------
160 //
161 // AUXILIARY FUNCTIONS THAT DIRECTLY USE LAPACK ROUTINES: Not to be used outside of this file!
162 //
163 // -------------------------------------------------------------------------------------------------
164 // -------------------------------------------------------------------------------------------------
165 
166 extern "C"
167 {
168  /*!***********************************************************************************************
169  * \brief This function is not (never) to be used.
170  *
171  * This function is \b not to be used in regular code. It only / solely is defined to allow the
172  * use of functions \c lapack_solve that will be implemented below.
173  ************************************************************************************************/
174  void daxpy_(int* n, double* alpha, double* dx, int* incx, double* dy, int* incy);
175  /*!***********************************************************************************************
176  * \brief This function is not (never) to be used.
177  *
178  * This function is \b not to be used in regular code. It only / solely is defined to allow the
179  * use of functions \c lapack_solve that will be implemented below.
180  ************************************************************************************************/
181  double dnrm2_(int* n, double* x, int* incx);
182  /*!***********************************************************************************************
183  * \brief This function is not (never) to be used.
184  *
185  * This function is \b not to be used in regular code. It only / solely is defined to allow the
186  * use of functions \c lapack_solve that will be implemented below.
187  ************************************************************************************************/
188  void dgetrf_(int* M, int* N, double* A, int* lda, int* IPIV, int* INFO);
189  /*!***********************************************************************************************
190  * \brief This function is not (never) to be used.
191  *
192  * This function is \b not to be used in regular code. It only / solely is defined to allow the
193  * use of functions \c lapack_solve that will be implemented below.
194  ************************************************************************************************/
195  void dgetrs_(char* C,
196  int* N,
197  int* NRHS,
198  double* A,
199  int* LDA,
200  int* IPIV,
201  double* B,
202  int* LDB,
203  int* INFO);
204  /*!***********************************************************************************************
205  * \brief This function is not (never) to be used.
206  *
207  * This function is \b not to be used in regular code. It only / solely is defined to allow the
208  * use of functions \c lapack_solve that will be implemented below.
209  ************************************************************************************************/
210  void dgesv_(int* n, int* nrhs, double* a, int* lda, int* ipiv, double* b, int* ldb, int* info);
211  /*!***********************************************************************************************
212  * \brief This function is not (never) to be used.
213  *
214  * This function is \b not to be used in regular code. It only / solely is defined to allow the
215  * use of functions \c lapack_solve that will be implemented below.
216  ************************************************************************************************/
217  void dgeqr2_(int* m, int* n, double* a, int* lda, double* tau, double* work, int* info);
218 
219  /*!***********************************************************************************************
220  * \brief This function is not (never) to be used.
221  *
222  * This function is \b not to be used in regular code. It only / solely is defined to allow the
223  * use of functions \c lapack_solve that will be implemented below.
224  ************************************************************************************************/
225  void saxpy_(int* n, float* alpha, float* dx, int* incx, float* dy, int* incy);
226  /*!***********************************************************************************************
227  * \brief This function is not (never) to be used.
228  *
229  * This function is \b not to be used in regular code. It only / solely is defined to allow the
230  * use of functions \c lapack_solve that will be implemented below.
231  ************************************************************************************************/
232  double snrm2_(int* n, float* x, int* incx);
233  /*!***********************************************************************************************
234  * \brief This function is not (never) to be used.
235  *
236  * This function is \b not to be used in regular code. It only / solely is defined to allow the
237  * use of functions \c lapack_solve that will be implemented below.
238  ************************************************************************************************/
239  void sgetrf_(int* M, int* N, float* A, int* lda, int* IPIV, int* INFO);
240  /*!***********************************************************************************************
241  * \brief This function is not (never) to be used.
242  *
243  * This function is \b not to be used in regular code. It only / solely is defined to allow the
244  * use of functions \c lapack_solve that will be implemented below.
245  ************************************************************************************************/
246  void
247  sgetrs_(char* C, int* N, int* NRHS, float* A, int* LDA, int* IPIV, float* B, int* LDB, int* INFO);
248  /*!***********************************************************************************************
249  * \brief This function is not (never) to be used.
250  *
251  * This function is \b not to be used in regular code. It only / solely is defined to allow the
252  * use of functions \c lapack_solve that will be implemented below.
253  ************************************************************************************************/
254  void sgesv_(int* n, int* nrhs, float* a, int* lda, int* ipiv, float* b, int* ldb, int* info);
255  /*!***********************************************************************************************
256  * \brief This function is not (never) to be used.
257  *
258  * This function is \b not to be used in regular code. It only / solely is defined to allow the
259  * use of functions \c lapack_solve that will be implemented below.
260  ************************************************************************************************/
261  void sgeqr2_(int* m, int* n, float* a, int* lda, float* tau, float* work, int* info);
262 } // end of extern "C"
263 
264 // -------------------------------------------------------------------------------------------------
266 // -------------------------------------------------------------------------------------------------
267 
268 /*!*************************************************************************************************
269  * \brief Solve local system of equations with \c double floating point numbers --- DO NOT USE.
270  *
271  * Solve linear (dense) system of equations \f$Ax=b\f$, where \$A\$ is an \f$n \times n\f$ square
272  * matrix, which enters as a pointer to \c double, \f$n\f$ is provided via \c system_size, and the
273  * \c double pointer \c rhs_b is both, input (i.e., \f$b\f$) and output (i.e., \f$x\f$) of the
274  * function.
275  *
276  * Independent of \c const expressions of functions using \c lapack_solve one should not use
277  * \c mat_a after calling this function. The input that has been within \c rhs_b will have been
278  * replaced by the solution of the system of equations.
279  *
280  * \param system_size Size of the system of equations.
281  * \param n_rhs_cols Number of right-hand sides the system is solved for.
282  * \param mat_a Pointer to the matrix describing the linear system of equations.
283  * \param rhs_b Pointer to the right-hand side of the system.
284  * \retval rhs_b Pointer to the solution of the system of equations.
285  **************************************************************************************************/
286 inline void lapack_solve(int system_size, int n_rhs_cols, double* mat_a, double* rhs_b)
287 {
288  int info = -1;
289  int* ipiv = new int[system_size];
290  dgesv_(&system_size, &n_rhs_cols, mat_a, &system_size, ipiv, rhs_b, &system_size, &info);
291  delete[] ipiv;
292  if (info != 0)
293  throw LAPACKexception();
294 }
295 /*!*************************************************************************************************
296  * \brief Solve local system of equations with \c float floating point numbers --- DO NOT USE.
297  *
298  * Solve linear (dense) system of equations \f$Ax=b\f$, where \$A\$ is an \f$n \times n\f$ square
299  * matrix, which enters as a pointer to \c float, \f$n\f$ is provided via \c system_size, and the
300  * \c float pointer \c rhs_b is both, input (i.e., \f$b\f$) and output (i.e., \f$x\f$) of the
301  * function.
302  *
303  * Independent of \c const expressions of functions using \c lapack_solve one should not use
304  * \c mat_a after calling this function. The input that has been within \c rhs_b will have been
305  * replaced by the solution of the system of equations.
306  *
307  * \param system_size Size of the system of equations.
308  * \param n_rhs_cols Number of right-hand sides the system is solved for.
309  * \param mat_a Pointer to the matrix describing the linear system of equations.
310  * \param rhs_b Pointer to the right-hand side of the system.
311  * \retval rhs_b Pointer to the solution of the system of equations.
312  **************************************************************************************************/
313 inline void lapack_solve(int system_size, int n_rhs_cols, float* mat_a, float* rhs_b)
314 {
315  int info = -1;
316  int* ipiv = new int[system_size];
317  sgesv_(&system_size, &n_rhs_cols, mat_a, &system_size, ipiv, rhs_b, &system_size, &info);
318  delete[] ipiv;
319  if (info != 0)
320  throw LAPACKexception();
321 }
322 /*!*************************************************************************************************
323  * \brief QR decomposition in \c double floating point arithmetic --- DO NOT USE.
324  *
325  * Perform QR decomposition of a given matrix using Householder transformations. The matrices Q and
326  * R are stored according to the LAPACK function in an encoded way in matrix mat_a and vector tau.
327  *
328  * \param n_rows Number of rows of matrix.
329  * \param n_cols Number of columns of matrix.
330  * \param mat_a Matrix that is to be decomposed.
331  * \param tau Space for return value, i.e. vector with auxiliary darta.
332  * \retval mat_a Encoded QR decomposition of the matrix.
333  * \retval tau Auxiliary parameters needed to reconstruct matrix Q.
334  **************************************************************************************************/
335 inline void lapack_qr(int n_rows, int n_cols, double* mat_a, double* tau)
336 {
337  int info = -1;
338  double* work = new double[n_cols];
339  dgeqr2_(&n_rows, &n_cols, mat_a, &n_rows, tau, work, &info);
340  delete[] work;
341  if (info != 0)
342  throw LAPACKexception();
343 }
344 /*!*************************************************************************************************
345  * \brief QR decomposition in \c float floating point arithmetic --- DO NOT USE.
346  *
347  * Perform QR decomposition of a given matrix using Householder transformations. The matrices Q and
348  * R are stored according to the LAPACK function in an encoded way in matrix mat_a and vector tau.
349  *
350  * \param n_rows Number of rows of matrix.
351  * \param n_cols Number of columns of matrix.
352  * \param mat_a Matrix that is to be decomposed.
353  * \param tau Space for return value, i.e. vector with auxiliary darta.
354  * \retval mat_a Encoded QR decomposition of the matrix.
355  * \retval tau Auxiliary parameters needed to reconstruct matrix Q.
356  **************************************************************************************************/
357 inline void lapack_qr(int n_rows, int n_cols, float* mat_a, float* tau)
358 {
359  int info = -1;
360  float* work = new float[n_cols];
361  sgeqr2_(&n_rows, &n_cols, mat_a, &n_rows, tau, work, &info);
362  delete[] work;
363  if (info != 0)
364  throw LAPACKexception();
365 }
366 
367 // -------------------------------------------------------------------------------------------------
368 // -------------------------------------------------------------------------------------------------
369 //
370 // IMPLEMENTATION OF MAIN FUNCTIONS THAT DO NOT NEED A DENSE LINEAR ALGEBRA IMPLEMENTATION:
371 // Functions have been declared above and are only implemented here!
372 //
373 // -------------------------------------------------------------------------------------------------
374 // -------------------------------------------------------------------------------------------------
375 
376 // -------------------------------------------------------------------------------------------------
377 // lapack_solve
378 // -------------------------------------------------------------------------------------------------
379 
380 template <unsigned int system_size, unsigned int n_rhs_cols, typename lapack_float_t>
381 std::array<lapack_float_t, system_size * n_rhs_cols> lapack_solve(
382  std::array<lapack_float_t, system_size * system_size>& dense_mat,
383  std::array<lapack_float_t, system_size * n_rhs_cols>& rhs)
384 {
385  lapack_solve(system_size, n_rhs_cols, dense_mat.data(), rhs.data());
386  return rhs;
387 }
388 
389 // -------------------------------------------------------------------------------------------------
390 // lapack_det
391 // -------------------------------------------------------------------------------------------------
392 
393 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
394 lapack_float_t lapack_det(std::array<lapack_float_t, n_rows * n_cols>& dense_mat)
395 {
396  constexpr unsigned int rank = std::min(n_rows, n_cols);
397  std::array<lapack_float_t, rank> tau;
398  lapack_qr(n_rows, n_cols, dense_mat.data(), tau.data());
399 
400  lapack_float_t determinant = 1.;
401  for (unsigned int i = 0; i < rank; ++i)
402  determinant *= -dense_mat[i * (n_rows + 1)];
403  return determinant;
404 }
405 
406 } // end of namespace Wrapper
407 
408 // -------------------------------------------------------------------------------------------------
409 // -------------------------------------------------------------------------------------------------
410 //
411 // IMPLEMENTATION OF AUXILIARY FUNCTIONS THAT REQUIRE A DENSE LA IMPLEMENTATION:
412 // Do not use these functions outside of this file.
413 //
414 // -------------------------------------------------------------------------------------------------
415 // -------------------------------------------------------------------------------------------------
416 
417 #include <HyperHDG/dense_la.hxx> // Dense linear algebra that is utilized in the following.
418 
419 namespace Wrapper
420 {
421 /*!*************************************************************************************************
422  * \brief Matrix Q of QR decomposition --- DO NOT USE.
423  *
424  * Return matrix Q of the Householder QR decomposition of the matrix. The matrix is supposed to
425  * contain the Householder QR decomposition which is encoded in the style of LAPACK. The function
426  * returns matrix Q in standard encoding.
427  *
428  * \tparam n_rows Number of rows of the matrix whose determinant should be calculated.
429  * \tparam n_cols Number of columns of the matrix whose determinant should be calculated.
430  * \tparam float_t Floating type which this function should be executed with. Only \c float and
431  * \c double are supported.
432  * \param dense_mat Array comprising the QR decomposed matrix.
433  * \param tau Auxiliary data from LAPACK to encode matrix Q.
434  * \retval mat_q Matrix Q of Householder QR decomposition.
435  **************************************************************************************************/
436 template <unsigned int n_rows, unsigned int n_cols, unsigned int rank, typename lapack_float_t>
437 inline std::array<lapack_float_t, n_rows * n_rows> get_q_from_lapack_qr_result(
438  const std::array<lapack_float_t, n_rows * n_cols>& dense_mat,
439  const std::array<lapack_float_t, rank>& tau)
440 {
441  SmallMat unity = diagonal<n_rows, n_rows, lapack_float_t>(1.), matQ = unity;
442  SmallVec<n_rows, lapack_float_t> vec;
443 
444  if constexpr (rank > 0)
445  for (unsigned int i = 0; i < rank; ++i) // i is column index column index, here!
446  {
447  for (unsigned int j = 0; j < n_rows; ++j) // j is row index, here!
448  if (j < i)
449  vec[j] = 0.;
450  else if (j == i)
451  vec[j] = 1.;
452  else
453  vec[j] = dense_mat[i * n_rows + j];
454  matQ = matQ * (unity - tau[i] * dyadic_product(vec, vec));
455  }
456 
457  return matQ.data();
458 }
459 /*!*************************************************************************************************
460  * \brief Matrix Q of QR decomposition --- DO NOT USE.
461  *
462  * Return matrix Q of the Householder QR decomposition of the matrix. The matrix is supposed to
463  * contain the Householder QR decomposition which is encoded in the style of LAPACK. The function
464  * returns matrix Q in standard encoding.
465  *
466  * \tparam n_rows Number of rows of the matrix whose determinant should be calculated.
467  * \tparam n_cols Number of columns of the matrix whose determinant should be calculated.
468  * \tparam float_t Floating type which this function should be executed with. Only \c float and
469  * \c double are supported.
470  * \param dense_mat Array comprising the QR decomposed matrix.
471  * \param tau Auxiliary data from LAPACK to encode matrix Q.
472  * \param mat_q Matrix to be filled with entries of the return value.
473  * \retval mat_q Matrix Q of Householder QR decomposition.
474  **************************************************************************************************/
475 template <unsigned int n_rows, unsigned int n_cols, unsigned int rank, typename lapack_float_t>
476 inline void get_q_from_lapack_qr_result(
477  const std::array<lapack_float_t, n_rows * n_cols>& dense_mat,
478  const std::array<lapack_float_t, rank>& tau,
479  std::array<lapack_float_t, n_rows * n_rows>& mat_q)
480 {
481  SmallMat unity = diagonal<n_rows, n_rows, lapack_float_t>(1.);
482  SmallMat<n_rows, n_rows, lapack_float_t> matQ(std::move(mat_q));
483  SmallVec<n_rows, lapack_float_t> vec;
484  matQ = unity;
485 
486  for (unsigned int i = 0; i < rank; ++i) // i is column index column index, here!
487  {
488  for (unsigned int j = 0; j < n_rows; ++j) // j is row index, here!
489  if (j < i)
490  vec[j] = 0.;
491  else if (j == i)
492  vec[j] = 1.;
493  else
494  vec[j] = dense_mat[i * n_rows + j];
495  matQ = matQ * (unity - tau[i] * dyadic_product(vec, vec));
496  }
497 
498  mat_q = std::move(matQ.data());
499 }
500 /*!*************************************************************************************************
501  * \brief Matrix R of QR decomposition --- DO NOT USE.
502  *
503  * Return matrix R of the Householder QR decomposition of the matrix. The matrix is supposed to
504  * contain the Householder QR decomposition which is encoded in the style of LAPACK. The function
505  * returns matrix R in standard encoding.
506  *
507  * The function also manipulates the input matrix to be R.
508  *
509  * \tparam n_rows Number of rows of the matrix whose determinant should be calculated.
510  * \tparam n_cols Number of columns of the matrix whose determinant should be calculated.
511  * \tparam float_t Floating type which this function should be executed with. Only \c float and
512  * \c double are supported.
513  * \param dense_mat Array comprising the QR decomposed matrix.
514  * \retval mat_r Matrix R of Householder QR decomposition.
515  **************************************************************************************************/
516 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
517 inline std::array<lapack_float_t, n_rows * n_cols>& get_r_from_lapack_qr_result(
518  std::array<lapack_float_t, n_rows * n_cols>& dense_mat)
519 {
520  for (unsigned int i = 0; i < n_cols; ++i) // i is column index, here!
521  for (unsigned int j = 0; j < n_rows; ++j) // j is row index, here!
522  if (j > i)
523  dense_mat[i * n_rows + j] = 0.;
524  return dense_mat;
525 }
526 /*!*************************************************************************************************
527  * \brief Matrix R of QR decomposition --- DO NOT USE.
528  *
529  * Return matrix R of the Householder QR decomposition of the matrix. The matrix is supposed to
530  * contain the Householder QR decomposition which is encoded in the style of LAPACK. The function
531  * returns matrix R in standard encoding, but with dimensions \c n_cols times \c n_cols.
532  *
533  * \tparam n_rows Number of rows of the matrix whose determinant should be calculated.
534  * \tparam n_cols Number of columns of the matrix whose determinant should be calculated.
535  * \tparam float_t Floating type which this function should be executed with. Only \c float and
536  * \c double are supported.
537  * \param lapack_mat Array comprising the QR decomposed matrix.
538  * \param mat_r Matrix to be filled with the entries of the return value.
539  * \retval mat_r Matrix R of Householder QR decomposition.
540  **************************************************************************************************/
541 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
542 inline void get_r_from_lapack_qr_result(
543  const std::array<lapack_float_t, n_rows * n_cols>& lapack_mat,
544  std::array<lapack_float_t, n_cols * n_cols>& mat_r)
545 {
546  static_assert(n_rows >= n_cols, "Function only defined for these matrices!");
547  for (unsigned int i = 0; i < n_cols; ++i) // i is column index, here!
548  for (unsigned int j = 0; j < n_cols; ++j) // j is row index, here!
549  if (j > i)
550  mat_r[i * n_rows + j] = 0.;
551  else
552  mat_r[i * n_cols + j] = lapack_mat[i * n_rows + j];
553 }
554 
555 // -------------------------------------------------------------------------------------------------
556 // -------------------------------------------------------------------------------------------------
557 //
558 // IMPLEMENTATION OF MAIN FUNCTIONS THAT NEED A DENSE LINEAR ALGEBRA IMPLEMENTATION:
559 // Functions have been declared above and are only implemented here!
560 //
561 // -------------------------------------------------------------------------------------------------
562 // -------------------------------------------------------------------------------------------------
563 
564 // -------------------------------------------------------------------------------------------------
565 // lapack_qr_decomp_q
566 // -------------------------------------------------------------------------------------------------
567 
568 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
569 std::array<lapack_float_t, n_rows * n_rows> lapack_qr_decomp_q(
570  std::array<lapack_float_t, n_rows * n_cols>& dense_mat)
571 {
572  constexpr unsigned int rank = std::min(n_rows, n_cols);
573  std::array<lapack_float_t, rank> tau;
574  lapack_qr(n_rows, n_cols, dense_mat.data(), tau.data());
575  return get_q_from_lapack_qr_result<n_rows, n_cols, rank, lapack_float_t>(dense_mat, tau);
576 }
577 
578 // -------------------------------------------------------------------------------------------------
579 // lapack_qr_decomp_r
580 // -------------------------------------------------------------------------------------------------
581 
582 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
583 std::array<lapack_float_t, n_rows * n_cols>& lapack_qr_decomp_r(
584  std::array<lapack_float_t, n_rows * n_cols>& dense_mat)
585 {
586  constexpr unsigned int rank = std::min(n_rows, n_cols);
587  std::array<lapack_float_t, rank> tau;
588  lapack_qr(n_rows, n_cols, dense_mat.data(), tau.data());
589  return get_r_from_lapack_qr_result<n_rows, n_cols, rank, lapack_float_t>(dense_mat);
590 }
591 
592 // -------------------------------------------------------------------------------------------------
593 // lapack_qr_decomp
594 // -------------------------------------------------------------------------------------------------
595 
596 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
597 void lapack_qr_decomp(std::array<lapack_float_t, n_rows * n_cols>& dense_mat,
598  std::array<lapack_float_t, n_rows * n_rows>& mat_q)
599 {
600  constexpr unsigned int rank = std::min(n_rows, n_cols);
601  std::array<lapack_float_t, rank> tau;
602  lapack_qr(n_rows, n_cols, dense_mat.data(), tau.data());
603 
604  get_q_from_lapack_qr_result<n_rows, n_cols, rank, lapack_float_t>(dense_mat, tau, mat_q);
605  get_r_from_lapack_qr_result<n_rows, n_cols, lapack_float_t>(dense_mat);
606 }
607 
608 // -------------------------------------------------------------------------------------------------
609 // lapack_qr_decomp
610 // -------------------------------------------------------------------------------------------------
611 
612 template <unsigned int n_rows, unsigned int n_cols, typename lapack_float_t>
613 void lapack_qr_decomp(std::array<lapack_float_t, n_rows * n_cols>& dense_mat,
614  std::array<lapack_float_t, n_rows * n_rows>& mat_q,
615  std::array<lapack_float_t, n_cols * n_cols>& mat_r)
616 {
617  constexpr unsigned int rank = std::min(n_rows, n_cols);
618  std::array<lapack_float_t, rank> tau;
619  lapack_qr(n_rows, n_cols, dense_mat.data(), tau.data());
620 
621  get_q_from_lapack_qr_result<n_rows, n_cols, rank, lapack_float_t>(dense_mat, tau, mat_q);
622  get_r_from_lapack_qr_result<n_rows, n_cols, lapack_float_t>(dense_mat, mat_r);
623  get_r_from_lapack_qr_result<n_rows, n_cols, lapack_float_t>(dense_mat);
624 }
625 
626 } // end of namespace Wrapper
Wrapper::get_r_from_lapack_qr_result
std::array< lapack_float_t, n_rows *n_cols > & get_r_from_lapack_qr_result(std::array< lapack_float_t, n_rows *n_cols > &dense_mat)
Matrix R of QR decomposition — DO NOT USE.
Definition: lapack.hxx:517
Wrapper
A namespace containing the different wrapper functions.
Definition: lapack.hxx:20
Wrapper::lapack_qr_decomp_r
std::array< lapack_float_t, n_rows *n_cols > & lapack_qr_decomp_r(std::array< lapack_float_t, n_rows *n_cols > &dense_mat)
Matrix R of QR decomposition.
Definition: lapack.hxx:583
Wrapper::get_q_from_lapack_qr_result
std::array< lapack_float_t, n_rows *n_rows > get_q_from_lapack_qr_result(const std::array< lapack_float_t, n_rows *n_cols > &dense_mat, const std::array< lapack_float_t, rank > &tau)
Matrix Q of QR decomposition — DO NOT USE.
Definition: lapack.hxx:437
Wrapper::LAPACKexception::what
const char * what() const
Definition: lapack.hxx:27
dyadic_product
SmallMat< n_rows, n_cols, mat_entry_t > dyadic_product(const SmallMat< n_rows, 1, mat_entry_t > &left, const SmallMat< n_cols, 1, mat_entry_t > &right)
Create dyadic product of two small vectors.
Definition: dense_la.hxx:538
Wrapper::lapack_det
lapack_float_t lapack_det(std::array< lapack_float_t, n_rows *n_cols > &dense_mat)
Determinant of a rectangular system.
Definition: lapack.hxx:394
work
the intent is to exercise the right to control the distribution of derivative or collective works based on the Library In mere aggregation of another work not based on the Library with the you must alter all the notices that refer to this so that they refer to the ordinary GNU General Public instead of to this it is irreversible for that so the ordinary GNU General Public License applies to all subsequent copies and derivative works made from that copy This option is useful when you wish to copy part of the code of the Library into a program that is not a library You may copy and distribute the which must be distributed under the terms of Sections and above on a medium customarily used for software interchange If distribution of object code is made by offering access to copy from a designated then offering equivalent access to copy the source code from the same place satisfies the requirement to distribute the source even though third parties are not compelled to copy the source along with the object code A program that contains no derivative of any portion of the but is designed to work with the Library by being compiled or linked with is called a work that uses the Library Such a work
Definition: License.txt:243
Wrapper::lapack_solve
std::array< lapack_float_t, system_size *n_rhs_cols > lapack_solve(std::array< lapack_float_t, system_size *system_size > &dense_mat, std::array< lapack_float_t, system_size *n_rhs_cols > &rhs)
Solve local system of equations.
Definition: lapack.hxx:381
Wrapper::lapack_qr_decomp_q
std::array< lapack_float_t, n_rows *n_rows > lapack_qr_decomp_q(std::array< lapack_float_t, n_rows *n_cols > &dense_mat)
Matrix Q of QR decomposition.
Definition: lapack.hxx:569
Wrapper::lapack_qr
void lapack_qr(int n_rows, int n_cols, double *mat_a, double *tau)
QR decomposition in double floating point arithmetic — DO NOT USE.
Definition: lapack.hxx:335
exception
if the work is an executable linked with the with the complete machine readable work that uses the as object code and or source so that the user can modify the Library and then relink to produce a modified executable containing the modified rather than copying library functions into the if the user installs as long as the modified version is interface compatible with the version that the work was made with c Accompany the work with a written valid for at least three to give the same user the materials specified in for a charge no more than the cost of performing this distribution d If distribution of the work is made by offering access to copy from a designated offer equivalent access to copy the above specified materials from the same place e Verify that the user has already received a copy of these materials or that you have already sent this user a copy For an the required form of the work that uses the Library must include any data and utility programs needed for reproducing the executable from it as a special exception
Definition: License.txt:320
SmallMat::data
std::array< mat_entry_t, size()> & data()
Return data array that allows to manipulate the SmallMat.
Definition: dense_la.hxx:190
diffusion_uniform.system_size
system_size
Definition: diffusion_uniform.py:38
Wrapper::LAPACKexception
Exception to be thrown if LAPACK's solve fails.
Definition: lapack.hxx:25
Wrapper::lapack_qr_decomp
void lapack_qr_decomp(std::array< lapack_float_t, n_rows *n_cols > &dense_mat, std::array< lapack_float_t, n_rows *n_rows > &mat_q)
Matrices Q and R of QR decomposition.
Definition: lapack.hxx:597
dense_la.hxx
diffusion_uniform.A
A
Definition: diffusion_uniform.py:39
determinant
mat_entry_t determinant(SmallMat< n_rows, n_cols, mat_entry_t > &mat)
Determinant of a rectangular system.
Definition: dense_la.hxx:1121
SmallMat
This class implements a small/dense matrix.
Definition: dense_la.hxx:34
Generated by   doxygen 1.8.17