Nektar++
RiemannSolver.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: RiemannSolver.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Abstract base class for Riemann solvers with factory.
32//
33///////////////////////////////////////////////////////////////////////////////
34
37
38#define EPSILON 0.000001
39
40#define CROSS(dest, v1, v2) \
41 { \
42 dest[0] = v1[1] * v2[2] - v1[2] * v2[1]; \
43 dest[1] = v1[2] * v2[0] - v1[0] * v2[2]; \
44 dest[2] = v1[0] * v2[1] - v1[1] * v2[0]; \
45 }
46
47#define DOT(v1, v2) (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2])
48
49#define SUB(dest, v1, v2) \
50 { \
51 dest[0] = v1[0] - v2[0]; \
52 dest[1] = v1[1] - v2[1]; \
53 dest[2] = v1[2] - v2[2]; \
54 }
55
56namespace Nektar::SolverUtils
57{
58/**
59 * Retrieves the singleton instance of the Riemann solver factory.
60 */
62{
63 static RiemannSolverFactory instance;
64 return instance;
65}
66
67/**
68 * @class RiemannSolver
69 *
70 * @brief The RiemannSolver class provides an abstract interface under
71 * which solvers for various Riemann problems can be implemented.
72 */
73
74RiemannSolver::RiemannSolver() : m_requiresRotation(false), m_rotStorage(3)
75{
76}
77
79 [[maybe_unused]] const LibUtilities::SessionReaderSharedPtr &pSession)
80 : m_requiresRotation(false), m_rotStorage(3)
81{
82}
83
84/**
85 * @brief Perform the Riemann solve given the forwards and backwards
86 * spaces.
87 *
88 * This routine calls the virtual function #v_Solve to perform the
89 * Riemann solve. If the flag #m_requiresRotation is set, then the
90 * velocity field is rotated to the normal direction to perform
91 * dimensional splitting, and the resulting fluxes are rotated back to
92 * the Cartesian directions before being returned. For the Rotation to
93 * work, the normal vectors "N" and the location of the vector
94 * components in Fwd "vecLocs"must be set via the SetAuxVec() method.
95 *
96 * @param Fwd Forwards trace space.
97 * @param Bwd Backwards trace space.
98 * @param flux Resultant flux along trace space.
99 */
100void RiemannSolver::Solve(const int nDim,
101 const Array<OneD, const Array<OneD, NekDouble>> &Fwd,
102 const Array<OneD, const Array<OneD, NekDouble>> &Bwd,
104{
106 {
107 ASSERTL1(CheckVectors("N"), "N not defined.");
108 ASSERTL1(CheckAuxVec("vecLocs"), "vecLocs not defined.");
110 m_vectors["N"]();
112 m_auxVec["vecLocs"]();
113
114 int nFields = Fwd.size();
115 int nPts = Fwd[0].size();
116
117 if (m_rotStorage[0].size() != nFields ||
118 m_rotStorage[0][0].size() != nPts)
119 {
120 for (int i = 0; i < 3; ++i)
121 {
123 for (int j = 0; j < nFields; ++j)
124 {
126 }
127 }
128 }
129
130 rotateToNormal(Fwd, normals, vecLocs, m_rotStorage[0]);
131 rotateToNormal(Bwd, normals, vecLocs, m_rotStorage[1]);
132
134 rotateFromNormal(m_rotStorage[2], normals, vecLocs, flux);
135
136 // Attempt to subtract (\vec{U}\vec{vg})\dot n for ALE
137 if (m_ALESolver)
138 {
139 auto vgt = m_vectors["vgt"]();
140 auto N = m_vectors["N"]();
141 for (int i = 0; i < flux.size(); ++i)
142 {
143 for (int j = 0; j < flux[i].size(); ++j)
144 {
145 NekDouble tmp = 0;
146 for (int k = 0; k < nDim; ++k)
147 {
148 tmp += N[k][j] * vgt[k][j];
149 }
150 flux[i][j] -= 0.5 * (Fwd[i][j] + Bwd[i][j]) * tmp;
151 }
152 }
153 }
154 }
155 else
156 {
157 v_Solve(nDim, Fwd, Bwd, flux);
158 }
159}
160
161/**
162 * @brief Rotate a vector field to trace normal.
163 *
164 * This function performs a rotation of a vector so that the first
165 * component aligns with the trace normal direction.
166 *
167 * The vectors components are stored in inarray. Their locations must
168 * be specified in the "vecLocs" array. vecLocs[0] contains the locations
169 * of the first vectors components, vecLocs[1] those of the second and
170 * so on.
171 *
172 * In 2D, this is accomplished through the transform:
173 *
174 * \f[ (u_x, u_y) = (n_x u_x + n_y u_y, -n_x v_x + n_y v_y) \f]
175 *
176 * In 3D, we generate a (non-unique) transformation using
177 * RiemannSolver::fromToRotation.
178 *
179 */
181 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
182 const Array<OneD, const Array<OneD, NekDouble>> &normals,
183 const Array<OneD, const Array<OneD, NekDouble>> &vecLocs,
185{
186 for (int i = 0; i < inarray.size(); ++i)
187 {
188 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
189 }
190
191 for (int i = 0; i < vecLocs.size(); i++)
192 {
193 ASSERTL1(vecLocs[i].size() == normals.size(),
194 "vecLocs[i] element count mismatch");
195
196 switch (normals.size())
197 {
198 case 1:
199 { // do nothing
200 const int nq = inarray[0].size();
201 const int vx = (int)vecLocs[i][0];
202 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
203 break;
204 }
205 case 2:
206 {
207 const int nq = inarray[0].size();
208 const int vx = (int)vecLocs[i][0];
209 const int vy = (int)vecLocs[i][1];
210
211 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
212 Vmath::Vvtvp(nq, inarray[vy], 1, normals[1], 1, outarray[vx], 1,
213 outarray[vx], 1);
214 Vmath::Vmul(nq, inarray[vx], 1, normals[1], 1, outarray[vy], 1);
215 Vmath::Vvtvm(nq, inarray[vy], 1, normals[0], 1, outarray[vy], 1,
216 outarray[vy], 1);
217 break;
218 }
219
220 case 3:
221 {
222 const int nq = inarray[0].size();
223 const int vx = (int)vecLocs[i][0];
224 const int vy = (int)vecLocs[i][1];
225 const int vz = (int)vecLocs[i][2];
226
227 // Generate matrices if they don't already exist.
228 if (m_rotMat.size() == 0)
229 {
231 }
232
233 // Apply rotation matrices.
234 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[0], 1, inarray[vy],
235 1, m_rotMat[1], 1, outarray[vx], 1);
236 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[2], 1, outarray[vx],
237 1, outarray[vx], 1);
238 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[3], 1, inarray[vy],
239 1, m_rotMat[4], 1, outarray[vy], 1);
240 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[5], 1, outarray[vy],
241 1, outarray[vy], 1);
242 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[6], 1, inarray[vy],
243 1, m_rotMat[7], 1, outarray[vz], 1);
244 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[8], 1, outarray[vz],
245 1, outarray[vz], 1);
246 break;
247 }
248
249 default:
250 ASSERTL1(false, "Invalid space dimension.");
251 break;
252 }
253 }
254}
255
256/**
257 * @brief Rotate a vector field from trace normal.
258 *
259 * This function performs a rotation of the triad of vector components
260 * provided in inarray so that the first component aligns with the
261 * Cartesian components; it performs the inverse operation of
262 * RiemannSolver::rotateToNormal.
263 */
265 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
266 const Array<OneD, const Array<OneD, NekDouble>> &normals,
267 const Array<OneD, const Array<OneD, NekDouble>> &vecLocs,
269{
270 for (int i = 0; i < inarray.size(); ++i)
271 {
272 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
273 }
274
275 for (int i = 0; i < vecLocs.size(); i++)
276 {
277 ASSERTL1(vecLocs[i].size() == normals.size(),
278 "vecLocs[i] element count mismatch");
279
280 switch (normals.size())
281 {
282 case 1:
283 { // do nothing
284 const int nq = normals[0].size();
285 const int vx = (int)vecLocs[i][0];
286 Vmath::Vmul(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1);
287 break;
288 }
289 case 2:
290 {
291 const int nq = normals[0].size();
292 const int vx = (int)vecLocs[i][0];
293 const int vy = (int)vecLocs[i][1];
294
295 Vmath::Vmul(nq, inarray[vy], 1, normals[1], 1, outarray[vx], 1);
296 Vmath::Vvtvm(nq, inarray[vx], 1, normals[0], 1, outarray[vx], 1,
297 outarray[vx], 1);
298 Vmath::Vmul(nq, inarray[vx], 1, normals[1], 1, outarray[vy], 1);
299 Vmath::Vvtvp(nq, inarray[vy], 1, normals[0], 1, outarray[vy], 1,
300 outarray[vy], 1);
301 break;
302 }
303
304 case 3:
305 {
306 const int nq = normals[0].size();
307 const int vx = (int)vecLocs[i][0];
308 const int vy = (int)vecLocs[i][1];
309 const int vz = (int)vecLocs[i][2];
310
311 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[0], 1, inarray[vy],
312 1, m_rotMat[3], 1, outarray[vx], 1);
313 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[6], 1, outarray[vx],
314 1, outarray[vx], 1);
315 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[1], 1, inarray[vy],
316 1, m_rotMat[4], 1, outarray[vy], 1);
317 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[7], 1, outarray[vy],
318 1, outarray[vy], 1);
319 Vmath::Vvtvvtp(nq, inarray[vx], 1, m_rotMat[2], 1, inarray[vy],
320 1, m_rotMat[5], 1, outarray[vz], 1);
321 Vmath::Vvtvp(nq, inarray[vz], 1, m_rotMat[8], 1, outarray[vz],
322 1, outarray[vz], 1);
323 break;
324 }
325
326 default:
327 ASSERTL1(false, "Invalid space dimension.");
328 break;
329 }
330 }
331}
332
333/**
334 * @brief Determine whether a scalar has been defined in #m_scalars.
335 *
336 * @param name Scalar name.
337 */
339{
340 return m_scalars.find(name) != m_scalars.end();
341}
342
343/**
344 * @brief Determine whether a vector has been defined in #m_vectors.
345 *
346 * @param name Vector name.
347 */
349{
350 return m_vectors.find(name) != m_vectors.end();
351}
352
353/**
354 * @brief Determine whether a parameter has been defined in #m_params.
355 *
356 * @param name Parameter name.
357 */
359{
360 return m_params.find(name) != m_params.end();
361}
362
363/**
364 * @brief Determine whether a scalar has been defined in #m_auxScal.
365 *
366 * @param name Scalar name.
367 */
369{
370 return m_auxScal.find(name) != m_auxScal.end();
371}
372
373/**
374 * @brief Determine whether a vector has been defined in #m_auxVec.
375 *
376 * @param name Vector name.
377 */
379{
380 return m_auxVec.find(name) != m_auxVec.end();
381}
382
383/**
384 * @brief Generate rotation matrices for 3D expansions.
385 */
387 const Array<OneD, const Array<OneD, NekDouble>> &normals)
388{
389 Array<OneD, NekDouble> xdir(3, 0.0);
391 NekDouble tmp[9];
392 const int nq = normals[0].size();
393 int i, j;
394 xdir[0] = 1.0;
395
396 // Allocate storage for rotation matrices.
398
399 for (i = 0; i < 9; ++i)
400 {
402 }
403 for (i = 0; i < normals[0].size(); ++i)
404 {
405 // Generate matrix which takes us from (1,0,0) vector to trace
406 // normal.
407 tn[0] = normals[0][i];
408 tn[1] = normals[1][i];
409 tn[2] = normals[2][i];
410 FromToRotation(tn, xdir, tmp);
411
412 for (j = 0; j < 9; ++j)
413 {
414 m_rotMat[j][i] = tmp[j];
415 }
416 }
417}
418
419/**
420 * @brief A function for creating a rotation matrix that rotates a
421 * vector @a from into another vector @a to.
422 *
423 * Authors: Tomas Möller, John Hughes
424 * "Efficiently Building a Matrix to Rotate One Vector to
425 * Another" Journal of Graphics Tools, 4(4):1-4, 1999
426 *
427 * @param from Normalised 3-vector to rotate from.
428 * @param to Normalised 3-vector to rotate to.
429 * @param out Resulting 3x3 rotation matrix (row-major order).
430 */
433 NekDouble *mat)
434{
435 NekDouble v[3];
436 NekDouble e, h, f;
437
438 CROSS(v, from, to);
439 e = DOT(from, to);
440 f = (e < 0) ? -e : e;
441 if (f > 1.0 - EPSILON)
442 {
443 NekDouble u[3], v[3];
444 NekDouble x[3];
445 NekDouble c1, c2, c3;
446 int i, j;
447
448 x[0] = (from[0] > 0.0) ? from[0] : -from[0];
449 x[1] = (from[1] > 0.0) ? from[1] : -from[1];
450 x[2] = (from[2] > 0.0) ? from[2] : -from[2];
451
452 if (x[0] < x[1])
453 {
454 if (x[0] < x[2])
455 {
456 x[0] = 1.0;
457 x[1] = x[2] = 0.0;
458 }
459 else
460 {
461 x[2] = 1.0;
462 x[0] = x[1] = 0.0;
463 }
464 }
465 else
466 {
467 if (x[1] < x[2])
468 {
469 x[1] = 1.0;
470 x[0] = x[2] = 0.0;
471 }
472 else
473 {
474 x[2] = 1.0;
475 x[0] = x[1] = 0.0;
476 }
477 }
478
479 u[0] = x[0] - from[0];
480 u[1] = x[1] - from[1];
481 u[2] = x[2] - from[2];
482 v[0] = x[0] - to[0];
483 v[1] = x[1] - to[1];
484 v[2] = x[2] - to[2];
485
486 c1 = 2.0 / DOT(u, u);
487 c2 = 2.0 / DOT(v, v);
488 c3 = c1 * c2 * DOT(u, v);
489
490 for (i = 0; i < 3; i++)
491 {
492 for (j = 0; j < 3; j++)
493 {
494 mat[3 * i + j] =
495 -c1 * u[i] * u[j] - c2 * v[i] * v[j] + c3 * v[i] * u[j];
496 }
497 mat[i + 3 * i] += 1.0;
498 }
499 }
500 else
501 {
502 NekDouble hvx, hvz, hvxy, hvxz, hvyz;
503 h = 1.0 / (1.0 + e);
504 hvx = h * v[0];
505 hvz = h * v[2];
506 hvxy = hvx * v[1];
507 hvxz = hvx * v[2];
508 hvyz = hvz * v[1];
509 mat[0] = e + hvx * v[0];
510 mat[1] = hvxy - v[2];
511 mat[2] = hvxz + v[1];
512 mat[3] = hvxy + v[2];
513 mat[4] = e + h * v[1] * v[1];
514 mat[5] = hvyz - v[0];
515 mat[6] = hvxz - v[1];
516 mat[7] = hvyz + v[0];
517 mat[8] = e + hvz * v[2];
518 }
519}
520
521/**
522 * @brief Calculate the flux jacobian of Fwd and Bwd
523 *
524 * @param Fwd Forwards trace space.
525 * @param Bwd Backwards trace space.
526 * @param flux Resultant flux along trace space.
527 */
529 const int nDim, const Array<OneD, const Array<OneD, NekDouble>> &Fwd,
530 const Array<OneD, const Array<OneD, NekDouble>> &Bwd,
532{
533 int nPts = Fwd[0].size();
534
536 {
537 ASSERTL1(CheckVectors("N"), "N not defined.");
538 ASSERTL1(CheckAuxVec("vecLocs"), "vecLocs not defined.");
540 m_vectors["N"]();
542 m_auxVec["vecLocs"]();
543
544 v_CalcFluxJacobian(nDim, Fwd, Bwd, normals, FJac, BJac);
545 }
546 else
547 {
549 for (int i = 0; i < nDim; i++)
550 {
551 normals[i] = Array<OneD, NekDouble>(nPts, 0.0);
552 }
553 Vmath::Fill(nPts, 1.0, normals[0], 1);
554
555 v_CalcFluxJacobian(nDim, Fwd, Bwd, normals, FJac, BJac);
556 }
557}
558
560 [[maybe_unused]] const int nDim,
561 [[maybe_unused]] const Array<OneD, const Array<OneD, NekDouble>> &Fwd,
562 [[maybe_unused]] const Array<OneD, const Array<OneD, NekDouble>> &Bwd,
563 [[maybe_unused]] const Array<OneD, const Array<OneD, NekDouble>> &normals,
564 [[maybe_unused]] DNekBlkMatSharedPtr &FJac,
565 [[maybe_unused]] DNekBlkMatSharedPtr &BJac)
566{
567 NEKERROR(ErrorUtil::efatal, "v_CalcFluxJacobian not specified.");
568}
569
570} // namespace Nektar::SolverUtils
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
#define EPSILON
#define CROSS(dest, v1, v2)
#define DOT(v1, v2)
Provides a generic Factory class.
SOLVER_UTILS_EXPORT bool CheckAuxVec(std::string name)
Determine whether a vector has been defined in m_auxVec.
SOLVER_UTILS_EXPORT void Solve(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, Array< OneD, Array< OneD, NekDouble > > &flux)
Perform the Riemann solve given the forwards and backwards spaces.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_rotStorage
Rotation storage.
bool m_requiresRotation
Indicates whether the Riemann solver requires a rotation to be applied to the velocity fields.
virtual SOLVER_UTILS_EXPORT void v_CalcFluxJacobian(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, const Array< OneD, const Array< OneD, NekDouble > > &normals, DNekBlkMatSharedPtr &FJac, DNekBlkMatSharedPtr &BJac)
SOLVER_UTILS_EXPORT bool CheckParams(std::string name)
Determine whether a parameter has been defined in m_params.
std::map< std::string, RSScalarFuncType > m_auxScal
Map of auxiliary scalar function types.
SOLVER_UTILS_EXPORT void rotateToNormal(const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)
Rotate a vector field to trace normal.
virtual void v_Solve(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, Array< OneD, Array< OneD, NekDouble > > &flux)=0
SOLVER_UTILS_EXPORT void rotateFromNormal(const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)
Rotate a vector field from trace normal.
Array< OneD, Array< OneD, NekDouble > > m_rotMat
Rotation matrices for each trace quadrature point.
SOLVER_UTILS_EXPORT bool CheckScalars(std::string name)
Determine whether a scalar has been defined in m_scalars.
SOLVER_UTILS_EXPORT void CalcFluxJacobian(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, DNekBlkMatSharedPtr &FJac, DNekBlkMatSharedPtr &BJac)
Calculate the flux jacobian of Fwd and Bwd.
std::map< std::string, RSVecFuncType > m_vectors
Map of vector function types.
std::map< std::string, RSScalarFuncType > m_scalars
Map of scalar function types.
SOLVER_UTILS_EXPORT bool CheckVectors(std::string name)
Determine whether a vector has been defined in m_vectors.
void FromToRotation(Array< OneD, const NekDouble > &from, Array< OneD, const NekDouble > &to, NekDouble *mat)
A function for creating a rotation matrix that rotates a vector from into another vector to.
SOLVER_UTILS_EXPORT void GenerateRotationMatrices(const Array< OneD, const Array< OneD, NekDouble > > &normals)
Generate rotation matrices for 3D expansions.
SOLVER_UTILS_EXPORT bool CheckAuxScal(std::string name)
Determine whether a scalar has been defined in m_auxScal.
std::map< std::string, RSParamFuncType > m_params
Map of parameter function types.
SOLVER_UTILS_EXPORT RiemannSolver()
bool m_ALESolver
Flag if using the ALE formulation.
std::map< std::string, RSVecFuncType > m_auxVec
Map of auxiliary vector function types.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:77
double NekDouble
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector minus vector): z = w*x - y
Definition: Vmath.hpp:381
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.hpp:439
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825