You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

499 lines
18 KiB

  1. /* -*- c++ -*- (enables emacs c++ mode) */
  2. /*===========================================================================
  3. Copyright (C) 2009-2015 Yves Renard
  4. This file is a part of GETFEM++
  5. Getfem++ is free software; you can redistribute it and/or modify it
  6. under the terms of the GNU Lesser General Public License as published
  7. by the Free Software Foundation; either version 3 of the License, or
  8. (at your option) any later version along with the GCC Runtime Library
  9. Exception either version 3.1 or (at your option) any later version.
  10. This program is distributed in the hope that it will be useful, but
  11. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12. or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
  13. License and GCC Runtime Library Exception for more details.
  14. You should have received a copy of the GNU Lesser General Public License
  15. along with this program; if not, write to the Free Software Foundation,
  16. Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
  17. As a special exception, you may use this file as it is a part of a free
  18. software library without restriction. Specifically, if other files
  19. instantiate templates or use macros or inline functions from this file,
  20. or you compile this file and link it with other files to produce an
  21. executable, this file does not by itself cause the resulting executable
  22. to be covered by the GNU Lesser General Public License. This exception
  23. does not however invalidate any other reasons why the executable file
  24. might be covered by the GNU Lesser General Public License.
  25. ===========================================================================*/
  26. /**@file gmm_range_basis.h
  27. @author Yves Renard <Yves.Renard@insa-lyon.fr>
  28. @date March 10, 2009.
  29. @brief Extract a basis of the range of a (large sparse) matrix from the
  30. columns of this matrix.
  31. */
  32. #ifndef GMM_RANGE_BASIS_H
  33. #define GMM_RANGE_BASIS_H
  34. #include "gmm_dense_qr.h"
  35. #include "gmm_dense_lu.h"
  36. #include "gmm_kernel.h"
  37. #include "gmm_iter.h"
  38. #include <set>
  39. #include <list>
  40. namespace gmm {
  41. template <typename T, typename VECT, typename MAT1>
  42. void tridiag_qr_algorithm
  43. (std::vector<typename number_traits<T>::magnitude_type> diag,
  44. std::vector<T> sdiag, const VECT &eigval_, const MAT1 &eigvect_,
  45. bool compvect, tol_type_for_qr tol = default_tol_for_qr) {
  46. VECT &eigval = const_cast<VECT &>(eigval_);
  47. MAT1 &eigvect = const_cast<MAT1 &>(eigvect_);
  48. typedef typename number_traits<T>::magnitude_type R;
  49. if (compvect) gmm::copy(identity_matrix(), eigvect);
  50. size_type n = diag.size(), q = 0, p, ite = 0;
  51. if (n == 0) return;
  52. if (n == 1) { eigval[0] = gmm::real(diag[0]); return; }
  53. symmetric_qr_stop_criterion(diag, sdiag, p, q, tol);
  54. while (q < n) {
  55. sub_interval SUBI(p, n-p-q), SUBJ(0, mat_ncols(eigvect)), SUBK(p, n-p-q);
  56. if (!compvect) SUBK = sub_interval(0,0);
  57. symmetric_Wilkinson_qr_step(sub_vector(diag, SUBI),
  58. sub_vector(sdiag, SUBI),
  59. sub_matrix(eigvect, SUBJ, SUBK), compvect);
  60. symmetric_qr_stop_criterion(diag, sdiag, p, q, tol*R(3));
  61. ++ite;
  62. GMM_ASSERT1(ite < n*100, "QR algorithm failed.");
  63. }
  64. gmm::copy(diag, eigval);
  65. }
  66. // Range basis with a restarted Lanczos method
  67. template <typename Mat>
  68. void range_basis_eff_Lanczos(const Mat &BB, std::set<size_type> &columns,
  69. double EPS=1E-12) {
  70. typedef std::set<size_type> TAB;
  71. typedef typename linalg_traits<Mat>::value_type T;
  72. typedef typename number_traits<T>::magnitude_type R;
  73. size_type nc_r = columns.size(), k;
  74. col_matrix< rsvector<T> > B(mat_nrows(BB), mat_ncols(BB));
  75. k = 0;
  76. for (TAB::iterator it = columns.begin(); it!=columns.end(); ++it, ++k){
  77. gmm::copy(scaled(mat_col(BB, *it), T(1)/vect_norm2(mat_col(BB, *it))),
  78. mat_col(B, *it));
  79. }
  80. std::vector<T> w(mat_nrows(B));
  81. size_type restart = 120;
  82. std::vector<T> sdiag(restart);
  83. std::vector<R> eigval(restart), diag(restart);
  84. dense_matrix<T> eigvect(restart, restart);
  85. R rho = R(-1), rho2;
  86. while (nc_r) {
  87. std::vector<T> v(nc_r), v0(nc_r), wl(nc_r);
  88. dense_matrix<T> lv(nc_r, restart);
  89. if (rho < R(0)) { // Estimate of the spectral radius of B^* B
  90. gmm::fill_random(v);
  91. for (size_type i = 0; i < 100; ++i) {
  92. gmm::scale(v, T(1)/vect_norm2(v));
  93. gmm::copy(v, v0);
  94. k = 0; gmm::clear(w);
  95. for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++k)
  96. add(scaled(mat_col(B, *it), v[k]), w);
  97. k = 0;
  98. for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++k)
  99. v[k] = vect_hp(w, mat_col(B, *it));
  100. rho = gmm::abs(vect_hp(v, v0) / vect_hp(v0, v0));
  101. }
  102. rho *= R(2);
  103. }
  104. // Computing vectors of the null space of de B^* B with restarted Lanczos
  105. rho2 = 0;
  106. gmm::fill_random(v);
  107. size_type iter = 0;
  108. for(;;++iter) {
  109. R rho_old = rho2;
  110. R beta = R(0), alpha;
  111. gmm::scale(v, T(1)/vect_norm2(v));
  112. size_type eff_restart = restart;
  113. if (sdiag.size() != restart) {
  114. sdiag.resize(restart); eigval.resize(restart); diag.resize(restart); gmm::resize(eigvect, restart, restart);
  115. gmm::resize(lv, nc_r, restart);
  116. }
  117. for (size_type i = 0; i < restart; ++i) { // Lanczos iterations
  118. gmm::copy(v, mat_col(lv, i));
  119. gmm::clear(w);
  120. k = 0;
  121. for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++k)
  122. add(scaled(mat_col(B, *it), v[k]), w);
  123. k = 0;
  124. for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++k)
  125. wl[k] = v[k]*rho - vect_hp(w, mat_col(B, *it)) - beta*v0[k];
  126. alpha = gmm::real(vect_hp(wl, v));
  127. diag[i] = alpha;
  128. gmm::add(gmm::scaled(v, -alpha), wl);
  129. sdiag[i] = beta = vect_norm2(wl);
  130. gmm::copy(v, v0);
  131. if (beta < EPS) { eff_restart = i+1; break; }
  132. gmm::copy(gmm::scaled(wl, T(1) / beta), v);
  133. }
  134. if (eff_restart != restart) {
  135. sdiag.resize(eff_restart); eigval.resize(eff_restart); diag.resize(eff_restart);
  136. gmm::resize(eigvect, eff_restart, eff_restart); gmm::resize(lv, nc_r, eff_restart);
  137. }
  138. tridiag_qr_algorithm(diag, sdiag, eigval, eigvect, true);
  139. size_type num = size_type(-1);
  140. rho2 = R(0);
  141. for (size_type j = 0; j < eff_restart; ++j)
  142. { R nvp=gmm::abs(eigval[j]); if (nvp > rho2) { rho2=nvp; num=j; }}
  143. GMM_ASSERT1(num != size_type(-1), "Internal error");
  144. gmm::mult(lv, mat_col(eigvect, num), v);
  145. if (gmm::abs(rho2-rho_old) < rho_old*R(EPS)) break;
  146. // if (gmm::abs(rho-rho2) < rho*R(gmm::sqrt(EPS))) break;
  147. if (gmm::abs(rho-rho2) < rho*R(EPS)*R(100)) break;
  148. }
  149. if (gmm::abs(rho-rho2) < rho*R(EPS*10.)) {
  150. size_type j_max = size_type(-1), j = 0;
  151. R val_max = R(0);
  152. for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it, ++j)
  153. if (gmm::abs(v[j]) > val_max)
  154. { val_max = gmm::abs(v[j]); j_max = *it; }
  155. columns.erase(j_max); nc_r = columns.size();
  156. }
  157. else break;
  158. }
  159. }
  160. // Range basis with LU decomposition. Not stable from a numerical viewpoint.
  161. // Complex version not verified
  162. template <typename Mat>
  163. void range_basis_eff_lu(const Mat &B, std::set<size_type> &columns,
  164. std::vector<bool> &c_ortho, double EPS) {
  165. typedef std::set<size_type> TAB;
  166. typedef typename linalg_traits<Mat>::value_type T;
  167. typedef typename number_traits<T>::magnitude_type R;
  168. size_type nc_r = 0, nc_o = 0, nc = mat_ncols(B), nr = mat_nrows(B), i, j;
  169. for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it)
  170. if (!(c_ortho[*it])) ++nc_r; else nc_o++;
  171. if (nc_r > 0) {
  172. gmm::row_matrix< gmm::rsvector<T> > Hr(nc, nc_r), Ho(nc, nc_o);
  173. gmm::row_matrix< gmm::rsvector<T> > BBr(nr, nc_r), BBo(nr, nc_o);
  174. i = j = 0;
  175. for (TAB::iterator it=columns.begin(); it!=columns.end(); ++it)
  176. if (!(c_ortho[*it]))
  177. { Hr(*it, i) = T(1)/ vect_norminf(mat_col(B, *it)); ++i; }
  178. else
  179. { Ho(*it, j) = T(1)/ vect_norm2(mat_col(B, *it)); ++j; }
  180. gmm::mult(B, Hr, BBr);
  181. gmm::mult(B, Ho, BBo);
  182. gmm::dense_matrix<T> M(nc_r, nc_r), BBB(nc_r, nc_o), MM(nc_r, nc_r);
  183. gmm::mult(gmm::conjugated(BBr), BBr, M);
  184. gmm::mult(gmm::conjugated(BBr), BBo, BBB);
  185. gmm::mult(BBB, gmm::conjugated(BBB), MM);
  186. gmm::add(gmm::scaled(MM, T(-1)), M);
  187. std::vector<int> ipvt(nc_r);
  188. gmm::lu_factor(M, ipvt);
  189. R emax = R(0);
  190. for (i = 0; i < nc_r; ++i) emax = std::max(emax, gmm::abs(M(i,i)));
  191. i = 0;
  192. std::set<size_type> c = columns;
  193. for (TAB::iterator it = c.begin(); it != c.end(); ++it)
  194. if (!(c_ortho[*it])) {
  195. if (gmm::abs(M(i,i)) <= R(EPS)*emax) columns.erase(*it);
  196. ++i;
  197. }
  198. }
  199. }
  200. // Range basis with Gram-Schmidt orthogonalization (sparse version)
  201. // The sparse version is better when the sparsity is high and less efficient
  202. // than the dense version for high degree elements (P3, P4 ...)
  203. // Complex version not verified
  204. template <typename Mat>
  205. void range_basis_eff_Gram_Schmidt_sparse(const Mat &BB,
  206. std::set<size_type> &columns,
  207. std::vector<bool> &c_ortho,
  208. double EPS) {
  209. typedef std::set<size_type> TAB;
  210. typedef typename linalg_traits<Mat>::value_type T;
  211. typedef typename number_traits<T>::magnitude_type R;
  212. size_type nc = mat_ncols(BB), nr = mat_nrows(BB);
  213. std::set<size_type> c = columns, rc = columns;
  214. gmm::col_matrix< rsvector<T> > B(nr, nc);
  215. for (std::set<size_type>::iterator it = columns.begin();
  216. it != columns.end(); ++it) {
  217. gmm::copy(mat_col(BB, *it), mat_col(B, *it));
  218. gmm::scale(mat_col(B, *it), T(1)/vect_norm2(mat_col(B, *it)));
  219. }
  220. for (std::set<size_type>::iterator it = c.begin(); it != c.end(); ++it)
  221. if (c_ortho[*it]) {
  222. for (std::set<size_type>::iterator it2 = rc.begin();
  223. it2 != rc.end(); ++it2)
  224. if (!(c_ortho[*it2])) {
  225. T r = -vect_hp(mat_col(B, *it2), mat_col(B, *it));
  226. if (r != T(0)) add(scaled(mat_col(B, *it), r), mat_col(B, *it2));
  227. }
  228. rc.erase(*it);
  229. }
  230. while (rc.size()) {
  231. R nmax = R(0); size_type cmax = size_type(-1);
  232. for (std::set<size_type>::iterator it=rc.begin(); it != rc.end();) {
  233. TAB::iterator itnext = it; ++itnext;
  234. R n = vect_norm2(mat_col(B, *it));
  235. if (nmax < n) { nmax = n; cmax = *it; }
  236. if (n < R(EPS)) { columns.erase(*it); rc.erase(*it); }
  237. it = itnext;
  238. }
  239. if (nmax < R(EPS)) break;
  240. gmm::scale(mat_col(B, cmax), T(1)/vect_norm2(mat_col(B, cmax)));
  241. rc.erase(cmax);
  242. for (std::set<size_type>::iterator it=rc.begin(); it!=rc.end(); ++it) {
  243. T r = -vect_hp(mat_col(B, *it), mat_col(B, cmax));
  244. if (r != T(0)) add(scaled(mat_col(B, cmax), r), mat_col(B, *it));
  245. }
  246. }
  247. for (std::set<size_type>::iterator it=rc.begin(); it!=rc.end(); ++it)
  248. columns.erase(*it);
  249. }
  250. // Range basis with Gram-Schmidt orthogonalization (dense version)
  251. template <typename Mat>
  252. void range_basis_eff_Gram_Schmidt_dense(const Mat &B,
  253. std::set<size_type> &columns,
  254. std::vector<bool> &c_ortho,
  255. double EPS) {
  256. typedef std::set<size_type> TAB;
  257. typedef typename linalg_traits<Mat>::value_type T;
  258. typedef typename number_traits<T>::magnitude_type R;
  259. size_type nc_r = columns.size(), nc = mat_ncols(B), nr = mat_nrows(B), i;
  260. std::set<size_type> rc;
  261. row_matrix< gmm::rsvector<T> > H(nc, nc_r), BB(nr, nc_r);
  262. std::vector<T> v(nc_r);
  263. std::vector<size_type> ind(nc_r);
  264. i = 0;
  265. for (TAB::iterator it = columns.begin(); it != columns.end(); ++it, ++i)
  266. H(*it, i) = T(1) / vect_norm2(mat_col(B, *it));
  267. mult(B, H, BB);
  268. dense_matrix<T> M(nc_r, nc_r);
  269. mult(gmm::conjugated(BB), BB, M);
  270. i = 0;
  271. for (TAB::iterator it = columns.begin(); it != columns.end(); ++it, ++i)
  272. if (c_ortho[*it]) {
  273. gmm::copy(mat_row(M, i), v);
  274. rank_one_update(M, scaled(v, T(-1)), v);
  275. M(i, i) = T(1);
  276. }
  277. else { rc.insert(i); ind[i] = *it; }
  278. while (rc.size() > 0) {
  279. // Next pivot
  280. R nmax = R(0); size_type imax = size_type(-1);
  281. for (TAB::iterator it = rc.begin(); it != rc.end();) {
  282. TAB::iterator itnext = it; ++itnext;
  283. R a = gmm::abs(M(*it, *it));
  284. if (a > nmax) { nmax = a; imax = *it; }
  285. if (a < R(EPS)) { columns.erase(ind[*it]); rc.erase(*it); }
  286. it = itnext;
  287. }
  288. if (nmax < R(EPS)) break;
  289. // Normalization
  290. gmm::scale(mat_row(M, imax), T(1) / sqrt(nmax));
  291. gmm::scale(mat_col(M, imax), T(1) / sqrt(nmax));
  292. // orthogonalization
  293. copy(mat_row(M, imax), v);
  294. rank_one_update(M, scaled(v, T(-1)), v);
  295. M(imax, imax) = T(1);
  296. rc.erase(imax);
  297. }
  298. for (std::set<size_type>::iterator it=rc.begin(); it!=rc.end(); ++it)
  299. columns.erase(ind[*it]);
  300. }
  301. template <typename L> size_type nnz_eps(const L& l, double eps) {
  302. typename linalg_traits<L>::const_iterator it = vect_const_begin(l),
  303. ite = vect_const_end(l);
  304. size_type res(0);
  305. for (; it != ite; ++it) if (gmm::abs(*it) >= eps) ++res;
  306. return res;
  307. }
  308. template <typename L>
  309. bool reserve__rb(const L& l, std::vector<bool> &b, double eps) {
  310. typename linalg_traits<L>::const_iterator it = vect_const_begin(l),
  311. ite = vect_const_end(l);
  312. bool ok = true;
  313. for (; it != ite; ++it)
  314. if (gmm::abs(*it) >= eps && b[it.index()]) ok = false;
  315. if (ok) {
  316. for (it = vect_const_begin(l); it != ite; ++it)
  317. if (gmm::abs(*it) >= eps) b[it.index()] = true;
  318. }
  319. return ok;
  320. }
  321. template <typename Mat>
  322. void range_basis(const Mat &B, std::set<size_type> &columns,
  323. double EPS, col_major, bool skip_init=false) {
  324. typedef typename linalg_traits<Mat>::value_type T;
  325. typedef typename number_traits<T>::magnitude_type R;
  326. size_type nc = mat_ncols(B), nr = mat_nrows(B);
  327. std::vector<R> norms(nc);
  328. std::vector<bool> c_ortho(nc), booked(nr);
  329. std::vector< std::set<size_type> > nnzs(mat_nrows(B));
  330. if (!skip_init) {
  331. R norm_max = R(0);
  332. for (size_type i = 0; i < nc; ++i) {
  333. norms[i] = vect_norminf(mat_col(B, i));
  334. norm_max = std::max(norm_max, norms[i]);
  335. }
  336. columns.clear();
  337. for (size_type i = 0; i < nc; ++i)
  338. if (norms[i] > norm_max*R(EPS)) {
  339. columns.insert(i);
  340. nnzs[nnz_eps(mat_col(B, i), R(EPS) * norms[i])].insert(i);
  341. }
  342. for (size_type i = 1; i < nr; ++i)
  343. for (std::set<size_type>::iterator it = nnzs[i].begin();
  344. it != nnzs[i].end(); ++it)
  345. if (reserve__rb(mat_col(B, *it), booked, R(EPS) * norms[*it]))
  346. c_ortho[*it] = true;
  347. }
  348. size_type sizesm[7] = {125, 200, 350, 550, 800, 1100, 1500}, actsize;
  349. for (int k = 0; k < 7; ++k) {
  350. size_type nc_r = columns.size();
  351. std::set<size_type> c1, cres;
  352. actsize = sizesm[k];
  353. for (std::set<size_type>::iterator it = columns.begin();
  354. it != columns.end(); ++it) {
  355. c1.insert(*it);
  356. if (c1.size() >= actsize) {
  357. range_basis_eff_Gram_Schmidt_dense(B, c1, c_ortho, EPS);
  358. for (std::set<size_type>::iterator it2=c1.begin(); it2 != c1.end();
  359. ++it2) cres.insert(*it2);
  360. c1.clear();
  361. }
  362. }
  363. if (c1.size() > 1)
  364. range_basis_eff_Gram_Schmidt_dense(B, c1, c_ortho, EPS);
  365. for (std::set<size_type>::iterator it = c1.begin(); it != c1.end(); ++it)
  366. cres.insert(*it);
  367. columns = cres;
  368. if (nc_r <= actsize) return;
  369. if (columns.size() == nc_r) break;
  370. if (sizesm[k] >= 350 && columns.size() > (nc_r*19)/20) break;
  371. }
  372. if (columns.size() > std::max(size_type(10), actsize))
  373. range_basis_eff_Lanczos(B, columns, EPS);
  374. else
  375. range_basis_eff_Gram_Schmidt_dense(B, columns, c_ortho, EPS);
  376. }
  377. template <typename Mat>
  378. void range_basis(const Mat &B, std::set<size_type> &columns,
  379. double EPS, row_major) {
  380. typedef typename linalg_traits<Mat>::value_type T;
  381. gmm::col_matrix< rsvector<T> > BB(mat_nrows(B), mat_ncols(B));
  382. GMM_WARNING3("A copy of a row matrix is done into a column matrix "
  383. "for range basis algorithm.");
  384. gmm::copy(B, BB);
  385. range_basis(BB, columns, EPS);
  386. }
  387. /** Range Basis :
  388. Extract a basis of the range of a (large sparse) matrix selecting some
  389. column vectors of this matrix. This is in particular useful to select
  390. an independent set of linear constraints.
  391. The algorithm is optimized for two cases :
  392. - when the (non trivial) kernel is small. An iterativ algorithm
  393. based on Lanczos method is applied
  394. - when the (non trivial) kernel is large and most of the dependencies
  395. can be detected locally. A block Gram-Schmidt is applied first then
  396. a restarted Lanczos method when the remaining kernel is greatly
  397. smaller.
  398. The restarted Lanczos method could be improved or replaced by a block
  399. Lanczos method, a block Wiedelann method (in order to be parallelized for
  400. instance) or simply could compute more than one vector of the null
  401. space at each iteration.
  402. The LU decomposition has been tested for local elimination but gives bad
  403. results : the algorithm is unstable and do not permit to give the right
  404. number of vector at the end of the process. Moreover, the number of final
  405. vectors depends greatly on the number of vectors in a block of the local
  406. analysis.
  407. */
  408. template <typename Mat>
  409. void range_basis(const Mat &B, std::set<size_type> &columns,
  410. double EPS=1E-12) {
  411. range_basis(B, columns, EPS,
  412. typename principal_orientation_type
  413. <typename linalg_traits<Mat>::sub_orientation>::potype());
  414. }
  415. }
  416. #endif