Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/quad.cc @ 10154:40dfc0c99116

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
DLD-FUNCTIONS/*.cc: untabify
author John W. Eaton <jwe@octave.org>
date Wed, 20 Jan 2010 17:33:41 -0500
parents 2cd940306a06
children d0ce5e973937
comparison
equal deleted inserted replaced
10153:2c28f9d0360f 10154:40dfc0c99116
69 if (quad_fcn) 69 if (quad_fcn)
70 { 70 {
71 octave_value_list tmp = quad_fcn->do_multi_index_op (1, args); 71 octave_value_list tmp = quad_fcn->do_multi_index_op (1, args);
72 72
73 if (error_state) 73 if (error_state)
74 { 74 {
75 quad_integration_error = 1; // FIXME 75 quad_integration_error = 1; // FIXME
76 gripe_user_supplied_eval ("quad"); 76 gripe_user_supplied_eval ("quad");
77 return retval; 77 return retval;
78 } 78 }
79 79
80 if (tmp.length () && tmp(0).is_defined ()) 80 if (tmp.length () && tmp(0).is_defined ())
81 { 81 {
82 if (! warned_imaginary && tmp(0).is_complex_type ()) 82 if (! warned_imaginary && tmp(0).is_complex_type ())
83 { 83 {
84 warning ("quad: ignoring imaginary part returned from user-supplied function"); 84 warning ("quad: ignoring imaginary part returned from user-supplied function");
85 warned_imaginary = true; 85 warned_imaginary = true;
86 } 86 }
87 87
88 retval = tmp(0).double_value (); 88 retval = tmp(0).double_value ();
89 89
90 if (error_state) 90 if (error_state)
91 { 91 {
92 quad_integration_error = 1; // FIXME 92 quad_integration_error = 1; // FIXME
93 gripe_user_supplied_eval ("quad"); 93 gripe_user_supplied_eval ("quad");
94 } 94 }
95 } 95 }
96 else 96 else
97 { 97 {
98 quad_integration_error = 1; // FIXME 98 quad_integration_error = 1; // FIXME
99 gripe_user_supplied_eval ("quad"); 99 gripe_user_supplied_eval ("quad");
100 } 100 }
101 } 101 }
102 102
103 return retval; 103 return retval;
104 } 104 }
105 105
114 if (quad_fcn) 114 if (quad_fcn)
115 { 115 {
116 octave_value_list tmp = quad_fcn->do_multi_index_op (1, args); 116 octave_value_list tmp = quad_fcn->do_multi_index_op (1, args);
117 117
118 if (error_state) 118 if (error_state)
119 { 119 {
120 quad_integration_error = 1; // FIXME 120 quad_integration_error = 1; // FIXME
121 gripe_user_supplied_eval ("quad"); 121 gripe_user_supplied_eval ("quad");
122 return retval; 122 return retval;
123 } 123 }
124 124
125 if (tmp.length () && tmp(0).is_defined ()) 125 if (tmp.length () && tmp(0).is_defined ())
126 { 126 {
127 if (! warned_imaginary && tmp(0).is_complex_type ()) 127 if (! warned_imaginary && tmp(0).is_complex_type ())
128 { 128 {
129 warning ("quad: ignoring imaginary part returned from user-supplied function"); 129 warning ("quad: ignoring imaginary part returned from user-supplied function");
130 warned_imaginary = true; 130 warned_imaginary = true;
131 } 131 }
132 132
133 retval = tmp(0).float_value (); 133 retval = tmp(0).float_value ();
134 134
135 if (error_state) 135 if (error_state)
136 { 136 {
137 quad_integration_error = 1; // FIXME 137 quad_integration_error = 1; // FIXME
138 gripe_user_supplied_eval ("quad"); 138 gripe_user_supplied_eval ("quad");
139 } 139 }
140 } 140 }
141 else 141 else
142 { 142 {
143 quad_integration_error = 1; // FIXME 143 quad_integration_error = 1; // FIXME
144 gripe_user_supplied_eval ("quad"); 144 gripe_user_supplied_eval ("quad");
145 } 145 }
146 } 146 }
147 147
148 return retval; 148 return retval;
149 } 149 }
150 150
151 #define QUAD_ABORT() \ 151 #define QUAD_ABORT() \
152 do \ 152 do \
153 { \ 153 { \
154 if (fcn_name.length()) \ 154 if (fcn_name.length()) \
155 clear_function (fcn_name); \ 155 clear_function (fcn_name); \
156 return retval; \ 156 return retval; \
157 } \ 157 } \
158 while (0) 158 while (0)
159 159
160 #define QUAD_ABORT1(msg) \ 160 #define QUAD_ABORT1(msg) \
231 int nargin = args.length (); 231 int nargin = args.length ();
232 232
233 if (nargin > 2 && nargin < 6 && nargout < 5) 233 if (nargin > 2 && nargin < 6 && nargout < 5)
234 { 234 {
235 if (args(0).is_function_handle () || args(0).is_inline_function ()) 235 if (args(0).is_function_handle () || args(0).is_inline_function ())
236 quad_fcn = args(0).function_value (); 236 quad_fcn = args(0).function_value ();
237 else 237 else
238 { 238 {
239 fcn_name = unique_symbol_name ("__quad_fcn_"); 239 fcn_name = unique_symbol_name ("__quad_fcn_");
240 std::string fname = "function y = "; 240 std::string fname = "function y = ";
241 fname.append (fcn_name); 241 fname.append (fcn_name);
242 fname.append ("(x) y = "); 242 fname.append ("(x) y = ");
243 quad_fcn = extract_function (args(0), "quad", fcn_name, fname, 243 quad_fcn = extract_function (args(0), "quad", fcn_name, fname,
244 "; endfunction"); 244 "; endfunction");
245 } 245 }
246 246
247 if (! quad_fcn) 247 if (! quad_fcn)
248 QUAD_ABORT (); 248 QUAD_ABORT ();
249 249
250 if (args(1).is_single_type () || args(2).is_single_type ()) 250 if (args(1).is_single_type () || args(2).is_single_type ())
251 { 251 {
252 float a = args(1).float_value (); 252 float a = args(1).float_value ();
253 253
254 if (error_state) 254 if (error_state)
255 QUAD_ABORT1 ("expecting second argument to be a scalar"); 255 QUAD_ABORT1 ("expecting second argument to be a scalar");
256 256
257 float b = args(2).float_value (); 257 float b = args(2).float_value ();
258 258
259 if (error_state) 259 if (error_state)
260 QUAD_ABORT1 ("expecting third argument to be a scalar"); 260 QUAD_ABORT1 ("expecting third argument to be a scalar");
261 261
262 int indefinite = 0; 262 int indefinite = 0;
263 FloatIndefQuad::IntegralType indef_type = FloatIndefQuad::doubly_infinite; 263 FloatIndefQuad::IntegralType indef_type = FloatIndefQuad::doubly_infinite;
264 float bound = 0.0; 264 float bound = 0.0;
265 if (xisinf (a) && xisinf (b)) 265 if (xisinf (a) && xisinf (b))
266 { 266 {
267 indefinite = 1; 267 indefinite = 1;
268 indef_type = FloatIndefQuad::doubly_infinite; 268 indef_type = FloatIndefQuad::doubly_infinite;
269 } 269 }
270 else if (xisinf (a)) 270 else if (xisinf (a))
271 { 271 {
272 indefinite = 1; 272 indefinite = 1;
273 bound = b; 273 bound = b;
274 indef_type = FloatIndefQuad::neg_inf_to_bound; 274 indef_type = FloatIndefQuad::neg_inf_to_bound;
275 } 275 }
276 else if (xisinf (b)) 276 else if (xisinf (b))
277 { 277 {
278 indefinite = 1; 278 indefinite = 1;
279 bound = a; 279 bound = a;
280 indef_type = FloatIndefQuad::bound_to_inf; 280 indef_type = FloatIndefQuad::bound_to_inf;
281 } 281 }
282 282
283 octave_idx_type ier = 0; 283 octave_idx_type ier = 0;
284 octave_idx_type nfun = 0; 284 octave_idx_type nfun = 0;
285 float abserr = 0.0; 285 float abserr = 0.0;
286 float val = 0.0; 286 float val = 0.0;
287 bool have_sing = false; 287 bool have_sing = false;
288 FloatColumnVector sing; 288 FloatColumnVector sing;
289 FloatColumnVector tol; 289 FloatColumnVector tol;
290 290
291 switch (nargin) 291 switch (nargin)
292 { 292 {
293 case 5: 293 case 5:
294 if (indefinite) 294 if (indefinite)
295 QUAD_ABORT1 ("singularities not allowed on infinite intervals"); 295 QUAD_ABORT1 ("singularities not allowed on infinite intervals");
296 296
297 have_sing = true; 297 have_sing = true;
298 298
299 sing = FloatColumnVector (args(4).float_vector_value ()); 299 sing = FloatColumnVector (args(4).float_vector_value ());
300 300
301 if (error_state) 301 if (error_state)
302 QUAD_ABORT1 ("expecting vector of singularities as fourth argument"); 302 QUAD_ABORT1 ("expecting vector of singularities as fourth argument");
303 303
304 case 4: 304 case 4:
305 tol = FloatColumnVector (args(3).float_vector_value ()); 305 tol = FloatColumnVector (args(3).float_vector_value ());
306 306
307 if (error_state) 307 if (error_state)
308 QUAD_ABORT1 ("expecting vector of tolerances as fifth argument"); 308 QUAD_ABORT1 ("expecting vector of tolerances as fifth argument");
309 309
310 switch (tol.capacity ()) 310 switch (tol.capacity ())
311 { 311 {
312 case 2: 312 case 2:
313 quad_opts.set_single_precision_relative_tolerance (tol (1)); 313 quad_opts.set_single_precision_relative_tolerance (tol (1));
314 314
315 case 1: 315 case 1:
316 quad_opts.set_single_precision_absolute_tolerance (tol (0)); 316 quad_opts.set_single_precision_absolute_tolerance (tol (0));
317 break; 317 break;
318 318
319 default: 319 default:
320 QUAD_ABORT1 ("expecting tol to contain no more than two values"); 320 QUAD_ABORT1 ("expecting tol to contain no more than two values");
321 } 321 }
322 322
323 case 3: 323 case 3:
324 if (indefinite) 324 if (indefinite)
325 { 325 {
326 FloatIndefQuad iq (quad_float_user_function, bound, 326 FloatIndefQuad iq (quad_float_user_function, bound,
327 indef_type); 327 indef_type);
328 iq.set_options (quad_opts); 328 iq.set_options (quad_opts);
329 val = iq.float_integrate (ier, nfun, abserr); 329 val = iq.float_integrate (ier, nfun, abserr);
330 } 330 }
331 else 331 else
332 { 332 {
333 if (have_sing) 333 if (have_sing)
334 { 334 {
335 FloatDefQuad dq (quad_float_user_function, a, b, sing); 335 FloatDefQuad dq (quad_float_user_function, a, b, sing);
336 dq.set_options (quad_opts); 336 dq.set_options (quad_opts);
337 val = dq.float_integrate (ier, nfun, abserr); 337 val = dq.float_integrate (ier, nfun, abserr);
338 } 338 }
339 else 339 else
340 { 340 {
341 FloatDefQuad dq (quad_float_user_function, a, b); 341 FloatDefQuad dq (quad_float_user_function, a, b);
342 dq.set_options (quad_opts); 342 dq.set_options (quad_opts);
343 val = dq.float_integrate (ier, nfun, abserr); 343 val = dq.float_integrate (ier, nfun, abserr);
344 } 344 }
345 } 345 }
346 break; 346 break;
347 347
348 default: 348 default:
349 panic_impossible (); 349 panic_impossible ();
350 break; 350 break;
351 } 351 }
352 352
353 retval(3) = abserr; 353 retval(3) = abserr;
354 retval(2) = nfun; 354 retval(2) = nfun;
355 retval(1) = ier; 355 retval(1) = ier;
356 retval(0) = val; 356 retval(0) = val;
357 357
358 } 358 }
359 else 359 else
360 { 360 {
361 double a = args(1).double_value (); 361 double a = args(1).double_value ();
362 362
363 if (error_state) 363 if (error_state)
364 QUAD_ABORT1 ("expecting second argument to be a scalar"); 364 QUAD_ABORT1 ("expecting second argument to be a scalar");
365 365
366 double b = args(2).double_value (); 366 double b = args(2).double_value ();
367 367
368 if (error_state) 368 if (error_state)
369 QUAD_ABORT1 ("expecting third argument to be a scalar"); 369 QUAD_ABORT1 ("expecting third argument to be a scalar");
370 370
371 int indefinite = 0; 371 int indefinite = 0;
372 IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite; 372 IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite;
373 double bound = 0.0; 373 double bound = 0.0;
374 if (xisinf (a) && xisinf (b)) 374 if (xisinf (a) && xisinf (b))
375 { 375 {
376 indefinite = 1; 376 indefinite = 1;
377 indef_type = IndefQuad::doubly_infinite; 377 indef_type = IndefQuad::doubly_infinite;
378 } 378 }
379 else if (xisinf (a)) 379 else if (xisinf (a))
380 { 380 {
381 indefinite = 1; 381 indefinite = 1;
382 bound = b; 382 bound = b;
383 indef_type = IndefQuad::neg_inf_to_bound; 383 indef_type = IndefQuad::neg_inf_to_bound;
384 } 384 }
385 else if (xisinf (b)) 385 else if (xisinf (b))
386 { 386 {
387 indefinite = 1; 387 indefinite = 1;
388 bound = a; 388 bound = a;
389 indef_type = IndefQuad::bound_to_inf; 389 indef_type = IndefQuad::bound_to_inf;
390 } 390 }
391 391
392 octave_idx_type ier = 0; 392 octave_idx_type ier = 0;
393 octave_idx_type nfun = 0; 393 octave_idx_type nfun = 0;
394 double abserr = 0.0; 394 double abserr = 0.0;
395 double val = 0.0; 395 double val = 0.0;
396 bool have_sing = false; 396 bool have_sing = false;
397 ColumnVector sing; 397 ColumnVector sing;
398 ColumnVector tol; 398 ColumnVector tol;
399 399
400 switch (nargin) 400 switch (nargin)
401 { 401 {
402 case 5: 402 case 5:
403 if (indefinite) 403 if (indefinite)
404 QUAD_ABORT1 ("singularities not allowed on infinite intervals"); 404 QUAD_ABORT1 ("singularities not allowed on infinite intervals");
405 405
406 have_sing = true; 406 have_sing = true;
407 407
408 sing = ColumnVector (args(4).vector_value ()); 408 sing = ColumnVector (args(4).vector_value ());
409 409
410 if (error_state) 410 if (error_state)
411 QUAD_ABORT1 ("expecting vector of singularities as fourth argument"); 411 QUAD_ABORT1 ("expecting vector of singularities as fourth argument");
412 412
413 case 4: 413 case 4:
414 tol = ColumnVector (args(3).vector_value ()); 414 tol = ColumnVector (args(3).vector_value ());
415 415
416 if (error_state) 416 if (error_state)
417 QUAD_ABORT1 ("expecting vector of tolerances as fifth argument"); 417 QUAD_ABORT1 ("expecting vector of tolerances as fifth argument");
418 418
419 switch (tol.capacity ()) 419 switch (tol.capacity ())
420 { 420 {
421 case 2: 421 case 2:
422 quad_opts.set_relative_tolerance (tol (1)); 422 quad_opts.set_relative_tolerance (tol (1));
423 423
424 case 1: 424 case 1:
425 quad_opts.set_absolute_tolerance (tol (0)); 425 quad_opts.set_absolute_tolerance (tol (0));
426 break; 426 break;
427 427
428 default: 428 default:
429 QUAD_ABORT1 ("expecting tol to contain no more than two values"); 429 QUAD_ABORT1 ("expecting tol to contain no more than two values");
430 } 430 }
431 431
432 case 3: 432 case 3:
433 if (indefinite) 433 if (indefinite)
434 { 434 {
435 IndefQuad iq (quad_user_function, bound, indef_type); 435 IndefQuad iq (quad_user_function, bound, indef_type);
436 iq.set_options (quad_opts); 436 iq.set_options (quad_opts);
437 val = iq.integrate (ier, nfun, abserr); 437 val = iq.integrate (ier, nfun, abserr);
438 } 438 }
439 else 439 else
440 { 440 {
441 if (have_sing) 441 if (have_sing)
442 { 442 {
443 DefQuad dq (quad_user_function, a, b, sing); 443 DefQuad dq (quad_user_function, a, b, sing);
444 dq.set_options (quad_opts); 444 dq.set_options (quad_opts);
445 val = dq.integrate (ier, nfun, abserr); 445 val = dq.integrate (ier, nfun, abserr);
446 } 446 }
447 else 447 else
448 { 448 {
449 DefQuad dq (quad_user_function, a, b); 449 DefQuad dq (quad_user_function, a, b);
450 dq.set_options (quad_opts); 450 dq.set_options (quad_opts);
451 val = dq.integrate (ier, nfun, abserr); 451 val = dq.integrate (ier, nfun, abserr);
452 } 452 }
453 } 453 }
454 break; 454 break;
455 455
456 default: 456 default:
457 panic_impossible (); 457 panic_impossible ();
458 break; 458 break;
459 } 459 }
460 460
461 retval(3) = abserr; 461 retval(3) = abserr;
462 retval(2) = nfun; 462 retval(2) = nfun;
463 retval(1) = ier; 463 retval(1) = ier;
464 retval(0) = val; 464 retval(0) = val;
465 } 465 }
466 466
467 if (fcn_name.length()) 467 if (fcn_name.length())
468 clear_function (fcn_name); 468 clear_function (fcn_name);
469 } 469 }
470 else 470 else
471 print_usage (); 471 print_usage ();
472 472
473 return retval; 473 return retval;