Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_DerivOrthogPolyExpansionImp.hpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#include "Teuchos_Assert.hpp"
45#include "Stokhos_DynamicArrayTraits.hpp"
46
47template <typename ordinal_type, typename value_type>
48Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
49DerivOrthogPolyExpansion(
50 const Teuchos::RCP<const Stokhos::DerivBasis<ordinal_type, value_type> >& basis_,
51 const Teuchos::RCP<const Teuchos::SerialDenseMatrix<ordinal_type, value_type> >& Bij_,
52 const Teuchos::RCP<const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk_,
53 const Teuchos::RCP<const Stokhos::Dense3Tensor<ordinal_type, value_type> >& Dijk_) :
54 basis(basis_),
55 Bij(Bij_),
56 Cijk(Cijk_),
57 Dijk(Dijk_),
58 sz(basis->size()),
59 A(2*sz,2*sz),
60 B(2*sz,2),
61 piv(2*sz),
62 lapack()
63{
64}
65
66extern "C" {
67 double dlange_(char*, int*, int*, double*, int*, double*);
68}
69
70template <typename ordinal_type, typename value_type>
71ordinal_type
72Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
73solve(ordinal_type s, ordinal_type nrhs)
74{
75 if (s == 0 || nrhs == 0)
76 return 0;
77
78 ordinal_type info;
79// lapack.GESV(s, nrhs, A.values(), A.numRows(), &(piv[0]), b.values(),
80// b.numRows(), &info);
81 lapack.GETRF(s, s, A.values(), A.numRows(), &(piv[0]), &info);
82 value_type norm, rcond;
83 Teuchos::Array<ordinal_type> iwork(4*s);
84 Teuchos::Array<value_type> work(4*s);
85 norm = 1.0;
86 ordinal_type n = A.numRows();
87 char t = '1';
88 norm = dlange_(&t, &s, &s, A.values(), &n, &work[0]);
89 lapack.GECON('1', s, A.values(), A.numRows(), norm, &rcond, &work[0],
90 &iwork[0], &info);
91 //std::cout << "condition number = " << 1.0/rcond << std::endl;
92 lapack.GETRS('N', s, nrhs, A.values(), A.numRows(), &(piv[0]), B.values(),
93 B.numRows(), &info);
94 return info;
95}
96
97template <typename ordinal_type, typename value_type>
98void
99Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
100unaryMinus(
101 Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
102 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
103{
104 ordinal_type pc = a.size();
105 if (c.size() != pc)
106 c.resize(pc);
107
108 value_type* cc = c.coeff();
109 const value_type* ca = a.coeff();
110
111 for (ordinal_type i=0; i<pc; i++)
112 cc[i] = -ca[i];
113}
114
115template <typename ordinal_type, typename value_type>
116void
123
124template <typename ordinal_type, typename value_type>
125void
132
133template <typename ordinal_type, typename value_type>
134void
144
145template <typename ordinal_type, typename value_type>
146void
156
157template <typename ordinal_type, typename value_type>
158void
159Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
160plusEqual(
161 Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
162 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& x)
163{
164 ordinal_type xp = x.size();
165 if (c.size() < xp)
166 c.resize(xp);
167
168 value_type* cc = c.coeff();
169 const value_type* xc = x.coeff();
170 for (ordinal_type i=0; i<xp; i++)
171 cc[i] += xc[i];
172}
173
174template <typename ordinal_type, typename value_type>
175void
176Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
177minusEqual(
178 Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
179 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& x)
180{
181 ordinal_type xp = x.size();
182 if (c.size() < xp)
183 c.resize(xp);
184
185 value_type* cc = c.coeff();
186 const value_type* xc = x.coeff();
187 for (ordinal_type i=0; i<xp; i++)
188 cc[i] -= xc[i];
189}
190
191template <typename ordinal_type, typename value_type>
192void
193Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
194timesEqual(
195 Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
196 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& x)
197{
198 ordinal_type p = c.size();
199 ordinal_type xp = x.size();
200 ordinal_type pc;
201 if (p > 1 && xp > 1)
202 pc = sz;
203 else
204 pc = p*xp;
205 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
206 "Stokhos::DerivOrthogPolyExpansion::timesEqual()" <<
207 ": Expansion size (" << sz <<
208 ") is too small for computation.");
209 if (c.size() != pc)
210 c.resize(pc);
211
212 value_type* cc = c.coeff();
213 const value_type* xc = x.coeff();
214
215 if (p > 1 && xp > 1) {
216 // Copy c coefficients into temporary array
217 value_type* tc = Stokhos::ds_array<value_type>::get_and_fill(cc,p);
218 value_type tmp, cijk;
219 ordinal_type i,j;
220 for (ordinal_type k=0; k<pc; k++) {
221 tmp = value_type(0.0);
222 ordinal_type n = Cijk->num_values(k);
223 for (ordinal_type l=0; l<n; l++) {
224 Cijk->value(k,l,i,j,cijk);
225 if (i < p && j < xp)
226 tmp += cijk*tc[i]*xc[j];
227 }
228 cc[k] = tmp / basis->norm_squared(k);
229 }
230 }
231 else if (p > 1) {
232 for (ordinal_type i=0; i<p; i++)
233 cc[i] *= xc[0];
234 }
235 else if (xp > 1) {
236 for (ordinal_type i=1; i<xp; i++)
237 cc[i] = cc[0]*xc[i];
238 cc[0] *= xc[0];
239 }
240 else {
241 cc[0] *= xc[0];
242 }
243}
244
245template <typename ordinal_type, typename value_type>
246void
247Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
248divideEqual(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
249 const OrthogPolyApprox<ordinal_type, value_type, node_type>& x)
250{
251 const char* func = "Stokhos::DerivOrthogPolyExpansion::divide()";
252 ordinal_type p = c.size();
253 ordinal_type xp = x.size();
254 ordinal_type pc;
255 if (xp > 1)
256 pc = sz;
257 else
258 pc = p;
259 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
260 "Stokhos::DerivOrthogPolyExpansion::divideEqual()" <<
261 ": Expansion size (" << sz <<
262 ") is too small for computation.");
263 if (c.size() != pc)
264 c.resize(pc);
265
266 value_type* cc = c.coeff();
267 const value_type* xc = x.coeff();
268
269 if (xp > 1) {
270
271 // Fill A
272 value_type cijk;
273 ordinal_type i,j;
274 for (ordinal_type i=0; i<pc; i++)
275 for (ordinal_type j=0; j<pc; j++)
276 A(i,j) = 0.0;
277 for (ordinal_type k=0; k<xp; k++) {
278 ordinal_type n = Cijk->num_values(k);
279 for (ordinal_type l=0; l<n; l++) {
280 Cijk->value(k,l,i,j,cijk);
281 A(i,j) += cijk*xc[k]/basis->norm_squared(i);
282 }
283 }
284
285 // Fill B
286 for (ordinal_type i=0; i<p; i++)
287 B(i,0) = cc[i];
288 for (ordinal_type i=p; i<pc; i++)
289 B(i,0) = value_type(0.0);
290
291 // Solve system
292 int info = solve(pc, 1);
293
294 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
295 func << ": Argument " << info
296 << " for solve had illegal value");
297 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
298 func << ": Diagonal entry " << info
299 << " in LU factorization is exactly zero");
300
301 // Get coefficients
302 for (ordinal_type i=0; i<pc; i++)
303 cc[i] = B(i,0);
304 }
305 else {
306 for (ordinal_type i=0; i<p; i++)
307 cc[i] /= xc[0];
308 }
309}
310
311template <typename ordinal_type, typename value_type>
312void
313Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
314plus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
315 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
316 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
317{
318 ordinal_type pa = a.size();
319 ordinal_type pb = b.size();
320 ordinal_type pc = pa > pb ? pa : pb;
321 if (c.size() != pc)
322 c.resize(pc);
323
324 const value_type* ca = a.coeff();
325 const value_type* cb = b.coeff();
326 value_type* cc = c.coeff();
327
328 if (pa > pb) {
329 for (ordinal_type i=0; i<pb; i++)
330 cc[i] = ca[i] + cb[i];
331 for (ordinal_type i=pb; i<pc; i++)
332 cc[i] = ca[i];
333 }
334 else {
335 for (ordinal_type i=0; i<pa; i++)
336 cc[i] = ca[i] + cb[i];
337 for (ordinal_type i=pa; i<pc; i++)
338 cc[i] = cb[i];
339 }
340}
341
342template <typename ordinal_type, typename value_type>
343void
344Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
345plus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
346 const value_type& a,
347 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
348{
349 ordinal_type pc = b.size();
350 if (c.size() != pc)
351 c.resize(pc);
352
353 const value_type* cb = b.coeff();
354 value_type* cc = c.coeff();
355
356 cc[0] = a + cb[0];
357 for (ordinal_type i=1; i<pc; i++)
358 cc[i] = cb[i];
359}
360
361template <typename ordinal_type, typename value_type>
362void
363Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
364plus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
365 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
366 const value_type& b)
367{
368 ordinal_type pc = a.size();
369 if (c.size() != pc)
370 c.resize(pc);
371
372 const value_type* ca = a.coeff();
373 value_type* cc = c.coeff();
374
375 cc[0] = ca[0] + b;
376 for (ordinal_type i=1; i<pc; i++)
377 cc[i] = ca[i];
378}
379
380template <typename ordinal_type, typename value_type>
381void
382Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
383minus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
384 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
385 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
386{
387 ordinal_type pa = a.size();
388 ordinal_type pb = b.size();
389 ordinal_type pc = pa > pb ? pa : pb;
390 if (c.size() != pc)
391 c.resize(pc);
392
393 const value_type* ca = a.coeff();
394 const value_type* cb = b.coeff();
395 value_type* cc = c.coeff();
396
397 if (pa > pb) {
398 for (ordinal_type i=0; i<pb; i++)
399 cc[i] = ca[i] - cb[i];
400 for (ordinal_type i=pb; i<pc; i++)
401 cc[i] = ca[i];
402 }
403 else {
404 for (ordinal_type i=0; i<pa; i++)
405 cc[i] = ca[i] - cb[i];
406 for (ordinal_type i=pa; i<pc; i++)
407 cc[i] = -cb[i];
408 }
409}
410
411template <typename ordinal_type, typename value_type>
412void
413Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
414minus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
415 const value_type& a,
416 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
417{
418 ordinal_type pc = b.size();
419 if (c.size() != pc)
420 c.resize(pc);
421
422 const value_type* cb = b.coeff();
423 value_type* cc = c.coeff();
424
425 cc[0] = a - cb[0];
426 for (ordinal_type i=1; i<pc; i++)
427 cc[i] = -cb[i];
428}
429
430template <typename ordinal_type, typename value_type>
431void
432Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
433minus(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
434 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
435 const value_type& b)
436{
437 ordinal_type pc = a.size();
438 if (c.size() != pc)
439 c.resize(pc);
440
441 const value_type* ca = a.coeff();
442 value_type* cc = c.coeff();
443
444 cc[0] = ca[0] - b;
445 for (ordinal_type i=1; i<pc; i++)
446 cc[i] = ca[i];
447}
448
449template <typename ordinal_type, typename value_type>
450void
451Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
452times(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
453 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
454 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
455{
456 ordinal_type pa = a.size();
457 ordinal_type pb = b.size();
458 ordinal_type pc;
459 if (pa > 1 && pb > 1)
460 pc = sz;
461 else
462 pc = pa*pb;
463 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
464 "Stokhos::DerivOrthogPolyExpansion::times()" <<
465 ": Expansion size (" << sz <<
466 ") is too small for computation.");
467 if (c.size() != pc)
468 c.resize(pc);
469
470 const value_type* ca = a.coeff();
471 const value_type* cb = b.coeff();
472 value_type* cc = c.coeff();
473
474 if (pa > 1 && pb > 1) {
475 value_type tmp, cijk;
476 ordinal_type i,j;
477 for (ordinal_type k=0; k<pc; k++) {
478 tmp = value_type(0.0);
479 ordinal_type n = Cijk->num_values(k);
480 for (ordinal_type l=0; l<n; l++) {
481 Cijk->value(k,l,i,j,cijk);
482 if (i < pa && j < pb)
483 tmp += cijk*ca[i]*cb[j];
484 }
485 cc[k] = tmp / basis->norm_squared(k);
486 }
487 }
488 else if (pa > 1) {
489 for (ordinal_type i=0; i<pc; i++)
490 cc[i] = ca[i]*cb[0];
491 }
492 else if (pb > 1) {
493 for (ordinal_type i=0; i<pc; i++)
494 cc[i] = ca[0]*cb[i];
495 }
496 else {
497 cc[0] = ca[0]*cb[0];
498 }
499}
500
501template <typename ordinal_type, typename value_type>
502void
503Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
504times(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
505 const value_type& a,
506 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
507{
508 ordinal_type pc = b.size();
509 if (c.size() != pc)
510 c.resize(pc);
511
512 const value_type* cb = b.coeff();
513 value_type* cc = c.coeff();
514
515 for (ordinal_type i=0; i<pc; i++)
516 cc[i] = a*cb[i];
517}
518
519template <typename ordinal_type, typename value_type>
520void
521Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
522times(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
523 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
524 const value_type& b)
525{
526 ordinal_type pc = a.size();
527 if (c.size() != pc)
528 c.resize(pc);
529
530 const value_type* ca = a.coeff();
531 value_type* cc = c.coeff();
532
533 for (ordinal_type i=0; i<pc; i++)
534 cc[i] = ca[i]*b;
535}
536
537template <typename ordinal_type, typename value_type>
538void
539Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
540divide(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
541 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
542 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
543{
544 const char* func = "Stokhos::DerivOrthogPolyExpansion::divide()";
545 ordinal_type pa = a.size();
546 ordinal_type pb = b.size();
547 ordinal_type pc;
548 if (pb > 1)
549 pc = sz;
550 else
551 pc = pa;
552 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
553 "Stokhos::DerivOrthogPolyExpansion::divide()" <<
554 ": Expansion size (" << sz <<
555 ") is too small for computation.");
556 if (c.size() != pc)
557 c.resize(pc);
558
559 const value_type* ca = a.coeff();
560 const value_type* cb = b.coeff();
561 value_type* cc = c.coeff();
562
563 if (pb > 1) {
564
565 // Fill A
566 value_type cijk;
567 ordinal_type i,j;
568 for (ordinal_type i=0; i<pc; i++)
569 for (ordinal_type j=0; j<pc; j++)
570 A(i,j) = 0.0;
571
572 for (ordinal_type k=0; k<pb; k++) {
573 ordinal_type n = Cijk->num_values(k);
574 for (ordinal_type l=0; l<n; l++) {
575 Cijk->value(k,l,i,j,cijk);
576 A(i,j) += cijk*cb[k] / basis->norm_squared(i);
577 }
578 }
579
580 // Fill B
581 for (ordinal_type i=0; i<pa; i++)
582 B(i,0) = ca[i];
583 for (ordinal_type i=pa; i<pc; i++)
584 B(i,0) = value_type(0.0);
585
586 // Solve system
587 int info = solve(pc, 1);
588
589 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
590 func << ": Argument " << info
591 << " for solve had illegal value");
592 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
593 func << ": Diagonal entry " << info
594 << " in LU factorization is exactly zero");
595
596 // Get coefficients
597 for (ordinal_type i=0; i<pc; i++)
598 cc[i] = B(i,0);
599 }
600 else {
601 for (ordinal_type i=0; i<pa; i++)
602 cc[i] = ca[i]/cb[0];
603 }
604}
605
606template <typename ordinal_type, typename value_type>
607void
608Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
609divide(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
610 const value_type& a,
611 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
612{
613 const char* func = "Stokhos::DerivOrthogPolyExpansion::divide()";
614 ordinal_type pb = b.size();
615 ordinal_type pc;
616 if (pb > 1)
617 pc = sz;
618 else
619 pc = 1;
620 if (c.size() != pc)
621 c.resize(pc);
622
623 const value_type* cb = b.coeff();
624 value_type* cc = c.coeff();
625
626 if (pb > 1) {
627
628 // Fill A
629 value_type cijk;
630 ordinal_type i,j;
631 for (ordinal_type i=0; i<pc; i++)
632 for (ordinal_type j=0; j<pc; j++)
633 A(i,j) = 0.0;
634
635 for (ordinal_type k=0; k<pb; k++) {
636 ordinal_type n = Cijk->num_values(k);
637 for (ordinal_type l=0; l<n; l++) {
638 Cijk->value(k,l,i,j,cijk);
639 A(i,j) += cijk*cb[k] / basis->norm_squared(i);
640 }
641 }
642
643 // Fill B
644 B(0,0) = a;
645 for (ordinal_type i=1; i<pc; i++)
646 B(i,0) = value_type(0.0);
647
648 // Solve system
649 int info = solve(pc, 1);
650
651 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
652 func << ": Argument " << info
653 << " for solve had illegal value");
654 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
655 func << ": Diagonal entry " << info
656 << " in LU factorization is exactly zero");
657
658 // Get coefficients
659 for (ordinal_type i=0; i<pc; i++)
660 cc[i] = B(i,0);
661 }
662 else
663 cc[0] = a / cb[0];
664}
665
666template <typename ordinal_type, typename value_type>
667void
668Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
669divide(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
670 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
671 const value_type& b)
672{
673 ordinal_type pc = a.size();
674 if (c.size() != pc)
675 c.resize(pc);
676
677 const value_type* ca = a.coeff();
678 value_type* cc = c.coeff();
679
680 for (ordinal_type i=0; i<pc; i++)
681 cc[i] = ca[i]/b;
682}
683
684template <typename ordinal_type, typename value_type>
685void
686Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
687exp(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
688 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
689{
690 const char* func = "Stokhos::DerivOrthogPolyExpansion::exp()";
691 ordinal_type pa = a.size();
692 ordinal_type pc;
693 if (pa > 1)
694 pc = sz;
695 else
696 pc = 1;
697 if (c.size() != pc)
698 c.resize(pc);
699
700 const value_type* ca = a.coeff();
701 value_type* cc = c.coeff();
702
703 if (pa > 1) {
704 value_type psi_0 = basis->evaluateZero(0);
705
706 // Fill A and B
707 for (ordinal_type i=1; i<pc; i++) {
708 B(i-1,0) = 0.0;
709 for (ordinal_type j=1; j<pc; j++) {
710 A(i-1,j-1) = (*Bij)(i-1,j);
711 for (ordinal_type k=1; k<pa; k++)
712 A(i-1,j-1) -= ca[k]*(*Dijk)(i-1,j,k);
713 B(i-1,0) += ca[j]*(*Bij)(i-1,j);
714 }
715 B(i-1,0) *= psi_0;
716 }
717
718 // Solve system
719 int info = solve(pc-1, 1);
720
721 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
722 func << ": Argument " << info
723 << " for solve had illegal value");
724 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
725 func << ": Diagonal entry " << info
726 << " in LU factorization is exactly zero");
727
728 // Compute order-0 coefficient
729 value_type s = psi_0 * ca[0];
730 value_type t = psi_0;
731 for (ordinal_type i=1; i<pc; i++) {
732 s += basis->evaluateZero(i) * ca[i];
733 t += basis->evaluateZero(i) * B(i-1,0);
734 }
735 s = std::exp(s);
736 cc[0] = (s/t);
737
738 // Compute remaining coefficients
739 for (ordinal_type i=1; i<pc; i++)
740 cc[i] = B(i-1,0) * cc[0];
741 }
742 else
743 cc[0] = std::exp(ca[0]);
744}
745
746template <typename ordinal_type, typename value_type>
747void
748Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
749log(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
750 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
751{
752 const char* func = "Stokhos::DerivOrthogPolyExpansion::log()";
753 ordinal_type pa = a.size();
754 ordinal_type pc;
755 if (pa > 1)
756 pc = sz;
757 else
758 pc = 1;
759 if (c.size() != pc)
760 c.resize(pc);
761
762 const value_type* ca = a.coeff();
763 value_type* cc = c.coeff();
764
765 if (pa > 1) {
766 value_type psi_0 = basis->evaluateZero(0);
767
768 // Fill A and B
769 for (ordinal_type i=1; i<pc; i++) {
770 B(i-1,0) = 0.0;
771 for (ordinal_type j=1; j<pc; j++) {
772 A(i-1,j-1) = 0.0;
773 for (ordinal_type k=0; k<pa; k++)
774 A(i-1,j-1) += ca[k]*(*Dijk)(i-1,k,j);
775 B(i-1,0) += ca[j]*(*Bij)(i-1,j);
776 }
777 }
778
779 // Solve system
780 int info = solve(pc-1, 1);
781
782 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
783 func << ": Argument " << info
784 << " for solve had illegal value");
785 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
786 func << ": Diagonal entry " << info
787 << " in LU factorization is exactly zero");
788
789 // Compute order-0 coefficient
790 value_type s = psi_0 * ca[0];
791 value_type t = value_type(0.0);
792 for (ordinal_type i=1; i<pc; i++) {
793 s += basis->evaluateZero(i) * ca[i];
794 t += basis->evaluateZero(i) * B(i-1,0);
795 }
796 cc[0] = (std::log(s) - t) / psi_0;
797
798 // Compute remaining coefficients
799 for (ordinal_type i=1; i<pc; i++)
800 cc[i] = B(i-1,0);
801 }
802 else
803 cc[0] = std::log(ca[0]);
804}
805
806template <typename ordinal_type, typename value_type>
807void
808Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
809log10(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
810 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
811{
812 if (a.size() > 1) {
813 log(c,a);
814 divide(c,c,std::log(10.0));
815 }
816 else {
817 if (c.size() != 1)
818 c.resize(1);
819 c[0] = std::log10(a[0]);
820 }
821}
822
823template <typename ordinal_type, typename value_type>
824void
825Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
826sqrt(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
827 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
828{
829 if (a.size() > 1) {
830 log(c,a);
831 timesEqual(c,value_type(0.5));
832 exp(c,c);
833 }
834 else {
835 if (c.size() != 1)
836 c.resize(1);
837 c[0] = std::sqrt(a[0]);
838 }
839}
840
841template <typename ordinal_type, typename value_type>
842void
843Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
844cbrt(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
845 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
846{
847 if (a.size() > 1) {
848 log(c,a);
849 timesEqual(c,value_type(1.0/3.0));
850 exp(c,c);
851 }
852 else {
853 if (c.size() != 1)
854 c.resize(1);
855 c[0] = std::cbrt(a[0]);
856 }
857}
858
859template <typename ordinal_type, typename value_type>
860void
861Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
862pow(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
863 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
864 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
865{
866 if (a.size() > 1 || b.size() > 1) {
867 log(c,a);
868 timesEqual(c,b);
869 exp(c,c);
870 }
871 else {
872 if (c.size() != 1)
873 c.resize(1);
874 c[0] = std::pow(a[0], b[0]);
875 }
876}
877
878template <typename ordinal_type, typename value_type>
879void
880Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
881pow(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
882 const value_type& a,
883 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
884{
885 if (b.size() > 1) {
886 times(c,std::log(a),b);
887 exp(c,c);
888 }
889 else {
890 if (c.size() != 1)
891 c.resize(1);
892 c[0] = std::pow(a, b[0]);
893 }
894}
895
896template <typename ordinal_type, typename value_type>
897void
898Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
899pow(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
900 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
901 const value_type& b)
902{
903 if (a.size() > 1) {
904 log(c,a);
905 timesEqual(c,b);
906 exp(c,c);
907 }
908 else {
909 if (c.size() != 1)
910 c.resize(1);
911 c[0] = std::pow(a[0], b);
912 }
913}
914
915template <typename ordinal_type, typename value_type>
916void
917Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
918sincos(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& s,
919 Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
920 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
921{
922 const char* func = "Stokhos::DerivOrthogPolyExpansion::sincos()";
923 ordinal_type pa = a.size();
924 ordinal_type pc;
925 if (pa > 1)
926 pc = sz;
927 else
928 pc = 1;
929 if (c.size() != pc)
930 c.resize(pc);
931 if (s.size() != pc)
932 s.resize(pc);
933
934 const value_type* ca = a.coeff();
935 value_type* cs = s.coeff();
936 value_type* cc = c.coeff();
937
938 if (pa > 1) {
939 value_type psi_0 = basis->evaluateZero(0);
940 ordinal_type offset = pc-1;
941
942 // Fill A and b
943 B.putScalar(value_type(0.0));
944 value_type tmp, tmp2;
945 for (ordinal_type i=1; i<pc; i++) {
946 tmp2 = value_type(0.0);
947 for (ordinal_type j=1; j<pc; j++) {
948 tmp = (*Bij)(i-1,j);
949 A(i-1,j-1) = tmp;
950 A(i-1+offset,j-1+offset) = tmp;
951 tmp = value_type(0.0);
952 for (ordinal_type k=1; k<pa; k++)
953 tmp += ca[k]*(*Dijk)(i-1,j,k);
954 A(i-1+offset,j-1) = tmp;
955 A(i-1,j-1+offset) = -tmp;
956 tmp2 += ca[j]*(*Bij)(i-1,j);
957 }
958 B(i-1,0) = tmp2*psi_0;
959 B(i-1+offset,1) = tmp2*psi_0;
960 }
961
962 // Solve system
963 int info = solve(2*pc-2, 2);
964
965 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
966 func << ": Argument " << info
967 << " for solve had illegal value");
968 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
969 func << ": Diagonal entry " << info
970 << " in LU factorization is exactly zero");
971
972 // Compute order-0 coefficients
973 value_type t = psi_0 * ca[0];
974 value_type a00 = psi_0;
975 value_type a01 = value_type(0.0);
976 value_type a10 = value_type(0.0);
977 value_type a11 = psi_0;
978 value_type b0 = B(0,0);
979 value_type b1 = B(1,0);
980 for (ordinal_type i=1; i<pc; i++) {
981 t += basis->evaluateZero(i) * ca[i];
982 a00 -= basis->evaluateZero(i) * B(i-1,1);
983 a01 += basis->evaluateZero(i) * B(i-1,0);
984 a10 -= basis->evaluateZero(i) * B(i-1+offset,1);
985 a11 += basis->evaluateZero(i) * B(i-1+offset,0);
986 }
987 A(0,0) = a00;
988 A(0,1) = a01;
989 A(1,0) = a10;
990 A(1,1) = a11;
991 B(0,0) = std::sin(t);
992 B(1,0) = std::cos(t);
993
994 info = solve(2, 1);
995
996 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
997 func << ": Argument " << info
998 << " for (2x2) solve had illegal value");
999 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1000 func << ": Diagonal entry " << info
1001 << " in (2x2) LU factorization is exactly zero");
1002 cs[0] = B(0,0);
1003 cc[0] = B(1,0);
1004
1005 // Compute remaining coefficients
1006 B(0,0) = b0;
1007 B(1,0) = b1;
1008 for (ordinal_type i=1; i<pc; i++) {
1009 cs[i] = cc[0]*B(i-1,0) - cs[0]*B(i-1,1);
1010 cc[i] = cc[0]*B(i-1+offset,0) - cs[0]*B(i-1+offset,1);
1011 }
1012 }
1013 else {
1014 cs[0] = std::sin(ca[0]);
1015 cc[0] = std::cos(ca[0]);
1016 }
1017}
1018
1019template <typename ordinal_type, typename value_type>
1020void
1035
1036template <typename ordinal_type, typename value_type>
1037void
1052
1053template <typename ordinal_type, typename value_type>
1054void
1055Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1056tan(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& t,
1057 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1058{
1059 if (a.size() > 1) {
1060 OrthogPolyApprox<ordinal_type, value_type, node_type> c(t);
1061 sincos(t, c, a);
1062 divideEqual(t,c);
1063 }
1064 else {
1065 if (t.size() != 1)
1066 t.resize(1);
1067 t[0] = std::tan(a[0]);
1068 }
1069}
1070
1071template <typename ordinal_type, typename value_type>
1072void
1073Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1074sinhcosh(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& s,
1075 Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
1076 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1077{
1078 const char* func = "Stokhos::DerivOrthogPolyExpansion::sinhcosh()";
1079 ordinal_type pa = a.size();
1080 ordinal_type pc;
1081 if (pa > 1)
1082 pc = sz;
1083 else
1084 pc = 1;
1085 if (c.size() != pc)
1086 c.resize(pc);
1087 if (s.size() != pc)
1088 s.resize(pc);
1089
1090 const value_type* ca = a.coeff();
1091 value_type* cs = s.coeff();
1092 value_type* cc = c.coeff();
1093
1094 if (pa > 1) {
1095 value_type psi_0 = basis->evaluateZero(0);
1096 ordinal_type offset = pc-1;
1097
1098 // Fill A and b
1099 B.putScalar(value_type(0.0));
1100 value_type tmp, tmp2;
1101 for (ordinal_type i=1; i<pc; i++) {
1102 tmp2 = value_type(0.0);
1103 for (ordinal_type j=1; j<pc; j++) {
1104 tmp = (*Bij)(i-1,j);
1105 A(i-1,j-1) = tmp;
1106 A(i-1+offset,j-1+offset) = tmp;
1107 tmp = value_type(0.0);
1108 for (ordinal_type k=1; k<pa; k++)
1109 tmp += ca[k]*(*Dijk)(i-1,j,k);
1110 A(i-1+offset,j-1) = -tmp;
1111 A(i-1,j-1+offset) = -tmp;
1112 tmp2 += ca[j]*(*Bij)(i-1,j);
1113 }
1114 B(i-1,0) = tmp2*psi_0;
1115 B(i-1+offset,1) = tmp2*psi_0;
1116 }
1117
1118 // Solve system
1119 int info = solve(2*pc-2, 2);
1120
1121 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
1122 func << ": Argument " << info
1123 << " for solve had illegal value");
1124 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1125 func << ": Diagonal entry " << info
1126 << " in LU factorization is exactly zero");
1127
1128 // Compute order-0 coefficients
1129 value_type t = psi_0 * ca[0];
1130 value_type a00 = psi_0;
1131 value_type a01 = value_type(0.0);
1132 value_type a10 = value_type(0.0);
1133 value_type a11 = psi_0;
1134 value_type b0 = B(0,0);
1135 value_type b1 = B(1,0);
1136 for (ordinal_type i=1; i<pc; i++) {
1137 t += basis->evaluateZero(i) * ca[i];
1138 a00 += basis->evaluateZero(i) * B(i-1,1);
1139 a01 += basis->evaluateZero(i) * B(i-1,0);
1140 a10 += basis->evaluateZero(i) * B(i-1+offset,1);
1141 a11 += basis->evaluateZero(i) * B(i-1+offset,0);
1142 }
1143 A(0,0) = a00;
1144 A(0,1) = a01;
1145 A(1,0) = a10;
1146 A(1,1) = a11;
1147 B(0,0) = std::sinh(t);
1148 B(1,0) = std::cosh(t);
1149 info = solve(2, 1);
1150 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
1151 func << ": Argument " << info
1152 << " for (2x2) solve had illegal value");
1153 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1154 func << ": Diagonal entry " << info
1155 << " in (2x2) LU factorization is exactly zero");
1156 cs[0] = B(0,0);
1157 cc[0] = B(1,0);
1158
1159 // Compute remaining coefficients
1160 B(0,0) = b0;
1161 B(1,0) = b1;
1162 for (ordinal_type i=1; i<pc; i++) {
1163 cs[i] = cc[0]*B(i-1,0) + cs[0]*B(i-1,1);
1164 cc[i] = cc[0]*B(i-1+offset,0) + cs[0]*B(i-1+offset,1);
1165 }
1166 }
1167 else {
1168 cs[0] = std::sinh(ca[0]);
1169 cc[0] = std::cosh(ca[0]);
1170 }
1171}
1172
1173template <typename ordinal_type, typename value_type>
1174void
1189
1190template <typename ordinal_type, typename value_type>
1191void
1206
1207template <typename ordinal_type, typename value_type>
1208void
1209Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1210tanh(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& t,
1211 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1212{
1213 if (a.size() > 1) {
1214 OrthogPolyApprox<ordinal_type, value_type, node_type> c(t);
1215 sinhcosh(t, c, a);
1216 divideEqual(t,c);
1217 }
1218 else {
1219 if (t.size() != 1)
1220 t.resize(1);
1221 t[0] = std::tanh(a[0]);
1222 }
1223}
1224
1225template <typename ordinal_type, typename value_type>
1226template <typename OpT>
1227void
1228Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1229quad(const OpT& quad_func,
1230 Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
1231 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a,
1232 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& b)
1233{
1234 const char* func = "Stokhos::DerivOrthogPolyExpansion::quad()";
1235 ordinal_type pa = a.size();
1236 ordinal_type pb = b.size();
1237 ordinal_type pc = sz;
1238 if (c.size() != pc)
1239 c.resize(pc);
1240
1241 const value_type* ca = a.coeff();
1242 const value_type* cb = b.coeff();
1243 value_type* cc = c.coeff();
1244
1245 value_type psi_0 = basis->evaluateZero(0);
1246
1247 // Fill A and B
1248 for (ordinal_type i=1; i<pc; i++) {
1249 B(i-1,0) = 0.0;
1250 for (ordinal_type j=1; j<pc; j++) {
1251 A(i-1,j-1) = 0.0;
1252 for (ordinal_type k=0; k<pb; k++)
1253 A(i-1,j-1) += cb[k]*(*Dijk)(i-1,k,j);
1254 }
1255 for (ordinal_type j=1; j<pa; j++) {
1256 B(i-1,0) += ca[j]*(*Bij)(i-1,j);
1257 }
1258 }
1259
1260 // Solve system
1261 int info = solve(pc-1, 1);
1262
1263 TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
1264 func << ": Argument " << info
1265 << " for solve had illegal value");
1266 TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1267 func << ": Diagonal entry " << info
1268 << " in LU factorization is exactly zero");
1269
1270 // Compute order-0 coefficient
1271 value_type s = psi_0 * ca[0];
1272 value_type t = value_type(0.0);
1273 for (ordinal_type i=1; i<pc; i++) {
1274 s += basis->evaluateZero(i) * ca[i];
1275 t += basis->evaluateZero(i) * B(i-1,0);
1276 }
1277 cc[0] = (quad_func(s) - t) / psi_0;
1278
1279 // Get coefficients
1280 for (ordinal_type i=1; i<pc; i++)
1281 cc[i] = B(i-1,0);
1282}
1283
1284template <typename ordinal_type, typename value_type>
1285void
1286Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1287acos(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
1288 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1289{
1290 if (a.size() > 1) {
1291 times(c,a,a);
1292 minus(c,value_type(1.0),c);
1293 sqrt(c,c);
1294 timesEqual(c,value_type(-1.0));
1295 quad(acos_quad_func(), c, a, c);
1296 }
1297 else {
1298 if (c.size() != 1)
1299 c.resize(1);
1300 c[0] = std::acos(a[0]);
1301 }
1302}
1303
1304template <typename ordinal_type, typename value_type>
1305void
1306Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1307asin(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
1308 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1309{
1310 if (a.size() > 1) {
1311 times(c,a,a);
1312 minus(c,value_type(1.0),c);
1313 sqrt(c,c);
1314 quad(asin_quad_func(), c, a, c);
1315 }
1316 else {
1317 if (c.size() != 1)
1318 c.resize(1);
1319 c[0] = std::asin(a[0]);
1320 }
1321}
1322
1323template <typename ordinal_type, typename value_type>
1324void
1325Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1326atan(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
1327 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1328{
1329 if (a.size() > 1) {
1330 times(c,a,a);
1331 plusEqual(c,value_type(1.0));
1332 quad(atan_quad_func(), c, a, c);
1333 }
1334 else {
1335 if (c.size() != 1)
1336 c.resize(1);
1337 c[0] = std::atan(a[0]);
1338 }
1339}
1340
1341// template <typename ordinal_type, typename value_type>
1342// Hermite<ordinal_type, value_type>
1343// atan2(const Hermite<ordinal_type, value_type>& a,
1344// const Hermite<ordinal_type, value_type>& b)
1345// {
1346// Hermite<ordinal_type, value_type> c = atan(a/b);
1347// c.fastAccessCoeff(0) = atan2(a.coeff(0),b.coeff(0));
1348// }
1349
1350// template <typename ordinal_type, typename value_type>
1351// Hermite<ordinal_type, value_type>
1352// atan2(const T& a,
1353// const Hermite<ordinal_type, value_type>& b)
1354// {
1355// Hermite<ordinal_type, value_type> c = atan(a/b);
1356// c.fastAccessCoeff(0) = atan2(a,b.coeff(0));
1357// }
1358
1359// template <typename ordinal_type, typename value_type>
1360// Hermite<ordinal_type, value_type>
1361// atan2(const Hermite<ordinal_type, value_type>& a,
1362// const T& b)
1363// {
1364// Hermite<ordinal_type, value_type> c = atan(a/b);
1365// c.fastAccessCoeff(0) = atan2(a.coeff(0),b);
1366// }
1367
1368template <typename ordinal_type, typename value_type>
1369void
1370Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1371acosh(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
1372 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1373{
1374 if (a.size() > 1) {
1375 times(c,a,a);
1376 minusEqual(c,value_type(1.0));
1377 sqrt(c,c);
1378 quad(acosh_quad_func(), c, a, c);
1379 }
1380 else {
1381 if (c.size() != 1)
1382 c.resize(1);
1383 c[0] = std::log(a[0]+std::sqrt(a[0]*a[0]-value_type(1.0)));
1384 }
1385}
1386
1387template <typename ordinal_type, typename value_type>
1388void
1389Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1390asinh(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
1391 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1392{
1393 if (a.size() > 1) {
1394 times(c,a,a);
1395 plusEqual(c,value_type(1.0));
1396 sqrt(c,c);
1397 quad(asinh_quad_func(), c, a, c);
1398 }
1399 else {
1400 if (c.size() != 1)
1401 c.resize(1);
1402 c[0] = std::log(a[0]+std::sqrt(a[0]*a[0]+value_type(1.0)));
1403 }
1404}
1405
1406template <typename ordinal_type, typename value_type>
1407void
1408Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1409atanh(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
1410 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1411{
1412 if (a.size() > 1) {
1413 times(c,a,a);
1414 minus(c,value_type(1.0),c);
1415 quad(atanh_quad_func(), c, a, c);
1416 }
1417 else {
1418 if (c.size() != 1)
1419 c.resize(1);
1420 c[0] = 0.5*std::log((value_type(1.0)+a[0])/(value_type(1.0)-a[0]));
1421 }
1422}
1423
1424template <typename ordinal_type, typename value_type>
1425void
1435
1436template <typename ordinal_type, typename value_type>
1437void
1447
1448template <typename ordinal_type, typename value_type>
1449void
1460
1461template <typename ordinal_type, typename value_type>
1462void
1475
1476template <typename ordinal_type, typename value_type>
1477void
1490
1491template <typename ordinal_type, typename value_type>
1492void
1503
1504template <typename ordinal_type, typename value_type>
1505void
1518
1519template <typename ordinal_type, typename value_type>
1520void
1533
1534template <typename ordinal_type, typename value_type>
1535void
1536Stokhos::DerivOrthogPolyExpansion<ordinal_type, value_type>::
1537derivative(Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& c,
1538 const Stokhos::OrthogPolyApprox<ordinal_type, value_type, node_type>& a)
1539{
1540 ordinal_type pc = a.size();
1541
1542 const value_type* ca = a.coeff();
1543 value_type* cc = c.coeff();
1544
1545 for (ordinal_type i=0; i<pc; i++) {
1546 cc[i] = 0.0;
1547 for (ordinal_type j=0; j<pc; j++)
1548 cc[i] += ca[j]*(*Bij)(i,j);
1549 cc[i] /= basis->norm_squared(i);
1550 }
1551}
exp
exp(expr.val())
j
j
Definition Sacado_Fad_Exp_MP_Vector.hpp:527
sqrt
sqrt(expr.val())
val
expr val()
dlange_
double dlange_(char *, int *, int *, double *, int *, double *)
Stokhos_DynamicArrayTraits.hpp
Stokhos::Dense3Tensor
Data structure storing a dense 3-tensor C(i,j,k).
Definition Stokhos_Dense3Tensor.hpp:57
Stokhos::DerivBasis
Abstract base class for multivariate orthogonal polynomials that support computing double and triple ...
Definition Stokhos_DerivBasis.hpp:58
Stokhos::DerivOrthogPolyExpansion::B
Teuchos::SerialDenseMatrix< ordinal_type, value_type > B
RHS.
Definition Stokhos_DerivOrthogPolyExpansion.hpp:267
Stokhos::DerivOrthogPolyExpansion::tan
void tan(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1056
Stokhos::DerivOrthogPolyExpansion::minus
void minus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:383
Stokhos::DerivOrthogPolyExpansion::cbrt
void cbrt(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:844
Stokhos::DerivOrthogPolyExpansion::A
Teuchos::SerialDenseMatrix< ordinal_type, value_type > A
Matrix.
Definition Stokhos_DerivOrthogPolyExpansion.hpp:264
Stokhos::DerivOrthogPolyExpansion::times
void times(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:452
Stokhos::DerivOrthogPolyExpansion::exp
void exp(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:687
Stokhos::DerivOrthogPolyExpansion::basis
Teuchos::RCP< const Stokhos::DerivBasis< ordinal_type, value_type > > basis
Basis.
Definition Stokhos_DerivOrthogPolyExpansion.hpp:249
Stokhos::DerivOrthogPolyExpansion::abs
void abs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1439
Stokhos::DerivOrthogPolyExpansion::Dijk
Teuchos::RCP< const Stokhos::Dense3Tensor< ordinal_type, value_type > > Dijk
Derivative Triple-product tensor.
Definition Stokhos_DerivOrthogPolyExpansion.hpp:258
Stokhos::DerivOrthogPolyExpansion::quad
void quad(const OpT &quad_func, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1229
Stokhos::DerivOrthogPolyExpansion::cosh
void cosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1193
Stokhos::DerivOrthogPolyExpansion::max
void max(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1451
Stokhos::DerivOrthogPolyExpansion::unaryMinus
void unaryMinus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:100
Stokhos::DerivOrthogPolyExpansion::plusEqual
void plusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:118
Stokhos::DerivOrthogPolyExpansion::asin
void asin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1307
Stokhos::DerivOrthogPolyExpansion::Cijk
Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > Cijk
Triple-product tensor.
Definition Stokhos_DerivOrthogPolyExpansion.hpp:255
Stokhos::DerivOrthogPolyExpansion::sinhcosh
void sinhcosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &s, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1074
Stokhos::DerivOrthogPolyExpansion::atan
void atan(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1326
Stokhos::DerivOrthogPolyExpansion::cos
void cos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1039
Stokhos::DerivOrthogPolyExpansion::acos
void acos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1287
Stokhos::DerivOrthogPolyExpansion::atanh
void atanh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1409
Stokhos::DerivOrthogPolyExpansion::derivative
void derivative(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1537
Stokhos::DerivOrthogPolyExpansion::divide
void divide(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:540
Stokhos::DerivOrthogPolyExpansion::minusEqual
void minusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:127
Stokhos::DerivOrthogPolyExpansion::timesEqual
void timesEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:136
Stokhos::DerivOrthogPolyExpansion::min
void min(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1494
Stokhos::DerivOrthogPolyExpansion::log10
void log10(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:809
Stokhos::DerivOrthogPolyExpansion::pow
void pow(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:862
Stokhos::DerivOrthogPolyExpansion::sin
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1022
Stokhos::DerivOrthogPolyExpansion::solve
ordinal_type solve(ordinal_type s, ordinal_type nrhs)
Solve linear system.
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:73
Stokhos::DerivOrthogPolyExpansion::sinh
void sinh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1176
Stokhos::DerivOrthogPolyExpansion::DerivOrthogPolyExpansion
DerivOrthogPolyExpansion(const Teuchos::RCP< const DerivBasis< ordinal_type, value_type > > &basis, const Teuchos::RCP< const Teuchos::SerialDenseMatrix< ordinal_type, value_type > > &Bij, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk, const Teuchos::RCP< const Stokhos::Dense3Tensor< ordinal_type, value_type > > &Dijk)
Constructor.
Stokhos::DerivOrthogPolyExpansion::plus
void plus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:314
Stokhos::DerivOrthogPolyExpansion::divideEqual
void divideEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:148
Stokhos::DerivOrthogPolyExpansion::fabs
void fabs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1427
Stokhos::DerivOrthogPolyExpansion::asinh
void asinh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1390
Stokhos::DerivOrthogPolyExpansion::piv
Teuchos::Array< ordinal_type > piv
Pivot array.
Definition Stokhos_DerivOrthogPolyExpansion.hpp:270
Stokhos::DerivOrthogPolyExpansion::acosh
void acosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1371
Stokhos::DerivOrthogPolyExpansion::sqrt
void sqrt(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:826
Stokhos::DerivOrthogPolyExpansion::Bij
Teuchos::RCP< const Teuchos::SerialDenseMatrix< ordinal_type, value_type > > Bij
Derivative double-product tensor.
Definition Stokhos_DerivOrthogPolyExpansion.hpp:252
Stokhos::DerivOrthogPolyExpansion::tanh
void tanh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:1210
Stokhos::DerivOrthogPolyExpansion::lapack
Teuchos::LAPACK< ordinal_type, value_type > lapack
LAPACK wrappers.
Definition Stokhos_DerivOrthogPolyExpansion.hpp:273
Stokhos::DerivOrthogPolyExpansion::log
void log(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:749
Stokhos::DerivOrthogPolyExpansion::sincos
void sincos(OrthogPolyApprox< ordinal_type, value_type, node_type > &s, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition Stokhos_DerivOrthogPolyExpansionImp.hpp:918
Stokhos::OrthogPolyApprox
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Definition Stokhos_OrthogPolyApprox.hpp:63
Stokhos::OrthogPolyApprox::resize
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
Definition Stokhos_OrthogPolyApproxImp.hpp:152
Stokhos::OrthogPolyApprox::coeff
pointer coeff()
Return coefficient array.
Definition Stokhos_OrthogPolyApproxImp.hpp:168
Stokhos::OrthogPolyApprox::size
ordinal_type size() const
Return size.
Definition Stokhos_OrthogPolyApproxImp.hpp:160
Stokhos::Sparse3Tensor
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Definition Stokhos_Sparse3Tensor.hpp:56
Stokhos::DerivOrthogPolyExpansion::acos_quad_func
Definition Stokhos_DerivOrthogPolyExpansion.hpp:280
Stokhos::DerivOrthogPolyExpansion::acosh_quad_func
Definition Stokhos_DerivOrthogPolyExpansion.hpp:298
Stokhos::DerivOrthogPolyExpansion::asin_quad_func
Definition Stokhos_DerivOrthogPolyExpansion.hpp:286
Stokhos::DerivOrthogPolyExpansion::asinh_quad_func
Definition Stokhos_DerivOrthogPolyExpansion.hpp:304
Stokhos::DerivOrthogPolyExpansion::atan_quad_func
Definition Stokhos_DerivOrthogPolyExpansion.hpp:292
Stokhos::DerivOrthogPolyExpansion::atanh_quad_func
Definition Stokhos_DerivOrthogPolyExpansion.hpp:310
Stokhos::ds_array::get_and_fill
static T * get_and_fill(int sz)
Get memory for new array of length sz and fill with zeros.
Definition Stokhos_DynamicArrayTraits.hpp:82
Generated by doxygen 1.13.2