GetFEM  5.4.2
gmm_condition_number.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-2020 Yves Renard, Julien Pommier
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_condition_number.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>, Julien Pommier <Julien.Pommier@insa-toulouse.fr>
34  @date August 27, 2003.
35  @brief computation of the condition number of dense matrices.
36 */
37 #ifndef GMM_CONDITION_NUMBER_H__
38 #define GMM_CONDITION_NUMBER_H__
39 
40 #include "gmm_dense_qr.h"
41 
42 namespace gmm {
43 
44  /** computation of the condition number of dense matrices using SVD.
45 
46  Uses symmetric_qr_algorithm => dense matrices only.
47 
48  @param M a matrix.
49  @param emin smallest (in magnitude) eigenvalue
50  @param emax largest eigenvalue.
51  */
52  template <typename MAT>
53  typename number_traits<typename
54  linalg_traits<MAT>::value_type>::magnitude_type
55  condition_number(const MAT& M,
56  typename number_traits<typename
57  linalg_traits<MAT>::value_type>::magnitude_type& emin,
58  typename number_traits<typename
59  linalg_traits<MAT>::value_type>::magnitude_type& emax) {
60  typedef typename linalg_traits<MAT>::value_type T;
61  typedef typename number_traits<T>::magnitude_type R;
62 
63  // Added because of errors in complex with zero det
64  if (sizeof(T) != sizeof(R) && gmm::abs(gmm::lu_det(M)) == R(0))
65  return gmm::default_max(R());
66 
67  size_type m = mat_nrows(M), n = mat_ncols(M);
68  emax = emin = R(0);
69  std::vector<R> eig(m+n);
70 
71  if (m+n == 0) return R(0);
72  if (is_hermitian(M)) {
73  eig.resize(m);
74  gmm::symmetric_qr_algorithm(M, eig);
75  }
76  else {
77  dense_matrix<T> B(m+n, m+n); // not very efficient ??
78  gmm::copy(conjugated(M), sub_matrix(B, sub_interval(m, n), sub_interval(0, m)));
79  gmm::copy(M, sub_matrix(B, sub_interval(0, m),
80  sub_interval(m, n)));
81  gmm::symmetric_qr_algorithm(B, eig);
82  }
83  emin = emax = gmm::abs(eig[0]);
84  for (size_type i = 1; i < eig.size(); ++i) {
85  R e = gmm::abs(eig[i]);
86  emin = std::min(emin, e);
87  emax = std::max(emax, e);
88  }
89  // cout << "emin = " << emin << " emax = " << emax << endl;
90  if (emin == R(0)) return gmm::default_max(R());
91  return emax / emin;
92  }
93 
94  template <typename MAT>
95  typename number_traits<typename
96  linalg_traits<MAT>::value_type>::magnitude_type
97  condition_number(const MAT& M) {
98  typename number_traits<typename
99  linalg_traits<MAT>::value_type>::magnitude_type emax, emin;
100  return condition_number(M, emin, emax);
101  }
102 
103  template <typename MAT>
104  typename number_traits<typename
105  linalg_traits<MAT>::value_type>::magnitude_type
106  Frobenius_condition_number_sqr(const MAT& M) {
107  typedef typename linalg_traits<MAT>::value_type T;
108  typedef typename number_traits<T>::magnitude_type R;
109  size_type m = mat_nrows(M), n = mat_ncols(M);
110  dense_matrix<T> B(std::min(m,n), std::min(m,n));
111  if (m < n) mult(M,gmm::conjugated(M),B);
112  else mult(gmm::conjugated(M),M,B);
113  R trB = abs(mat_trace(B));
114  lu_inverse(B);
115  return trB*abs(mat_trace(B));
116  }
117 
118  template <typename MAT>
119  typename number_traits<typename
120  linalg_traits<MAT>::value_type>::magnitude_type
121  Frobenius_condition_number(const MAT& M)
122  { return sqrt(Frobenius_condition_number_sqr(M)); }
123 
124  /** estimation of the condition number (TO BE DONE...)
125  */
126  template <typename MAT>
127  typename number_traits<typename
128  linalg_traits<MAT>::value_type>::magnitude_type
129  condest(const MAT& M,
130  typename number_traits<typename
131  linalg_traits<MAT>::value_type>::magnitude_type& emin,
132  typename number_traits<typename
133  linalg_traits<MAT>::value_type>::magnitude_type& emax) {
134  return condition_number(M, emin, emax);
135  }
136 
137  template <typename MAT>
138  typename number_traits<typename
139  linalg_traits<MAT>::value_type>::magnitude_type
140  condest(const MAT& M) {
141  typename number_traits<typename
142  linalg_traits<MAT>::value_type>::magnitude_type emax, emin;
143  return condest(M, emin, emax);
144  }
145 }
146 
147 #endif
gmm::is_hermitian
bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol=magnitude_of_linalg(MAT)(-1))
*‍/
Definition: gmm_blas.h:2229
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::conjugated
conjugated_return< L >::return_type conjugated(const L &v)
return a conjugated view of the input matrix or vector.
Definition: gmm_conjugated.h:294
gmm::condest
number_traits< typename linalg_traits< MAT >::value_type >::magnitude_type condest(const MAT &M, typename number_traits< typename linalg_traits< MAT >::value_type >::magnitude_type &emin, typename number_traits< typename linalg_traits< MAT >::value_type >::magnitude_type &emax)
estimation of the condition number (TO BE DONE...)
Definition: gmm_condition_number.h:129
gmm::mat_trace
linalg_traits< M >::value_type mat_trace(const M &m)
Trace of a matrix.
Definition: gmm_blas.h:528
gmm_dense_qr.h
Dense QR factorization.
gmm::condition_number
number_traits< typename linalg_traits< MAT >::value_type >::magnitude_type condition_number(const MAT &M, typename number_traits< typename linalg_traits< MAT >::value_type >::magnitude_type &emin, typename number_traits< typename linalg_traits< MAT >::value_type >::magnitude_type &emax)
computation of the condition number of dense matrices using SVD.
Definition: gmm_condition_number.h:55