Nektar++
Static Public Member Functions | List of all members
Nektar::LinearSystemSolver Struct Reference

#include <NekLinSys.hpp>

Static Public Member Functions

template<typename BVectorType , typename XVectorType >
static void Solve (const BVectorType &b, XVectorType &x, MatrixStorage m_matrixType, const Array< OneD, const int > &m_ipivot, unsigned int n, const Array< OneD, const double > &A, char m_transposeFlag, unsigned int m_numberOfSubDiagonals, unsigned int m_numberOfSuperDiagonals)
 
template<typename BVectorType , typename XVectorType >
static void SolveTranspose (const BVectorType &b, XVectorType &x, MatrixStorage m_matrixType, const Array< OneD, const int > &m_ipivot, unsigned int n, const Array< OneD, const double > &A, char m_transposeFlag, unsigned int m_numberOfSubDiagonals, unsigned int m_numberOfSuperDiagonals)
 

Detailed Description

Definition at line 72 of file NekLinSys.hpp.

Member Function Documentation

◆ Solve()

template<typename BVectorType , typename XVectorType >
static void Nektar::LinearSystemSolver::Solve ( const BVectorType &  b,
XVectorType &  x,
MatrixStorage  m_matrixType,
const Array< OneD, const int > &  m_ipivot,
unsigned int  n,
const Array< OneD, const double > &  A,
char  m_transposeFlag,
unsigned int  m_numberOfSubDiagonals,
unsigned int  m_numberOfSuperDiagonals 
)
inlinestatic

Definition at line 75 of file NekLinSys.hpp.

81 {
82 switch (m_matrixType)
83 {
84 case eFULL:
85 {
86 x = b;
87 int info = 0;
88 Lapack::Dgetrs('N', n, 1, A.get(), n, (int *)m_ipivot.get(),
89 x.GetRawPtr(), n, info);
90 if (info < 0)
91 {
92 std::string message =
93 "ERROR: The " + std::to_string(-info) +
94 "th parameter had an illegal parameter for dgetrs";
95 ASSERTL0(false, message.c_str());
96 }
97 }
98 break;
99 case eDIAGONAL:
100 for (unsigned int i = 0; i < A.size(); ++i)
101 {
102 x[i] = b[i] * A[i];
103 }
104 break;
105 case eUPPER_TRIANGULAR:
106 {
107 x = b;
108 int info = 0;
109 Lapack::Dtptrs('U', m_transposeFlag, 'N', n, 1, A.get(),
110 x.GetRawPtr(), n, info);
111
112 if (info < 0)
113 {
114 std::string message =
115 "ERROR: The " + std::to_string(-info) +
116 "th parameter had an illegal parameter for dtptrs";
117 ASSERTL0(false, message.c_str());
118 }
119 else if (info > 0)
120 {
121 std::string message =
122 "ERROR: The " + std::to_string(-info) +
123 "th diagonal element of A is 0 for dtptrs";
124 ASSERTL0(false, message.c_str());
125 }
126 }
127 break;
128 case eLOWER_TRIANGULAR:
129 {
130 x = b;
131 int info = 0;
132 Lapack::Dtptrs('L', m_transposeFlag, 'N', n, 1, A.get(),
133 x.GetRawPtr(), n, info);
134
135 if (info < 0)
136 {
137 std::string message =
138 "ERROR: The " + std::to_string(-info) +
139 "th parameter had an illegal parameter for dtptrs";
140 ASSERTL0(false, message.c_str());
141 }
142 else if (info > 0)
143 {
144 std::string message =
145 "ERROR: The " + std::to_string(-info) +
146 "th diagonal element of A is 0 for dtptrs";
147 ASSERTL0(false, message.c_str());
148 }
149 }
150 break;
151 case eSYMMETRIC:
152 {
153 x = b;
154 int info = 0;
155 Lapack::Dsptrs('U', n, 1, A.get(), m_ipivot.get(),
156 x.GetRawPtr(), x.GetRows(), info);
157 if (info < 0)
158 {
159 std::string message =
160 "ERROR: The " + std::to_string(-info) +
161 "th parameter had an illegal parameter for dsptrs";
162 ASSERTL0(false, message.c_str());
163 }
164 }
165 break;
166 case ePOSITIVE_DEFINITE_SYMMETRIC:
167 {
168 x = b;
169 int info = 0;
170 Lapack::Dpptrs('U', n, 1, A.get(), x.GetRawPtr(), x.GetRows(),
171 info);
172 if (info < 0)
173 {
174 std::string message =
175 "ERROR: The " + std::to_string(-info) +
176 "th parameter had an illegal parameter for dpptrs";
177 ASSERTL0(false, message.c_str());
178 }
179 }
180 break;
181 case eBANDED:
182 {
183 x = b;
184 int KL = m_numberOfSubDiagonals;
185 int KU = m_numberOfSuperDiagonals;
186 int info = 0;
187
188 Lapack::Dgbtrs(m_transposeFlag, n, KL, KU, 1, A.get(),
189 2 * KL + KU + 1, m_ipivot.get(), x.GetRawPtr(),
190 n, info);
191
192 if (info < 0)
193 {
194 std::string message =
195 "ERROR: The " + std::to_string(-info) +
196 "th parameter had an illegal parameter for dgbtrs";
197 ASSERTL0(false, message.c_str());
198 }
199 }
200 break;
201 case ePOSITIVE_DEFINITE_SYMMETRIC_BANDED:
202 {
203 x = b;
204 int KU = m_numberOfSuperDiagonals;
205 int info = 0;
206
207 Lapack::Dpbtrs('U', n, KU, 1, A.get(), KU + 1, x.GetRawPtr(), n,
208 info);
209
210 if (info < 0)
211 {
212 std::string message =
213 "ERROR: The " + std::to_string(-info) +
214 "th parameter had an illegal parameter for dpbtrs";
215 ASSERTL0(false, message.c_str());
216 }
217 }
218 break;
219 case eSYMMETRIC_BANDED:
220 NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
221 break;
222 case eUPPER_TRIANGULAR_BANDED:
223 NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
224 break;
225 case eLOWER_TRIANGULAR_BANDED:
226 NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
227 break;
228
229 default:
230 NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
231 }
232 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
Lapack::Dpptrs
static void Dpptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, double *b, const int &ldb, int &info)
Solve a real positive definite symmetric matrix problem using Cholesky factorization.
Definition: Lapack.hpp:206
Lapack::Dsptrs
static void Dsptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, const int *ipiv, double *b, const int &ldb, int &info)
Solve a real packed-symmetric matrix problem using Bunch-Kaufman pivoting.
Definition: Lapack.hpp:154
Lapack::Dgetrs
static void Dgetrs(const char &trans, const int &n, const int &nrhs, const double *a, const int &lda, int *ipiv, double *b, const int &ldb, int &info)
General matrix LU backsolve.
Definition: Lapack.hpp:269
Lapack::Dgbtrs
static void Dgbtrs(const char &trans, const int &n, const int &kl, const int &ku, const int &nrhs, const double *a, const int &lda, const int *ipiv, double *b, const int &ldb, int &info)
Solve general banded matrix using LU factorisation.
Definition: Lapack.hpp:239
Lapack::Dtptrs
static void Dtptrs(const char &uplo, const char &trans, const char &diag, const int &n, const int &nrhs, const double *a, double *b, const int &ldb, int &info)
Solve a triangular system.
Definition: Lapack.hpp:184
Lapack::Dpbtrs
static void Dpbtrs(const char &uplo, const int &n, const int &kd, const int &nrhs, const double *ab, const int &ldab, double *b, const int &ldb, int &info)
Solve a real, positive definite banded-symmetric matrix problem using Cholesky factorization.
Definition: Lapack.hpp:223
Nektar::eLOWER_TRIANGULAR_BANDED
@ eLOWER_TRIANGULAR_BANDED
Definition: MatrixStorageType.h:53
Nektar::ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
Definition: MatrixStorageType.h:51
Nektar::ePOSITIVE_DEFINITE_SYMMETRIC
@ ePOSITIVE_DEFINITE_SYMMETRIC
Definition: MatrixStorageType.h:48
Nektar::eUPPER_TRIANGULAR_BANDED
@ eUPPER_TRIANGULAR_BANDED
Definition: MatrixStorageType.h:52

