GetFEM  5.4.2
gmm_dense_matrix_functions.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2014-2020 Konstantinos Poulios
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file gmm_dense_matrix_functions.h
33  @author Konstantinos Poulios <poulios.konstantinos@gmail.com>
34  @date December 10, 2014.
35  @brief Common matrix functions for dense matrices.
36 */
37 #ifndef GMM_DENSE_MATRIX_FUNCTIONS_H
38 #define GMM_DENSE_MATRIX_FUNCTIONS_H
39 
40 
41 namespace gmm {
42 
43 
44  /**
45  Matrix square root for upper triangular matrices (from GNU Octave).
46  */
47  template <typename T>
48  void sqrtm_utri_inplace(dense_matrix<T>& A)
49  {
50  typedef typename number_traits<T>::magnitude_type R;
51  bool singular = false;
52 
53  // The following code is equivalent to this triple loop:
54  //
55  // n = rows (A);
56  // for j = 1:n
57  // A(j,j) = sqrt (A(j,j));
58  // for i = j-1:-1:1
59  // A(i,j) /= (A(i,i) + A(j,j));
60  // k = 1:i-1;
61  // t storing a A(k,j) -= A(k,i) * A(i,j);
62  // endfor
63  // endfor
64 
65  R tol = R(0); // default_tol(R()) * gmm::mat_maxnorm(A);
66 
67  const size_type n = mat_nrows(A);
68  for (int j=0; j < int(n); j++) {
69  typename dense_matrix<T>::iterator colj = A.begin() + j*n;
70  if (gmm::abs(colj[j]) > tol)
71  colj[j] = gmm::sqrt(colj[j]);
72  else
73  singular = true;
74 
75  for (int i=j-1; i >= 0; i--) {
76  typename dense_matrix<T>::const_iterator coli = A.begin() + i*n;
77  T colji = colj[i] = safe_divide(colj[i], (coli[i] + colj[j]));
78  for (int k = 0; k < i; k++)
79  colj[k] -= coli[k] * colji;
80  }
81  }
82 
83  if (singular)
84  GMM_WARNING1("Matrix is singular, may not have a square root");
85  }
86 
87 
88  template <typename T>
89  void sqrtm(const dense_matrix<std::complex<T> >& A,
90  dense_matrix<std::complex<T> >& SQRTMA)
91  {
92  GMM_ASSERT1(gmm::mat_nrows(A) == gmm::mat_ncols(A),
93  "Matrix square root requires a square matrix");
94  gmm::resize(SQRTMA, gmm::mat_nrows(A), gmm::mat_ncols(A));
95  dense_matrix<std::complex<T> > S(A), Q(A), TMP(A);
96  #if defined(GMM_USES_LAPACK)
97  schur(TMP, S, Q);
98  #else
99  GMM_ASSERT1(false, "Please recompile with lapack and blas librairies "
100  "to use sqrtm matrix function.");
101  #endif
102  sqrtm_utri_inplace(S);
103  gmm::mult(Q, S, TMP);
104  gmm::mult(TMP, gmm::transposed(Q), SQRTMA);
105  }
106 
107  template <typename T>
108  void sqrtm(const dense_matrix<T>& A,
109  dense_matrix<std::complex<T> >& SQRTMA)
110  {
111  dense_matrix<std::complex<T> > cA(mat_nrows(A), mat_ncols(A));
112  gmm::copy(A, gmm::real_part(cA));
113  sqrtm(cA, SQRTMA);
114  }
115 
116  template <typename T>
117  void sqrtm(const dense_matrix<T>& A, dense_matrix<T>& SQRTMA)
118  {
119  dense_matrix<std::complex<T> > cA(mat_nrows(A), mat_ncols(A));
120  gmm::copy(A, gmm::real_part(cA));
121  dense_matrix<std::complex<T> > cSQRTMA(cA);
122  sqrtm(cA, cSQRTMA);
123  gmm::resize(SQRTMA, gmm::mat_nrows(A), gmm::mat_ncols(A));
124  gmm::copy(gmm::real_part(cSQRTMA), SQRTMA);
125 // dense_matrix<std::complex<T1> >::const_reference
126 // it = cSQRTMA.begin(), ite = cSQRTMA.end();
127 // dense_matrix<std::complex<T1> >::reference
128 // rit = SQRTMA.begin();
129 // for (; it != ite; ++it, ++rit) *rit = it->real();
130  }
131 
132 
133  /**
134  Matrix logarithm for upper triangular matrices (from GNU/Octave)
135  */
136  template <typename T>
137  void logm_utri_inplace(dense_matrix<T>& S)
138  {
139  typedef typename number_traits<T>::magnitude_type R;
140 
141  size_type n = gmm::mat_nrows(S);
142  GMM_ASSERT1(n == gmm::mat_ncols(S),
143  "Matrix logarithm is not defined for non-square matrices");
144  for (size_type i=0; i < n-1; ++i)
145  if (gmm::abs(S(i+1,i)) > default_tol(T())) {
146  GMM_ASSERT1(false, "An upper triangular matrix is expected");
147  break;
148  }
149  for (size_type i=0; i < n-1; ++i)
150  if (gmm::real(S(i,i)) <= -default_tol(R()) &&
151  gmm::abs(gmm::imag(S(i,i))) <= default_tol(R())) {
152  GMM_ASSERT1(false, "Principal matrix logarithm is not defined "
153  "for matrices with negative eigenvalues");
154  break;
155  }
156 
157  // Algorithm 11.9 in "Function of matrices", by N. Higham
158  R theta[] = { R(0),R(0),R(1.61e-2),R(5.38e-2),R(1.13e-1),R(1.86e-1),R(2.6429608311114350e-1) };
159 
160  R scaling(1);
161  size_type p(0), m(6), opt_iters(100);
162  for (size_type k=0; k < opt_iters; ++k, scaling *= R(2)) {
163  dense_matrix<T> auxS(S);
164  for (size_type i = 0; i < n; ++i) auxS(i,i) -= R(1);
165  R tau = gmm::mat_norm1(auxS);
166  if (tau <= theta[6]) {
167  ++p;
168  size_type j1(6), j2(6);
169  for (size_type j=0; j < 6; ++j)
170  if (tau <= theta[j]) { j1 = j; break; }
171  for (size_type j=0; j < j1; ++j)
172  if (tau <= 2*theta[j]) { j2 = j; break; }
173  if (j1 - j2 <= 1 || p == 2) { m = j1; break; }
174  }
175  sqrtm_utri_inplace(S);
176  if (k == opt_iters-1)
177  GMM_WARNING1 ("Maximum number of square roots exceeded; "
178  "the calculated matrix logarithm may still be accurate");
179  }
180 
181  for (size_type i = 0; i < n; ++i) S(i,i) -= R(1);
182 
183  if (m > 0) {
184 
185  std::vector<R> nodes, wts;
186  switch(m) {
187  case 0: {
188  R nodes_[] = { R(0.5) };
189  R wts_[] = { R(1) };
190  nodes.assign(nodes_, nodes_+m+1);
191  wts.assign(wts_, wts_+m+1);
192  } break;
193  case 1: {
194  R nodes_[] = { R(0.211324865405187),R(0.788675134594813) };
195  R wts_[] = { R(0.5),R(0.5) };
196  nodes.assign(nodes_, nodes_+m+1);
197  wts.assign(wts_, wts_+m+1);
198  } break;
199  case 2: {
200  R nodes_[] = { R(0.112701665379258),R(0.500000000000000),R(0.887298334620742) };
201  R wts_[] = { R(0.277777777777778),R(0.444444444444444),R(0.277777777777778) };
202  nodes.assign(nodes_, nodes_+m+1);
203  wts.assign(wts_, wts_+m+1);
204  } break;
205  case 3: {
206  R nodes_[] = { R(0.0694318442029737),R(0.3300094782075718),R(0.6699905217924281),R(0.9305681557970263) };
207  R wts_[] = { R(0.173927422568727),R(0.326072577431273),R(0.326072577431273),R(0.173927422568727) };
208  nodes.assign(nodes_, nodes_+m+1);
209  wts.assign(wts_, wts_+m+1);
210  } break;
211  case 4: {
212  R nodes_[] = { R(0.0469100770306681),R(0.2307653449471584),R(0.5000000000000000),R(0.7692346550528415),R(0.9530899229693319) };
213  R wts_[] = { R(0.118463442528095),R(0.239314335249683),R(0.284444444444444),R(0.239314335249683),R(0.118463442528094) };
214  nodes.assign(nodes_, nodes_+m+1);
215  wts.assign(wts_, wts_+m+1);
216  } break;
217  case 5: {
218  R nodes_[] = { R(0.0337652428984240),R(0.1693953067668678),R(0.3806904069584015),R(0.6193095930415985),R(0.8306046932331322),R(0.9662347571015761) };
219  R wts_[] = { R(0.0856622461895853),R(0.1803807865240693),R(0.2339569672863452),R(0.2339569672863459),R(0.1803807865240693),R(0.0856622461895852) };
220  nodes.assign(nodes_, nodes_+m+1);
221  wts.assign(wts_, wts_+m+1);
222  } break;
223  case 6: {
224  R nodes_[] = { R(0.0254460438286208),R(0.1292344072003028),R(0.2970774243113015),R(0.4999999999999999),R(0.7029225756886985),R(0.8707655927996973),R(0.9745539561713792) };
225  R wts_[] = { R(0.0647424830844348),R(0.1398526957446384),R(0.1909150252525594),R(0.2089795918367343),R(0.1909150252525595),R(0.1398526957446383),R(0.0647424830844349) };
226  nodes.assign(nodes_, nodes_+m+1);
227  wts.assign(wts_, wts_+m+1);
228  } break;
229  }
230 
231  dense_matrix<T> auxS1(S), auxS2(S);
232  std::vector<T> auxvec(n);
233  gmm::clear(S);
234  for (size_type j=0; j <= m; ++j) {
235  gmm::copy(gmm::scaled(auxS1, nodes[j]), auxS2);
236  gmm::add(gmm::identity_matrix(), auxS2);
237  // S += wts[i] * auxS1 * inv(auxS2)
238  for (size_type i=0; i < n; ++i) {
239  gmm::copy(gmm::mat_row(auxS1, i), auxvec);
240  gmm::lower_tri_solve(gmm::transposed(auxS2), auxvec, false);
241  gmm::add(gmm::scaled(auxvec, wts[j]), gmm::mat_row(S, i));
242  }
243  }
244  }
245  gmm::scale(S, scaling);
246  }
247 
248  /**
249  Matrix logarithm (from GNU/Octave)
250  */
251  template <typename T>
252  void logm(const dense_matrix<T>& A, dense_matrix<T>& LOGMA)
253  {
254  typedef typename number_traits<T>::magnitude_type R;
255  size_type n = gmm::mat_nrows(A);
256  GMM_ASSERT1(n == gmm::mat_ncols(A),
257  "Matrix logarithm is not defined for non-square matrices");
258  dense_matrix<T> S(A), Q(A);
259  #if defined(GMM_USES_LAPACK)
260  schur(A, S, Q); // A = Q * S * Q^T
261  #else
262  GMM_ASSERT1(false, "Please recompile with lapack and blas librairies "
263  "to use logm matrix function.");
264  #endif
265 
266  bool convert_to_complex(false);
267  if (!is_complex(T()))
268  for (size_type i=0; i < n-1; ++i)
269  if (gmm::abs(S(i+1,i)) > default_tol(T())) {
270  convert_to_complex = true;
271  break;
272  }
273 
274  gmm::resize(LOGMA, n, n);
275  if (convert_to_complex) {
276  dense_matrix<std::complex<R> > cS(n,n), cQ(n,n), auxmat(n,n);
277  gmm::copy(gmm::real_part(S), gmm::real_part(cS));
278  gmm::copy(gmm::real_part(Q), gmm::real_part(cQ));
279  block2x2_reduction(cS, cQ, default_tol(R())*R(3));
280  for (size_type j=0; j < n-1; ++j)
281  for (size_type i=j+1; i < n; ++i)
282  cS(i,j) = T(0);
283  logm_utri_inplace(cS);
284  gmm::mult(cQ, cS, auxmat);
285  gmm::mult(auxmat, gmm::transposed(cQ), cS);
286  // Remove small complex values which may have entered calculation
287  gmm::copy(gmm::real_part(cS), LOGMA);
288 // GMM_ASSERT1(gmm::mat_norm1(gmm::imag_part(cS)) < n*default_tol(T()),
289 // "Internal error, imag part should be zero");
290  } else {
291  dense_matrix<T> auxmat(n,n);
292  logm_utri_inplace(S);
293  gmm::mult(Q, S, auxmat);
294  gmm::mult(auxmat, gmm::transposed(Q), LOGMA);
295  }
296 
297  }
298 
299 }
300 
301 #endif
302 
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::mat_norm1
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_norm1(const M &m)
*‍/
Definition: gmm_blas.h:782
gmm::sqrtm_utri_inplace
void sqrtm_utri_inplace(dense_matrix< T > &A)
Matrix square root for upper triangular matrices (from GNU Octave).
Definition: gmm_dense_matrix_functions.h:48
gmm::logm_utri_inplace
void logm_utri_inplace(dense_matrix< T > &S)
Matrix logarithm for upper triangular matrices (from GNU/Octave)
Definition: gmm_dense_matrix_functions.h:137
gmm::logm
void logm(const dense_matrix< T > &A, dense_matrix< T > &LOGMA)
Matrix logarithm (from GNU/Octave)
Definition: gmm_dense_matrix_functions.h:252