Mercurial > octave
annotate liboctave/numeric/CollocWt.cc @ 21724:aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
* configure.ac, Doxyfile.in, make_int.cc, Backend.cc, BaseControl.cc,
ButtonControl.cc, Canvas.cc, CheckBoxControl.cc, Container.cc,
ContextMenu.cc, EditControl.cc, Figure.cc, FigureWindow.cc,
GLCanvas.cc, KeyMap.cc, ListBoxControl.cc, Logger.cc, Menu.cc,
MouseModeActionGroup.cc, Object.cc, ObjectFactory.cc,
ObjectProxy.cc, Panel.cc, PopupMenuControl.cc, PushButtonControl.cc,
PushTool.cc, QtHandlesUtils.cc, RadioButtonControl.cc,
SliderControl.cc, TextControl.cc, TextEdit.cc,
ToggleButtonControl.cc, ToggleTool.cc, ToolBar.cc, __init_qt__.cc,
annotation-dialog.cc, gl-select.cc, module.mk, color-picker.cc,
dialog.cc, documentation-dock-widget.cc, files-dock-widget.cc,
find-files-dialog.cc, find-files-model.cc, history-dock-widget.cc,
liboctgui-build-info.in.cc, file-editor-tab.cc, file-editor.cc,
find-dialog.cc, marker.cc, marker.h, octave-qscintilla.cc,
octave-txt-lexer.cc, main-window.cc, octave-cmd.cc,
octave-dock-widget.cc, octave-gui.cc, octave-interpreter.cc,
octave-qt-link.cc, parser.cc, webinfo.cc, resource-manager.cc,
settings-dialog.cc, shortcut-manager.cc, terminal-dock-widget.cc,
thread-manager.cc, welcome-wizard.cc, workspace-model.cc,
workspace-view.cc, build-env.in.cc, Cell.cc, __contourc__.cc,
__dispatch__.cc, __dsearchn__.cc, __ichol__.cc, __ilu__.cc,
__lin_interpn__.cc, __luinc__.cc, __magick_read__.cc,
__pchip_deriv__.cc, __qp__.cc, balance.cc, besselj.cc, betainc.cc,
bitfcns.cc, bsxfun.cc, c-file-ptr-stream.cc, c-file-ptr-stream.h,
cdisplay.c, cdisplay.h, cellfun.cc, coct-hdf5-types.c, colloc.cc,
comment-list.cc, conv2.cc, daspk.cc, dasrt.cc, dassl.cc, data.cc,
debug.cc, defaults.cc, defaults.in.h, defun.cc, det.cc, dirfns.cc,
display.cc, dlmread.cc, dot.cc, dynamic-ld.cc, eig.cc, ellipj.cc,
error.cc, errwarn.cc, event-queue.cc, fft.cc, fft2.cc, fftn.cc,
file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, gammainc.cc,
gcd.cc, getgrent.cc, getpwent.cc, getrusage.cc, givens.cc,
gl-render.cc, gl2ps-print.cc, graphics.cc, gripes.cc, hash.cc,
help.cc, hess.cc, hex2num.cc, input.cc, inv.cc, jit-ir.cc, jit-ir.h,
jit-typeinfo.cc, jit-typeinfo.h, jit-util.cc, jit-util.h, kron.cc,
load-path.cc, load-save.cc, lookup.cc, ls-ascii-helper.cc,
ls-hdf5.cc, ls-mat-ascii.cc, ls-mat4.cc, ls-mat5.cc,
ls-oct-binary.cc, ls-oct-text.cc, ls-oct-text.h, ls-utils.cc,
lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex.cc, mex.h,
mexproto.h, mgorth.cc, nproc.cc, oct-errno.in.cc, oct-fstrm.cc,
oct-hdf5-types.cc, oct-hdf5-types.h, oct-hist.cc, oct-iostrm.cc,
oct-lvalue.cc, oct-map.cc, oct-prcstrm.cc, oct-procbuf.cc,
oct-stdstrm.h, oct-stream.cc, oct-strstrm.cc, oct-tex-lexer.in.ll,
oct-tex-parser.in.yy, octave-link.cc, ordschur.cc, pager.cc,
pinv.cc, pr-output.cc, procstream.cc, profiler.cc, psi.cc,
pt-jit.cc, pt-jit.h, quad.cc, quadcc.cc, qz.cc, rand.cc, rcond.cc,
regexp.cc, schur.cc, sighandlers.cc, sighandlers.h, siglist.c,
siglist.h, sparse-xdiv.cc, sparse-xpow.cc, sparse.cc, spparms.cc,
sqrtm.cc, str2double.cc, strfind.cc, strfns.cc, sub2ind.cc, svd.cc,
sylvester.cc, symtab.cc, syscalls.cc, sysdep.cc, text-renderer.cc,
time.cc, toplev.cc, tril.cc, tsearch.cc, txt-eng.cc, typecast.cc,
urlwrite.cc, utils.cc, variables.cc, xdiv.cc, xgl2ps.c, xnorm.cc,
xpow.cc, zfstream.cc, zfstream.h, __delaunayn__.cc, __eigs__.cc,
__fltk_uigetfile__.cc, __glpk__.cc, __init_fltk__.cc,
__init_gnuplot__.cc, __osmesa_print__.cc, __voronoi__.cc, amd.cc,
audiodevinfo.cc, audioread.cc, ccolamd.cc, chol.cc, colamd.cc,
convhulln.cc, dmperm.cc, fftw.cc, qr.cc, symbfact.cc, symrcm.cc,
liboctinterp-build-info.in.cc, mkbuiltins, mkops, ov-base.cc,
ov-bool-mat.cc, ov-bool-sparse.cc, ov-bool.cc, ov-builtin.cc,
ov-cell.cc, ov-ch-mat.cc, ov-class.cc, ov-classdef.cc, ov-colon.cc,
ov-complex.cc, ov-cs-list.cc, ov-cx-diag.cc, ov-cx-mat.cc,
ov-cx-sparse.cc, ov-dld-fcn.cc, ov-fcn-handle.cc, ov-fcn-inline.cc,
ov-fcn.cc, ov-float.cc, ov-flt-complex.cc, ov-flt-cx-diag.cc,
ov-flt-cx-mat.cc, ov-flt-re-diag.cc, ov-flt-re-mat.cc, ov-int16.cc,
ov-int32.cc, ov-int64.cc, ov-int8.cc, ov-java.cc, ov-lazy-idx.cc,
ov-mex-fcn.cc, ov-null-mat.cc, ov-oncleanup.cc, ov-perm.cc,
ov-range.cc, ov-re-diag.cc, ov-re-mat.cc, ov-re-sparse.cc,
ov-scalar.cc, ov-str-mat.cc, ov-struct.cc, ov-typeinfo.cc,
ov-uint16.cc, ov-uint32.cc, ov-uint64.cc, ov-uint8.cc,
ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ovl.cc, octave.cc, octave.h,
op-b-b.cc, op-b-bm.cc, op-b-sbm.cc, op-bm-b.cc, op-bm-bm.cc,
op-bm-sbm.cc, op-cdm-cdm.cc, op-cdm-cm.cc, op-cdm-dm.cc,
op-cdm-m.cc, op-cell.cc, op-chm.cc, op-class.cc, op-cm-cdm.cc,
op-cm-cm.cc, op-cm-cs.cc, op-cm-dm.cc, op-cm-m.cc, op-cm-s.cc,
op-cm-scm.cc, op-cm-sm.cc, op-cs-cm.cc, op-cs-cs.cc, op-cs-m.cc,
op-cs-s.cc, op-cs-scm.cc, op-cs-sm.cc, op-dm-cdm.cc, op-dm-cm.cc,
op-dm-dm.cc, op-dm-m.cc, op-dm-scm.cc, op-dm-sm.cc,
op-dm-template.cc, op-dms-template.cc, op-double-conv.cc,
op-fcdm-fcdm.cc, op-fcdm-fcm.cc, op-fcdm-fdm.cc, op-fcdm-fm.cc,
op-fcm-fcdm.cc, op-fcm-fcm.cc, op-fcm-fcs.cc, op-fcm-fdm.cc,
op-fcm-fm.cc, op-fcm-fs.cc, op-fcn.cc, op-fcs-fcm.cc, op-fcs-fcs.cc,
op-fcs-fm.cc, op-fcs-fs.cc, op-fdm-fcdm.cc, op-fdm-fcm.cc,
op-fdm-fdm.cc, op-fdm-fm.cc, op-float-conv.cc, op-fm-fcdm.cc,
op-fm-fcm.cc, op-fm-fcs.cc, op-fm-fdm.cc, op-fm-fm.cc, op-fm-fs.cc,
op-fs-fcm.cc, op-fs-fcs.cc, op-fs-fm.cc, op-fs-fs.cc, op-i16-i16.cc,
op-i32-i32.cc, op-i64-i64.cc, op-i8-i8.cc, op-int-concat.cc,
op-int-conv.cc, op-m-cdm.cc, op-m-cm.cc, op-m-cs.cc, op-m-dm.cc,
op-m-m.cc, op-m-s.cc, op-m-scm.cc, op-m-sm.cc, op-pm-pm.cc,
op-pm-scm.cc, op-pm-sm.cc, op-pm-template.cc, op-range.cc,
op-s-cm.cc, op-s-cs.cc, op-s-m.cc, op-s-s.cc, op-s-scm.cc,
op-s-sm.cc, op-sbm-b.cc, op-sbm-bm.cc, op-sbm-sbm.cc, op-scm-cm.cc,
op-scm-cs.cc, op-scm-m.cc, op-scm-s.cc, op-scm-scm.cc, op-scm-sm.cc,
op-sm-cm.cc, op-sm-cs.cc, op-sm-m.cc, op-sm-s.cc, op-sm-scm.cc,
op-sm-sm.cc, op-str-m.cc, op-str-s.cc, op-str-str.cc, op-struct.cc,
op-ui16-ui16.cc, op-ui32-ui32.cc, op-ui64-ui64.cc, op-ui8-ui8.cc,
lex.ll, oct-parse.in.yy, pt-arg-list.cc, pt-array-list.cc,
pt-assign.cc, pt-binop.cc, pt-bp.cc, pt-cbinop.cc, pt-cell.cc,
pt-check.cc, pt-classdef.cc, pt-cmd.cc, pt-colon.cc, pt-const.cc,
pt-decl.cc, pt-eval.cc, pt-except.cc, pt-exp.cc, pt-fcn-handle.cc,
pt-funcall.cc, pt-id.cc, pt-idx.cc, pt-jump.cc, pt-loop.cc,
pt-loop.h, pt-mat.cc, pt-misc.cc, pt-pr-code.cc, pt-select.cc,
pt-stmt.cc, pt-unop.cc, pt.cc, token.cc, Array-jit.cc, Array-tc.cc,
version.cc, Array-C.cc, Array-b.cc, Array-ch.cc, Array-d.cc,
Array-f.cc, Array-fC.cc, Array-i.cc, Array-idx-vec.cc, Array-s.cc,
Array-str.cc, Array-util.cc, Array-voidp.cc, Array.cc,
CColVector.cc, CDiagMatrix.cc, CMatrix.cc, CNDArray.cc,
CRowVector.cc, CSparse.cc, MArray-C.cc, MArray-d.cc, MArray-f.cc,
MArray-fC.cc, MArray-i.cc, MArray-s.cc, MSparse-C.cc, MSparse-d.cc,
MatrixType.cc, PermMatrix.cc, Range.cc, Sparse-C.cc, Sparse-b.cc,
Sparse-d.cc, boolMatrix.cc, boolNDArray.cc, boolSparse.cc,
chMatrix.cc, chNDArray.cc, dColVector.cc, dDiagMatrix.cc,
dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, dim-vector.cc,
dim-vector.h, fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc,
fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc,
fMatrix.cc, fNDArray.cc, fRowVector.cc, idx-vector.cc,
int16NDArray.cc, int32NDArray.cc, int64NDArray.cc, int8NDArray.cc,
uint16NDArray.cc, uint32NDArray.cc, uint64NDArray.cc,
uint8NDArray.cc, Faddeeva.cc, blaswrap.c, cquit.c, f77-extern.cc,
f77-fcn.c, f77-fcn.h, lo-error.c, lo-error.h, quit.cc, quit.h,
liboctave-build-info.in.cc, CollocWt.cc, DASPK.cc, DASRT.cc,
DASSL.cc, EIG.cc, LSODE.cc, ODES.cc, Quad.cc, aepbalance.cc,
chol.cc, eigs-base.cc, fEIG.cc, gepbalance.cc, hess.cc,
lo-mappers.cc, lo-specfun.cc, lu.cc, oct-convn.cc, oct-fftw.cc,
oct-norm.cc, oct-rand.cc, oct-spparms.cc, qr.cc, qrp.cc,
randmtzig.cc, randpoisson.cc, schur.cc, sparse-chol.cc,
sparse-dmsolve.cc, sparse-lu.cc, sparse-qr.cc, svd.cc, mk-ops.awk,
mx-defs.h, dir-ops.cc, file-ops.cc, file-stat.cc, lo-sysdep.cc,
mach-info.cc, oct-env.cc, oct-group.cc, oct-passwd.cc,
oct-syscalls.cc, oct-time.cc, oct-uname.cc, pathlen.h, syswait.h,
cmd-edit.cc, cmd-hist.cc, data-conv.cc, f2c-main.c, glob-match.cc,
kpse.cc, lo-array-errwarn.cc, lo-array-gripes.cc, lo-cutils.c,
lo-cutils.h, lo-ieee.cc, lo-ieee.h, lo-regexp.cc, lo-utils.cc,
oct-base64.cc, oct-glob.cc, oct-inttypes.cc, oct-inttypes.h,
oct-locbuf.cc, oct-mutex.cc, oct-rl-edit.c, oct-rl-edit.h,
oct-rl-hist.c, oct-rl-hist.h, oct-shlib.cc, oct-sort.cc,
pathsearch.cc, singleton-cleanup.cc, sparse-sort.cc, sparse-util.cc,
statdefs.h, str-vec.cc, unwind-prot.cc, url-transfer.cc,
acinclude.m4, display-available.c, display-available.h, main-cli.cc,
main-gui.cc, main.in.cc, mkoctfile.in.cc, octave-build-info.in.cc,
octave-config.in.cc, shared-fcns.h:
Use "#if ..." consistently instead of "#ifdef" and "#ifndef".
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 17 May 2016 12:09:30 -0400 |
parents | 3d60ed163b70 |
children | b571fc85953f |
rev | line source |
---|---|
3 | 1 /* |
2 | |
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
17744
diff
changeset
|
3 Copyright (C) 1993-2015 John W. Eaton |
3 | 4 |
5 This file is part of Octave. | |
6 | |
7 Octave is free software; you can redistribute it and/or modify it | |
8 under the terms of the GNU General Public License as published by the | |
7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
10 option) any later version. | |
3 | 11 |
12 Octave is distributed in the hope that it will be useful, but WITHOUT | |
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
15 for more details. | |
16 | |
17 You should have received a copy of the GNU General Public License | |
7016 | 18 along with Octave; see the file COPYING. If not, see |
19 <http://www.gnu.org/licenses/>. | |
3 | 20 |
21 */ | |
22 | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21568
diff
changeset
|
23 #if defined (HAVE_CONFIG_H) |
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21202
diff
changeset
|
24 # include "config.h" |
3 | 25 #endif |
26 | |
3503 | 27 #include <iostream> |
238 | 28 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
29 #include <cfloat> |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
30 |
3 | 31 #include "CollocWt.h" |
1847 | 32 #include "f77-fcn.h" |
227 | 33 #include "lo-error.h" |
3 | 34 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
35 // The following routines jcobi, dif, and dfopr are based on the code |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
36 // found in Villadsen, J. and M. L. Michelsen, Solution of Differential |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
37 // Equation Models by Polynomial Approximation, Prentice-Hall (1978) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
38 // pages 418-420. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
39 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
40 // Translated to C++ by jwe. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
41 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
42 // Compute the first three derivatives of the node polynomial. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
43 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
44 // n0 (alpha,beta) n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
45 // p (x) = (x) * p (x) * (1 - x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
46 // nt n |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
47 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
48 // at the interpolation points. Each of the parameters n0 and n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
49 // may be given the value 0 or 1. The total number of points |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
50 // nt = n + n0 + n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
51 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
52 // The values of root must be known before a call to dif is possible. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
53 // They may be computed using jcobi. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
54 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
55 static void |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
56 dif (octave_idx_type nt, double *root, double *dif1, double *dif2, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
57 double *dif3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
58 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
59 // Evaluate derivatives of node polynomial using recursion formulas. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
60 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
61 for (octave_idx_type i = 0; i < nt; i++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
62 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
63 double x = root[i]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
64 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
65 dif1[i] = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
66 dif2[i] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
67 dif3[i] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
68 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
69 for (octave_idx_type j = 0; j < nt; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
70 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
71 if (j != i) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
72 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
73 double y = x - root[j]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
74 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
75 dif3[i] = y * dif3[i] + 3.0 * dif2[i]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
76 dif2[i] = y * dif2[i] + 2.0 * dif1[i]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
77 dif1[i] = y * dif1[i]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
78 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
79 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
80 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
81 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
82 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
83 // Compute the zeros of the Jacobi polynomial. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
84 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
85 // (alpha,beta) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
86 // p (x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
87 // n |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
88 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
89 // Use dif to compute the derivatives of the node |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
90 // polynomial |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
91 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
92 // n0 (alpha,beta) n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
93 // p (x) = (x) * p (x) * (1 - x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
94 // nt n |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
95 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
96 // at the interpolation points. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
97 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
98 // See Villadsen and Michelsen, pages 131-132 and 418. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
99 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
100 // Input parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
101 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
102 // nd : the dimension of the vectors dif1, dif2, dif3, and root |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
103 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
104 // n : the degree of the jacobi polynomial, (i.e. the number |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
105 // of interior interpolation points) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
106 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
107 // n0 : determines whether x = 0 is included as an |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
108 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
109 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
110 // n0 = 0 ==> x = 0 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
111 // n0 = 1 ==> x = 0 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
112 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
113 // n1 : determines whether x = 1 is included as an |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
114 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
115 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
116 // n1 = 0 ==> x = 1 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
117 // n1 = 1 ==> x = 1 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
118 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
119 // alpha : the value of alpha in the description of the jacobi |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
120 // polynomial |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
121 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
122 // beta : the value of beta in the description of the jacobi |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
123 // polynomial |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
124 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
125 // For a more complete explanation of alpha an beta, see Villadsen |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
126 // and Michelsen, pages 57 to 59. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
127 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
128 // Output parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
129 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
130 // root : one dimensional vector containing on exit the |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
131 // n + n0 + n1 zeros of the node polynomial used in the |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
132 // interpolation routine |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
133 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
134 // dif1 : one dimensional vector containing the first derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
135 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
136 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
137 // dif2 : one dimensional vector containing the second derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
138 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
139 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
140 // dif3 : one dimensional vector containing the third derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
141 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
142 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
143 static bool |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
144 jcobi (octave_idx_type n, octave_idx_type n0, octave_idx_type n1, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
145 double alpha, double beta, double *dif1, double *dif2, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
146 double *dif3, double *root) |
3 | 147 { |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
148 assert (n0 == 0 || n0 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
149 assert (n1 == 0 || n1 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
150 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
151 octave_idx_type nt = n + n0 + n1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
152 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
153 assert (nt > 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
154 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
155 // -- first evaluation of coefficients in recursion formulas. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
156 // -- recursion coefficients are stored in dif1 and dif2. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
157 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
158 double ab = alpha + beta; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
159 double ad = beta - alpha; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
160 double ap = beta * alpha; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
161 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
162 dif1[0] = (ad / (ab + 2.0) + 1.0) / 2.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
163 dif2[0] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
164 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
165 if (n >= 2) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
166 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
167 for (octave_idx_type i = 1; i < n; i++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
168 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
169 double z1 = i; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
170 double z = ab + 2 * z1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
171 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
172 dif1[i] = (ab * ad / z / (z + 2.0) + 1.0) / 2.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
173 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
174 if (i == 1) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
175 dif2[i] = (ab + ap + z1) / z / z / (z + 1.0); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
176 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
177 { |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
178 z *= z; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
179 double y = z1 * (ab + z1); |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
180 y *= (ap + y); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
181 dif2[i] = y / z / (z - 1.0); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
182 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
183 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
184 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
185 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
186 // Root determination by Newton method with suppression of previously |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
187 // determined roots. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
188 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
189 double x = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
190 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
191 for (octave_idx_type i = 0; i < n; i++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
192 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
193 bool done = false; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
194 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
195 int k = 0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
196 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
197 while (! done) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
198 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
199 double xd = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
200 double xn = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
201 double xd1 = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
202 double xn1 = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
203 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
204 for (octave_idx_type j = 0; j < n; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
205 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
206 double xp = (dif1[j] - x) * xn - dif2[j] * xd; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
207 double xp1 = (dif1[j] - x) * xn1 - dif2[j] * xd1 - xn; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
208 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
209 xd = xn; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
210 xd1 = xn1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
211 xn = xp; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
212 xn1 = xp1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
213 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
214 |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
215 double zc = 1.0; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
216 double z = xn / xn1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
217 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
218 if (i != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
219 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
220 for (octave_idx_type j = 1; j <= i; j++) |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
221 zc -= z / (x - root[j-1]); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
222 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
223 |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
224 z /= zc; |
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
225 x -= z; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
226 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
227 // Famous last words: 100 iterations should be more than |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
228 // enough in all cases. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
229 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
230 if (++k > 100 || xisnan (z)) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
231 return false; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
232 |
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15018
diff
changeset
|
233 if (std::abs (z) <= 100 * std::numeric_limits<double>::epsilon ()) |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
234 done = true; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
235 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
236 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
237 root[i] = x; |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
238 x += sqrt (std::numeric_limits<double>::epsilon ()); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
239 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
240 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
241 // Add interpolation points at x = 0 and/or x = 1. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
242 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
243 if (n0 != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
244 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
245 for (octave_idx_type i = n; i > 0; i--) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
246 root[i] = root[i-1]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
247 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
248 root[0] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
249 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
250 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
251 if (n1 != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
252 root[nt-1] = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
253 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
254 dif (nt, root, dif1, dif2, dif3); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
255 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
256 return true; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
257 } |
3 | 258 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
259 // Compute derivative weights for orthogonal collocation. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
260 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
261 // See Villadsen and Michelsen, pages 133-134, 419. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
262 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
263 // Input parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
264 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
265 // nd : the dimension of the vectors dif1, dif2, dif3, and root |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
266 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
267 // n : the degree of the jacobi polynomial, (i.e. the number |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
268 // of interior interpolation points) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
269 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
270 // n0 : determines whether x = 0 is included as an |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
271 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
272 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
273 // n0 = 0 ==> x = 0 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
274 // n0 = 1 ==> x = 0 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
275 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
276 // n1 : determines whether x = 1 is included as an |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
277 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
278 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
279 // n1 = 0 ==> x = 1 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
280 // n1 = 1 ==> x = 1 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
281 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
282 // i : the index of the node for which the weights are to be |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
283 // calculated |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
284 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
285 // id : indicator |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
286 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
287 // id = 1 ==> first derivative weights are computed |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
288 // id = 2 ==> second derivative weights are computed |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
289 // id = 3 ==> gaussian weights are computed (in this |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
290 // case, the value of i is irrelevant) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
291 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
292 // Output parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
293 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
294 // dif1 : one dimensional vector containing the first derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
295 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
296 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
297 // dif2 : one dimensional vector containing the second derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
298 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
299 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
300 // dif3 : one dimensional vector containing the third derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
301 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
302 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
303 // vect : one dimensional vector of computed weights |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
304 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
305 static void |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
306 dfopr (octave_idx_type n, octave_idx_type n0, octave_idx_type n1, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
307 octave_idx_type i, octave_idx_type id, double *dif1, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
308 double *dif2, double *dif3, double *root, double *vect) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
309 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
310 assert (n0 == 0 || n0 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
311 assert (n1 == 0 || n1 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
312 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
313 octave_idx_type nt = n + n0 + n1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
314 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
315 assert (nt > 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
316 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
317 assert (id == 1 || id == 2 || id == 3); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
318 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
319 if (id != 3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
320 assert (i >= 0 && i < nt); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
321 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
322 // Evaluate discretization matrices and Gaussian quadrature weights. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
323 // Quadrature weights are normalized to sum to one. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
324 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
325 if (id != 3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
326 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
327 for (octave_idx_type j = 0; j < nt; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
328 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
329 if (j == i) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
330 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
331 if (id == 1) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
332 vect[i] = dif2[i] / dif1[i] / 2.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
333 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
334 vect[i] = dif3[i] / dif1[i] / 3.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
335 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
336 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
337 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
338 double y = root[i] - root[j]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
339 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
340 vect[j] = dif1[i] / dif1[j] / y; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
341 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
342 if (id == 2) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
343 vect[j] = vect[j] * (dif2[i] / dif1[i] - 2.0 / y); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
344 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
345 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
346 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
347 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
348 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
349 double y = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
350 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
351 for (octave_idx_type j = 0; j < nt; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
352 { |
21568
3d60ed163b70
maint: Eliminate bad spacing around '='.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
353 double x = root[j]; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
354 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
355 double ax = x * (1.0 - x); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
356 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
357 if (n0 == 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
358 ax = ax / x / x; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
359 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
360 if (n1 == 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
361 ax = ax / (1.0 - x) / (1.0 - x); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
362 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
363 vect[j] = ax / (dif1[j] * dif1[j]); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
364 |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
365 y += vect[j]; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
366 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
367 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
368 for (octave_idx_type j = 0; j < nt; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
369 vect[j] = vect[j] / y; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
370 } |
3 | 371 } |
372 | |
373 // Error handling. | |
374 | |
375 void | |
376 CollocWt::error (const char* msg) | |
377 { | |
20726
25d676f9619c
Preface error() messages with name of function when possible.
Rik <rik@octave.org>
parents:
20230
diff
changeset
|
378 (*current_liboctave_error_handler) ("CollocWt: fatal error '%s'", msg); |
3 | 379 } |
380 | |
381 CollocWt& | |
382 CollocWt::set_left (double val) | |
383 { | |
384 if (val >= rb) | |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
385 error ("CollocWt: left bound greater than right bound"); |
3 | 386 |
387 lb = val; | |
388 initialized = 0; | |
389 return *this; | |
390 } | |
391 | |
392 CollocWt& | |
393 CollocWt::set_right (double val) | |
394 { | |
395 if (val <= lb) | |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
396 error ("CollocWt: right bound less than left bound"); |
3 | 397 |
398 rb = val; | |
399 initialized = 0; | |
400 return *this; | |
401 } | |
402 | |
403 void | |
404 CollocWt::init (void) | |
405 { | |
1360 | 406 // Check for possible errors. |
3 | 407 |
408 double wid = rb - lb; | |
409 if (wid <= 0.0) | |
227 | 410 { |
20726
25d676f9619c
Preface error() messages with name of function when possible.
Rik <rik@octave.org>
parents:
20230
diff
changeset
|
411 error ("CollocWt: width less than or equal to zero"); |
227 | 412 return; |
413 } | |
3 | 414 |
5275 | 415 octave_idx_type nt = n + inc_left + inc_right; |
1870 | 416 |
3 | 417 if (nt < 0) |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
418 error ("CollocWt: total number of collocation points less than zero"); |
3 | 419 else if (nt == 0) |
420 return; | |
421 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
422 Array<double> dif1 (dim_vector (nt, 1)); |
1938 | 423 double *pdif1 = dif1.fortran_vec (); |
424 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
425 Array<double> dif2 (dim_vector (nt, 1)); |
1938 | 426 double *pdif2 = dif2.fortran_vec (); |
427 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
428 Array<double> dif3 (dim_vector (nt, 1)); |
1938 | 429 double *pdif3 = dif3.fortran_vec (); |
430 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
431 Array<double> vect (dim_vector (nt, 1)); |
1938 | 432 double *pvect = vect.fortran_vec (); |
3 | 433 |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
434 r.resize (nt, 1); |
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
435 q.resize (nt, 1); |
3 | 436 A.resize (nt, nt); |
437 B.resize (nt, nt); | |
438 | |
439 double *pr = r.fortran_vec (); | |
440 | |
1360 | 441 // Compute roots. |
3 | 442 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
443 if (! jcobi (n, inc_left, inc_right, Alpha, Beta, pdif1, pdif2, pdif3, pr)) |
21056
d48fdf3a8c0c
maint: one more addition to cset 5e00ed38a58b.
Rik <rik@octave.org>
parents:
21055
diff
changeset
|
444 error ("jcobi: newton iteration failed"); |
3 | 445 |
5275 | 446 octave_idx_type id; |
3 | 447 |
1360 | 448 // First derivative weights. |
3 | 449 |
450 id = 1; | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
451 for (octave_idx_type i = 0; i < nt; i++) |
3 | 452 { |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
453 dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); |
3 | 454 |
5275 | 455 for (octave_idx_type j = 0; j < nt; j++) |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
456 A(i,j) = vect(j); |
3 | 457 } |
458 | |
1360 | 459 // Second derivative weights. |
3 | 460 |
461 id = 2; | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
462 for (octave_idx_type i = 0; i < nt; i++) |
3 | 463 { |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
464 dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); |
3 | 465 |
5275 | 466 for (octave_idx_type j = 0; j < nt; j++) |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
467 B(i,j) = vect(j); |
3 | 468 } |
469 | |
1360 | 470 // Gaussian quadrature weights. |
3 | 471 |
472 id = 3; | |
473 double *pq = q.fortran_vec (); | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
474 dfopr (n, inc_left, inc_right, id, id, pdif1, pdif2, pdif3, pr, pq); |
3 | 475 |
476 initialized = 1; | |
477 } | |
478 | |
3504 | 479 std::ostream& |
480 operator << (std::ostream& os, const CollocWt& a) | |
3 | 481 { |
482 if (a.left_included ()) | |
483 os << "left boundary is included\n"; | |
484 else | |
485 os << "left boundary is not included\n"; | |
486 | |
487 if (a.right_included ()) | |
488 os << "right boundary is included\n"; | |
489 else | |
490 os << "right boundary is not included\n"; | |
491 | |
492 os << "\n"; | |
493 | |
494 os << a.Alpha << " " << a.Beta << "\n\n" | |
495 << a.r << "\n\n" | |
496 << a.q << "\n\n" | |
497 << a.A << "\n" | |
498 << a.B << "\n"; | |
499 | |
500 return os; | |
501 } |