References ASSERTL0, Lapack::Dgbtrs(), Lapack::Dgetrs(), Lapack::Dpbtrs(), Lapack::Dpptrs(), Lapack::Dsptrs(), Lapack::Dtptrs(), Nektar::eBANDED, Nektar::eDIAGONAL, Nektar::ErrorUtil::efatal, Nektar::eFULL, Nektar::eLOWER_TRIANGULAR, Nektar::eLOWER_TRIANGULAR_BANDED, Nektar::ePOSITIVE_DEFINITE_SYMMETRIC, Nektar::ePOSITIVE_DEFINITE_SYMMETRIC_BANDED, Nektar::eSYMMETRIC, Nektar::eSYMMETRIC_BANDED, Nektar::eUPPER_TRIANGULAR, Nektar::eUPPER_TRIANGULAR_BANDED, and NEKERROR.

Referenced by Nektar::LinearSystem::Solve(), and SolveTranspose().

◆ SolveTranspose()

template<typename BVectorType , typename XVectorType >
static void Nektar::LinearSystemSolver::SolveTranspose ( const BVectorType &  b,
XVectorType &  x,
MatrixStorage  m_matrixType,
const Array< OneD, const int > &  m_ipivot,
unsigned int  n,
const Array< OneD, const double > &  A,
char  m_transposeFlag,
unsigned int  m_numberOfSubDiagonals,
unsigned int  m_numberOfSuperDiagonals 
)
inlinestatic

Definition at line 235 of file NekLinSys.hpp.

