GetFEM  5.5
gmm_tri_solve.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-2026 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, see https://www.gnu.org/licenses/.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 
31 /**@file gmm_tri_solve.h
32  @author Yves Renard
33  @date October 13, 2002.
34  @brief Solve triangular linear system for dense matrices.
35 */
36 
37 #ifndef GMM_TRI_SOLVE_H__
38 #define GMM_TRI_SOLVE_H__
39 
40 #include "gmm_interface.h"
41 
42 namespace gmm {
43 
44  template <typename TriMatrix, typename VecX>
45  void upper_tri_solve__(const TriMatrix& T, VecX& x, size_t k,
46  col_major, abstract_sparse, bool is_unit) {
47  typename linalg_traits<TriMatrix>::value_type x_j;
48  for (int j = int(k) - 1; j >= 0; --j) {
49  typedef typename linalg_traits<TriMatrix>::const_sub_col_type COL;
50  COL c = mat_const_col(T, j);
51  typename linalg_traits<typename org_type<COL>::t>::const_iterator
52  it = vect_const_begin(c), ite = vect_const_end(c);
53  if (!is_unit) x[j] /= c[j];
54  for (x_j = x[j]; it != ite ; ++it)
55  if (int(it.index()) < j) x[it.index()] -= x_j * (*it);
56  }
57  }
58 
59  template <typename TriMatrix, typename VecX>
60  void upper_tri_solve__(const TriMatrix& T, VecX& x, size_t k,
61  col_major, abstract_dense, bool is_unit) {
62  typename linalg_traits<TriMatrix>::value_type x_j;
63  for (int j = int(k) - 1; j >= 0; --j) {
64  typedef typename linalg_traits<TriMatrix>::const_sub_col_type COL;
65  COL c = mat_const_col(T, j);
66  typename linalg_traits<typename org_type<COL>::t>::const_iterator
67  it = vect_const_begin(c), ite = it + j;
68  typename linalg_traits<VecX>::iterator itx = vect_begin(x);
69  if (!is_unit) x[j] /= c[j];
70  for (x_j = x[j]; it != ite ; ++it, ++itx) *itx -= x_j * (*it);
71  }
72  }
73 
74  template <typename TriMatrix, typename VecX>
75  void lower_tri_solve__(const TriMatrix& T, VecX& x, size_t k,
76  col_major, abstract_sparse, bool is_unit) {
77  typename linalg_traits<TriMatrix>::value_type x_j;
78  // cout << "(lower col)The Tri Matrix = " << T << endl;
79  // cout << "k = " << endl;
80  for (int j = 0; j < int(k); ++j) {
81  typedef typename linalg_traits<TriMatrix>::const_sub_col_type COL;
82  COL c = mat_const_col(T, j);
83  typename linalg_traits<typename org_type<COL>::t>::const_iterator
84  it = vect_const_begin(c), ite = vect_const_end(c);
85  if (!is_unit) x[j] /= c[j];
86  for (x_j = x[j]; it != ite ; ++it)
87  if (int(it.index()) > j && it.index() < k) x[it.index()] -= x_j*(*it);
88  }
89  }
90 
91  template <typename TriMatrix, typename VecX>
92  void lower_tri_solve__(const TriMatrix& T, VecX& x, size_t k,
93  col_major, abstract_dense, bool is_unit) {
94  typename linalg_traits<TriMatrix>::value_type x_j;
95  for (int j = 0; j < int(k); ++j) {
96  typedef typename linalg_traits<TriMatrix>::const_sub_col_type COL;
97  COL c = mat_const_col(T, j);
98  typename linalg_traits<typename org_type<COL>::t>::const_iterator
99  it = vect_const_begin(c) + (j+1), ite = vect_const_begin(c) + k;
100  typename linalg_traits<VecX>::iterator itx = vect_begin(x) + (j+1);
101  if (!is_unit) x[j] /= c[j];
102  for (x_j = x[j]; it != ite ; ++it, ++itx) *itx -= x_j * (*it);
103  }
104  }
105 
106 
107  template <typename TriMatrix, typename VecX>
108  void upper_tri_solve__(const TriMatrix& T, VecX& x, size_t k,
109  row_major, abstract_sparse, bool is_unit) {
110  typedef typename linalg_traits<TriMatrix>::const_sub_row_type ROW;
111  typename linalg_traits<TriMatrix>::value_type t;
112  typename linalg_traits<TriMatrix>::const_row_iterator
113  itr = mat_row_const_end(T);
114  for (int i = int(k) - 1; i >= 0; --i) {
115  --itr;
116  ROW c = linalg_traits<TriMatrix>::row(itr);
117  typename linalg_traits<typename org_type<ROW>::t>::const_iterator
118  it = vect_const_begin(c), ite = vect_const_end(c);
119  for (t = x[i]; it != ite; ++it)
120  if (int(it.index()) > i && it.index() < k) t -= (*it) * x[it.index()];
121  if (!is_unit) x[i] = t / c[i]; else x[i] = t;
122  }
123  }
124 
125  template <typename TriMatrix, typename VecX>
126  void upper_tri_solve__(const TriMatrix& T, VecX& x, size_t k,
127  row_major, abstract_dense, bool is_unit) {
128  typename linalg_traits<TriMatrix>::value_type t;
129 
130  for (int i = int(k) - 1; i >= 0; --i) {
131  typedef typename linalg_traits<TriMatrix>::const_sub_row_type ROW;
132  ROW c = mat_const_row(T, i);
133  typename linalg_traits<typename org_type<ROW>::t>::const_iterator
134  it = vect_const_begin(c) + (i + 1), ite = vect_const_begin(c) + k;
135  typename linalg_traits<VecX>::iterator itx = vect_begin(x) + (i+1);
136 
137  for (t = x[i]; it != ite; ++it, ++itx) t -= (*it) * (*itx);
138  if (!is_unit) x[i] = t / c[i]; else x[i] = t;
139  }
140  }
141 
142  template <typename TriMatrix, typename VecX>
143  void lower_tri_solve__(const TriMatrix& T, VecX& x, size_t k,
144  row_major, abstract_sparse, bool is_unit) {
145  typename linalg_traits<TriMatrix>::value_type t;
146 
147  for (int i = 0; i < int(k); ++i) {
148  typedef typename linalg_traits<TriMatrix>::const_sub_row_type ROW;
149  ROW c = mat_const_row(T, i);
150  typename linalg_traits<typename org_type<ROW>::t>::const_iterator
151  it = vect_const_begin(c), ite = vect_const_end(c);
152 
153  for (t = x[i]; it != ite; ++it)
154  if (int(it.index()) < i) t -= (*it) * x[it.index()];
155  if (!is_unit) x[i] = t / c[i]; else x[i] = t;
156  }
157  }
158 
159  template <typename TriMatrix, typename VecX>
160  void lower_tri_solve__(const TriMatrix& T, VecX& x, size_t k,
161  row_major, abstract_dense, bool is_unit) {
162  typename linalg_traits<TriMatrix>::value_type t;
163 
164  for (int i = 0; i < int(k); ++i) {
165  typedef typename linalg_traits<TriMatrix>::const_sub_row_type ROW;
166  ROW c = mat_const_row(T, i);
167  typename linalg_traits<typename org_type<ROW>::t>::const_iterator
168  it = vect_const_begin(c), ite = it + i;
169  typename linalg_traits<VecX>::iterator itx = vect_begin(x);
170 
171  for (t = x[i]; it != ite; ++it, ++itx) t -= (*it) * (*itx);
172  if (!is_unit) x[i] = t / c[i]; else x[i] = t;
173  }
174  }
175 
176 
177 // Triangular Solve: x <-- T^{-1} * x
178 
179  template <typename TriMatrix, typename VecX> inline
180  void upper_tri_solve(const TriMatrix& T, VecX &x_, bool is_unit = false)
181  { upper_tri_solve(T, x_, mat_nrows(T), is_unit); }
182 
183  template <typename TriMatrix, typename VecX> inline
184  void lower_tri_solve(const TriMatrix& T, VecX &x_, bool is_unit = false)
185  { lower_tri_solve(T, x_, mat_nrows(T), is_unit); }
186 
187  template <typename TriMatrix, typename VecX> inline
188  void upper_tri_solve(const TriMatrix& T, VecX &x_, size_t k,
189  bool is_unit) {
190  VecX& x = const_cast<VecX&>(x_);
191  GMM_ASSERT2(mat_nrows(T) >= k && vect_size(x) >= k
192  && mat_ncols(T) >= k && !is_sparse(x_), "dimensions mismatch");
193  upper_tri_solve__(T, x, k,
194  typename principal_orientation_type<typename
195  linalg_traits<TriMatrix>::sub_orientation>::potype(),
196  typename linalg_traits<TriMatrix>::storage_type(),
197  is_unit);
198  }
199 
200  template <typename TriMatrix, typename VecX> inline
201  void lower_tri_solve(const TriMatrix& T, VecX &x_, size_t k,
202  bool is_unit) {
203  VecX& x = const_cast<VecX&>(x_);
204  GMM_ASSERT2(mat_nrows(T) >= k && vect_size(x) >= k
205  && mat_ncols(T) >= k && !is_sparse(x_), "dimensions mismatch");
206  lower_tri_solve__(T, x, k,
207  typename principal_orientation_type<typename
208  linalg_traits<TriMatrix>::sub_orientation>::potype(),
209  typename linalg_traits<TriMatrix>::storage_type(),
210  is_unit);
211  }
212 
213 
214 
215 
216 
217 
218 }
219 
220 
221 #endif // GMM_TRI_SOLVE_H__
gmm_interface.h
gmm interface for STL vectors.