GetFEM  5.4.2
gmm_range_basis.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2009-2020 Yves Renard
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_range_basis.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>
34  @date March 10, 2009.
35  @brief Extract a basis of the range of a (large sparse) matrix from the
36  columns of this matrix.
37 */
38 #ifndef GMM_RANGE_BASIS_H
39 #define GMM_RANGE_BASIS_H
40 #include "gmm_dense_qr.h"
41 #include "gmm_dense_lu.h"
42 
43 #include "gmm_kernel.h"
44 #include "gmm_iter.h"
45 #include <set>
46 #include <list>
47 
48 
49 namespace gmm {
50 
51 
52  template <typename T, typename VECT, typename MAT1>
53  void tridiag_qr_algorithm
54  (std::vector<typename number_traits<T>::magnitude_type> diag,
55  std::vector<T> sdiag, const VECT &eigval_, const MAT1 &eigvect_,
56  bool compvect, tol_type_for_qr tol = default_tol_for_qr) {
57  VECT &eigval = const_cast<VECT &>(eigval_);
58  MAT1 &eigvect = const_cast<MAT1 &>(eigvect_);
59  typedef typename number_traits<T>::magnitude_type R;
60 
61  if (compvect) gmm::copy(identity_matrix(), eigvect);
62 
63  size_type n = diag.size(), q = 0, p, ite = 0;
64  if (n == 0) return;
65  if (n == 1) { eigval[0] = gmm::real(diag[0]); return; }
66 
67  symmetric_qr_stop_criterion(diag, sdiag, p, q, tol);
68 
69  while (q < n) {
70  sub_interval SUBI(p, n-p-q), SUBJ(0, mat_ncols(eigvect)), SUBK(p, n-p-q);
71  if (!compvect) SUBK = sub_interval(0,0);
72 
73  symmetric_Wilkinson_qr_step(sub_vector(diag, SUBI),
74  sub_vector(sdiag, SUBI),
75  sub_matrix(eigvect, SUBJ, SUBK), compvect);
76 
77  symmetric_qr_stop_criterion(diag, sdiag, p, q, tol*R(3));
78  ++ite;
79  GMM_ASSERT1(ite < n*100, "QR algorithm failed.");
80  }
81 
82  gmm::copy(diag, eigval);
83  }
84 
85  // Range basis with a restarted Lanczos method
86  template <typename Mat>
87  void range_basis_eff_Lanczos(const Mat &BB, std::set<size_type> &columns,
88  double EPS=1E-12) {
89  typedef std::set<size_type> TAB;
90  typedef typename linalg_traits<Mat>::value_type T;
91  typedef typename number_traits<T>::magnitude_type R;
92 
93  size_type nc_r = columns.size(), k;
94  col_matrix< rsvector<T> > B(mat_nrows(BB), mat_ncols(BB));
95 
96  k = 0;
97  for (TAB::iterator it = columns.begin(); it!=columns.end(); ++it, ++k){
98  gmm::copy(scaled(mat_col(BB, *it), T(1)/vect_norm2(mat_col(BB, *it))),
99  mat_col(B, *it));
100  }
101  std::vector<T> w(mat_nrows(B));
102  size_type restart = 120;
103  std::vector<T> sdiag(restart);
104  std::vector<R> eigval(restart), diag(restart);
105  dense_matrix<T> eigvect(restart, restart);
106 
107  R rho = R(-1), rho2;
108  while (nc_r) {
109 
110  std::vector<T> v(nc_r), v0(nc_r), wl(nc_r);
111  dense_matrix<T> lv(nc_r, restart);
112 
113  if (rho < R(0)) { // Estimate of the spectral radius of B^* B
114  gmm::fill_random(v);
115  for (size_type i = 0; i < 100; ++i) {
116  gmm::scale(v, T(1)/vect_norm2(v));
117  gmm::copy(v, v0);
118  k = 0; gmm::clear(w);
119  for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++k)
120  add(scaled(mat_col(B, *it), v[k]), w);
121  k = 0;
122  for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++k)
123  v[k] = vect_hp(w, mat_col(B, *it));
124  rho = gmm::abs(vect_hp(v, v0) / vect_hp(v0, v0));
125  }
126  rho *= R(2);
127  }
128 
129  // Computing vectors of the null space of de B^* B with restarted Lanczos
130  rho2 = 0;
131  gmm::fill_random(v);
132  size_type iter = 0;
133  for(;;++iter) {
134  R rho_old = rho2;
135  R beta = R(0), alpha;
136  gmm::scale(v, T(1)/vect_norm2(v));
137  size_type eff_restart = restart;
138  if (sdiag.size() != restart) {
139  sdiag.resize(restart); eigval.resize(restart); diag.resize(restart); gmm::resize(eigvect, restart, restart);
140  gmm::resize(lv, nc_r, restart);
141  }
142 
143  for (size_type i = 0; i < restart; ++i) { // Lanczos iterations
144  gmm::copy(v, mat_col(lv, i));
145  gmm::clear(w);
146  k = 0;
147  for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++k)
148  add(scaled(mat_col(B, *it), v[k]), w);
149 
150  k = 0;
151  for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++k)
152  wl[k] = v[k]*rho - vect_hp(w, mat_col(B, *it)) - beta*v0[k];
153  alpha = gmm::real(vect_hp(wl, v));
154  diag[i] = alpha;
155  gmm::add(gmm::scaled(v, -alpha), wl);
156  sdiag[i] = beta = vect_norm2(wl);
157  gmm::copy(v, v0);
158  if (beta < EPS) { eff_restart = i+1; break; }
159  gmm::copy(gmm::scaled(wl, T(1) / beta), v);
160  }
161  if (eff_restart != restart) {
162  sdiag.resize(eff_restart); eigval.resize(eff_restart); diag.resize(eff_restart);
163  gmm::resize(eigvect, eff_restart, eff_restart); gmm::resize(lv, nc_r, eff_restart);
164  }
165  tridiag_qr_algorithm(diag, sdiag, eigval, eigvect, true);
166 
167  size_type num = size_type(-1);
168  rho2 = R(0);
169  for (size_type j = 0; j < eff_restart; ++j)
170  { R nvp=gmm::abs(eigval[j]); if (nvp > rho2) { rho2=nvp; num=j; }}
171 
172  GMM_ASSERT1(num != size_type(-1), "Internal error");
173 
174  gmm::mult(lv, mat_col(eigvect, num), v);
175 
176  if (gmm::abs(rho2-rho_old) < rho_old*R(EPS)) break;
177  // if (gmm::abs(rho-rho2) < rho*R(gmm::sqrt(EPS))) break;
178  if (gmm::abs(rho-rho2) < rho*R(EPS)*R(100)) break;
179  }
180 
181  if (gmm::abs(rho-rho2) < rho*R(EPS*10.)) {
182  size_type j_max = size_type(-1), j = 0;
183  R val_max = R(0);
184  for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++j)
185  if (gmm::abs(v[j]) > val_max)
186  { val_max = gmm::abs(v[j]); j_max = *it; }
187  columns.erase(j_max); nc_r = columns.size();
188  }
189  else break;
190  }
191  }
192 
193  // Range basis with LU decomposition. Not stable from a numerical viewpoint.
194  // Complex version not verified
195  template <typename Mat>
196  void range_basis_eff_lu(const Mat &B, std::set<size_type> &columns,
197  std::vector<bool> &c_ortho, double EPS) {
198 
199  typedef std::set<size_type> TAB;
200  typedef typename linalg_traits<Mat>::value_type T;
201  typedef typename number_traits<T>::magnitude_type R;
202 
203  size_type nc_r = 0, nc_o = 0, nc = mat_ncols(B), nr = mat_nrows(B), i, j;
204 
205  for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it)
206  if (!(c_ortho[*it])) ++nc_r; else nc_o++;
207 
208  if (nc_r > 0) {
209 
210  gmm::row_matrix< gmm::rsvector<T> > Hr(nc, nc_r), Ho(nc, nc_o);
211  gmm::row_matrix< gmm::rsvector<T> > BBr(nr, nc_r), BBo(nr, nc_o);
212 
213  i = j = 0;
214  for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it)
215  if (!(c_ortho[*it]))
216  { Hr(*it, i) = T(1)/ vect_norminf(mat_col(B, *it)); ++i; }
217  else
218  { Ho(*it, j) = T(1)/ vect_norm2(mat_col(B, *it)); ++j; }
219 
220  gmm::mult(B, Hr, BBr);
221  gmm::mult(B, Ho, BBo);
222  gmm::dense_matrix<T> M(nc_r, nc_r), BBB(nc_r, nc_o), MM(nc_r, nc_r);
223  gmm::mult(gmm::conjugated(BBr), BBr, M);
224  gmm::mult(gmm::conjugated(BBr), BBo, BBB);
225  gmm::mult(BBB, gmm::conjugated(BBB), MM);
226  gmm::add(gmm::scaled(MM, T(-1)), M);
227 
228  std::vector<int> ipvt(nc_r);
229  gmm::lu_factor(M, ipvt);
230 
231  R emax = R(0);
232  for (i = 0; i < nc_r; ++i) emax = std::max(emax, gmm::abs(M(i,i)));
233 
234  i = 0;
235  std::set<size_type> c = columns;
236  for (TAB::iterator it = c.begin(); it != c.end(); ++it)
237  if (!(c_ortho[*it])) {
238  if (gmm::abs(M(i,i)) <= R(EPS)*emax) columns.erase(*it);
239  ++i;
240  }
241  }
242  }
243 
244 
245  // Range basis with Gram-Schmidt orthogonalization (sparse version)
246  // The sparse version is better when the sparsity is high and less efficient
247  // than the dense version for high degree elements (P3, P4 ...)
248  // Complex version not verified
249  template <typename Mat>
250  void range_basis_eff_Gram_Schmidt_sparse(const Mat &BB,
251  std::set<size_type> &columns,
252  std::vector<bool> &c_ortho,
253  double EPS) {
254 
255  typedef std::set<size_type> TAB;
256  typedef typename linalg_traits<Mat>::value_type T;
257  typedef typename number_traits<T>::magnitude_type R;
258 
259  size_type nc = mat_ncols(BB), nr = mat_nrows(BB);
260  std::set<size_type> c = columns, rc = columns;
261 
262  gmm::col_matrix< rsvector<T> > B(nr, nc);
263  for (std::set<size_type>::iterator it = columns.begin();
264  it != columns.end(); ++it) {
265  gmm::copy(mat_col(BB, *it), mat_col(B, *it));
266  gmm::scale(mat_col(B, *it), T(1)/vect_norm2(mat_col(B, *it)));
267  }
268 
269  for (std::set<size_type>::iterator it = c.begin(); it != c.end(); ++it)
270  if (c_ortho[*it]) {
271  for (std::set<size_type>::iterator it2 = rc.begin();
272  it2 != rc.end(); ++it2)
273  if (!(c_ortho[*it2])) {
274  T r = -vect_hp(mat_col(B, *it2), mat_col(B, *it));
275  if (r != T(0)) add(scaled(mat_col(B, *it), r), mat_col(B, *it2));
276  }
277  rc.erase(*it);
278  }
279 
280  while (rc.size()) {
281  R nmax = R(0); size_type cmax = size_type(-1);
282  for (std::set<size_type>::iterator it=rc.begin(); it != rc.end();) {
283  TAB::iterator itnext = it; ++itnext;
284  R n = vect_norm2(mat_col(B, *it));
285  if (nmax < n) { nmax = n; cmax = *it; }
286  if (n < R(EPS)) { columns.erase(*it); rc.erase(*it); }
287  it = itnext;
288  }
289 
290  if (nmax < R(EPS)) break;
291 
292  gmm::scale(mat_col(B, cmax), T(1)/vect_norm2(mat_col(B, cmax)));
293  rc.erase(cmax);
294  for (std::set<size_type>::iterator it=rc.begin(); it!=rc.end(); ++it) {
295  T r = -vect_hp(mat_col(B, *it), mat_col(B, cmax));
296  if (r != T(0)) add(scaled(mat_col(B, cmax), r), mat_col(B, *it));
297  }
298  }
299  for (std::set<size_type>::iterator it=rc.begin(); it!=rc.end(); ++it)
300  columns.erase(*it);
301  }
302 
303 
304  // Range basis with Gram-Schmidt orthogonalization (dense version)
305  template <typename Mat>
306  void range_basis_eff_Gram_Schmidt_dense(const Mat &B,
307  std::set<size_type> &columns,
308  std::vector<bool> &c_ortho,
309  double EPS) {
310 
311  typedef std::set<size_type> TAB;
312  typedef typename linalg_traits<Mat>::value_type T;
313  typedef typename number_traits<T>::magnitude_type R;
314 
315  size_type nc_r = columns.size(), nc = mat_ncols(B), nr = mat_nrows(B), i;
316  std::set<size_type> rc;
317 
318  row_matrix< gmm::rsvector<T> > H(nc, nc_r), BB(nr, nc_r);
319  std::vector<T> v(nc_r);
320  std::vector<size_type> ind(nc_r);
321 
322  i = 0;
323  for (TAB::iterator it = columns.begin(); it != columns.end(); ++it, ++i)
324  H(*it, i) = T(1) / vect_norm2(mat_col(B, *it));
325 
326  mult(B, H, BB);
327  dense_matrix<T> M(nc_r, nc_r);
328  mult(gmm::conjugated(BB), BB, M);
329 
330  i = 0;
331  for (TAB::iterator it = columns.begin(); it != columns.end(); ++it, ++i)
332  if (c_ortho[*it]) {
333  gmm::copy(mat_row(M, i), v);
334  rank_one_update(M, scaled(v, T(-1)), v);
335  M(i, i) = T(1);
336  }
337  else { rc.insert(i); ind[i] = *it; }
338 
339  while (rc.size() > 0) {
340 
341  // Next pivot
342  R nmax = R(0); size_type imax = size_type(-1);
343  for (TAB::iterator it = rc.begin(); it != rc.end();) {
344  TAB::iterator itnext = it; ++itnext;
345  R a = gmm::abs(M(*it, *it));
346  if (a > nmax) { nmax = a; imax = *it; }
347  if (a < R(EPS)) { columns.erase(ind[*it]); rc.erase(*it); }
348  it = itnext;
349  }
350 
351  if (nmax < R(EPS)) break;
352 
353  // Normalization
354  gmm::scale(mat_row(M, imax), T(1) / sqrt(nmax));
355  gmm::scale(mat_col(M, imax), T(1) / sqrt(nmax));
356 
357  // orthogonalization
358  copy(mat_row(M, imax), v);
359  rank_one_update(M, scaled(v, T(-1)), v);
360  M(imax, imax) = T(1);
361 
362  rc.erase(imax);
363  }
364  for (std::set<size_type>::iterator it=rc.begin(); it!=rc.end(); ++it)
365  columns.erase(ind[*it]);
366  }
367 
368  template <typename L> size_type nnz_eps(const L& l, double eps) {
369  typename linalg_traits<L>::const_iterator it = vect_const_begin(l),
370  ite = vect_const_end(l);
371  size_type res(0);
372  for (; it != ite; ++it) if (gmm::abs(*it) >= eps) ++res;
373  return res;
374  }
375 
376  template <typename L>
377  bool reserve__rb(const L& l, std::vector<bool> &b, double eps) {
378  typename linalg_traits<L>::const_iterator it = vect_const_begin(l),
379  ite = vect_const_end(l);
380  bool ok = true;
381  for (; it != ite; ++it)
382  if (gmm::abs(*it) >= eps && b[it.index()]) ok = false;
383  if (ok) {
384  for (it = vect_const_begin(l); it != ite; ++it)
385  if (gmm::abs(*it) >= eps) b[it.index()] = true;
386  }
387  return ok;
388  }
389 
390  template <typename Mat>
391  void range_basis(const Mat &B, std::set<size_type> &columns,
392  double EPS, col_major, bool skip_init=false) {
393 
394  typedef typename linalg_traits<Mat>::value_type T;
395  typedef typename number_traits<T>::magnitude_type R;
396 
397  size_type nc = mat_ncols(B), nr = mat_nrows(B);
398 
399  std::vector<R> norms(nc);
400  std::vector<bool> c_ortho(nc), booked(nr);
401  std::vector< std::set<size_type> > nnzs(mat_nrows(B));
402 
403  if (!skip_init) {
404 
405  R norm_max = R(0);
406  for (size_type i = 0; i < nc; ++i) {
407  norms[i] = vect_norminf(mat_col(B, i));
408  norm_max = std::max(norm_max, norms[i]);
409  }
410 
411  columns.clear();
412  for (size_type i = 0; i < nc; ++i)
413  if (norms[i] > norm_max*R(EPS)) {
414  columns.insert(i);
415  nnzs[nnz_eps(mat_col(B, i), R(EPS) * norms[i])].insert(i);
416  }
417 
418  for (size_type i = 1; i < nr; ++i)
419  for (std::set<size_type>::iterator it = nnzs[i].begin();
420  it != nnzs[i].end(); ++it)
421  if (reserve__rb(mat_col(B, *it), booked, R(EPS) * norms[*it]))
422  c_ortho[*it] = true;
423  }
424 
425  size_type sizesm[7] = {125, 200, 350, 550, 800, 1100, 1500}, actsize;
426  for (int k = 0; k < 7; ++k) {
427  size_type nc_r = columns.size();
428  std::set<size_type> c1, cres;
429  actsize = sizesm[k];
430  for (std::set<size_type>::iterator it = columns.begin();
431  it != columns.end(); ++it) {
432  c1.insert(*it);
433  if (c1.size() >= actsize) {
434  range_basis_eff_Gram_Schmidt_dense(B, c1, c_ortho, EPS);
435  for (std::set<size_type>::iterator it2=c1.begin(); it2 != c1.end();
436  ++it2) cres.insert(*it2);
437  c1.clear();
438  }
439  }
440  if (c1.size() > 1)
441  range_basis_eff_Gram_Schmidt_dense(B, c1, c_ortho, EPS);
442  for (std::set<size_type>::iterator it = c1.begin(); it != c1.end(); ++it)
443  cres.insert(*it);
444  columns = cres;
445  if (nc_r <= actsize) return;
446  if (columns.size() == nc_r) break;
447  if (sizesm[k] >= 350 && columns.size() > (nc_r*19)/20) break;
448  }
449  if (columns.size() > std::max(size_type(10), actsize))
450  range_basis_eff_Lanczos(B, columns, EPS);
451  else
452  range_basis_eff_Gram_Schmidt_dense(B, columns, c_ortho, EPS);
453  }
454 
455 
456  template <typename Mat>
457  void range_basis(const Mat &B, std::set<size_type> &columns,
458  double EPS, row_major) {
459  typedef typename linalg_traits<Mat>::value_type T;
460  gmm::col_matrix< rsvector<T> > BB(mat_nrows(B), mat_ncols(B));
461  GMM_WARNING3("A copy of a row matrix is done into a column matrix "
462  "for range basis algorithm.");
463  gmm::copy(B, BB);
464  range_basis(BB, columns, EPS);
465  }
466 
467  /** Range Basis :
468  Extract a basis of the range of a (large sparse) matrix selecting some
469  column vectors of this matrix. This is in particular useful to select
470  an independent set of linear constraints.
471 
472  The algorithm is optimized for two cases :
473  - when the (non trivial) kernel is small. An iterativ algorithm
474  based on Lanczos method is applied
475  - when the (non trivial) kernel is large and most of the dependencies
476  can be detected locally. A block Gram-Schmidt is applied first then
477  a restarted Lanczos method when the remaining kernel is greatly
478  smaller.
479  The restarted Lanczos method could be improved or replaced by a block
480  Lanczos method, a block Wiedelann method (in order to be parallelized for
481  instance) or simply could compute more than one vector of the null
482  space at each iteration.
483  The LU decomposition has been tested for local elimination but gives bad
484  results : the algorithm is unstable and do not permit to give the right
485  number of vector at the end of the process. Moreover, the number of final
486  vectors depends greatly on the number of vectors in a block of the local
487  analysis.
488  */
489  template <typename Mat>
490  void range_basis(const Mat &B, std::set<size_type> &columns,
491  double EPS=1E-12) {
492  range_basis(B, columns, EPS,
493  typename principal_orientation_type
494  <typename linalg_traits<Mat>::sub_orientation>::potype());
495 }
496 
497 }
498 
499 #endif
gmm::lu_factor
size_type lu_factor(DenseMatrix &A, lapack_ipvt &ipvt)
LU Factorization of a general (dense) matrix (real or complex).
Definition: gmm_dense_lu.h:129
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::vect_hp
strongest_value_type< V1, V2 >::value_type vect_hp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:511
gmm_dense_lu.h
LU factorizations and determinant computation for dense matrices.
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
gmm::fill_random
void fill_random(L &l, double cfill)
*‍/
Definition: gmm_blas.h:154
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
gmm_kernel.h
Include the base gmm files.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:977
gmm_dense_qr.h
Dense QR factorization.
gmm::vect_norminf
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norminf(const V &v)
Infinity norm of a vector.
Definition: gmm_blas.h:693
gmm_iter.h
Iteration object.
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1268