243 {
244 switch (m_matrixType)
245 {
246 case eFULL:
247 {
248 x = b;
249 int info = 0;
250 Lapack::Dgetrs('T', n, 1, A.get(), n, (int *)m_ipivot.get(),
251 x.GetRawPtr(), n, info);
252
253 if (info < 0)
254 {
255 std::string message =
256 "ERROR: The " + std::to_string(-info) +
257 "th parameter had an illegal parameter for dgetrs";
258 ASSERTL0(false, message.c_str());
259 }
260 }
261
262 break;
263 case eDIAGONAL:
264 Solve(b, x, m_matrixType, m_ipivot, n, A, m_transposeFlag,
265 m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
266 break;
267 case eUPPER_TRIANGULAR:
268 {
269 char trans = m_transposeFlag;
270 if (trans == 'N')
271 {
272 trans = 'T';
273 }
274 else
275 {
276 trans = 'N';
277 }
278
279 x = b;
280 int info = 0;
281 Lapack::Dtptrs('U', trans, 'N', n, 1, A.get(), x.GetRawPtr(), n,
282 info);
283
284 if (info < 0)
285 {
286 std::string message =
287 "ERROR: The " + std::to_string(-info) +
288 "th parameter had an illegal parameter for dtptrs";
289 ASSERTL0(false, message.c_str());
290 }
291 else if (info > 0)
292 {
293 std::string message =
294 "ERROR: The " + std::to_string(-info) +
295 "th diagonal element of A is 0 for dtptrs";
296 ASSERTL0(false, message.c_str());
297 }
298 }
299
300 break;
301 case eLOWER_TRIANGULAR:
302 {
303 char trans = m_transposeFlag;
304 if (trans == 'N')
305 {
306 trans = 'T';
307 }
308 else
309 {
310 trans = 'N';
311 }
312
313 x = b;
314 int info = 0;
315 Lapack::Dtptrs('L', trans, 'N', n, 1, A.get(), x.GetRawPtr(), n,
316 info);
317
318 if (info < 0)
319 {
320 std::string message =
321 "ERROR: The " + std::to_string(-info) +
322 "th parameter had an illegal parameter for dtptrs";
323 ASSERTL0(false, message.c_str());
324 }
325 else if (info > 0)
326 {
327 std::string message =
328 "ERROR: The " + std::to_string(-info) +
329 "th diagonal element of A is 0 for dtptrs";
330 ASSERTL0(false, message.c_str());
331 }
332 }
333 break;
334 case eSYMMETRIC:
335 case ePOSITIVE_DEFINITE_SYMMETRIC:
336 case ePOSITIVE_DEFINITE_SYMMETRIC_BANDED:
337 Solve(b, x, m_matrixType, m_ipivot, n, A, m_transposeFlag,
338 m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
339 break;
340 case eBANDED:
341 {
342 x = b;
343 int KL = m_numberOfSubDiagonals;
344 int KU = m_numberOfSuperDiagonals;
345 int info = 0;
346
347 Lapack::Dgbtrs(m_transposeFlag, n, KL, KU, 1, A.get(),
348 2 * KL + KU + 1, m_ipivot.get(), x.GetRawPtr(),
349 n, info);
350
351 if (info < 0)
352 {
353 std::string message =
354 "ERROR: The " + std::to_string(-info) +
355 "th parameter had an illegal parameter for dgbtrs";
356 ASSERTL0(false, message.c_str());
357 }
358 }
359 break;
360 case eSYMMETRIC_BANDED:
361 NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
362 break;
363 case eUPPER_TRIANGULAR_BANDED:
364 NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
365 break;
366 case eLOWER_TRIANGULAR_BANDED:
367 NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
368 break;
369
370 default:
371 NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
372 }
373 }
Nektar::LinearSystemSolver::Solve
static void Solve(const BVectorType &b, XVectorType &x, MatrixStorage m_matrixType, const Array< OneD, const int > &m_ipivot, unsigned int n, const Array< OneD, const double > &A, char m_transposeFlag, unsigned int m_numberOfSubDiagonals, unsigned int m_numberOfSuperDiagonals)
Definition: NekLinSys.hpp:75

References ASSERTL0, Lapack::Dgbtrs(), Lapack::Dgetrs(), Lapack::Dtptrs(), Nektar::eBANDED, Nektar::eDIAGONAL, Nektar::ErrorUtil::efatal, Nektar::eFULL, Nektar::eLOWER_TRIANGULAR, Nektar::eLOWER_TRIANGULAR_BANDED, Nektar::ePOSITIVE_DEFINITE_SYMMETRIC, Nektar::ePOSITIVE_DEFINITE_SYMMETRIC_BANDED, Nektar::eSYMMETRIC, Nektar::eSYMMETRIC_BANDED, Nektar::eUPPER_TRIANGULAR, Nektar::eUPPER_TRIANGULAR_BANDED, NEKERROR, and Solve().

Referenced by Nektar::LinearSystem::SolveTranspose().

Generated on Wed Jan 31 2024 09:36:00 for Nektar++ by doxygen 1.9.4