1 // Copyright (c) 2015-2016, Massachusetts Institute of Technology
2 // Copyright (c) 2016-2017 Sandia Corporation
3 // Copyright (c) 2017 NTESS, LLC.
4
5 // This file is part of the Compressed Continuous Computation (C3) Library
6 // Author: Alex A. Gorodetsky
7 // Contact: alex@alexgorodetsky.com
8
9 // All rights reserved.
10
11 // Redistribution and use in source and binary forms, with or without modification,
12 // are permitted provided that the following conditions are met:
13
14 // 1. Redistributions of source code must retain the above copyright notice,
15 // this list of conditions and the following disclaimer.
16
17 // 2. Redistributions in binary form must reproduce the above copyright notice,
18 // this list of conditions and the following disclaimer in the documentation
19 // and/or other materials provided with the distribution.
20
21 // 3. Neither the name of the copyright holder nor the names of its contributors
22 // may be used to endorse or promote products derived from this software
23 // without specific prior written permission.
24
25 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
28 // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
29 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
30 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
31 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
32 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
33 // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
34 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35
36 //Code
37
38
39 /** \file polynomials.c
40 * Provides routines for manipulating orthogonal polynomials
41 */
42
43 #include <stdlib.h>
44 #include <stdio.h>
45 #include <math.h>
46 #include <string.h>
47 #include <float.h>
48 #include <assert.h>
49 #include <complex.h>
50
51 #include "futil.h"
52
53 #ifndef ZEROTHRESH
54 #define ZEROTHRESH 1e0 * DBL_EPSILON
55 #endif
56
57 #include "stringmanip.h"
58 #include "array.h"
59 #include "fft.h"
60 #include "lib_quadrature.h"
61 #include "linalg.h"
62 #include "legtens.h"
63 #include "fourier.h"
64
fourierortho(size_t n)65 inline static double fourierortho(size_t n){
66 (void) n;
67 return 2.0*M_PI;
68 }
69
70 /********************************************************//**
71 * Initialize a Fourier Basis
72 *
73 * \return p - polynomial
74 *************************************************************/
init_fourier_poly()75 struct OrthPoly * init_fourier_poly(){
76
77 struct OrthPoly * p;
78 if ( NULL == (p = malloc(sizeof(struct OrthPoly)))){
79 fprintf(stderr, "failed to allocate memory for poly exp.\n");
80 exit(1);
81 }
82 p->ptype = FOURIER;
83 p->an = NULL;
84 p->bn = NULL;
85 p->cn = NULL;
86
87 p->lower = 0.0;
88 p->upper = 2.0 * M_PI;
89
90 p->const_term = 0.0;
91 p->lin_coeff = 0.0;
92 p->lin_const = 0.0;
93
94 p->norm = fourierortho;
95
96 return p;
97 }
98
99
100 /********************************************************//**
101 * Evaluate a fourier expansion
102 *
103 * \param[in] poly - polynomial expansion
104 * \param[in] x - location at which to evaluate
105 *
106 * \return out - polynomial value
107 *************************************************************/
fourier_expansion_eval(const struct OrthPolyExpansion * poly,double x)108 double fourier_expansion_eval(const struct OrthPolyExpansion * poly, double x)
109 {
110 assert (poly != NULL);
111 assert (poly->kristoffel_eval == 0);
112
113
114 double x_norm = space_mapping_map(poly->space_transform, x);
115
116 double complex val = cexp((double _Complex)I*x_norm);
117 double complex out = 0.0;
118 double complex basis = 1.0;
119 for (size_t ii = 0; ii < poly->num_poly; ii++){
120 out = out + poly->ccoeff[ii]*basis;
121 basis = basis * val;
122 }
123
124 basis = 1.0 / val;
125 for (size_t ii = 1; ii < poly->num_poly; ii++){
126 out = out + conj(poly->ccoeff[ii])*basis;
127 basis = basis / val;
128 }
129
130 return creal(out);
131 }
132
133 /********************************************************//**
134 * Evaluate the derivative of orth normal polynomial expansion
135 *
136 * \param[in] poly - pointer to orth poly expansion
137 * \param[in] x - location at which to evaluate
138 *
139 *
140 * \return out - value of derivative
141 *************************************************************/
fourier_expansion_deriv_eval(const struct OrthPolyExpansion * poly,double x)142 double fourier_expansion_deriv_eval(const struct OrthPolyExpansion * poly, double x)
143 {
144 assert (poly != NULL);
145 assert (poly->kristoffel_eval == 0);
146
147 double x_norm = space_mapping_map(poly->space_transform,x);
148
149 double complex val = cexp((double _Complex)I*x_norm);
150 double complex out = 0.0;
151 double complex basis = val;
152 for (size_t ii = 1; ii < poly->num_poly; ii++){
153 out = out + poly->ccoeff[ii]*basis * ii * (double _Complex)I;
154 basis = basis * val;
155 }
156
157 basis = 1.0 / val;
158 for (size_t ii = 1; ii < poly->num_poly; ii++){
159 out = out - conj(poly->ccoeff[ii])*basis * ii * (double _Complex)I;
160 basis = basis / val;
161 }
162
163 double rout = creal(out);
164 rout *= space_mapping_map_deriv(poly->space_transform,x);
165 return rout;
166 }
167
168 /********************************************************//**
169 Compute an expansion for the derivtive
170
171 \param[in] p - orthogonal polynomial expansion
172
173 \return derivative
174 *************************************************************/
fourier_expansion_deriv(const struct OrthPolyExpansion * p)175 struct OrthPolyExpansion * fourier_expansion_deriv(const struct OrthPolyExpansion * p)
176 {
177 if (p == NULL) return NULL;
178 assert (p->kristoffel_eval == 0);
179
180 double dx = space_mapping_map_deriv(p->space_transform, 0.0);
181 struct OrthPolyExpansion * out = orth_poly_expansion_copy(p);
182 out->ccoeff[0] = 0.0;
183 for (size_t ii = 1; ii < p->num_poly; ii++){
184 out->ccoeff[ii] *= (double _Complex)I*ii*dx;
185 }
186 return out;
187 }
188
189 /********************************************************//**
190 Compute an expansion for the second derivative
191
192 \param[in] p - orthogonal polynomial expansion
193
194 \return derivative
195 *************************************************************/
fourier_expansion_dderiv(const struct OrthPolyExpansion * p)196 struct OrthPolyExpansion * fourier_expansion_dderiv(const struct OrthPolyExpansion * p)
197 {
198 if (p == NULL) return NULL;
199 assert (p->kristoffel_eval == 0);
200
201 double dx = space_mapping_map_deriv(p->space_transform, 0.0);
202 struct OrthPolyExpansion * out = orth_poly_expansion_copy(p);
203 out->ccoeff[0] = 0.0;
204 for (size_t ii = 1; ii < p->num_poly; ii++){
205 out->ccoeff[ii] *= -1.0*(double)(ii*ii)*dx*dx;
206 }
207 return out;
208 }
209
210 /********************************************************//**
211 * Approximating a function that can take a vector of points as
212 * input
213 *
214 * \param[in,out] poly - orthogonal polynomial expansion
215 * \param[in] f - wrapped function
216 * \param[in] opts - approximation options
217 *
218 * \return 0 - no problems, > 0 problem
219 *
220 * \note Maximum quadrature limited to 200 nodes
221 *************************************************************/
222 int
fourier_expansion_approx_vec(struct OrthPolyExpansion * poly,struct Fwrap * f,const struct OpeOpts * opts)223 fourier_expansion_approx_vec(struct OrthPolyExpansion * poly,
224 struct Fwrap * f,
225 const struct OpeOpts * opts)
226 {
227 (void) opts;
228 int return_val = 1;
229 /* printf("I am here!\n"); */
230
231
232 /* size_t N = poly->num_poly; */
233 size_t N = (poly->num_poly-1)*2;
234 /* size_t N = 8; */
235 size_t nquad = N;
236 double frac = 2.0*M_PI/(nquad);
237
238
239 double * quad_pts = calloc_double(N);
240 for (size_t ii = 0; ii < nquad; ii++){
241 quad_pts[ii] = frac * ii;
242 }
243
244 /* printf("what!\n"); */
245 double * pts = calloc_double(nquad);
246 double * fvals = calloc_double(nquad);
247 for (size_t ii = 0; ii < nquad; ii++){
248 pts[ii] = space_mapping_map_inverse(poly->space_transform,
249 quad_pts[ii]);
250 /* pts[ii] = -M_PI + ii * frac; */
251 }
252
253
254 // Evaluate functions
255
256 return_val = fwrap_eval(nquad,pts,fvals,f);
257 if (return_val != 0){
258 return return_val;
259 }
260
261 /* printf("pts = "); dprint(nquad, pts); */
262 /* printf("fvals = "); dprint(nquad, fvals); */
263
264 double complex * coeff = malloc(nquad * sizeof(double complex));
265 double complex * fvc = malloc(nquad * sizeof(double complex));
266 for (size_t ii = 0; ii < nquad; ii++){
267 fvc[ii] = fvals[ii];
268 }
269
270 int fft_res = fft(N, fvc, 1, coeff, 1);
271 if (fft_res != 0){
272 fprintf(stderr, "fft is not successfull\n");
273 return 1;
274 }
275
276 /* printf("\n"); */
277 /* for (size_t ii = 0; ii < N; ii++){ */
278 /* fprintf(stdout, "%3.5G %3.5G\n", creal(coeff[ii]), cimag(coeff[ii])); */
279 /* } */
280
281
282 /* printf("copy coeffs %zu= \n", poly->num_poly); */
283 /* poly->num_poly = N/2; */
284 poly->ccoeff[0] = coeff[0]/N;
285 for (size_t ii = 1; ii < poly->num_poly; ii++){
286 poly->ccoeff[ii] = coeff[ii]/N;
287 /* printf("ii = %zu %3.5G %3.5G\n", ii, creal(poly->ccoeff[ii]), */
288 /* cimag(poly->ccoeff[ii])); */
289 }
290
291 /* poly->num_poly = N/2; */
292 /* for (size_t ii = 0; ii < poly->num_poly/2+1; ii++){ */
293 /* /\* printf("ii = %zu\n", ii); *\/ */
294 /* poly->ccoeff[ii] = coeff[ii]/N; */
295 /* } */
296 /* for (size_t ii = 1; ii < poly->num_poly/2; ii++){ */
297 /* poly->ccoeff[ii+poly->num_poly/2] = coeff[poly->num_poly-ii]/N; */
298 /* } */
299
300 free(pts); pts = NULL;
301 free(quad_pts); quad_pts = NULL;
302 free(fvals); fvals = NULL;
303 free(fvc); fvc = NULL;
304 free(coeff); coeff = NULL;
305 return return_val;
306 }
307
308 /********************************************************//**
309 * Integrate a fourier expansion
310 *
311 * \param[in] poly - polynomial to integrate
312 *
313 * \return out - integral of approximation
314 *************************************************************/
fourier_integrate(const struct OrthPolyExpansion * poly)315 double fourier_integrate(const struct OrthPolyExpansion * poly)
316 {
317 double out = 0.0;
318
319 double m = space_mapping_map_inverse_deriv(poly->space_transform, 0.0);
320 out = creal(poly->ccoeff[0] * 2.0 * M_PI);
321 out = out * m;
322 return out;
323 }
324
325 typedef struct Pair
326 {
327 const struct OrthPolyExpansion * a;
328 const struct OrthPolyExpansion * b;
329 } pair_t;
330
prod_eval(size_t n,const double * x,double * out,void * args)331 static int prod_eval(size_t n, const double * x, double * out, void * args)
332 {
333 pair_t * pairs = args;
334 double e1, e2;
335 for (size_t ii = 0; ii < n; ii++){
336 e1 = orth_poly_expansion_eval(pairs->a, x[ii]);
337 e2 = orth_poly_expansion_eval(pairs->b, x[ii]);
338 out[ii] = e1 * e2;
339 }
340 return 0;
341 }
342
343 /********************************************************//**
344 Multiply two fourier series'
345
346 \param[in] a - first fs
347 \param[in] b - second fs
348
349 \return c - product
350 *************************************************************/
351 struct OrthPolyExpansion *
fourier_expansion_prod(const struct OrthPolyExpansion * a,const struct OrthPolyExpansion * b)352 fourier_expansion_prod(const struct OrthPolyExpansion * a,
353 const struct OrthPolyExpansion * b)
354 {
355
356 /* (void) a; */
357 /* (void) b; */
358 /* fprintf(stderr, "HAVE NOT FINISHED IMPLEMENTING PRODUCT of fourier\n"); */
359 /* exit(1); */
360 double lb = a->lower_bound;
361 double ub = a->upper_bound;
362
363 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
364 ope_opts_set_lb(opts, lb);
365 ope_opts_set_ub(opts, ub);
366 size_t n = a->num_poly-1;
367 /* ope_opts_set_start(opts, 2*n+1); */
368 ope_opts_set_start(opts, n+1);
369
370 /* printf("n = %zu\n", n); */
371 /* struct OrthPolyExpansion * fourier = */
372 /* orth_poly_expansion_init_from_opts(opts, 2*n+1); */
373 struct OrthPolyExpansion * fourier =
374 orth_poly_expansion_init_from_opts(opts, n+1);
375
376 if (n > 0){
377 struct Fwrap * fw = fwrap_create(1, "general-vec");
378 pair_t pairs = {a, b};
379 fwrap_set_fvec(fw, prod_eval, &pairs);
380 int res = orth_poly_expansion_approx_vec(fourier, fw, opts);
381 assert (res == 0);
382 fwrap_destroy(fw); fw = NULL;
383 }
384 ope_opts_free(opts); opts = NULL;
385 return fourier;
386 }
387
388 /********************************************************//**
389 * Multiply by scalar and add two expansions
390 *
391 * \param[in] a - scaling factor for first polynomial
392 * \param[in] x - first polynomial
393 * \param[in] y - second polynomial
394 *
395 * \return 0 if successfull 1 if error with allocating more space for y
396 *
397 * \note
398 * Computes z=ax+by, where x and y are polynomial expansionx
399 * Requires both polynomials to have the same upper
400 * and lower bounds
401 *
402 **************************************************************/
fourier_expansion_axpy(double a,const struct OrthPolyExpansion * x,struct OrthPolyExpansion * y)403 int fourier_expansion_axpy(double a, const struct OrthPolyExpansion * x,
404 struct OrthPolyExpansion * y)
405 {
406
407 if (x->num_poly <= y->num_poly){
408 for (size_t ii = 0; ii < x->num_poly; ii++){
409 y->ccoeff[ii] += a * x->ccoeff[ii];
410 if (cabs(y->ccoeff[ii]) < ZEROTHRESH){
411 y->ccoeff[ii] = 0.0;
412 }
413 }
414 }
415 else{
416 if (x->num_poly > y->nalloc){
417 //printf("hereee\n");
418 y->nalloc = x->num_poly+10;
419 double complex * temp = realloc(y->ccoeff, (y->nalloc)*sizeof(double complex));
420 if (temp == NULL){
421 return 0;
422 }
423 else{
424 y->ccoeff = temp;
425 for (size_t ii = y->num_poly; ii < y->nalloc; ii++){
426 y->ccoeff[ii] = 0.0;
427 }
428 }
429 //printf("finished\n");
430 }
431 for (size_t ii = y->num_poly; ii < x->num_poly; ii++){
432 y->ccoeff[ii] = a * x->ccoeff[ii];
433 if (cabs(y->ccoeff[ii]) < ZEROTHRESH){
434 y->ccoeff[ii] = 0.0;
435 }
436 }
437 for (size_t ii = 0; ii < y->num_poly; ii++){
438 y->ccoeff[ii] += a * x->ccoeff[ii];
439 if (cabs(y->ccoeff[ii]) < ZEROTHRESH){
440 y->ccoeff[ii] = 0.0;
441 }
442 }
443 y->num_poly = x->num_poly;
444 size_t nround = y->num_poly;
445 for (size_t ii = 0; ii < y->num_poly-1;ii++){
446 if (cabs(y->ccoeff[y->num_poly-1-ii]) > ZEROTHRESH){
447 break;
448 }
449 else{
450 nround = nround-1;
451 }
452 }
453 y->num_poly = nround;
454 }
455
456 return 0;
457 }
458
459
460 /* /\********************************************************\//\** */
461 /* * Multiply by scalar and add two orthgonal */
462 /* * expansions of the same family together */
463 /* * */
464 /* * \param[in] a - scaling factor for first polynomial */
465 /* * \param[in] x - first polynomial */
466 /* * \param[in] b - scaling factor for second polynomial */
467 /* * \param[in] y second polynomial */
468 /* * */
469 /* * \return p - orthogonal polynomial expansion */
470 /* * */
471 /* * \note */
472 /* * Computes z=ax+by, where x and y are polynomial expansionx */
473 /* * Requires both polynomials to have the same upper */
474 /* * and lower bounds */
475 /* * */
476 /* *************************************************************\/ */
477 /* struct OrthPolyExpansion * */
478 /* orth_poly_expansion_daxpby(double a, struct OrthPolyExpansion * x, */
479 /* double b, struct OrthPolyExpansion * y) */
480 /* { */
481 /* /\* */
482 /* printf("a=%G b=%G \n",a,b); */
483 /* printf("x=\n"); */
484 /* print_orth_poly_expansion(x,0,NULL); */
485 /* printf("y=\n"); */
486 /* print_orth_poly_expansion(y,0,NULL); */
487 /* *\/ */
488
489 /* //double diffa = cabs(a-ZEROTHRESH); */
490 /* //double diffb = cabs(b-ZEROTHRESH); */
491 /* size_t ii; */
492 /* struct OrthPolyExpansion * p ; */
493 /* //if ( (x == NULL && y != NULL) || ((diffa <= ZEROTHRESH) && (y != NULL))){ */
494 /* if ( (x == NULL && y != NULL)){ */
495 /* //printf("b = %G\n",b); */
496 /* //if (diffb <= ZEROTHRESH){ */
497 /* // p = orth_poly_expansion_init(y->p->ptype,1,y->lower_bound, y->upper_bound); */
498 /* // } */
499 /* // else{ */
500 /* p = orth_poly_expansion_init(y->p->ptype, */
501 /* y->num_poly, y->lower_bound, y->upper_bound); */
502 /* space_mapping_free(p->space_transform); */
503 /* p->space_transform = space_mapping_copy(y->space_transform); */
504 /* for (ii = 0; ii < y->num_poly; ii++){ */
505 /* p->coeff[ii] = y->coeff[ii] * b; */
506 /* } */
507 /* //} */
508 /* orth_poly_expansion_round(&p); */
509 /* return p; */
510 /* } */
511 /* //if ( (y == NULL && x != NULL) || ((diffb <= ZEROTHRESH) && (x != NULL))){ */
512 /* if ( (y == NULL && x != NULL)){ */
513 /* //if (a <= ZEROTHRESH){ */
514 /* // p = orth_poly_expansion_init(x->p->ptype,1, x->lower_bound, x->upper_bound); */
515 /* // } */
516 /* //else{ */
517 /* p = orth_poly_expansion_init(x->p->ptype, */
518 /* x->num_poly, x->lower_bound, x->upper_bound); */
519 /* space_mapping_free(p->space_transform); */
520 /* p->space_transform = space_mapping_copy(x->space_transform); */
521 /* for (ii = 0; ii < x->num_poly; ii++){ */
522 /* p->coeff[ii] = x->coeff[ii] * a; */
523 /* } */
524 /* //} */
525 /* orth_poly_expansion_round(&p); */
526 /* return p; */
527 /* } */
528
529 /* size_t N = x->num_poly > y->num_poly ? x->num_poly : y->num_poly; */
530
531 /* p = orth_poly_expansion_init(x->p->ptype, */
532 /* N, x->lower_bound, x->upper_bound); */
533 /* space_mapping_free(p->space_transform); */
534 /* p->space_transform = space_mapping_copy(x->space_transform); */
535
536 /* size_t xN = x->num_poly; */
537 /* size_t yN = y->num_poly; */
538
539 /* //printf("diffa = %G, x==NULL %d\n",diffa,x==NULL); */
540 /* //printf("diffb = %G, y==NULL %d\n",diffb,y==NULL); */
541 /* // assert(diffa > ZEROTHRESH); */
542 /* // assert(diffb > ZEROTHRESH); */
543 /* if (xN > yN){ */
544 /* for (ii = 0; ii < yN; ii++){ */
545 /* p->coeff[ii] = x->coeff[ii]*a + y->coeff[ii]*b; */
546 /* //if ( cabs(p->coeff[ii]) < ZEROTHRESH){ */
547 /* // p->coeff[ii] = 0.0; */
548 /* // } */
549 /* } */
550 /* for (ii = yN; ii < xN; ii++){ */
551 /* p->coeff[ii] = x->coeff[ii]*a; */
552 /* //if ( cabs(p->coeff[ii]) < ZEROTHRESH){ */
553 /* // p->coeff[ii] = 0.0; */
554 /* // } */
555 /* } */
556 /* } */
557 /* else{ */
558 /* for (ii = 0; ii < xN; ii++){ */
559 /* p->coeff[ii] = x->coeff[ii]*a + y->coeff[ii]*b; */
560 /* //if ( cabs(p->coeff[ii]) < ZEROTHRESH){ */
561 /* // p->coeff[ii] = 0.0; */
562 /* // } */
563 /* } */
564 /* for (ii = xN; ii < yN; ii++){ */
565 /* p->coeff[ii] = y->coeff[ii]*b; */
566 /* //if ( cabs(p->coeff[ii]) < ZEROTHRESH){ */
567 /* // p->coeff[ii] = 0.0; */
568 /* //} */
569 /* } */
570 /* } */
571
572 /* orth_poly_expansion_round(&p); */
573 /* return p; */
574 /* } */
575
576 /* //////////////////////////////////////////////////////////////////////////// */
577 /* // Algorithms */
578
579 /* /\********************************************************\//\** */
580 /* * Obtain the real roots of a standard polynomial */
581 /* * */
582 /* * \param[in] p - standard polynomial */
583 /* * \param[in,out] nkeep - returns how many real roots tehre are */
584 /* * */
585 /* * \return real_roots - real roots of a standard polynomial */
586 /* * */
587 /* * \note */
588 /* * Only roots within the bounds are returned */
589 /* *************************************************************\/ */
590 /* double * */
591 /* standard_poly_real_roots(struct StandardPoly * p, size_t * nkeep) */
592 /* { */
593 /* if (p->num_poly == 1) // constant function */
594 /* { */
595 /* double * real_roots = NULL; */
596 /* *nkeep = 0; */
597 /* return real_roots; */
598 /* } */
599 /* else if (p->num_poly == 2){ // linear */
600 /* double root = -p->coeff[0] / p->coeff[1]; */
601
602 /* if ((root > p->lower_bound) && (root < p->upper_bound)){ */
603 /* *nkeep = 1; */
604 /* } */
605 /* else{ */
606 /* *nkeep = 0; */
607 /* } */
608 /* double * real_roots = NULL; */
609 /* if (*nkeep == 1){ */
610 /* real_roots = calloc_double(1); */
611 /* real_roots[0] = root; */
612 /* } */
613 /* return real_roots; */
614 /* } */
615
616 /* size_t nrows = p->num_poly-1; */
617 /* //printf("coeffs = \n"); */
618 /* //dprint(p->num_poly, p->coeff); */
619 /* while (cabs(p->coeff[nrows]) < ZEROTHRESH ){ */
620 /* //while (cabs(p->coeff[nrows]) < DBL_MIN){ */
621 /* nrows--; */
622 /* if (nrows == 1){ */
623 /* break; */
624 /* } */
625 /* } */
626
627 /* //printf("nrows left = %zu \n", nrows); */
628 /* if (nrows == 1) // linear */
629 /* { */
630 /* double root = -p->coeff[0] / p->coeff[1]; */
631 /* if ((root > p->lower_bound) && (root < p->upper_bound)){ */
632 /* *nkeep = 1; */
633 /* } */
634 /* else{ */
635 /* *nkeep = 0; */
636 /* } */
637 /* double * real_roots = NULL; */
638 /* if (*nkeep == 1){ */
639 /* real_roots = calloc_double(1); */
640 /* real_roots[0] = root; */
641 /* } */
642 /* return real_roots; */
643 /* } */
644 /* else if (nrows == 0) */
645 /* { */
646 /* double * real_roots = NULL; */
647 /* *nkeep = 0; */
648 /* return real_roots; */
649 /* } */
650
651 /* // transpose of the companion matrix */
652 /* double * t_companion = calloc_double((p->num_poly-1)*(p->num_poly-1)); */
653 /* size_t ii; */
654
655
656 /* // size_t m = nrows; */
657 /* t_companion[nrows-1] = -p->coeff[0]/p->coeff[nrows]; */
658 /* for (ii = 1; ii < nrows; ii++){ */
659 /* t_companion[ii * nrows + ii-1] = 1.0; */
660 /* t_companion[ii * nrows + nrows-1] = -p->coeff[ii]/p->coeff[nrows]; */
661 /* } */
662 /* double * real = calloc_double(nrows); */
663 /* double * img = calloc_double(nrows); */
664 /* int info; */
665 /* int lwork = 8 * nrows; */
666 /* double * iwork = calloc_double(8 * nrows); */
667 /* //double * vl; */
668 /* //double * vr; */
669 /* int n = nrows; */
670
671 /* //printf("hello! n=%d \n",n); */
672 /* dgeev_("N","N", &n, t_companion, &n, real, img, NULL, &n, */
673 /* NULL, &n, iwork, &lwork, &info); */
674
675 /* //printf("info = %d", info); */
676
677 /* free (iwork); */
678
679 /* int * keep = calloc_int(nrows); */
680 /* *nkeep = 0; */
681 /* // the 1e-10 is kinda hacky */
682 /* for (ii = 0; ii < nrows; ii++){ */
683 /* //printf("real[ii] - p->lower_bound = %G\n",real[ii]-p->lower_bound); */
684 /* //printf("real root = %3.15G, imag = %G \n",real[ii],img[ii]); */
685 /* //printf("lower thresh = %3.20G\n",p->lower_bound-1e-8); */
686 /* //printf("zero thresh = %3.20G\n",1e-8); */
687 /* //printf("upper thresh = %G\n",p->upper_bound+ZEROTHRESH); */
688 /* //printf("too low? %d \n", real[ii] < (p->lower_bound-1e-8)); */
689 /* if ((cabs(img[ii]) < 1e-7) && */
690 /* (real[ii] > (p->lower_bound-1e-8)) && */
691 /* //(real[ii] >= (p->lower_bound-1e-7)) && */
692 /* (real[ii] < (p->upper_bound+1e-8))) { */
693 /* //(real[ii] <= (p->upper_bound+1e-7))) { */
694
695 /* //\* */
696 /* if (real[ii] < p->lower_bound){ */
697 /* real[ii] = p->lower_bound; */
698 /* } */
699 /* if (real[ii] > p->upper_bound){ */
700 /* real[ii] = p->upper_bound; */
701 /* } */
702 /* //\*\/ */
703
704 /* keep[ii] = 1; */
705 /* *nkeep = *nkeep + 1; */
706 /* //printf("keep\n"); */
707 /* } */
708 /* else{ */
709 /* keep[ii] = 0; */
710 /* } */
711 /* } */
712
713 /* /\* */
714 /* printf("real portions roots = "); */
715 /* dprint(nrows, real); */
716 /* printf("imag portions roots = "); */
717 /* for (ii = 0; ii < nrows; ii++) printf("%E ",img[ii]); */
718 /* printf("\n"); */
719 /* //dprint(nrows, img); */
720 /* *\/ */
721
722 /* double * real_roots = calloc_double(*nkeep); */
723 /* size_t counter = 0; */
724 /* for (ii = 0; ii < nrows; ii++){ */
725 /* if (keep[ii] == 1){ */
726 /* real_roots[counter] = real[ii]; */
727 /* counter++; */
728 /* } */
729
730 /* } */
731
732 /* free(t_companion); */
733 /* free(real); */
734 /* free(img); */
735 /* free(keep); */
736
737 /* return real_roots; */
738 /* } */
739
740 /* static int dblcompare(const void * a, const void * b) */
741 /* { */
742 /* const double * aa = a; */
743 /* const double * bb = b; */
744 /* if ( *aa < *bb){ */
745 /* return -1; */
746 /* } */
747 /* return 1; */
748 /* } */
749
750 /* /\********************************************************\//\** */
751 /* * Obtain the real roots of a legendre polynomial expansion */
752 /* * */
753 /* * \param[in] p - orthogonal polynomial expansion */
754 /* * \param[in,out] nkeep - returns how many real roots tehre are */
755 /* * */
756 /* * \return real_roots - real roots of an orthonormal polynomial expansion */
757 /* * */
758 /* * \note */
759 /* * Only roots within the bounds are returned */
760 /* * Algorithm is based on eigenvalues of non-standard companion matrix from */
761 /* * Roots of Polynomials Expressed in terms of orthogonal polynomials */
762 /* * David Day and Louis Romero 2005 */
763 /* * */
764 /* * Multiplying by a factor of sqrt(2*N+1) because using orthonormal, */
765 /* * rather than orthogonal polynomials */
766 /* *************************************************************\/ */
767 /* double * */
768 /* legendre_expansion_real_roots(struct OrthPolyExpansion * p, size_t * nkeep) */
769 /* { */
770
771 /* double * real_roots = NULL; // output */
772 /* *nkeep = 0; */
773
774 /* double m = (p->upper_bound - p->lower_bound) / */
775 /* (p->p->upper - p->p->lower); */
776 /* double off = p->upper_bound - m * p->p->upper; */
777
778 /* orth_poly_expansion_round(&p); */
779 /* // print_orth_poly_expansion(p,3,NULL); */
780 /* //printf("last 2 = %G\n",p->coeff[p->num_poly-1]); */
781 /* size_t N = p->num_poly-1; */
782 /* //printf("N = %zu\n",N); */
783 /* if (N == 0){ */
784 /* return real_roots; */
785 /* } */
786 /* else if (N == 1){ */
787 /* if (cabs(p->coeff[N]) <= ZEROTHRESH){ */
788 /* return real_roots; */
789 /* } */
790 /* else{ */
791 /* double root = -p->coeff[0] / p->coeff[1]; */
792 /* if ( (root >= -1.0-ZEROTHRESH) && (root <= 1.0 - ZEROTHRESH)){ */
793 /* if (root <-1.0){ */
794 /* root = -1.0; */
795 /* } */
796 /* else if (root > 1.0){ */
797 /* root = 1.0; */
798 /* } */
799 /* *nkeep = 1; */
800 /* real_roots = calloc_double(1); */
801 /* real_roots[0] = m*root+off; */
802 /* } */
803 /* } */
804 /* } */
805 /* else{ */
806 /* /\* printf("I am here\n"); *\/ */
807 /* double * nscompanion = calloc_double(N*N); // nonstandard companion */
808 /* size_t ii; */
809 /* double hnn1 = - (double) (N) / (2.0 * (double) (N) - 1.0); */
810 /* /\* double hnn1 = - 1.0 / p->p->an(N); *\/ */
811 /* nscompanion[1] = 1.0; */
812 /* /\* nscompanion[(N-1)*N] += hnn1 * p->coeff[0] / p->coeff[N]; *\/ */
813 /* nscompanion[(N-1)*N] += hnn1 * p->coeff[0] / (p->coeff[N] * sqrt(2*N+1)); */
814 /* for (ii = 1; ii < N-1; ii++){ */
815 /* assert (cabs(p->p->bn(ii)) < 1e-14); */
816 /* double in = (double) ii; */
817 /* nscompanion[ii*N+ii-1] = in / ( 2.0 * in + 1.0); */
818 /* nscompanion[ii*N+ii+1] = (in + 1.0) / ( 2.0 * in + 1.0); */
819
820 /* /\* nscompanion[(N-1)*N + ii] += hnn1 * p->coeff[ii] / p->coeff[N]; *\/ */
821 /* nscompanion[(N-1)*N + ii] += hnn1 * p->coeff[ii] * sqrt(2*ii+1)/ p->coeff[N] / sqrt(2*N+1); */
822 /* } */
823 /* nscompanion[N*N-2] += (double) (N-1) / (2.0 * (double) (N-1) + 1.0); */
824
825
826 /* /\* nscompanion[N*N-1] += hnn1 * p->coeff[N-1] / p->coeff[N]; *\/ */
827 /* nscompanion[N*N-1] += hnn1 * p->coeff[N-1] * sqrt(2*(N-1)+1)/ p->coeff[N] / sqrt(2*N+1); */
828
829 /* //printf("good up to here!\n"); */
830 /* //dprint2d_col(N,N,nscompanion); */
831
832 /* int info; */
833 /* double * scale = calloc_double(N); */
834 /* //\* */
835 /* //Balance */
836 /* int ILO, IHI; */
837 /* //printf("am I here? N=%zu \n",N); */
838 /* //dprint(N*N,nscompanion); */
839 /* dgebal_("S", (int*)&N, nscompanion, (int *)&N,&ILO,&IHI,scale,&info); */
840 /* //printf("yep\n"); */
841 /* if (info < 0){ */
842 /* fprintf(stderr, "Calling dgebl had error in %d-th input in the legendre_expansion_real_roots function\n",info); */
843 /* exit(1); */
844 /* } */
845
846 /* //printf("balanced!\n"); */
847 /* //dprint2d_col(N,N,nscompanion); */
848
849 /* //IHI = M1; */
850 /* //printf("M1=%zu\n",M1); */
851 /* //printf("ilo=%zu\n",ILO); */
852 /* //printf("IHI=%zu\n",IHI); */
853 /* //\*\/ */
854
855 /* double * real = calloc_double(N); */
856 /* double * img = calloc_double(N); */
857 /* //printf("allocated eigs N = %zu\n",N); */
858 /* int lwork = 8 * (int)N; */
859 /* //printf("got lwork\n"); */
860 /* double * iwork = calloc_double(8*N); */
861 /* //printf("go here"); */
862
863 /* //dgeev_("N","N", &N, nscompanion, &N, real, img, NULL, &N, */
864 /* // NULL, &N, iwork, &lwork, &info); */
865 /* dhseqr_("E","N",(int*)&N,&ILO,&IHI,nscompanion,(int*)&N,real,img,NULL,(int*)&N,iwork,&lwork,&info); */
866 /* //printf("done here"); */
867
868 /* if (info < 0){ */
869 /* fprintf(stderr, "Calling dhesqr had error in %d-th input in the legendre_expansion_real_roots function\n",info); */
870 /* exit(1); */
871 /* } */
872 /* else if(info > 0){ */
873 /* //fprintf(stderr, "Eigenvalues are still uncovered in legendre_expansion_real_roots function\n"); */
874 /* // printf("coeffs are \n"); */
875 /* // dprint(p->num_poly, p->coeff); */
876 /* // printf("last 2 = %G\n",p->coeff[p->num_poly-1]); */
877 /* // exit(1); */
878 /* } */
879
880 /* // printf("eigenvalues \n"); */
881 /* size_t * keep = calloc_size_t(N); */
882 /* for (ii = 0; ii < N; ii++){ */
883 /* /\* printf("(%3.15G, %3.15G)\n",real[ii],img[ii]); *\/ */
884 /* if ((cabs(img[ii]) < 1e-6) && (real[ii] > -1.0-1e-12) && (real[ii] < 1.0+1e-12)){ */
885 /* if (real[ii] < -1.0){ */
886 /* real[ii] = -1.0; */
887 /* } */
888 /* else if (real[ii] > 1.0){ */
889 /* real[ii] = 1.0; */
890 /* } */
891 /* keep[ii] = 1; */
892 /* *nkeep = *nkeep + 1; */
893 /* } */
894 /* } */
895
896
897 /* if (*nkeep > 0){ */
898 /* real_roots = calloc_double(*nkeep); */
899 /* size_t counter = 0; */
900 /* for (ii = 0; ii < N; ii++){ */
901 /* if (keep[ii] == 1){ */
902 /* real_roots[counter] = real[ii]*m+off; */
903 /* counter++; */
904 /* } */
905 /* } */
906 /* } */
907
908
909 /* free(keep); keep = NULL; */
910 /* free(iwork); iwork = NULL; */
911 /* free(real); real = NULL; */
912 /* free(img); img = NULL; */
913 /* free(nscompanion); nscompanion = NULL; */
914 /* free(scale); scale = NULL; */
915 /* } */
916
917 /* if (*nkeep > 1){ */
918 /* qsort(real_roots, *nkeep, sizeof(double), dblcompare); */
919 /* } */
920 /* return real_roots; */
921 /* } */
922
923 /* /\********************************************************\//\** */
924 /* * Obtain the real roots of a chebyshev polynomial expansion */
925 /* * */
926 /* * \param[in] p - orthogonal polynomial expansion */
927 /* * \param[in,out] nkeep - returns how many real roots tehre are */
928 /* * */
929 /* * \return real_roots - real roots of an orthonormal polynomial expansion */
930 /* * */
931 /* * \note */
932 /* * Only roots within the bounds are returned */
933 /* * Algorithm is based on eigenvalues of non-standard companion matrix from */
934 /* * Roots of Polynomials Expressed in terms of orthogonal polynomials */
935 /* * David Day and Louis Romero 2005 */
936 /* * */
937 /* * Multiplying by a factor of sqrt(2*N+1) because using orthonormal, */
938 /* * rather than orthogonal polynomials */
939 /* *************************************************************\/ */
940 /* double * */
941 /* chebyshev_expansion_real_roots(struct OrthPolyExpansion * p, size_t * nkeep) */
942 /* { */
943 /* /\* fprintf(stderr, "Chebyshev real_roots not finished yet\n"); *\/ */
944 /* /\* exit(1); *\/ */
945 /* double * real_roots = NULL; // output */
946 /* *nkeep = 0; */
947
948 /* double m = (p->upper_bound - p->lower_bound) / (p->p->upper - p->p->lower); */
949 /* double off = p->upper_bound - m * p->p->upper; */
950
951
952 /* /\* printf("coeff pre truncate = "); dprint(p->num_, p->coeff); *\/ */
953 /* /\* for (size_t ii = 0; ii < p->num_poly; ii++){ *\/ */
954 /* /\* if (cabs(p->coeff[ii]) < 1e-13){ *\/ */
955 /* /\* p->coeff[ii] = 0.0; *\/ */
956 /* /\* } *\/ */
957 /* /\* } *\/ */
958 /* orth_poly_expansion_round(&p); */
959
960 /* size_t N = p->num_poly-1; */
961 /* if (N == 0){ */
962 /* return real_roots; */
963 /* } */
964 /* else if (N == 1){ */
965 /* if (cabs(p->coeff[N]) <= ZEROTHRESH){ */
966 /* return real_roots; */
967 /* } */
968 /* else{ */
969 /* double root = -p->coeff[0] / p->coeff[1]; */
970 /* if ( (root >= -1.0-ZEROTHRESH) && (root <= 1.0 - ZEROTHRESH)){ */
971 /* if (root <-1.0){ */
972 /* root = -1.0; */
973 /* } */
974 /* else if (root > 1.0){ */
975 /* root = 1.0; */
976 /* } */
977 /* *nkeep = 1; */
978 /* real_roots = calloc_double(1); */
979 /* real_roots[0] = m*root+off; */
980 /* } */
981 /* } */
982 /* } */
983 /* else{ */
984 /* /\* printf("I am heare\n"); *\/ */
985 /* /\* dprint(N+1, p->coeff); *\/ */
986 /* double * nscompanion = calloc_double(N*N); // nonstandard companion */
987 /* size_t ii; */
988
989 /* double hnn1 = 0.5; */
990 /* double gamma = p->coeff[N]; */
991
992 /* nscompanion[1] = 1.0; */
993 /* nscompanion[(N-1)*N] -= hnn1*p->coeff[0] / gamma; */
994 /* for (ii = 1; ii < N-1; ii++){ */
995 /* assert (cabs(p->p->bn(ii)) < 1e-14); */
996
997 /* nscompanion[ii*N+ii-1] = 0.5; // ii-th column */
998 /* nscompanion[ii*N+ii+1] = 0.5; */
999
1000 /* // update last column */
1001 /* nscompanion[(N-1)*N + ii] -= hnn1 * p->coeff[ii] / gamma; */
1002 /* } */
1003 /* nscompanion[N*N-2] += 0.5; */
1004 /* nscompanion[N*N-1] -= hnn1 * p->coeff[N-1] / gamma; */
1005
1006 /* //printf("good up to here!\n"); */
1007 /* /\* dprint2d_col(N,N,nscompanion); *\/ */
1008
1009 /* int info; */
1010 /* double * scale = calloc_double(N); */
1011 /* //\* */
1012 /* //Balance */
1013 /* int ILO, IHI; */
1014 /* //printf("am I here? N=%zu \n",N); */
1015 /* //dprint(N*N,nscompanion); */
1016 /* dgebal_("S", (int*)&N, nscompanion, (int *)&N,&ILO,&IHI,scale,&info); */
1017 /* //printf("yep\n"); */
1018 /* if (info < 0){ */
1019 /* fprintf(stderr, "Calling dgebl had error in %d-th input in the chebyshev_expansion_real_roots function\n",info); */
1020 /* exit(1); */
1021 /* } */
1022
1023 /* //printf("balanced!\n"); */
1024 /* //dprint2d_col(N,N,nscompanion); */
1025
1026 /* //IHI = M1; */
1027 /* //printf("M1=%zu\n",M1); */
1028 /* //printf("ilo=%zu\n",ILO); */
1029 /* //printf("IHI=%zu\n",IHI); */
1030 /* //\*\/ */
1031
1032 /* double * real = calloc_double(N); */
1033 /* double * img = calloc_double(N); */
1034 /* //printf("allocated eigs N = %zu\n",N); */
1035 /* int lwork = 8 * (int)N; */
1036 /* //printf("got lwork\n"); */
1037 /* double * iwork = calloc_double(8*N); */
1038 /* //printf("go here"); */
1039
1040 /* //dgeev_("N","N", &N, nscompanion, &N, real, img, NULL, &N, */
1041 /* // NULL, &N, iwork, &lwork, &info); */
1042 /* dhseqr_("E","N",(int*)&N,&ILO,&IHI,nscompanion,(int*)&N,real,img,NULL,(int*)&N,iwork,&lwork,&info); */
1043 /* //printf("done here"); */
1044
1045 /* if (info < 0){ */
1046 /* fprintf(stderr, "Calling dhesqr had error in %d-th input in the legendre_expansion_real_roots function\n",info); */
1047 /* exit(1); */
1048 /* } */
1049 /* else if(info > 0){ */
1050 /* //fprintf(stderr, "Eigenvalues are still uncovered in legendre_expansion_real_roots function\n"); */
1051 /* // printf("coeffs are \n"); */
1052 /* // dprint(p->num_poly, p->coeff); */
1053 /* // printf("last 2 = %G\n",p->coeff[p->num_poly-1]); */
1054 /* // exit(1); */
1055 /* } */
1056
1057 /* /\* printf("eigenvalues \n"); *\/ */
1058 /* size_t * keep = calloc_size_t(N); */
1059 /* for (ii = 0; ii < N; ii++){ */
1060 /* /\* printf("(%3.15G, %3.15G)\n",real[ii],img[ii]); *\/ */
1061 /* if ((cabs(img[ii]) < 1e-6) && (real[ii] > -1.0-1e-12) && (real[ii] < 1.0+1e-12)){ */
1062 /* /\* if ((real[ii] > -1.0-1e-12) && (real[ii] < 1.0+1e-12)){ *\/ */
1063 /* if (real[ii] < -1.0){ */
1064 /* real[ii] = -1.0; */
1065 /* } */
1066 /* else if (real[ii] > 1.0){ */
1067 /* real[ii] = 1.0; */
1068 /* } */
1069 /* keep[ii] = 1; */
1070 /* *nkeep = *nkeep + 1; */
1071 /* } */
1072 /* } */
1073
1074 /* /\* printf("nkeep = %zu\n", *nkeep); *\/ */
1075
1076 /* if (*nkeep > 0){ */
1077 /* real_roots = calloc_double(*nkeep); */
1078 /* size_t counter = 0; */
1079 /* for (ii = 0; ii < N; ii++){ */
1080 /* if (keep[ii] == 1){ */
1081 /* real_roots[counter] = real[ii]*m+off; */
1082 /* counter++; */
1083 /* } */
1084 /* } */
1085 /* } */
1086
1087
1088 /* free(keep); keep = NULL; */
1089 /* free(iwork); iwork = NULL; */
1090 /* free(real); real = NULL; */
1091 /* free(img); img = NULL; */
1092 /* free(nscompanion); nscompanion = NULL; */
1093 /* free(scale); scale = NULL; */
1094 /* } */
1095
1096 /* if (*nkeep > 1){ */
1097 /* qsort(real_roots, *nkeep, sizeof(double), dblcompare); */
1098 /* } */
1099 /* return real_roots; */
1100 /* } */
1101
1102 /* /\********************************************************\//\** */
1103 /* * Obtain the real roots of a orthogonal polynomial expansion */
1104 /* * */
1105 /* * \param[in] p - orthogonal polynomial expansion */
1106 /* * \param[in] nkeep - returns how many real roots tehre are */
1107 /* * */
1108 /* * \return real_roots - real roots of an orthonormal polynomial expansion */
1109 /* * */
1110 /* * \note */
1111 /* * Only roots within the bounds are returned */
1112 /* *************************************************************\/ */
1113 /* double * */
1114 /* orth_poly_expansion_real_roots(struct OrthPolyExpansion * p, size_t * nkeep) */
1115 /* { */
1116 /* double * real_roots = NULL; */
1117 /* enum poly_type ptype = p->p->ptype; */
1118 /* switch(ptype){ */
1119 /* case LEGENDRE: */
1120 /* real_roots = legendre_expansion_real_roots(p,nkeep); */
1121 /* break; */
1122 /* case STANDARD: */
1123 /* assert (1 == 0); */
1124 /* //x need to convert polynomial to standard polynomial first */
1125 /* //real_roots = standard_poly_real_roots(sp,nkeep); */
1126 /* //break; */
1127 /* case CHEBYSHEV: */
1128 /* real_roots = chebyshev_expansion_real_roots(p,nkeep); */
1129 /* break; */
1130 /* case HERMITE: */
1131 /* assert (1 == 0); */
1132 /* } */
1133 /* return real_roots; */
1134 /* } */
1135
1136 /* /\********************************************************\//\** */
1137 /* * Obtain the maximum of an orthogonal polynomial expansion */
1138 /* * */
1139 /* * \param[in] p - orthogonal polynomial expansion */
1140 /* * \param[in] x - location of maximum value */
1141 /* * */
1142 /* * \return maxval - maximum value */
1143 /* * */
1144 /* * \note */
1145 /* * if constant function then just returns the left most point */
1146 /* *************************************************************\/ */
1147 /* double orth_poly_expansion_max(struct OrthPolyExpansion * p, double * x) */
1148 /* { */
1149
1150 /* double maxval; */
1151 /* double tempval; */
1152
1153 /* maxval = orth_poly_expansion_eval(p,p->lower_bound); */
1154 /* *x = p->lower_bound; */
1155
1156 /* tempval = orth_poly_expansion_eval(p,p->upper_bound); */
1157 /* if (tempval > maxval){ */
1158 /* maxval = tempval; */
1159 /* *x = p->upper_bound; */
1160 /* } */
1161
1162 /* if (p->num_poly > 2){ */
1163 /* size_t nroots; */
1164 /* struct OrthPolyExpansion * deriv = orth_poly_expansion_deriv(p); */
1165 /* double * roots = orth_poly_expansion_real_roots(deriv,&nroots); */
1166 /* if (nroots > 0){ */
1167 /* size_t ii; */
1168 /* for (ii = 0; ii < nroots; ii++){ */
1169 /* tempval = orth_poly_expansion_eval(p, roots[ii]); */
1170 /* if (tempval > maxval){ */
1171 /* *x = roots[ii]; */
1172 /* maxval = tempval; */
1173 /* } */
1174 /* } */
1175 /* } */
1176
1177 /* free(roots); roots = NULL; */
1178 /* orth_poly_expansion_free(deriv); deriv = NULL; */
1179 /* } */
1180 /* return maxval; */
1181 /* } */
1182
1183 /* /\********************************************************\//\** */
1184 /* * Obtain the minimum of an orthogonal polynomial expansion */
1185 /* * */
1186 /* * \param[in] p - orthogonal polynomial expansion */
1187 /* * \param[in,out] x - location of minimum value */
1188 /* * */
1189 /* * \return minval - minimum value */
1190 /* *************************************************************\/ */
1191 /* double orth_poly_expansion_min(struct OrthPolyExpansion * p, double * x) */
1192 /* { */
1193
1194 /* double minval; */
1195 /* double tempval; */
1196
1197 /* minval = orth_poly_expansion_eval(p,p->lower_bound); */
1198 /* *x = p->lower_bound; */
1199
1200 /* tempval = orth_poly_expansion_eval(p,p->upper_bound); */
1201 /* if (tempval < minval){ */
1202 /* minval = tempval; */
1203 /* *x = p->upper_bound; */
1204 /* } */
1205
1206 /* if (p->num_poly > 2){ */
1207 /* size_t nroots; */
1208 /* struct OrthPolyExpansion * deriv = orth_poly_expansion_deriv(p); */
1209 /* double * roots = orth_poly_expansion_real_roots(deriv,&nroots); */
1210 /* if (nroots > 0){ */
1211 /* size_t ii; */
1212 /* for (ii = 0; ii < nroots; ii++){ */
1213 /* tempval = orth_poly_expansion_eval(p, roots[ii]); */
1214 /* if (tempval < minval){ */
1215 /* *x = roots[ii]; */
1216 /* minval = tempval; */
1217 /* } */
1218 /* } */
1219 /* } */
1220 /* free(roots); roots = NULL; */
1221 /* orth_poly_expansion_free(deriv); deriv = NULL; */
1222 /* } */
1223 /* return minval; */
1224 /* } */
1225
1226 /* /\********************************************************\//\** */
1227 /* * Obtain the maximum in absolute value of an orthogonal polynomial expansion */
1228 /* * */
1229 /* * \param[in] p - orthogonal polynomial expansion */
1230 /* * \param[in,out] x - location of maximum */
1231 /* * \param[in] oargs - optimization arguments */
1232 /* * required for HERMITE otherwise can set NULL */
1233 /* * */
1234 /* * \return maxval : maximum value (absolute value) */
1235 /* * */
1236 /* * \note */
1237 /* * if no roots then either lower or upper bound */
1238 /* *************************************************************\/ */
1239 /* double orth_poly_expansion_absmax( */
1240 /* struct OrthPolyExpansion * p, double * x, void * oargs) */
1241 /* { */
1242
1243 /* //printf("in absmax\n"); */
1244 /* // print_orth_poly_expansion(p,3,NULL); */
1245 /* //printf("%G\n", orth_poly_expansion_norm(p)); */
1246
1247 /* enum poly_type ptype = p->p->ptype; */
1248 /* if (oargs != NULL){ */
1249
1250 /* struct c3Vector * optnodes = oargs; */
1251 /* double mval = cabs(orth_poly_expansion_eval(p,optnodes->elem[0])); */
1252 /* *x = optnodes->elem[0]; */
1253 /* double cval = mval; */
1254 /* if (ptype == HERMITE){ */
1255 /* mval *= exp(-pow(optnodes->elem[0],2)/2.0); */
1256 /* } */
1257 /* *x = optnodes->elem[0]; */
1258 /* for (size_t ii = 0; ii < optnodes->size; ii++){ */
1259 /* double val = cabs(orth_poly_expansion_eval(p,optnodes->elem[ii])); */
1260 /* double tval = val; */
1261 /* if (ptype == HERMITE){ */
1262 /* val *= exp(-pow(optnodes->elem[ii],2)/2.0); */
1263 /* //printf("ii=%zu, x = %G. val=%G, tval=%G\n",ii,optnodes->elem[ii],val,tval); */
1264 /* } */
1265 /* if (val > mval){ */
1266 /* // printf("min achieved\n"); */
1267 /* mval = val; */
1268 /* cval = tval; */
1269 /* *x = optnodes->elem[ii]; */
1270 /* } */
1271 /* } */
1272 /* // printf("optloc=%G .... cval = %G\n",*x,cval); */
1273 /* return cval; */
1274 /* } */
1275 /* else if (ptype == HERMITE){ */
1276 /* fprintf(stderr,"Must specify optimizatino arguments\n"); */
1277 /* fprintf(stderr,"In the form of candidate points for \n"); */
1278 /* fprintf(stderr,"finding the absmax of hermite expansion\n"); */
1279 /* exit(1); */
1280
1281 /* } */
1282 /* double maxval; */
1283 /* double norm = orth_poly_expansion_norm(p); */
1284
1285 /* if (norm < ZEROTHRESH) { */
1286 /* *x = p->lower_bound; */
1287 /* maxval = 0.0; */
1288 /* } */
1289 /* else{ */
1290 /* //printf("nroots=%zu\n", nroots); */
1291 /* double tempval; */
1292
1293 /* maxval = cabs(orth_poly_expansion_eval(p,p->lower_bound)); */
1294 /* *x = p->lower_bound; */
1295
1296 /* tempval = cabs(orth_poly_expansion_eval(p,p->upper_bound)); */
1297 /* if (tempval > maxval){ */
1298 /* maxval = tempval; */
1299 /* *x = p->upper_bound; */
1300 /* } */
1301 /* if (p->num_poly > 2){ */
1302 /* size_t nroots; */
1303 /* struct OrthPolyExpansion * deriv = orth_poly_expansion_deriv(p); */
1304 /* double * roots = orth_poly_expansion_real_roots(deriv,&nroots); */
1305 /* if (nroots > 0){ */
1306 /* size_t ii; */
1307 /* for (ii = 0; ii < nroots; ii++){ */
1308 /* tempval = cabs(orth_poly_expansion_eval(p, roots[ii])); */
1309 /* if (tempval > maxval){ */
1310 /* *x = roots[ii]; */
1311 /* maxval = tempval; */
1312 /* } */
1313 /* } */
1314 /* } */
1315
1316 /* free(roots); roots = NULL; */
1317 /* orth_poly_expansion_free(deriv); deriv = NULL; */
1318 /* } */
1319 /* } */
1320 /* //printf("done\n"); */
1321 /* return maxval; */
1322 /* } */
1323
1324
1325
1326