Nektar++
StdExpansion2D.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: StdExpansion2D.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: Daughter of StdExpansion. This class contains routine
32// which are common to 2D expansion. Typically this inolves physiocal
33// space operations.
34//
35///////////////////////////////////////////////////////////////////////////////
36
39
40#ifdef max
41#undef max
42#endif
43
44namespace Nektar::StdRegions
45{
46
48 [[maybe_unused]] int numcoeffs,
49 [[maybe_unused]] const LibUtilities::BasisKey &Ba,
50 [[maybe_unused]] const LibUtilities::BasisKey &Bb)
51{
52}
53
54//----------------------------
55// Differentiation Methods
56//----------------------------
58 const Array<OneD, const NekDouble> &inarray,
59 Array<OneD, NekDouble> &outarray_d0, Array<OneD, NekDouble> &outarray_d1)
60{
61 int nquad0 = m_base[0]->GetNumPoints();
62 int nquad1 = m_base[1]->GetNumPoints();
63
64 if (outarray_d0.size() > 0) // calculate du/dx_0
65 {
66 DNekMatSharedPtr D0 = m_base[0]->GetD();
67 if (inarray.data() == outarray_d0.data())
68 {
69 Array<OneD, NekDouble> wsp(nquad0 * nquad1);
70 Vmath::Vcopy(nquad0 * nquad1, inarray.data(), 1, wsp.data(), 1);
71 Blas::Dgemm('N', 'N', nquad0, nquad1, nquad0, 1.0,
72 &(D0->GetPtr())[0], nquad0, &wsp[0], nquad0, 0.0,
73 &outarray_d0[0], nquad0);
74 }
75 else
76 {
77 Blas::Dgemm('N', 'N', nquad0, nquad1, nquad0, 1.0,
78 &(D0->GetPtr())[0], nquad0, &inarray[0], nquad0, 0.0,
79 &outarray_d0[0], nquad0);
80 }
81 }
82
83 if (outarray_d1.size() > 0) // calculate du/dx_1
84 {
85 DNekMatSharedPtr D1 = m_base[1]->GetD();
86 if (inarray.data() == outarray_d1.data())
87 {
88 Array<OneD, NekDouble> wsp(nquad0 * nquad1);
89 Vmath::Vcopy(nquad0 * nquad1, inarray.data(), 1, wsp.data(), 1);
90 Blas::Dgemm('N', 'T', nquad0, nquad1, nquad1, 1.0, &wsp[0], nquad0,
91 &(D1->GetPtr())[0], nquad1, 0.0, &outarray_d1[0],
92 nquad0);
93 }
94 else
95 {
96 Blas::Dgemm('N', 'T', nquad0, nquad1, nquad1, 1.0, &inarray[0],
97 nquad0, &(D1->GetPtr())[0], nquad1, 0.0,
98 &outarray_d1[0], nquad0);
99 }
100 }
101}
102
104 const Array<OneD, const NekDouble> &coords,
105 const Array<OneD, const NekDouble> &physvals)
106{
107 ASSERTL2(coords[0] > -1 - NekConstants::kNekZeroTol, "coord[0] < -1");
108 ASSERTL2(coords[0] < 1 + NekConstants::kNekZeroTol, "coord[0] > 1");
109 ASSERTL2(coords[1] > -1 - NekConstants::kNekZeroTol, "coord[1] < -1");
110 ASSERTL2(coords[1] < 1 + NekConstants::kNekZeroTol, "coord[1] > 1");
111
113 LocCoordToLocCollapsed(coords, coll);
114
115 const int nq0 = m_base[0]->GetNumPoints();
116 const int nq1 = m_base[1]->GetNumPoints();
117
118 Array<OneD, NekDouble> wsp(nq1);
119 for (int i = 0; i < nq1; ++i)
120 {
121 wsp[i] = StdExpansion::BaryEvaluate<0>(coll[0], &physvals[0] + i * nq0);
122 }
123
124 return StdExpansion::BaryEvaluate<1>(coll[1], &wsp[0]);
125}
126
129 const Array<OneD, const NekDouble> &physvals)
130{
131 NekDouble val;
132 int i;
133 int nq0 = m_base[0]->GetNumPoints();
134 int nq1 = m_base[1]->GetNumPoints();
135 Array<OneD, NekDouble> wsp1(nq1);
136
137 // interpolate first coordinate direction
138 for (i = 0; i < nq1; ++i)
139 {
140 wsp1[i] =
141 Blas::Ddot(nq0, &(I[0]->GetPtr())[0], 1, &physvals[i * nq0], 1);
142 }
143
144 // interpolate in second coordinate direction
145 val = Blas::Ddot(nq1, I[1]->GetPtr(), 1, wsp1, 1);
146
147 return val;
148}
149
150//////////////////////////////
151// Integration Methods
152//////////////////////////////
153
157{
158 int i;
159 NekDouble Int = 0.0;
160 int nquad0 = m_base[0]->GetNumPoints();
161 int nquad1 = m_base[1]->GetNumPoints();
162 Array<OneD, NekDouble> tmp(nquad0 * nquad1);
163
164 // multiply by integration constants
165 for (i = 0; i < nquad1; ++i)
166 {
167 Vmath::Vmul(nquad0, &inarray[0] + i * nquad0, 1, w0.data(), 1,
168 &tmp[0] + i * nquad0, 1);
169 }
170
171 for (i = 0; i < nquad0; ++i)
172 {
173 Vmath::Vmul(nquad1, &tmp[0] + i, nquad0, w1.data(), 1, &tmp[0] + i,
174 nquad0);
175 }
176 Int = Vmath::Vsum(nquad0 * nquad1, tmp, 1);
177
178 return Int;
179}
180
182 const Array<OneD, const NekDouble> &base0,
183 const Array<OneD, const NekDouble> &base1,
184 const Array<OneD, const NekDouble> &inarray,
186 bool doCheckCollDir0, bool doCheckCollDir1)
187{
188 v_BwdTrans_SumFacKernel(base0, base1, inarray, outarray, wsp,
189 doCheckCollDir0, doCheckCollDir1);
190}
191
193 const Array<OneD, const NekDouble> &base0,
194 const Array<OneD, const NekDouble> &base1,
195 const Array<OneD, const NekDouble> &inarray,
197 bool doCheckCollDir0, bool doCheckCollDir1)
198{
199 v_IProductWRTBase_SumFacKernel(base0, base1, inarray, outarray, wsp,
200 doCheckCollDir0, doCheckCollDir1);
201}
202
204{
205 ASSERTL1((dir == 0) || (dir == 1), "Invalid direction.");
206
207 int nquad0 = m_base[0]->GetNumPoints();
208 int nquad1 = m_base[1]->GetNumPoints();
209 int nqtot = nquad0 * nquad1;
210 int nmodes0 = m_base[0]->GetNumModes();
211
212 Array<OneD, NekDouble> tmp1(2 * nqtot + m_ncoeffs + nmodes0 * nquad1, 0.0);
213 Array<OneD, NekDouble> tmp3(tmp1 + 2 * nqtot);
214 Array<OneD, NekDouble> tmp4(tmp1 + 2 * nqtot + m_ncoeffs);
215
216 switch (dir)
217 {
218 case 0:
219 for (int i = 0; i < nqtot; i++)
220 {
221 tmp1[i] = 1.0;
222 IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
223 m_base[1]->GetBdata(), tmp1, tmp3,
224 tmp4, false, true);
225 tmp1[i] = 0.0;
226
227 for (int j = 0; j < m_ncoeffs; j++)
228 {
229 (*mat)(j, i) = tmp3[j];
230 }
231 }
232 break;
233 case 1:
234 for (int i = 0; i < nqtot; i++)
235 {
236 tmp1[i] = 1.0;
238 m_base[1]->GetDbdata(), tmp1, tmp3,
239 tmp4, true, false);
240 tmp1[i] = 0.0;
241
242 for (int j = 0; j < m_ncoeffs; j++)
243 {
244 (*mat)(j, i) = tmp3[j];
245 }
246 }
247 break;
248 default:
249 NEKERROR(ErrorUtil::efatal, "Not a 2D expansion.");
250 break;
251 }
252}
253
255 const Array<OneD, const NekDouble> &inarray,
257{
258 if (mkey.GetNVarCoeff() == 0 &&
261 {
262 using std::max;
263
264 // This implementation is only valid when there are no
265 // coefficients associated to the Laplacian operator
266 int nquad0 = m_base[0]->GetNumPoints();
267 int nquad1 = m_base[1]->GetNumPoints();
268 int nqtot = nquad0 * nquad1;
269 int nmodes0 = m_base[0]->GetNumModes();
270 int nmodes1 = m_base[1]->GetNumModes();
271 int wspsize =
272 max(max(max(nqtot, m_ncoeffs), nquad1 * nmodes0), nquad0 * nmodes1);
273
274 const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
275 const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
276
277 // Allocate temporary storage
278 Array<OneD, NekDouble> wsp0(4 * wspsize); // size wspsize
279 Array<OneD, NekDouble> wsp1(wsp0 + wspsize); // size 3*wspsize
280
281 if (!(m_base[0]->Collocation() && m_base[1]->Collocation()))
282 {
283 // LAPLACIAN MATRIX OPERATION
284 // wsp0 = u = B * u_hat
285 // wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
286 // wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
287 BwdTrans_SumFacKernel(base0, base1, inarray, wsp0, wsp1, true,
288 true);
289 LaplacianMatrixOp_MatFree_Kernel(wsp0, outarray, wsp1);
290 }
291 else
292 {
293 LaplacianMatrixOp_MatFree_Kernel(inarray, outarray, wsp1);
294 }
295 }
296 else
297 {
299 mkey);
300 }
301}
302
304 const Array<OneD, const NekDouble> &inarray,
306{
307 if (mkey.GetNVarCoeff() == 0 &&
310 {
311 using std::max;
312
313 int nquad0 = m_base[0]->GetNumPoints();
314 int nquad1 = m_base[1]->GetNumPoints();
315 int nqtot = nquad0 * nquad1;
316 int nmodes0 = m_base[0]->GetNumModes();
317 int nmodes1 = m_base[1]->GetNumModes();
318 int wspsize =
319 max(max(max(nqtot, m_ncoeffs), nquad1 * nmodes0), nquad0 * nmodes1);
321
322 const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
323 const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
324
325 // Allocate temporary storage
326 Array<OneD, NekDouble> wsp0(5 * wspsize); // size wspsize
327 Array<OneD, NekDouble> wsp1(wsp0 + wspsize); // size wspsize
328 Array<OneD, NekDouble> wsp2(wsp0 + 2 * wspsize); // size 3*wspsize
329
330 if (!(m_base[0]->Collocation() && m_base[1]->Collocation()))
331 {
332 // MASS MATRIX OPERATION
333 // The following is being calculated:
334 // wsp0 = B * u_hat = u
335 // wsp1 = W * wsp0
336 // outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
337 BwdTrans_SumFacKernel(base0, base1, inarray, wsp0, wsp2, true,
338 true);
339 MultiplyByQuadratureMetric(wsp0, wsp1);
340 IProductWRTBase_SumFacKernel(base0, base1, wsp1, outarray, wsp2,
341 true, true);
342
343 LaplacianMatrixOp_MatFree_Kernel(wsp0, wsp1, wsp2);
344 }
345 else
346 {
347 MultiplyByQuadratureMetric(inarray, outarray);
348 LaplacianMatrixOp_MatFree_Kernel(inarray, wsp1, wsp2);
349 }
350
351 // outarray = lambda * outarray + wsp1
352 // = (lambda * M + L ) * u_hat
353 Vmath::Svtvp(m_ncoeffs, lambda, &outarray[0], 1, &wsp1[0], 1,
354 &outarray[0], 1);
355 }
356 else
357 {
359 mkey);
360 }
361}
362
364 [[maybe_unused]] const unsigned int traceid,
365 [[maybe_unused]] Array<OneD, unsigned int> &maparray)
366{
367 ASSERTL0(false,
368 "This method must be defined at the individual shape level");
369}
370
371/** \brief Determine the mapping to re-orientate the
372 coefficients along the element trace (assumed to align
373 with the standard element) into the orientation of the
374 local trace given by edgeOrient.
375 */
377 const unsigned int eid, Array<OneD, unsigned int> &maparray,
378 Array<OneD, int> &signarray, Orientation edgeOrient, int P,
379 [[maybe_unused]] int Q)
380{
381 unsigned int i;
382
383 int dir;
384 // determine basis direction for edge.
386 {
387 dir = (eid == 0) ? 0 : 1;
388 }
389 else
390 {
391 dir = eid % 2;
392 }
393
394 int numModes = m_base[dir]->GetNumModes();
395
396 // P is the desired length of the map
397 P = (P == -1) ? numModes : P;
398
399 // decalare maparray
400 if (maparray.size() != P)
401 {
402 maparray = Array<OneD, unsigned int>(P);
403 }
404
405 // fill default mapping as increasing index
406 for (i = 0; i < P; ++i)
407 {
408 maparray[i] = i;
409 }
410
411 if (signarray.size() != P)
412 {
413 signarray = Array<OneD, int>(P, 1);
414 }
415 else
416 {
417 std::fill(signarray.data(), signarray.data() + P, 1);
418 }
419
420 // Zero signmap and set maparray to zero if
421 // elemental modes are not as large as trace modes
422 for (i = numModes; i < P; ++i)
423 {
424 signarray[i] = 0.0;
425 maparray[i] = maparray[0];
426 }
427
428 if (edgeOrient == eBackwards)
429 {
430 const LibUtilities::BasisType bType = GetBasisType(dir);
431
432 if ((bType == LibUtilities::eModified_A) ||
433 (bType == LibUtilities::eModified_B))
434 {
435 std::swap(maparray[0], maparray[1]);
436
437 for (i = 3; i < std::min(P, numModes); i += 2)
438 {
439 signarray[i] *= -1;
440 }
441 }
442 else if (bType == LibUtilities::eGLL_Lagrange ||
444 {
445 ASSERTL1(P == numModes, "Different trace space edge dimension "
446 "and element edge dimension not currently "
447 "possible for GLL-Lagrange bases");
448
449 std::reverse(maparray.data(), maparray.data() + P);
450 }
451 else
452 {
453 ASSERTL0(false, "Mapping not defined for this type of basis");
454 }
455 }
456}
457
460 Array<OneD, int> &signarray,
461 Orientation edgeOrient, int P,
462 int Q)
463{
464 Array<OneD, unsigned int> map1, map2;
465 v_GetTraceCoeffMap(eid, map1);
466 v_GetElmtTraceToTraceMap(eid, map2, signarray, edgeOrient, P, Q);
467
468 if (maparray.size() != map2.size())
469 {
470 maparray = Array<OneD, unsigned int>(map2.size());
471 }
472
473 for (int i = 0; i < map2.size(); ++i)
474 {
475 maparray[i] = map1[map2[i]];
476 }
477}
478
479void StdExpansion2D::v_PhysInterp(std::shared_ptr<StdExpansion> fromExp,
480 const Array<OneD, const NekDouble> &fromData,
482{
483
484 LibUtilities::Interp2D(fromExp->GetBasis(0)->GetPointsKey(),
485 fromExp->GetBasis(1)->GetPointsKey(), fromData,
486 m_base[0]->GetPointsKey(), m_base[1]->GetPointsKey(),
487 toData);
488}
489
490} // namespace Nektar::StdRegions
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#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 ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:265
Describes the specification for a Basis.
Definition: Basis.h:45
void v_GenStdMatBwdDeriv(const int dir, DNekMatSharedPtr &mat) override
virtual void v_IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1)=0
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d0, Array< OneD, NekDouble > &outarray_d1)
Calculate the 2D derivative in the local tensor/collapsed coordinate at the physical points.
void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_GetTraceToElementMap(const int eid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation edgeOrient=eForwards, int P=-1, int Q=-1) override
virtual void v_BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1)=0
NekDouble v_PhysEvaluateInterp(const Array< OneD, DNekMatSharedPtr > &I, const Array< OneD, const NekDouble > &physvals) override
NekDouble Integral(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &w0, const Array< OneD, const NekDouble > &w1)
void BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0=true, bool doCheckCollDir1=true)
void v_GetTraceCoeffMap(const unsigned int traceid, Array< OneD, unsigned int > &maparray) override
void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0=true, bool doCheckCollDir1=true)
void v_PhysInterp(std::shared_ptr< StdExpansion > fromExp, const Array< OneD, const NekDouble > &fromData, Array< OneD, NekDouble > &toData) override
NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals) override
This function evaluates the expansion at a single (arbitrary) point of the domain.
void v_GetElmtTraceToTraceMap(const unsigned int eid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation edgeOrient, int P, int Q) override
Determine the mapping to re-orientate the coefficients along the element trace (assumed to align with...
void LaplacianMatrixOp_MatFree_GenericImpl(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:156
void LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:370
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:732
void HelmholtzMatrixOp_MatFree_GenericImpl(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:124
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:133
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:163
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:383
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:101
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
@ P
Monomial polynomials .
Definition: BasisType.h:62
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:57
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:56
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
static const NekDouble kNekZeroTol
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
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 Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
Definition: Vmath.hpp:396
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825