Mercurial > hg > octave-nkf
annotate src/DLD-FUNCTIONS/__lin_interpn__.cc @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author | David Bateman <dbateman@free.fr> |
---|---|
date | Sun, 27 Apr 2008 22:34:17 +0200 |
parents | a938cd7869b2 |
children | 25bc2d31e1bf |
rev | line source |
---|---|
6702 | 1 /* |
2 | |
3 Copyright (C) 2007 Alexander Barth | |
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. | |
6702 | 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/>. | |
6702 | 20 |
21 */ | |
22 | |
23 #ifdef HAVE_CONFIG_H | |
24 #include <config.h> | |
25 #endif | |
26 | |
27 #include "dNDArray.h" | |
28 | |
29 #include "defun-dld.h" | |
30 #include "error.h" | |
31 #include "oct-obj.h" | |
32 | |
33 // equivalent to isvector.m | |
34 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
35 template <class T> |
6702 | 36 bool |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
37 isvector (const T& array) |
6702 | 38 { |
39 const dim_vector dv = array.dims (); | |
40 return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1); | |
41 } | |
42 | |
43 // lookup a value in a sorted table (lookup.m) | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
44 template <class T> |
6702 | 45 octave_idx_type |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
46 lookup (const T *x, octave_idx_type n, T y) |
6702 | 47 { |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
48 octave_idx_type j; |
6702 | 49 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
50 if (x[0] < x[n-1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
51 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
52 // increasing x |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
53 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
54 if (y > x[n-1] || y < x[0]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
55 return -1; |
6702 | 56 |
57 #ifdef EXHAUSTIF | |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
58 for (j = 0; j < n - 1; j++) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
59 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
60 if (x[j] <= y && y <= x[j+1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
61 return j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
62 } |
6702 | 63 #else |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
64 octave_idx_type j0 = 0; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
65 octave_idx_type j1 = n - 1; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
66 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
67 while (true) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
68 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
69 j = (j0+j1)/2; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
70 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
71 if (y <= x[j+1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
72 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
73 if (x[j] <= y) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
74 return j; |
6702 | 75 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
76 j1 = j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
77 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
78 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
79 if (x[j] <= y) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
80 j0 = j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
81 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
82 #endif |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
83 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
84 else |
6702 | 85 { |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
86 // decreasing x |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
87 // previous code with x -> -x and y -> -y |
6702 | 88 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
89 if (y > x[0] || y < x[n-1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
90 return -1; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
91 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
92 #ifdef EXHAUSTIF |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
93 for (j = 0; j < n - 1; j++) |
6702 | 94 { |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
95 if (x[j+1] <= y && y <= x[j]) |
6702 | 96 return j; |
97 } | |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
98 #else |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
99 octave_idx_type j0 = 0; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
100 octave_idx_type j1 = n - 1; |
6702 | 101 |
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
102 while (true) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
103 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
104 j = (j0+j1)/2; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
105 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
106 if (y >= x[j+1]) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
107 { |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
108 if (x[j] >= y) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
109 return j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
110 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
111 j1 = j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
112 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
113 |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
114 if (x[j] >= y) |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
115 j0 = j; |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
116 } |
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7016
diff
changeset
|
117 #endif |
6702 | 118 } |
119 } | |
120 | |
121 // n-dimensional linear interpolation | |
122 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
123 template <class T> |
6702 | 124 void |
125 lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale, | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
126 octave_idx_type Ni, T extrapval, const T **x, |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
127 const T *v, const T **y, T *vi) |
6702 | 128 { |
129 bool out = false; | |
130 int bit; | |
131 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
132 OCTAVE_LOCAL_BUFFER (T, coef, 2*n); |
6702 | 133 OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n); |
134 | |
135 // loop over all points | |
136 for (octave_idx_type m = 0; m < Ni; m++) | |
137 { | |
138 // loop over all dimensions | |
139 for (int i = 0; i < n; i++) | |
140 { | |
141 index[i] = lookup (x[i], size[i], y[i][m]); | |
142 out = index[i] == -1; | |
143 | |
144 if (out) | |
145 break; | |
146 else | |
147 { | |
148 octave_idx_type j = index[i]; | |
149 coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]); | |
150 coef[2*i] = 1 - coef[2*i+1]; | |
151 } | |
152 } | |
153 | |
154 | |
155 if (out) | |
156 vi[m] = extrapval; | |
157 else | |
158 { | |
159 vi[m] = 0; | |
160 | |
161 // loop over all corners of hypercube (1<<n = 2^n) | |
162 for (int i = 0; i < (1 << n); i++) | |
163 { | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
164 T c = 1; |
6702 | 165 octave_idx_type l = 0; |
166 | |
167 // loop over all dimensions | |
168 for (int j = 0; j < n; j++) | |
169 { | |
170 // test if the jth bit in i is set | |
171 bit = i >> j & 1; | |
172 l += scale[j] * (index[j] + bit); | |
173 c *= coef[2*j+bit]; | |
174 } | |
175 | |
176 vi[m] += c * v[l]; | |
177 } | |
178 } | |
179 } | |
180 } | |
181 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
182 template <class T, class M> |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
183 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
184 lin_interpn (int n, M *X, const M V, M *Y) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
185 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
186 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
187 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
188 M Vi = M (Y[0].dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
189 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
190 OCTAVE_LOCAL_BUFFER (const T *, y, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
191 OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
192 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
193 for (int i = 0; i < n; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
194 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
195 y[i] = Y[i].data (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
196 size[i] = V.dims()(i); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
197 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
198 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
199 OCTAVE_LOCAL_BUFFER (const T *, x, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
200 OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
201 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
202 const T *v = V.data (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
203 T *vi = Vi.fortran_vec (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
204 octave_idx_type Ni = Vi.numel (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
205 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
206 T extrapval = octave_NA; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
207 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
208 // offset in memory of each dimension |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
209 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
210 scale[0] = 1; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
211 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
212 for (int i = 1; i < n; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
213 scale[i] = scale[i-1] * size[i-1]; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
214 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
215 // tests if X[0] is a vector, if yes, assume that all elements of X are |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
216 // in the ndgrid format. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
217 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
218 if (! isvector (X[0])) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
219 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
220 for (int i = 0; i < n; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
221 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
222 if (X[i].dims () != V.dims ()) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
223 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
224 error ("interpn: incompatible size of argument number %d", i+1); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
225 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
226 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
227 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
228 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
229 M tmp = M (dim_vector (size[i], 1)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
230 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
231 for (octave_idx_type j = 0; j < size[i]; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
232 tmp(j) = X[i](scale[i]*j); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
233 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
234 X[i] = tmp; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
235 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
236 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
237 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
238 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
239 for (int i = 0; i < n; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
240 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
241 if (! isvector (X[i]) && X[i].numel () != size[i]) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
242 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
243 error ("interpn: incompatible size of argument number %d", i+1); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
244 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
245 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
246 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
247 x[i] = X[i].data (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
248 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
249 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
250 lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
251 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
252 retval = Vi; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
253 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
254 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
255 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
256 |
6945 | 257 // Perform @var{n}-dimensional interpolation. Each element of then |
258 // @var{n}-dimensional array @var{v} represents a value at a location | |
259 // given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters | |
260 // @var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional | |
261 // arrays of the same size as the array @var{v} in the \"ndgrid\" format | |
262 // or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are | |
263 // all @var{n}-dimensional arrays of the same size and represent the | |
264 // points at which the array @var{vi} is interpolated. | |
265 // | |
266 //This function only performs linear interpolation. | |
267 | |
6702 | 268 DEFUN_DLD (__lin_interpn__, args, , |
269 "-*- texinfo -*-\n\ | |
270 @deftypefn {Loadable Function} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\ | |
6945 | 271 Undocumented internal function.\n\ |
6702 | 272 @end deftypefn") |
273 { | |
274 octave_value retval; | |
275 | |
276 int nargin = args.length (); | |
277 | |
278 if (nargin < 2 || nargin % 2 == 0) | |
279 { | |
280 print_usage (); | |
281 return retval; | |
282 } | |
283 | |
284 // dimension of the problem | |
285 int n = (nargin-1)/2; | |
286 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
287 if (args(n).is_single_type()) |
6702 | 288 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
289 OCTAVE_LOCAL_BUFFER (FloatNDArray, X, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
290 OCTAVE_LOCAL_BUFFER (FloatNDArray, Y, n); |
6702 | 291 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
292 const FloatNDArray V = args(n).float_array_value (); |
6702 | 293 |
294 if (error_state) | |
295 { | |
296 print_usage (); | |
297 return retval; | |
298 } | |
299 | |
300 for (int i = 0; i < n; i++) | |
301 { | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
302 X[i] = args(i).float_array_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
303 Y[i] = args(n+i+1).float_array_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
304 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
305 if (error_state) |
6702 | 306 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
307 print_usage (); |
6702 | 308 return retval; |
309 } | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
310 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
311 if (Y[0].dims () != Y[i].dims ()) |
6702 | 312 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
313 error ("interpn: incompatible size of argument number %d", n+i+2); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
314 return retval; |
6702 | 315 } |
316 } | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
317 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
318 retval = lin_interpn<float, FloatNDArray> (n, X, V, Y); |
6702 | 319 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
320 else |
6702 | 321 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
322 OCTAVE_LOCAL_BUFFER (NDArray, X, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
323 OCTAVE_LOCAL_BUFFER (NDArray, Y, n); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
324 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
325 const NDArray V = args(n).array_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
326 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
327 if (error_state) |
6702 | 328 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
329 print_usage (); |
6702 | 330 return retval; |
331 } | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
332 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
333 for (int i = 0; i < n; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
334 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
335 X[i] = args(i).array_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
336 Y[i] = args(n+i+1).array_value (); |
6702 | 337 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
338 if (error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
339 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
340 print_usage (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
341 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
342 } |
6702 | 343 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
344 if (Y[0].dims () != Y[i].dims ()) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
345 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
346 error ("interpn: incompatible size of argument number %d", n+i+2); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
347 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
348 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
349 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
350 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
351 retval = lin_interpn<double, NDArray> (n, X, V, Y); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7561
diff
changeset
|
352 } |
6702 | 353 |
354 return retval; | |
355 } |