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