1 /***************************************************************************
2 * data.cpp is part of Math Graphic Library
3 * Copyright (C) 2007-2016 Alexey Balakin <mathgl.abalakin@gmail.ru> *
4 * *
5 * This program is free software; you can redistribute it and/or modify *
6 * it under the terms of the GNU Lesser General Public License as *
7 * published by the Free Software Foundation; either version 3 of the *
8 * License, or (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU Lesser General Public *
16 * License along with this program; if not, write to the *
17 * Free Software Foundation, Inc., *
18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19 ***************************************************************************/
20 #include <time.h>
21 #include "mgl2/data.h"
22 #include "mgl2/datac.h"
23 #include "mgl2/eval.h"
24 #include "mgl2/thread.h"
25 #include "interp.hpp"
26
27 MGL_EXPORT int mglNumThr=0;
28 MGL_EXPORT std::vector<mglDataA*> mglDataList;
29
30 //-----------------------------------------------------------------------------
31 #if MGL_HAVE_PTHREAD
32 #ifdef WIN32
33 #include <windows.h>
34 #include <process.h>
35 #elif defined(__APPLE__) || defined (__FreeBSD__)
36 #include <sys/sysctl.h>
37 #elif defined(unix) || defined(__unix) || defined(__unix__)
38 #include <sys/sysinfo.h>
39 #endif
mgl_set_num_thr(int n)40 void MGL_EXPORT mgl_set_num_thr(int n)
41 {
42 #ifdef WIN32
43 SYSTEM_INFO systemInfo;
44 GetSystemInfo(&systemInfo);
45 mglNumThr = n>0 ? n : systemInfo.dwNumberOfProcessors;
46 #elif defined (__APPLE__) || defined(__FreeBSD__)
47 int numProcessors = 1;
48 size_t size = sizeof(numProcessors);
49 sysctlbyname("hw.ncpu", &numProcessors, &size, NULL, 0);
50 mglNumThr = n>0 ? n : numProcessors;
51 #else
52 mglNumThr = n>0 ? n : get_nprocs_conf();
53 #endif
54 }
55 #else
mgl_set_num_thr(int)56 void MGL_EXPORT mgl_set_num_thr(int) { mglNumThr = 1; }
57 #endif
mgl_set_num_thr_(int * n)58 void MGL_EXPORT mgl_set_num_thr_(int *n) { mgl_set_num_thr(*n); }
59 //-----------------------------------------------------------------------------
mglStartThread(void * (* func)(void *),void (* post)(mglThreadD *,mreal *),long n,mreal * a,const mreal * b,const mreal * c,const long * p,const void * v,const mreal * d,const mreal * e,const char * s)60 void MGL_EXPORT mglStartThread(void *(*func)(void *), void (*post)(mglThreadD *,mreal *), long n,
61 mreal *a, const mreal *b, const mreal *c, const long *p,
62 const void *v, const mreal *d, const mreal *e, const char *s)
63 {
64 if(!func) return;
65 #if MGL_HAVE_PTHREAD
66 if(mglNumThr<1) mgl_set_num_thr(0);
67 if(mglNumThr>1)
68 {
69 pthread_t *tmp=new pthread_t[mglNumThr];
70 mglThreadD *par=new mglThreadD[mglNumThr];
71 for(long i=0;i<mglNumThr;i++) // put parameters into the structure
72 { par[i].n=n; par[i].a=a; par[i].b=b; par[i].c=c; par[i].d=d;
73 par[i].p=p; par[i].v=v; par[i].s=s; par[i].e=e; par[i].id=i; }
74 for(long i=0;i<mglNumThr;i++) pthread_create(tmp+i, 0, func, par+i);
75 for(long i=0;i<mglNumThr;i++) pthread_join(tmp[i], 0);
76 if(post) post(par,a);
77 delete []tmp; delete []par;
78 }
79 else
80 #endif
81 {
82 mglNumThr = 1;
83 mglThreadD par;
84 par.n=n; par.a=a; par.b=b; par.c=c; par.d=d;
85 par.p=p; par.v=v; par.s=s; par.e=e; par.id=0;
86 func(&par);
87 if(post) post(&par,a);
88 }
89 }
90 //-----------------------------------------------------------------------------
mglStartThreadV(void * (* func)(void *),long n,mreal * a,const void * b,const void * c,const long * p,const void * v,const mreal * d)91 void MGL_EXPORT mglStartThreadV(void *(*func)(void *), long n, mreal *a, const void *b,
92 const void *c, const long *p, const void *v, const mreal *d)
93 {
94 if(!func) return;
95 #if MGL_HAVE_PTHREAD
96 if(mglNumThr<1) mgl_set_num_thr(0);
97 if(mglNumThr>1)
98 {
99 pthread_t *tmp=new pthread_t[mglNumThr];
100 mglThreadV *par=new mglThreadV[mglNumThr];
101 for(long i=0;i<mglNumThr;i++) // put parameters into the structure
102 { par[i].n=n; par[i].a=a; par[i].b=b; par[i].c=c; par[i].d=d;
103 par[i].p=p; par[i].v=v; par[i].id=i;par[i].aa=0; }
104 for(long i=0;i<mglNumThr;i++) pthread_create(tmp+i, 0, func, par+i);
105 for(long i=0;i<mglNumThr;i++) pthread_join(tmp[i], 0);
106 delete []tmp; delete []par;
107 }
108 else
109 #endif
110 {
111 mglNumThr = 1;
112 mglThreadV par;
113 par.n=n; par.a=a; par.b=b; par.c=c; par.d=d;
114 par.p=p; par.v=v; par.id=0; par.aa=0;
115 func(&par);
116 }
117 }
118 //-----------------------------------------------------------------------------
mglSpline3s(const mreal * a,long nx,long ny,long nz,mreal x,mreal y,mreal z)119 mreal MGL_EXPORT_PURE mglSpline3s(const mreal *a, long nx, long ny, long nz, mreal x, mreal y, mreal z)
120 { return mglSpline3st<mreal>(a,nx,ny,nz,x,y,z); }
121 //-----------------------------------------------------------------------------
mglSpline3(const mreal * a,long nx,long ny,long nz,mreal x,mreal y,mreal z,mreal * dx,mreal * dy,mreal * dz)122 mreal MGL_EXPORT mglSpline3(const mreal *a, long nx, long ny, long nz, mreal x, mreal y, mreal z,mreal *dx, mreal *dy, mreal *dz)
123 { return mglSpline3t<mreal>(a,nx,ny,nz,x,y,z,dx,dy,dz); }
124 //-----------------------------------------------------------------------------
mglLinear(const mreal * a,long nx,long ny,long nz,mreal x,mreal y,mreal z)125 mreal MGL_EXPORT_PURE mglLinear(const mreal *a, long nx, long ny, long nz, mreal x, mreal y, mreal z)
126 { return mglLineart<mreal>(a,nx,ny,nz,x,y,z); }
127 //-----------------------------------------------------------------------------
mgl_ipow(double x,int n)128 double MGL_EXPORT_CONST mgl_ipow(double x,int n)
129 {
130 double t;
131 if(n==2) return x*x;
132 if(n==1) return x;
133 if(n<0) return 1./mgl_ipow(x,-n);
134 if(n==0) return 1;
135 t = mgl_ipow(x,n/2); t = t*t;
136 if(n%2==1) t *= x;
137 return t;
138 }
mgl_ipow_(mreal * x,int * n)139 double MGL_EXPORT_PURE mgl_ipow_(mreal *x,int *n) { return mgl_ipow(*x,*n); }
140 //-----------------------------------------------------------------------------
mgl_get_time(const char * time,const char * fmt)141 double mgl_get_time(const char *time, const char *fmt)
142 {
143 #if !defined(WIN32)
144 tm t;
145 strptime(time,fmt,&t);
146 return timegm(&t);
147 #else
148 return NAN;
149 #endif
150 }
mgl_get_time_(const char * time,const char * fmt,int l,int m)151 double mgl_get_time_(const char *time, const char *fmt,int l,int m)
152 { char *s=new char[l+1]; memcpy(s,time,l); s[l]=0;
153 char *f=new char[m+1]; memcpy(f,fmt,m); f[m]=0;
154 double t=mgl_get_time(s,f); delete []s; delete []f; return t; }
155 //-----------------------------------------------------------------------------
mgl_smth_x(void * par)156 static void *mgl_smth_x(void *par)
157 {
158 mglThreadD *t=(mglThreadD *)par;
159 long nx=t->p[0], kind=t->p[2];
160 mreal *b=t->a, delta=t->c[0];
161 const mreal *a=t->b;
162 if(kind>0)
163 #if !MGL_HAVE_PTHREAD
164 #pragma omp parallel for
165 #endif
166 for(long i=t->id;i<t->n;i+=mglNumThr)
167 {
168 long j = i%nx;
169 double s0=0,s1=0,s2=0,sf=0,sxf=0;
170 for(long k=-kind;k<=kind;k++)
171 if(j+k>=0 && j+k<nx && mgl_isnum(a[i+k]))
172 {
173 s0+=1; s1+=k; s2+=k*k;
174 sf+=a[i+k]; sxf+=a[i+k]*k;
175 }
176 double d = s0*s2-s1*s1;
177 b[i] = d?(s2*sf - s1*sxf)/d:a[i];
178 }
179 else if(kind<-1)
180 #if !MGL_HAVE_PTHREAD
181 #pragma omp parallel for
182 #endif
183 for(long i=t->id;i<t->n;i+=mglNumThr)
184 {
185 long j = i%nx;
186 double s0=0,s1=0,s2=0,s3=0,s4=0,sf=0,sxf=0,sxxf=0;
187 for(long k=kind;k<=-kind;k++)
188 if(j+k>=0 && j+k<nx && mgl_isnum(a[i+k]))
189 {
190 long k2 = k*k, i0 = i+k;
191 s0+=1; s1+=k; s2+=k2; s3+=k2*k; s4+=k2*k2;
192 sf+=a[i0]; sxf+=a[i0]*k; sxxf+=a[i0]*k2;
193 }
194 double d = (s0*s4-s2*s2)*s2 + (s2*s3-s4*s1)*s1 + (s1*s2-s0*s3)*s3;
195 b[i] = d?( (s1*s3-s2*s2)*sxxf + (s2*s3-s1*s4)*sxf + (s2*s4-s3*s3)*sf )/d:a[i];
196 }
197 else if(kind==0)
198 #if !MGL_HAVE_PTHREAD
199 #pragma omp parallel for
200 #endif
201 for(long i=t->id;i<t->n;i+=mglNumThr)
202 {
203 long j = i%nx;
204 mreal v = (j>0 && j<nx-1) ? (a[i-1]+a[i+1])/2 : NAN;
205 b[i] = v<a[i]?v:a[i];
206 }
207 else if(kind==-1)
208 #if !MGL_HAVE_PTHREAD
209 #pragma omp parallel for
210 #endif
211 for(long i=t->id;i<t->n;i+=mglNumThr)
212 {
213 long j = i%nx;
214 mreal v = (j>0 && j<nx-1) ? (a[i-1]+a[i+1])/2 : NAN;
215 b[i] = v>a[i]?v:a[i];
216 }
217 if(delta>0)
218 #if !MGL_HAVE_PTHREAD
219 #pragma omp parallel for
220 #endif
221 for(long i=t->id;i<t->n;i+=mglNumThr)
222 {
223 double ab = fabs(a[i]-b[i]);
224 if(ab>delta) b[i] = a[i]+(delta/ab)*(b[i]-a[i]);
225 }
226 return 0;
227 }
mgl_smth_y(void * par)228 static void *mgl_smth_y(void *par)
229 {
230 mglThreadD *t=(mglThreadD *)par;
231 long nx=t->p[0],ny=t->p[1], kind=t->p[2];
232 mreal *b=t->a, delta=t->c[0];
233 const mreal *a=t->b;
234 if(kind>0)
235 #if !MGL_HAVE_PTHREAD
236 #pragma omp parallel for
237 #endif
238 for(long i=t->id;i<t->n;i+=mglNumThr)
239 {
240 long j = (i/nx)%ny;
241 double s0=0,s1=0,s2=0,sf=0,sxf=0;
242 for(long k=-kind;k<=kind;k++)
243 if(j+k>=0 && j+k<ny && mgl_isnum(a[i+k*nx]))
244 {
245 s0+=1; s1+=k; s2+=k*k;
246 sf+=a[i+k*nx]; sxf+=a[i+k*nx]*k;
247 }
248 double d = s0*s2-s1*s1;
249 b[i] = d?(s2*sf - s1*sxf)/d:a[i];
250 }
251 else if(kind<-1)
252 #if !MGL_HAVE_PTHREAD
253 #pragma omp parallel for
254 #endif
255 for(long i=t->id;i<t->n;i+=mglNumThr)
256 {
257 long j = (i/nx)%ny;
258 double s0=0,s1=0,s2=0,s3=0,s4=0,sf=0,sxf=0,sxxf=0;
259 for(long k=kind;k<=-kind;k++)
260 if(j+k>=0 && j+k<ny && mgl_isnum(a[i+k*nx]))
261 {
262 long k2 = k*k, i0 = i+k*nx;
263 s0+=1; s1+=k; s2+=k2; s3+=k2*k; s4+=k2*k2;
264 sf+=a[i0]; sxf+=a[i0]*k; sxxf+=a[i0]*k2;
265 }
266 double d = (s0*s4-s2*s2)*s2 + (s2*s3-s4*s1)*s1 + (s1*s2-s0*s3)*s3;
267 b[i] = d?( (s1*s3-s2*s2)*sxxf + (s2*s3-s1*s4)*sxf + (s2*s4-s3*s3)*sf )/d:a[i];
268 }
269 else if(kind==0)
270 #if !MGL_HAVE_PTHREAD
271 #pragma omp parallel for
272 #endif
273 for(long i=t->id;i<t->n;i+=mglNumThr)
274 {
275 long j = (i/nx)%ny;
276 mreal v = (j>0 && j<ny-1) ? (a[i-nx]+a[i+nx])/2 : NAN;
277 b[i] = v<a[i]?v:a[i];
278 }
279 else if(kind==-1)
280 #if !MGL_HAVE_PTHREAD
281 #pragma omp parallel for
282 #endif
283 for(long i=t->id;i<t->n;i+=mglNumThr)
284 {
285 long j = (i/nx)%ny;
286 mreal v = (j>0 && j<ny-1) ? (a[i-nx]+a[i+nx])/2 : NAN;
287 b[i] = v>a[i]?v:a[i];
288 }
289 if(delta>0)
290 #if !MGL_HAVE_PTHREAD
291 #pragma omp parallel for
292 #endif
293 for(long i=t->id;i<t->n;i+=mglNumThr)
294 {
295 double ab = fabs(a[i]-b[i]);
296 if(ab>delta) b[i] = a[i]+(delta/ab)*(b[i]-a[i]);
297 }
298 return 0;
299 }
mgl_smth_z(void * par)300 static void *mgl_smth_z(void *par)
301 {
302 mglThreadD *t=(mglThreadD *)par;
303 long nn=t->p[0]*t->p[1], nz=t->n/nn, kind=t->p[2];
304 mreal *b=t->a, delta=t->c[0];
305 const mreal *a=t->b;
306 if(kind>1)
307 #if !MGL_HAVE_PTHREAD
308 #pragma omp parallel for
309 #endif
310 for(long i=t->id;i<t->n;i+=mglNumThr)
311 {
312 long j = i/nn;
313 double s0=0,s1=0,s2=0,sf=0,sxf=0;
314 for(long k=-kind;k<=kind;k++)
315 if(j+k>=0 && j+k<nz && mgl_isnum(a[i+k*nn]))
316 {
317 s0+=1; s1+=k; s2+=k*k;
318 sf+=a[i+k*nn]; sxf+=a[i+k*nn]*k;
319 }
320 double d = s0*s2-s1*s1;
321 b[i] = d?(s2*sf - s1*sxf)/d:a[i];
322 }
323 else if(kind<-1)
324 #if !MGL_HAVE_PTHREAD
325 #pragma omp parallel for
326 #endif
327 for(long i=t->id;i<t->n;i+=mglNumThr)
328 {
329 long j = i/nn;
330 double s0=0,s1=0,s2=0,s3=0,s4=0,sf=0,sxf=0,sxxf=0;
331 for(long k=kind;k<=-kind;k++)
332 if(j+k>=0 && j+k<nz && mgl_isnum(a[i+k*nn]))
333 {
334 long k2 = k*k, i0 = i+k*nn;
335 s0+=1; s1+=k; s2+=k2; s3+=k2*k; s4+=k2*k2;
336 sf+=a[i0]; sxf+=a[i0]*k; sxxf+=a[i0]*k2;
337 }
338 double d = (s0*s4-s2*s2)*s2 + (s2*s3-s4*s1)*s1 + (s1*s2-s0*s3)*s3;
339 b[i] = d?( (s1*s3-s2*s2)*sxxf + (s2*s3-s1*s4)*sxf + (s2*s4-s3*s3)*sf )/d:a[i];
340 }
341 else if(kind==0)
342 #if !MGL_HAVE_PTHREAD
343 #pragma omp parallel for
344 #endif
345 for(long i=t->id;i<t->n;i+=mglNumThr)
346 {
347 long j = i/nn;
348 mreal v = (j>0 && j<nz-1) ? (a[i-nn]+a[i+nn])/2 : NAN;
349 b[i] = v<a[i]?v:a[i];
350 }
351 else if(kind==-1)
352 #if !MGL_HAVE_PTHREAD
353 #pragma omp parallel for
354 #endif
355 for(long i=t->id;i<t->n;i+=mglNumThr)
356 {
357 long j = i/nn;
358 mreal v = (j>0 && j<nz-1) ? (a[i-nn]+a[i+nn])/2 : NAN;
359 b[i] = v>a[i]?v:a[i];
360 }
361 if(delta>0)
362 #if !MGL_HAVE_PTHREAD
363 #pragma omp parallel for
364 #endif
365 for(long i=t->id;i<t->n;i+=mglNumThr)
366 {
367 double ab = fabs(a[i]-b[i]);
368 if(ab>delta) b[i] = a[i]+(delta/ab)*(b[i]-a[i]);
369 }
370 return 0;
371 }
mgl_data_smooth(HMDT d,const char * dirs,mreal delta)372 void MGL_EXPORT mgl_data_smooth(HMDT d, const char *dirs, mreal delta)
373 {
374 long Type = -2;
375 if(mglchr(dirs,'0')) return;
376 bool xdir=mglchr(dirs,'x'), ydir=mglchr(dirs,'y'), zdir=mglchr(dirs,'z');
377 if(!xdir && !ydir && !zdir) xdir=ydir=zdir=true;
378 if(mglchr(dirs,'d'))
379 {
380 if(mglchr(dirs,'1')) Type = 1;
381 if(mglchr(dirs,'2')) Type = 2;
382 if(mglchr(dirs,'3')) Type = 3;
383 if(mglchr(dirs,'4')) Type = 4;
384 if(mglchr(dirs,'5')) Type = 5;
385 if(mglchr(dirs,'6')) Type = 6;
386 if(mglchr(dirs,'7')) Type = 7;
387 if(mglchr(dirs,'8')) Type = 8;
388 if(mglchr(dirs,'9')) Type = 9;
389 }
390 else if(mglchr(dirs,'p'))
391 {
392 if(mglchr(dirs,'2')) Type = -2;
393 if(mglchr(dirs,'3')) Type = -3;
394 if(mglchr(dirs,'4')) Type = -4;
395 if(mglchr(dirs,'5')) Type = -5;
396 if(mglchr(dirs,'6')) Type = -6;
397 if(mglchr(dirs,'7')) Type = -7;
398 if(mglchr(dirs,'8')) Type = -8;
399 if(mglchr(dirs,'9')) Type = -9;
400 }
401 else
402 {
403 if(mglchr(dirs,'1')) return;
404 if(mglchr(dirs,'3')) Type = 1;
405 if(mglchr(dirs,'5')) Type = 2;
406 }
407 if(mglchr(dirs,'_')) Type = 0;
408 if(mglchr(dirs,'^')) Type = -1;
409 long nx=d->nx,ny=d->ny,nz=d->nz;
410 // if(Type == SMOOTH_NONE) return;
411 long p[3]={nx,ny,Type};
412 mreal *b = new mreal[nx*ny*nz],dd=delta;
413 // ����������� �� x
414 memset(b,0,nx*ny*nz*sizeof(mreal));
415 if(nx>4 && xdir)
416 {
417 mglStartThread(mgl_smth_x,0,nx*ny*nz,b,d->a,&dd,p);
418 memcpy(d->a,b,nx*ny*nz*sizeof(mreal));
419 memset(b,0,nx*ny*nz*sizeof(mreal));
420 }
421 if(ny>4 && ydir)
422 {
423 mglStartThread(mgl_smth_y,0,nx*ny*nz,b,d->a,&dd,p);
424 memcpy(d->a,b,nx*ny*nz*sizeof(mreal));
425 memset(b,0,nx*ny*nz*sizeof(mreal));
426 }
427 if(nz>4 && zdir)
428 {
429 mglStartThread(mgl_smth_z,0,nx*ny*nz,b,d->a,&dd,p);
430 memcpy(d->a,b,nx*ny*nz*sizeof(mreal));
431 }
432 delete []b;
433 }
mgl_data_smooth_(uintptr_t * d,const char * dir,mreal * delta,int l)434 void MGL_EXPORT mgl_data_smooth_(uintptr_t *d, const char *dir, mreal *delta,int l)
435 { char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
436 mgl_data_smooth(_DT_,s,*delta); delete []s; }
437 //-----------------------------------------------------------------------------
mgl_csum_z(void * par)438 static void *mgl_csum_z(void *par)
439 {
440 mglThreadD *t=(mglThreadD *)par;
441 long nz=t->p[2], nn=t->n;
442 mreal *b=t->a;
443 const mreal *a=t->b;
444 #if !MGL_HAVE_PTHREAD
445 #pragma omp parallel for
446 #endif
447 for(long i=t->id;i<nn;i+=mglNumThr)
448 {
449 b[i] = a[i];
450 for(long j=1;j<nz;j++) b[i+j*nn] = b[i+j*nn-nn] + a[i+j*nn];
451 }
452 return 0;
453 }
mgl_csum_y(void * par)454 static void *mgl_csum_y(void *par)
455 {
456 mglThreadD *t=(mglThreadD *)par;
457 long nx=t->p[0], ny=t->p[1], nn=t->n;
458 mreal *b=t->a;
459 const mreal *a=t->b;
460 #if !MGL_HAVE_PTHREAD
461 #pragma omp parallel for
462 #endif
463 for(long i=t->id;i<nn;i+=mglNumThr)
464 {
465 long k = (i%nx)+nx*ny*(i/nx); b[k] = a[k];
466 for(long j=1;j<ny;j++) b[k+j*nx] = b[k+j*nx-nx] + a[k+nx*j];
467 }
468 return 0;
469 }
mgl_csum_x(void * par)470 static void *mgl_csum_x(void *par)
471 {
472 mglThreadD *t=(mglThreadD *)par;
473 long nx=t->p[0], nn=t->n;
474 mreal *b=t->a;
475 const mreal *a=t->b;
476 #if !MGL_HAVE_PTHREAD
477 #pragma omp parallel for
478 #endif
479 for(long i=t->id;i<nn;i+=mglNumThr)
480 {
481 long k = i*nx; b[k] = a[k];
482 for(long j=1;j<nx;j++) b[j+k] = b[j+k-1] + a[j+k];
483 }
484 return 0;
485 }
mgl_data_cumsum(HMDT d,const char * dir)486 void MGL_EXPORT mgl_data_cumsum(HMDT d, const char *dir)
487 {
488 if(!dir || *dir==0) return;
489 long nx=d->nx,ny=d->ny,nz=d->nz,nn=nx*ny*nz;
490 long p[3]={nx,ny,nz};
491 mreal *b = new mreal[nn];
492 memcpy(b,d->a,nn*sizeof(mreal));
493 if(strchr(dir,'z') && nz>1)
494 {
495 mglStartThread(mgl_csum_z,0,nx*ny,b,d->a,0,p);
496 memcpy(d->a,b,nn*sizeof(mreal));
497 }
498 if(strchr(dir,'y') && ny>1)
499 {
500 mglStartThread(mgl_csum_y,0,nx*nz,b,d->a,0,p);
501 memcpy(d->a,b,nn*sizeof(mreal));
502 }
503 if(strchr(dir,'x') && nx>1)
504 {
505 mglStartThread(mgl_csum_x,0,nz*ny,b,d->a,0,p);
506 memcpy(d->a,b,nn*sizeof(mreal));
507 }
508 delete []b;
509 }
mgl_data_cumsum_(uintptr_t * d,const char * dir,int l)510 void MGL_EXPORT mgl_data_cumsum_(uintptr_t *d, const char *dir,int l)
511 { char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
512 mgl_data_cumsum(_DT_,s); delete []s; }
513 //-----------------------------------------------------------------------------
mgl_int_z(void * par)514 static void *mgl_int_z(void *par)
515 {
516 mglThreadD *t=(mglThreadD *)par;
517 long nz=t->p[2], nn=t->n;
518 mreal *b=t->a, dd=0.5/nz;
519 const mreal *a=t->b;
520 #if !MGL_HAVE_PTHREAD
521 #pragma omp parallel for
522 #endif
523 for(long i=t->id;i<nn;i+=mglNumThr)
524 {
525 b[i] = 0;
526 for(long j=1;j<nz;j++) b[i+j*nn] = b[i+j*nn-nn] + (a[i+nn*j]+a[i+j*nn-nn])*dd;
527 }
528 return 0;
529 }
mgl_int_y(void * par)530 static void *mgl_int_y(void *par)
531 {
532 mglThreadD *t=(mglThreadD *)par;
533 long nx=t->p[0], ny=t->p[1], nn=t->n;
534 mreal *b=t->a, dd=0.5/ny;
535 const mreal *a=t->b;
536 #if !MGL_HAVE_PTHREAD
537 #pragma omp parallel for
538 #endif
539 for(long i=t->id;i<nn;i+=mglNumThr)
540 {
541 long k = (i%nx)+nx*ny*(i/nx); b[k] = 0;
542 for(long j=1;j<ny;j++) b[k+j*nx] = b[k+j*nx-nx] + (a[k+j*nx]+a[k+j*nx-nx])*dd;
543 }
544 return 0;
545 }
mgl_int_x(void * par)546 static void *mgl_int_x(void *par)
547 {
548 mglThreadD *t=(mglThreadD *)par;
549 long nx=t->p[0], nn=t->n;
550 mreal *b=t->a, dd=0.5/nx;
551 const mreal *a=t->b;
552 #if !MGL_HAVE_PTHREAD
553 #pragma omp parallel for
554 #endif
555 for(long i=t->id;i<nn;i+=mglNumThr)
556 {
557 long k = i*nx; b[k] = 0;
558 for(long j=1;j<nx;j++) b[j+k] = b[j+k-1] + (a[j+k]+a[j+k-1])*dd;
559 }
560 return 0;
561 }
mgl_data_integral(HMDT d,const char * dir)562 void MGL_EXPORT mgl_data_integral(HMDT d, const char *dir)
563 {
564 if(!dir || *dir==0) return;
565 long nx=d->nx,ny=d->ny,nz=d->nz,nn=nx*ny*nz;
566 long p[3]={nx,ny,nz};
567 mreal *b = new mreal[nn];
568 memcpy(b,d->a,nn*sizeof(mreal));
569 if(strchr(dir,'z') && nz>1)
570 {
571 mglStartThread(mgl_int_z,0,nx*ny,b,d->a,0,p);
572 memcpy(d->a,b,nn*sizeof(mreal));
573 }
574 if(strchr(dir,'y') && ny>1)
575 {
576 mglStartThread(mgl_int_y,0,nx*nz,b,d->a,0,p);
577 memcpy(d->a,b,nn*sizeof(mreal));
578 }
579 if(strchr(dir,'x') && nx>1)
580 {
581 mglStartThread(mgl_int_x,0,nz*ny,b,d->a,0,p);
582 memcpy(d->a,b,nn*sizeof(mreal));
583 }
584 delete []b;
585 }
mgl_data_integral_(uintptr_t * d,const char * dir,int l)586 void MGL_EXPORT mgl_data_integral_(uintptr_t *d, const char *dir,int l)
587 { char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
588 mgl_data_integral(_DT_,s); delete []s; }
589 //-----------------------------------------------------------------------------
mgl_dif_z(void * par)590 static void *mgl_dif_z(void *par)
591 {
592 mglThreadD *t=(mglThreadD *)par;
593 long nz=t->p[2], nn=t->n;
594 mreal *b=t->a, dd=0.5*nz;
595 const mreal *a=t->b;
596 #if !MGL_HAVE_PTHREAD
597 #pragma omp parallel for
598 #endif
599 for(long i=t->id;i<nn;i+=mglNumThr)
600 {
601 b[i] = -(3*a[i]-4*a[i+nn]+a[i+2*nn])*dd;
602 b[i+(nz-1)*nn] = (3*a[i+(nz-1)*nn]-4*a[i+(nz-2)*nn]+a[i+(nz-3)*nn])*dd;
603 for(long j=1;j<nz-1;j++) b[i+j*nn] = (a[i+j*nn+nn]-a[i+j*nn-nn])*dd;
604 }
605 return 0;
606 }
mgl_dif_y(void * par)607 static void *mgl_dif_y(void *par)
608 {
609 mglThreadD *t=(mglThreadD *)par;
610 long nx=t->p[0], ny=t->p[1], nn=t->n;
611 mreal *b=t->a, dd=0.5*ny;
612 const mreal *a=t->b;
613 #if !MGL_HAVE_PTHREAD
614 #pragma omp parallel for
615 #endif
616 for(long i=t->id;i<nn;i+=mglNumThr)
617 {
618 long k = (i%nx)+nx*ny*(i/nx);
619 b[k] = -(3*a[k]-4*a[k+nx]+a[k+2*nx])*dd;
620 b[k+(ny-1)*nx] = (3*a[k+(ny-1)*nx]-4*a[k+(ny-2)*nx]+a[k+(ny-3)*nx])*dd;
621 for(long j=1;j<ny-1;j++) b[k+j*nx] = (a[k+j*nx+nx]-a[k+j*nx-nx])*dd;
622 }
623 return 0;
624 }
mgl_dif_x(void * par)625 static void *mgl_dif_x(void *par)
626 {
627 mglThreadD *t=(mglThreadD *)par;
628 long nx=t->p[0], nn=t->n;
629 mreal *b=t->a, dd=0.5*nx;
630 const mreal *a=t->b;
631 #if !MGL_HAVE_PTHREAD
632 #pragma omp parallel for
633 #endif
634 for(long i=t->id;i<nn;i+=mglNumThr)
635 {
636 long k = i*nx;
637 b[k] = -(3*a[k]-4*a[k+1]+a[k+2])*dd;
638 b[k+nx-1] = (3*a[k+nx-1]-4*a[k+nx-2]+a[k+nx-3])*dd;
639 for(long j=1;j<nx-1;j++) b[j+k] = (a[j+k+1]-a[j+k-1])*dd;
640 }
641 return 0;
642 }
mgl_data_diff(HMDT d,const char * dir)643 void MGL_EXPORT mgl_data_diff(HMDT d, const char *dir)
644 {
645 if(!dir || *dir==0) return;
646 long nx=d->nx,ny=d->ny,nz=d->nz,nn=nx*ny*nz;
647 long p[3]={nx,ny,nz};
648 mreal *b = new mreal[nn];
649 if(strchr(dir,'z') && nz>1)
650 {
651 mglStartThread(mgl_dif_z,0,nx*ny,b,d->a,0,p);
652 memcpy(d->a,b,nn*sizeof(mreal));
653 }
654 if(strchr(dir,'y') && ny>1)
655 {
656 mglStartThread(mgl_dif_y,0,nx*nz,b,d->a,0,p);
657 memcpy(d->a,b,nn*sizeof(mreal));
658 }
659 if(strchr(dir,'x') && nx>1)
660 {
661 mglStartThread(mgl_dif_x,0,nz*ny,b,d->a,0,p);
662 memcpy(d->a,b,nn*sizeof(mreal));
663 }
664 delete []b;
665 }
mgl_data_diff_(uintptr_t * d,const char * dir,int l)666 void MGL_EXPORT mgl_data_diff_(uintptr_t *d, const char *dir,int l)
667 { char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
668 mgl_data_diff(_DT_,s); delete []s; }
669 //-----------------------------------------------------------------------------
mgl_dif2_z(void * par)670 static void *mgl_dif2_z(void *par)
671 {
672 mglThreadD *t=(mglThreadD *)par;
673 long nz=t->p[2], nn=t->n;
674 mreal *b=t->a, dd=nz*nz;
675 const mreal *a=t->b;
676 #if !MGL_HAVE_PTHREAD
677 #pragma omp parallel for
678 #endif
679 for(long i=t->id;i<nn;i+=mglNumThr)
680 {
681 b[i] = b[i+(nz-1)*nn] = 0;
682 for(long j=1;j<nz-1;j++) b[i+j*nn] = (a[i+j*nn+nn]+a[i+j*nn-nn]-2*a[i+j*nn])*dd;
683 }
684 return 0;
685 }
mgl_dif2_y(void * par)686 static void *mgl_dif2_y(void *par)
687 {
688 mglThreadD *t=(mglThreadD *)par;
689 long nx=t->p[0], ny=t->p[1], nn=t->n;
690 mreal *b=t->a, dd=ny*ny;
691 const mreal *a=t->b;
692 #if !MGL_HAVE_PTHREAD
693 #pragma omp parallel for
694 #endif
695 for(long i=t->id;i<nn;i+=mglNumThr)
696 {
697 long k = (i%nx)+nx*ny*(i/nx); b[k] = b[k+(ny-1)*nx] = 0;
698 for(long j=1;j<ny-1;j++) b[k+j*nx] = (a[k+j*nx+nx]+a[k+j*nx-nx]-2*a[k+j*nx])*dd;
699 }
700 return 0;
701 }
mgl_dif2_x(void * par)702 static void *mgl_dif2_x(void *par)
703 {
704 mglThreadD *t=(mglThreadD *)par;
705 long nx=t->p[0], nn=t->n;
706 mreal *b=t->a, dd=nx*nx;
707 const mreal *a=t->b;
708 #if !MGL_HAVE_PTHREAD
709 #pragma omp parallel for
710 #endif
711 for(long i=t->id;i<nn;i+=mglNumThr)
712 {
713 long k = i*nx; b[k] = b[k+nx-1] = 0;
714 for(long j=1;j<nx-1;j++) b[j+k] = (a[j+k+1]+a[j+k-1]-2*a[j+k])*dd;
715 }
716 return 0;
717 }
mgl_data_diff2(HMDT d,const char * dir)718 void MGL_EXPORT mgl_data_diff2(HMDT d, const char *dir)
719 {
720 if(!dir || *dir==0) return;
721 long nx=d->nx,ny=d->ny,nz=d->nz,nn=nx*ny*nz;
722 long p[3]={nx,ny,nz};
723 mreal *b = new mreal[nn];
724 if(strchr(dir,'z') && nz>1)
725 {
726 mglStartThread(mgl_dif2_z,0,nx*ny,b,d->a,0,p);
727 memcpy(d->a,b,nn*sizeof(mreal));
728 }
729 if(strchr(dir,'y') && ny>1)
730 {
731 mglStartThread(mgl_dif2_y,0,nx*nz,b,d->a,0,p);
732 memcpy(d->a,b,nn*sizeof(mreal));
733 }
734 if(strchr(dir,'x') && nx>1)
735 {
736 mglStartThread(mgl_dif2_x,0,nz*ny,b,d->a,0,p);
737 memcpy(d->a,b,nn*sizeof(mreal));
738 }
739 delete []b;
740 }
mgl_data_diff2_(uintptr_t * d,const char * dir,int l)741 void MGL_EXPORT mgl_data_diff2_(uintptr_t *d, const char *dir,int l)
742 { char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
743 mgl_data_diff2(_DT_,s); delete []s; }
744 //-----------------------------------------------------------------------------
mgl_data_swap(HMDT d,const char * dir)745 void MGL_EXPORT mgl_data_swap(HMDT d, const char *dir)
746 {
747 if(!dir || *dir==0) return;
748 if(strchr(dir,'z') && d->nz>1) mgl_data_roll(d,'z',d->nz/2);
749 if(strchr(dir,'y') && d->ny>1) mgl_data_roll(d,'y',d->ny/2);
750 if(strchr(dir,'x') && d->nx>1) mgl_data_roll(d,'x',d->nx/2);
751 }
mgl_data_swap_(uintptr_t * d,const char * dir,int l)752 void MGL_EXPORT mgl_data_swap_(uintptr_t *d, const char *dir,int l)
753 { char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
754 mgl_data_swap(_DT_,s); delete []s; }
755 //-----------------------------------------------------------------------------
mgl_data_roll(HMDT dd,char dir,long num)756 void MGL_EXPORT mgl_data_roll(HMDT dd, char dir, long num)
757 {
758 long nx=dd->nx,ny=dd->ny,nz=dd->nz, d;
759 mreal *b,*a=dd->a;
760 if(dir=='z' && nz>1)
761 {
762 d = num>0 ? num%nz : (num+nz*(1-num/nz))%nz;
763 if(d==0) return; // nothing to do
764 b = new mreal[nx*ny*nz];
765 memcpy(b,a+nx*ny*d,nx*ny*(nz-d)*sizeof(mreal));
766 memcpy(b+nx*ny*(nz-d),a,nx*ny*d*sizeof(mreal));
767 memcpy(a,b,nx*ny*nz*sizeof(mreal)); delete []b;
768 }
769 if(dir=='y' && ny>1)
770 {
771 d = num>0 ? num%ny : (num+ny*(1-num/ny))%ny;
772 if(d==0) return; // nothing to do
773 b = new mreal[nx*ny*nz];
774 memcpy(b,a+nx*d,(nx*ny*nz-nx*d)*sizeof(mreal));
775 #pragma omp parallel for
776 for(long i=0;i<nz;i++)
777 memcpy(b+nx*(ny-d)+nx*ny*i,a+nx*ny*i,nx*d*sizeof(mreal));
778 memcpy(a,b,nx*ny*nz*sizeof(mreal)); delete []b;
779 }
780 if(dir=='x' && nx>1)
781 {
782 d = num>0 ? num%nx : (num+nx*(1-num/nx))%nx;
783 if(d==0) return; // nothing to do
784 b = new mreal[nx*ny*nz];
785 memcpy(b,a+d,(nx*ny*nz-d)*sizeof(mreal));
786 #pragma omp parallel for
787 for(long i=0;i<nz*ny;i++)
788 memcpy(b+nx-d+nx*i,a+nx*i,d*sizeof(mreal));
789 memcpy(a,b,nx*ny*nz*sizeof(mreal)); delete []b;
790 }
791 }
mgl_data_roll_(uintptr_t * d,const char * dir,int * num,int)792 void MGL_EXPORT mgl_data_roll_(uintptr_t *d, const char *dir, int *num, int)
793 { mgl_data_roll(_DT_,*dir,*num); }
794 //-----------------------------------------------------------------------------
mgl_data_mirror(HMDT d,const char * dir)795 void MGL_EXPORT mgl_data_mirror(HMDT d, const char *dir)
796 {
797 if(!dir || *dir==0) return;
798 long nx=d->nx,ny=d->ny,nz=d->nz;
799 mreal *a=d->a;
800 if(strchr(dir,'z') && nz>1)
801 {
802 #pragma omp parallel for collapse(2)
803 for(long j=0;j<nz/2;j++) for(long i=0;i<nx*ny;i++)
804 {
805 long i0 = i+j*nx*ny, j0 = i+(nz-1-j)*nx*ny;
806 mreal b = a[i0]; a[i0] = a[j0]; a[j0] = b;
807 }
808 }
809 if(strchr(dir,'y') && ny>1)
810 {
811 #pragma omp parallel for
812 for(long i=0;i<nx*nz;i++)
813 {
814 long j0 = (i%nx)+nx*ny*(i/nx);
815 for(long j=0;j<ny/2;j++)
816 {
817 long i0 = j0+(ny-1-j)*nx;
818 mreal b = a[j0+j*nx]; a[j0+j*nx] = a[i0]; a[i0] = b;
819 }
820 }
821 }
822 if(strchr(dir,'x') && nx>1)
823 {
824 #pragma omp parallel for
825 for(long j=0;j<ny*nz;j++)
826 {
827 long j0 = j*nx;
828 for(long i=0;i<nx/2;i++)
829 {
830 long i0 = nx-1-i+j0;
831 mreal b = a[i+j0]; a[i+j0] = a[i0]; a[i0] = b;
832 }
833 }
834 }
835 }
mgl_data_mirror_(uintptr_t * d,const char * dir,int l)836 void MGL_EXPORT mgl_data_mirror_(uintptr_t *d, const char *dir,int l)
837 { char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
838 mgl_data_mirror(_DT_,s); delete []s; }
839 //-----------------------------------------------------------------------------
mgl_data_clean(HMDT d,long id)840 void MGL_EXPORT mgl_data_clean(HMDT d, long id)
841 {
842 if(id<0 || id+1>d->nx) return;
843 long i,j,n=d->nx,m=d->ny;
844 mreal *b = new mreal[m*n], *a=d->a;
845 for(i=j=0;i+1<m;i++)
846 {
847 if(a[id+n*i]!=a[id+n*i+n]) // this can be saved
848 {
849 #pragma omp parallel for
850 for(long k=0;k<n;k++) b[k+n*j]=a[k+n*i];
851 j++;
852 }
853 }
854 // always save last row
855 i=n*(m-1);
856 #pragma omp parallel for
857 for(long k=0;k<n;k++) b[k+n*j]=a[k+i];
858 j++;
859 memcpy(a,b,n*j*sizeof(mreal)); d->ny = j;
860 delete []b;
861 }
mgl_data_clean_(uintptr_t * d,int * id)862 void MGL_EXPORT mgl_data_clean_(uintptr_t *d, int *id) { mgl_data_clean(_DT_,*id); }
863 //-----------------------------------------------------------------------------
mgl_solve_x(void * par)864 static void *mgl_solve_x(void *par)
865 {
866 mglThreadD *t=(mglThreadD *)par;
867 HCDT d=(HCDT)t->v;
868 long nx=t->p[0], ny=t->p[1], nz=t->p[2], n1=t->p[3]?nx-1:1, nn=t->n;
869 const mreal *a=t->b, *ii=t->c;
870 mreal *b=t->a,val=t->d[0],da = 1e-5*(val?fabs(val):1);
871 #if !MGL_HAVE_PTHREAD
872 #pragma omp parallel for
873 #endif
874 for(long l=t->id;l<nn;l+=mglNumThr)
875 {
876 long j=l%ny, k=l/ny; b[l] = NAN;
877 if(ii && ii[l]!=ii[l]) continue;
878 long i0 = ii?(ii[l]*n1+1):0;
879 if(i0>nx-2) continue;
880 if(a) for(long i=i0+1;i<nx;i++)
881 {
882 mreal y1=a[i-1+nx*l], y2=a[i+nx*l];
883 if((y1-val)*(y2-val)<=0)
884 {
885 mreal x = i-1 + (val-y1)/(y2-y1), dx;
886 mreal v0=mglSpline3(a,nx,ny,nz,x,j,k, &dx,0,0), v=v0;
887 unsigned kk=0;
888 while(fabs(v-val)>da || dx==0)
889 {
890 x += (val-v)/dx; kk++;
891 v = mglSpline3(a,nx,ny,nz,x,j,k, &dx,0,0);
892 if(kk>=10)
893 {
894 b[l] = x = fabs(v-val)<fabs(v0-val) ? x:i-1 + (val-y1)/(y2-y1);
895 break;
896 }
897 }
898 b[l] = x; break;
899 }
900 }
901 else for(long i=i0+1;i<nx;i++)
902 {
903 mreal y1=d->v(i-1,j,k), y2=d->v(i,j,k);
904 if((y1-val)*(y2-val)<=0)
905 { b[l] = i-1 + (val-y1)/(y2-y1); break; }
906 }
907 b[l] /= n1;
908 }
909 return 0;
910 }
mgl_solve_y(void * par)911 static void *mgl_solve_y(void *par)
912 {
913 mglThreadD *t=(mglThreadD *)par;
914 HCDT d=(HCDT)t->v;
915 long nx=t->p[0], ny=t->p[1], nz=t->p[2], n1=t->p[3]?ny-1:1, nn=t->n;
916 const mreal *a=t->b, *ii=t->c;
917 mreal *b=t->a,val=t->d[0],da = 1e-5*(val?fabs(val):1);
918 #if !MGL_HAVE_PTHREAD
919 #pragma omp parallel for
920 #endif
921 for(long l=t->id;l<nn;l+=mglNumThr)
922 {
923 long j=l%nx, k=l/nx; b[l] = NAN;
924 if(ii && ii[l]!=ii[l]) continue;
925 long i0 = ii?(ii[l]*n1+1):0;
926 if(i0>ny-2) continue;
927 if(a) for(long i=i0+1;i<ny;i++)
928 {
929 mreal y1=a[j+nx*(i-1+ny*k)], y2=a[j+nx*(i+ny*k)];
930 if((y1-val)*(y2-val)<=0)
931 {
932 mreal x = i-1 + (val-y1)/(y2-y1), dy;
933 mreal v0=mglSpline3(a,nx,ny,nz,j,x,k, 0,&dy,0), v=v0;
934 unsigned kk=0;
935 while(fabs(v-val)>da || dy==0)
936 {
937 x += (val-v)/dy; kk++;
938 v = mglSpline3(a,nx,ny,nz,j,x,k, 0,&dy,0);
939 if(kk>=10)
940 {
941 b[l] = x = fabs(v-val)<fabs(v0-val) ? x:i-1 + (val-y1)/(y2-y1);
942 break;
943 }
944 }
945 b[l] = x; break;
946 }
947 }
948 else for(long i=i0+1;i<ny;i++)
949 {
950 mreal y1=d->v(j,i-1,k), y2=d->v(j,i,k);
951 if((y1-val)*(y2-val)<=0)
952 { b[l] = i-1 + (val-y1)/(y2-y1); break; }
953 }
954 b[l] /= n1;
955 }
956 return 0;
957 }
mgl_solve_z(void * par)958 static void *mgl_solve_z(void *par)
959 {
960 mglThreadD *t=(mglThreadD *)par;
961 HCDT d=(HCDT)t->v;
962 long nx=t->p[0], ny=t->p[1], nz=t->p[2], n1=t->p[3]?nz-1:1, nn=t->n;
963 const mreal *a=t->b, *ii=t->c;
964 mreal *b=t->a,val=t->d[0],da = 1e-5*(val?fabs(val):1);
965 #if !MGL_HAVE_PTHREAD
966 #pragma omp parallel for
967 #endif
968 for(long l=t->id;l<nn;l+=mglNumThr)
969 {
970 long j=l%nx, k=l/nx; b[l] = NAN;
971 if(ii && ii[l]!=ii[l]) continue;
972 long i0 = ii?(ii[l]*n1+1):0;
973 if(i0>nz-2) continue;
974 if(a) for(long i=i0+1;i<nz;i++)
975 {
976 mreal y1=a[nn*i-nn+l], y2=a[nn*i+l];
977 if((y1-val)*(y2-val)<=0)
978 {
979 mreal x = i-1 + (val-y1)/(y2-y1), dz;
980 mreal v0=mglSpline3(a,nx,ny,nz,j,k,x, 0,0,&dz), v=v0;
981 unsigned kk=0;
982 while(fabs(v-val)>da || dz==0)
983 {
984 x += (val-v)/dz; kk++;
985 v = mglSpline3(a,nx,ny,nz,j,k,x, 0,0,&dz);
986 if(kk>=10)
987 {
988 b[l] = x = fabs(v-val)<fabs(v0-val) ? x:i-1 + (val-y1)/(y2-y1);
989 break;
990 }
991 }
992 b[l] = x; break;
993 }
994 }
995 else for(long i=i0+1;i<nz;i++)
996 {
997 mreal y1=d->v(j,k,i-1), y2=d->v(j,k,i);
998 if((y1-val)*(y2-val)<=0)
999 { b[l] = i-1 + (val-y1)/(y2-y1); break; }
1000 }
1001 b[l] /= n1;
1002 }
1003 return 0;
1004 }
mgl_data_solve(HCDT dat,mreal val,char dir,HCDT i0,int norm)1005 HMDT MGL_EXPORT mgl_data_solve(HCDT dat, mreal val, char dir, HCDT i0, int norm)
1006 {
1007 const mglData *i = dynamic_cast<const mglData *>(i0);
1008 const mglData *d = dynamic_cast<const mglData *>(dat);
1009 long p[4]={dat->GetNx(), dat->GetNy(), dat->GetNz(), norm};
1010 const mreal *ii=0;
1011 mglData *r=new mglData, id0;
1012 if(i0 && !i) { id0.Set(i0); i=&id0; } // <-- slow but should work
1013 if(dir=='x' && p[0]>1)
1014 {
1015 r->Create(p[1],p[2]);
1016 ii = (i && i->nx*i->ny==p[1]*p[2])?i->a:0;
1017 mglStartThread(mgl_solve_x,0,p[1]*p[2],r->a,d?d->a:0,ii,p,dat,&val);
1018 }
1019 if(dir=='y' && p[1]>1)
1020 {
1021 r->Create(p[0],p[2]);
1022 ii = (i && i->nx*i->ny==p[0]*p[2])?i->a:0;
1023 mglStartThread(mgl_solve_y,0,p[0]*p[2],r->a,d?d->a:0,ii,p,dat,&val);
1024 }
1025 if(dir=='z' && p[2]>1)
1026 {
1027 r->Create(p[0],p[1]);
1028 ii = (i && i->nx*i->ny==p[0]*p[1])?i->a:0;
1029 mglStartThread(mgl_solve_z,0,p[0]*p[1],r->a,d?d->a:0,ii,p,dat,&val);
1030 }
1031 return r;
1032 }
1033 //-----------------------------------------------------------------------------
mgl_data_solve_1d(HCDT d,mreal val,int spl,long i0)1034 mreal MGL_EXPORT mgl_data_solve_1d(HCDT d, mreal val, int spl, long i0)
1035 {
1036 mreal x=0, y1, y2, a, a0, dx=0, da = 1e-5*(val?fabs(val):1);
1037 long nx = d->GetNx();
1038 if(i0<0 || i0>=nx) i0=0;
1039 if(val==d->v(i0+1)) return i0+1;
1040 const mglData *dd=dynamic_cast<const mglData *>(d);
1041 const mglDataC *dc=dynamic_cast<const mglDataC *>(d);
1042 if(dd) for(long i=i0+1;i<nx;i++)
1043 {
1044 y1=dd->a[i-1]; y2=dd->a[i];
1045 if((y1-val)*(y2-val)<=0)
1046 {
1047 x = i-1 + (val-y1)/(y2-y1);
1048 a0 = a = mglSpline1t<mreal>(dd->a,nx,x,&dx);
1049 if(spl) for(unsigned k=0;fabs(a-val)>da || dx==0;)
1050 {
1051 x += (val-a)/dx; k++;
1052 a = mglSpline1t<mreal>(dd->a,nx,x,&dx);
1053 if(k>=10)
1054 return fabs(a-val)<fabs(a0-val) ? x:i-1 + (val-y1)/(y2-y1);
1055 }
1056 return x;
1057 }
1058 }
1059 else if(dc) for(long i=i0+1;i<nx;i++)
1060 {
1061 y1=abs(dc->a[i-1]); y2=abs(dc->a[i]);
1062 if((y1-val)*(y2-val)<=0)
1063 {
1064 x = i-1 + (val-y1)/(y2-y1);
1065 dual cx, ca = mglSpline1t<dual>(dc->a,nx,x,&cx);
1066 a0 = a = abs(ca); dx = a?(cx.real()*ca.real()+cx.imag()*ca.imag())/a:0;
1067 if(spl) for(unsigned k=0;fabs(a-val)>da || dx==0;)
1068 {
1069 x += (val-a)/dx; k++;
1070 ca = mglSpline1t<dual>(dc->a,nx,x,&cx);
1071 a = abs(ca); dx = a?(cx.real()*ca.real()+cx.imag()*ca.imag())/a:0;
1072 if(k>=10)
1073 return fabs(a-val)<fabs(a0-val) ? x:i-1 + (val-y1)/(y2-y1);
1074 }
1075 return x;
1076 }
1077 }
1078 else for(long i=i0+1;i<nx;i++)
1079 {
1080 y1=d->v(i-1); y2=d->v(i);
1081 if((y1-val)*(y2-val)<=0)
1082 return i-1 + (val-y1)/(y2-y1);
1083 }
1084 return NAN;
1085 }
1086 //-----------------------------------------------------------------------------
mgl_data_linear_ext(HCDT d,mreal x,mreal y,mreal z,mreal * dx,mreal * dy,mreal * dz)1087 mreal MGL_EXPORT mgl_data_linear_ext(HCDT d, mreal x,mreal y,mreal z, mreal *dx,mreal *dy,mreal *dz)
1088 {
1089 if(!d) return NAN;
1090 long kx=long(x), ky=long(y), kz=long(z);
1091 mreal b0,b1;
1092 const mglData *dd=dynamic_cast<const mglData *>(d);
1093 if(dd)
1094 {
1095 long nx=dd->nx, ny=dd->ny, nz=dd->nz, dn=ny>1?nx:0;
1096 kx = kx>=0 ? kx:0; kx = kx<nx-1 ? kx:nx-2;
1097 ky = ky>=0 ? ky:0; ky = ky<ny-1 ? ky:ny-2;
1098 kz = kz>=0 ? kz:0; kz = kz<nz-1 ? kz:nz-2;
1099 x -= kx; y -= ky; z -= kz;
1100 const mreal *aa = dd->a, *bb;
1101 if(kz>=0)
1102 {
1103 aa=dd->a+kx+nx*(ky+ny*kz); bb = aa+nx*ny;
1104 b0 = aa[0]*(1-x-y+x*y) + x*(1-y)*aa[1] + y*(1-x)*aa[dn] + x*y*aa[1+dn];
1105 b1 = bb[0]*(1-x-y+x*y) + x*(1-y)*bb[1] + y*(1-x)*bb[dn] + x*y*bb[1+dn];
1106 }
1107 else
1108 {
1109 z=0;
1110 if(ky>=0)
1111 {
1112 aa=dd->a+kx+nx*ky;
1113 b0 = b1 = aa[0]*(1-x-y+x*y) + x*(1-y)*aa[1] + y*(1-x)*aa[dn] + x*y*aa[1+dn];
1114 }
1115 else if(kx>=0)
1116 {
1117 aa=dd->a+kx; b0 = b1 = aa[0]*(1-x) + x*aa[1];
1118 }
1119 else b0 = b1 = dd->a[0];
1120 }
1121 if(dx) *dx = kx>=0?aa[1]-aa[0]:0;
1122 if(dy) *dy = ky>=0?aa[dn]-aa[0]:0;
1123 if(dz) *dz = b1-b0;
1124 }
1125 else
1126 {
1127 long nx=d->GetNx(), ny=d->GetNy(), nz=d->GetNz();
1128 kx = kx>=0 ? kx:0; kx = kx<nx-1 ? kx:nx-2;
1129 ky = ky>=0 ? ky:0; ky = ky<ny-1 ? ky:ny-2;
1130 kz = kz>=0 ? kz:0; kz = kz<nz-1 ? kz:nz-2;
1131 x -= kx; y -= ky; z -= kz;
1132 mreal a0 = 0, a1 = 0, a2 = 0;
1133 if(kz>=0)
1134 {
1135 a0 = d->v(kx,ky,kz); a1 = d->v(kx+1,ky,kz); a2 = d->v(kx,ky+1,kz);
1136 b0 = a0*(1-x-y+x*y) + x*(1-y)*a1 + y*(1-x)*a2 + x*y*d->v(kx+1,ky+1,kz);;
1137 b1 = d->v(kx,ky,kz+1)*(1-x-y+x*y) + x*(1-y)*d->v(kx+1,ky,kz+1) + y*(1-x)*d->v(kx,ky+1,kz+1) + x*y*d->v(kx+1,ky+1,kz+1);
1138 }
1139 else
1140 {
1141 z=0;
1142 if(ky>=0)
1143 {
1144 a0 = d->v(kx,ky); a1 = d->v(kx+1,ky); a2 = d->v(kx,ky+1);
1145 b0 = b1 = a0*(1-x-y+x*y) + x*(1-y)*a1 + y*(1-x)*a2 + x*y*d->v(kx+1,ky+1);
1146 }
1147 else if(kx>=0)
1148 {
1149 a2=a0 = d->v(kx); a1 = d->v(kx+1); b0 = b1 = a0*(1-x) + x*a1;
1150 }
1151 else b0 = b1 = d->v(0);
1152 }
1153 if(dx) *dx = a1-a0;
1154 if(dy) *dy = a2-a0;
1155 if(dz) *dz = b1-b0;
1156 }
1157 return b0 + z*(b1-b0);
1158 }
1159 //-----------------------------------------------------------------------------
mgl_data_linear(HCDT d,mreal x,mreal y,mreal z)1160 mreal MGL_EXPORT mgl_data_linear(HCDT d, mreal x,mreal y,mreal z)
1161 { return mgl_data_linear_ext(d, x,y,z, 0,0,0); }
1162 //-----------------------------------------------------------------------------
mgl_data_spline(HCDT d,mreal x,mreal y,mreal z)1163 mreal MGL_EXPORT mgl_data_spline(HCDT d, mreal x,mreal y,mreal z)
1164 {
1165 if(mgl_isbad(x) || mgl_isbad(y) || mgl_isbad(z)) return NAN;
1166 return d->value(x,y,z);
1167 }
1168 //-----------------------------------------------------------------------------
mgl_data_spline_ext(HCDT d,mreal x,mreal y,mreal z,mreal * dx,mreal * dy,mreal * dz)1169 mreal MGL_EXPORT mgl_data_spline_ext(HCDT d, mreal x,mreal y,mreal z, mreal *dx,mreal *dy,mreal *dz)
1170 {
1171 if(mgl_isbad(x) || mgl_isbad(y) || mgl_isbad(z)) return NAN;
1172 return d->valueD(x,y,z,dx,dy,dz);
1173 }
1174 //-----------------------------------------------------------------------------
mgl_data_spline_(uintptr_t * d,mreal * x,mreal * y,mreal * z)1175 mreal MGL_EXPORT mgl_data_spline_(uintptr_t *d, mreal *x,mreal *y,mreal *z)
1176 { return mgl_data_spline(_DA_(d),*x,*y,*z); }
mgl_data_linear_(uintptr_t * d,mreal * x,mreal * y,mreal * z)1177 mreal MGL_EXPORT mgl_data_linear_(uintptr_t *d, mreal *x,mreal *y,mreal *z)
1178 { return mgl_data_linear(_DA_(d),*x,*y,*z); }
mgl_data_spline_ext_(uintptr_t * d,mreal * x,mreal * y,mreal * z,mreal * dx,mreal * dy,mreal * dz)1179 mreal MGL_EXPORT mgl_data_spline_ext_(uintptr_t *d, mreal *x,mreal *y,mreal *z, mreal *dx,mreal *dy,mreal *dz)
1180 { return mgl_data_spline_ext(_DA_(d),*x,*y,*z,dx,dy,dz); }
mgl_data_linear_ext_(uintptr_t * d,mreal * x,mreal * y,mreal * z,mreal * dx,mreal * dy,mreal * dz)1181 mreal MGL_EXPORT mgl_data_linear_ext_(uintptr_t *d, mreal *x,mreal *y,mreal *z, mreal *dx,mreal *dy,mreal *dz)
1182 { return mgl_data_linear_ext(_DA_(d),*x,*y,*z,dx,dy,dz); }
mgl_data_solve_1d_(uintptr_t * d,mreal * val,int * spl,int * i0)1183 mreal MGL_EXPORT mgl_data_solve_1d_(uintptr_t *d, mreal *val, int *spl, int *i0)
1184 { return mgl_data_solve_1d(_DA_(d),*val, *spl, *i0); }
mgl_data_solve_(uintptr_t * d,mreal * val,const char * dir,uintptr_t * i0,int * norm,int)1185 uintptr_t MGL_EXPORT mgl_data_solve_(uintptr_t *d, mreal *val, const char *dir, uintptr_t *i0, int *norm,int)
1186 { return uintptr_t(mgl_data_solve(_DA_(d),*val, *dir, _DA_(i0), *norm)); }
1187 //-----------------------------------------------------------------------------
int_pow(long x,long n)1188 long MGL_LOCAL_CONST int_pow(long x, long n)
1189 {
1190 if(n==2) return x*x;
1191 if(n==1) return x;
1192 if(n==0) return 1;
1193 if(n<0) return 0;
1194 long t = int_pow(x,n/2); t = t*t;
1195 if(n%2==1) t *= x;
1196 return t;
1197 }
mgl_powers(long N,const char * how)1198 long MGL_NO_EXPORT mgl_powers(long N, const char *how)
1199 {
1200 bool k2 = mglchr(how,'2'), k3 = mglchr(how,'3'), k5 = mglchr(how,'5');
1201 const double lN=log(N), l2=log(2), l3=log(3), l5=log(5);
1202 if(k2 && k3 && k5)
1203 {
1204 double dm=lN; long im=0, jm=0, km=0;
1205 for(long i=0;i<=lN/l2;i++) for(long j=0;j<=(lN-i*l2)/l3;j++) for(long k=0;k<=(lN-i*l2-j*l3)/l5;k++)
1206 {
1207 double d = lN-i*l2-j*l3-k*l5;
1208 if(d>0 && d<dm) { im=i; jm=j; km=k; dm=d; }
1209 }
1210 return int_pow(2,im)*int_pow(3,jm)*int_pow(5,km);
1211 }
1212 else if(k2 && !k3 && !k5) return int_pow(2,lN/l2);
1213 else if(k3 && !k2 && !k5) return int_pow(3,lN/l3);
1214 else if(k5 && !k3 && !k2) return int_pow(5,lN/l5);
1215 else if(k2 && k3 && !k5)
1216 {
1217 double dm=lN; long im=0, jm=0;
1218 for(long i=0;i<=lN/l2;i++) for(long j=0;j<=(lN-i*l2)/l3;j++)
1219 {
1220 double d = lN-i*l2-j*l3;
1221 if(d>0 && d<dm) { im=i; jm=j; dm=d; }
1222 }
1223 return int_pow(2,im)*int_pow(3,jm);
1224 }
1225 else if(k2 && k5 && !k3)
1226 {
1227 double dm=lN; long im=0, jm=0;
1228 for(long i=0;i<=lN/l2;i++) for(long j=0;j<=(lN-i*l2)/l5;j++)
1229 {
1230 double d = lN-i*l2-j*l5;
1231 if(d>0 && d<dm) { im=i; jm=j; dm=d; }
1232 }
1233 return int_pow(2,im)*int_pow(5,jm);
1234 }
1235 else if(k5 && k3 && !k2)
1236 {
1237 double dm=lN; long im=0, jm=0;
1238 for(long i=0;i<=lN/l5;i++) for(long j=0;j<=(lN-i*l5)/l3;j++)
1239 {
1240 double d = lN-i*l5-j*l3;
1241 if(d>0 && d<dm) { im=i; jm=j; dm=d; }
1242 }
1243 return int_pow(5,im)*int_pow(3,jm);
1244 }
1245 return 0;
1246 }
1247 //-----------------------------------------------------------------------------
mgl_data_crop_opt(HMDT d,const char * how)1248 void MGL_EXPORT mgl_data_crop_opt(HMDT d, const char *how)
1249 {
1250 const char *h = "235";
1251 if(mglchr(how,'2') || mglchr(how,'3') || mglchr(how,'5')) h = how;
1252 if(mglchr(how,'x')) mgl_data_crop(d, 0, mgl_powers(d->nx, h), 'x');
1253 if(mglchr(how,'y')) mgl_data_crop(d, 0, mgl_powers(d->ny, h), 'y');
1254 if(mglchr(how,'z')) mgl_data_crop(d, 0, mgl_powers(d->nz, h), 'z');
1255 }
mgl_data_crop_opt_(uintptr_t * d,const char * how,int l)1256 void MGL_EXPORT mgl_data_crop_opt_(uintptr_t *d, const char *how, int l)
1257 { char *s=new char[l+1]; memcpy(s,how,l); s[l]=0;
1258 mgl_data_crop_opt(_DT_,s); delete []s; }
1259 //-----------------------------------------------------------------------------
mgl_data_crop(HMDT d,long n1,long n2,char dir)1260 void MGL_EXPORT mgl_data_crop(HMDT d, long n1, long n2, char dir)
1261 {
1262 long nx=d->nx,ny=d->ny,nz=d->nz, nn;
1263 mreal *b;
1264 if(n1<0) n1=0;
1265 switch(dir)
1266 {
1267 case 'x':
1268 if(n1>=nx) break;
1269 n2 = n2>0 ? n2 : nx+n2;
1270 if(n2<0 || n2>=nx || n2<n1) n2 = nx;
1271 nn = n2-n1; b = new mreal[nn*ny*nz];
1272 #pragma omp parallel for
1273 for(long i=0;i<ny*nz;i++)
1274 memcpy(b+nn*i,d->a+nx*i+n1,nn*sizeof(mreal));
1275 d->nx = nn; if(!d->link) delete []d->a;
1276 d->a = b; d->link=false; d->NewId();
1277 break;
1278 case 'y':
1279 if(n1>=ny) break;
1280 n2 = n2>0 ? n2 : ny+n2;
1281 if(n2<0 || n2>=ny || n2<n1) n2 = ny;
1282 nn = n2-n1; b = new mreal[nn*nx*nz];
1283 #pragma omp parallel for collapse(2)
1284 for(long j=0;j<nz;j++) for(long i=0;i<nn;i++)
1285 memcpy(b+nx*(i+nn*j),d->a+nx*(n1+i+ny*j),nx*sizeof(mreal));
1286 d->ny = nn; if(!d->link) delete []d->a;
1287 d->a = b; d->link=false;
1288 break;
1289 case 'z':
1290 if(n1>=nz) break;
1291 n2 = n2>0 ? n2 : nz+n2;
1292 if(n2<0 || n2>=nz || n2<n1) n2 = nz;
1293 nn = n2-n1; b = new mreal[nn*nx*ny];
1294 memcpy(b,d->a+nx*ny*n1,nn*nx*ny*sizeof(mreal));
1295 d->nz = nn; if(!d->link) delete []d->a;
1296 d->a = b; d->link=false;
1297 break;
1298 }
1299 }
mgl_data_crop_(uintptr_t * d,int * n1,int * n2,const char * dir,int)1300 void MGL_EXPORT mgl_data_crop_(uintptr_t *d, int *n1, int *n2, const char *dir,int)
1301 { mgl_data_crop(_DT_,*n1,*n2,*dir); }
1302 //-----------------------------------------------------------------------------
mgl_data_last(HCDT d,const char * cond,long * i,long * j,long * k)1303 mreal MGL_EXPORT mgl_data_last(HCDT d, const char *cond, long *i, long *j, long *k)
1304 {
1305 long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
1306 if(!cond) cond = "u";
1307 mglFormula eq(cond);
1308 if(*i<0 || *i>=nx) *i=nx;
1309 if(*j<0 || *j>=ny) *j=ny-1;
1310 if(*k<0 || *k>=nz) *k=nz-1;
1311 long i0 = *i+nx*(*j+ny*(*k))-1;
1312 mreal x,y,z,dx=nx>1?1/(nx-1.):0,dy=ny>1?1/(ny-1.):0,dz=nz>1?1/(nz-1.):0;
1313 for(;i0>=0;i0--)
1314 {
1315 x = dx*(i0%nx); y = dy*((i0/nx)%ny); z = dz*(i0/(nx*ny));
1316 if(eq.Calc(x,y,z,d->vthr(i0))) break;
1317 }
1318 *i = i0%nx; *j = (i0/nx)%ny; *k = i0/(nx*ny);
1319 return i0>=0 ? d->vthr(i0) : NAN; // NOTE: Return NAN if false
1320 }
mgl_data_last_(uintptr_t * d,const char * cond,int * i,int * j,int * k,int l)1321 mreal MGL_EXPORT mgl_data_last_(uintptr_t *d, const char *cond, int *i, int *j, int *k, int l)
1322 { long ii=*i,jj=*j,kk=*k; char *s=new char[l+1]; memcpy(s,cond,l); s[l]=0;
1323 mreal res = mgl_data_last(_DT_,s,&ii,&jj,&kk); *i=ii; *j=jj; *k=kk;
1324 delete []s; return res; }
1325 //-----------------------------------------------------------------------------
mgl_data_first(HCDT d,const char * cond,long * i,long * j,long * k)1326 mreal MGL_EXPORT mgl_data_first(HCDT d, const char *cond, long *i, long *j, long *k)
1327 {
1328 long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
1329 if(!cond) cond = "u";
1330 mglFormula eq(cond);
1331 if(*i<0 || *i>=nx) *i=nx;
1332 if(*j<0 || *j>=ny) *j=ny-1;
1333 if(*k<0 || *k>=nz) *k=nz-1;
1334 long i0 = *i+nx*(*j+ny*(*k))-1;
1335 mreal x,y,z,dx=nx>1?1/(nx-1.):0,dy=ny>1?1/(ny-1.):0,dz=nz>1?1/(nz-1.):0;
1336 for(;i0<nx*ny*nz;i0--)
1337 {
1338 x = dx*(i0%nx); y = dy*((i0/nx)%ny); z = dz*(i0/(nx*ny));
1339 if(eq.Calc(x,y,z,d->vthr(i0))) break;
1340 }
1341 *i = i0%nx; *j = (i0/nx)%ny; *k = i0/(nx*ny);
1342 return i0<nx*ny*nz ? d->vthr(i0) : NAN; // NOTE: Return NAN if false
1343 }
mgl_data_first_(uintptr_t * d,const char * cond,int * i,int * j,int * k,int l)1344 mreal MGL_EXPORT mgl_data_first_(uintptr_t *d, const char *cond, int *i, int *j, int *k, int l)
1345 { long ii=*i,jj=*j,kk=*k; char *s=new char[l+1]; memcpy(s,cond,l); s[l]=0;
1346 mreal res = mgl_data_first(_DT_,s,&ii,&jj,&kk); *i=ii; *j=jj; *k=kk;
1347 delete []s; return res; }
1348 //-----------------------------------------------------------------------------
mgl_data_find(HCDT d,const char * cond,char dir,long i,long j,long k)1349 long MGL_EXPORT mgl_data_find(HCDT d, const char *cond, char dir, long i, long j, long k)
1350 {
1351 long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
1352 long m=-1;
1353 if(!cond) cond = "u";
1354 mglFormula eq(cond);
1355 mreal x=i/(nx-1.),y=j/(ny-1.),z=k/(nz-1.);
1356 if(dir=='x' && nx>1) for(m=i;m<nx;m++)
1357 if(eq.Calc(m/(nx-1.),y,z,d->v(m,j,k))) break;
1358 if(dir=='y' && ny>1) for(m=j;m<ny;m++)
1359 if(eq.Calc(x,m/(ny-1.),z,d->v(i,m,k))) break;
1360 if(dir=='z' && nz>1) for(m=k;m<nz;m++)
1361 if(eq.Calc(x,y,m/(nz-1.),d->v(i,j,m))) break;
1362 return m;
1363 }
mgl_data_find_(uintptr_t * d,const char * cond,char * dir,int * i,int * j,int * k,int l,int)1364 int MGL_EXPORT mgl_data_find_(uintptr_t *d, const char *cond, char *dir, int *i, int *j, int *k, int l, int)
1365 { char *s=new char[l+1]; memcpy(s,cond,l); s[l]=0;
1366 int res = mgl_data_find(_DT_,s,*dir,*i,*j,*k); delete []s; return res; }
1367 //-----------------------------------------------------------------------------
mgl_data_find_any(HCDT d,const char * cond)1368 int MGL_EXPORT mgl_data_find_any(HCDT d, const char *cond)
1369 {
1370 long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
1371 bool cc = false;
1372 if(!cond || *cond==0) cond = "u";
1373 mglFormula eq(cond);
1374 #pragma omp parallel for collapse(3)
1375 for(long k=0;k<nz;k++) for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
1376 {
1377 if(cc) continue;
1378 if(eq.Calc(i/(nx-1.),j/(ny-1.),k/(nz-1.),d->v(i,j,k))) cc = true;
1379 }
1380 return cc;
1381 }
mgl_data_find_any_(uintptr_t * d,const char * cond,int l)1382 int MGL_EXPORT mgl_data_find_any_(uintptr_t *d, const char *cond, int l)
1383 { char *s=new char[l+1]; memcpy(s,cond,l); s[l]=0;
1384 int res = mgl_data_find_any(_DT_,s); delete []s; return res; }
1385 //-----------------------------------------------------------------------------
mgl_data_momentum_val(HCDT dd,char dir,mreal * x,mreal * w,mreal * s,mreal * k)1386 mreal MGL_EXPORT mgl_data_momentum_val(HCDT dd, char dir, mreal *x, mreal *w, mreal *s, mreal *k)
1387 {
1388 long nx=dd->GetNx(),ny=dd->GetNy(),nz=dd->GetNz();
1389 mreal i0=0,i1=0,i2=0,i3=0,i4=0,is;
1390 switch(dir)
1391 {
1392 case 'x':
1393 #pragma omp parallel for reduction(+:i0,i1)
1394 for(long i=0;i<nx*ny*nz;i++)
1395 {
1396 mreal d = mreal(i%nx), v = dd->vthr(i);
1397 i0+= v; i1+= v*d;
1398 }
1399 is = i1/i0;
1400 #pragma omp parallel for reduction(+:i2,i3,i4)
1401 for(long i=0;i<nx*ny*nz;i++)
1402 {
1403 mreal d = mreal(i%nx)-is, t = d*d, v = dd->vthr(i);
1404 i2+= v*t; i3+= v*d*t; i4+= v*t*t;
1405 }
1406 break;
1407 case 'y':
1408 #pragma omp parallel for reduction(+:i0,i1)
1409 for(long i=0;i<nx*ny*nz;i++)
1410 {
1411 mreal d = mreal((i/nx)%ny), v = dd->vthr(i);
1412 i0+= v; i1+= v*d;
1413 }
1414 is = i1/i0;
1415 #pragma omp parallel for reduction(+:i2,i3,i4)
1416 for(long i=0;i<nx*ny*nz;i++)
1417 {
1418 mreal d = mreal((i/nx)%ny)-is, t = d*d, v = dd->vthr(i);
1419 i2+= v*t; i3+= v*d*t; i4+= v*t*t;
1420 }
1421 break;
1422 case 'z':
1423 #pragma omp parallel for reduction(+:i0,i1)
1424 for(long i=0;i<nx*ny*nz;i++)
1425 {
1426 mreal d = mreal(i/(nx*ny)), v = dd->vthr(i);
1427 i0+= v; i1+= v*d;
1428 }
1429 is = i1/i0;
1430 #pragma omp parallel for reduction(+:i2,i3,i4)
1431 for(long i=0;i<nx*ny*nz;i++)
1432 {
1433 mreal d = mreal(i/(nx*ny))-is, t = d*d, v = dd->vthr(i);
1434 i2+= v*t; i3+= v*d*t; i4+= v*t*t;
1435 }
1436 break;
1437 default: // "self-dispersion"
1438 i0 = nx*ny*nz;
1439 #pragma omp parallel for reduction(+:i1)
1440 for(long i=0;i<nx*ny*nz;i++) i1 += dd->vthr(i);
1441 is = i1/i0;
1442 #pragma omp parallel for reduction(+:i2,i3,i4)
1443 for(long i=0;i<nx*ny*nz;i++)
1444 {
1445 mreal v = dd->vthr(i)-is, t = v*v;
1446 i2+= t; i3+= v*t; i4+= t*t;
1447 }
1448 }
1449 if(i0==0) return 0;
1450 mreal w2=i2/i0, ww = sqrt(w2);
1451 if(x) *x=is;
1452 if(w) *w=ww;
1453 if(s) *s=i3/i0/ww/w2;
1454 if(k) *k=i4/(i0*3)/w2/w2;
1455 return i0;
1456 }
mgl_data_momentum_val_(uintptr_t * d,char * dir,mreal * m,mreal * w,mreal * s,mreal * k,int)1457 mreal MGL_EXPORT mgl_data_momentum_val_(uintptr_t *d, char *dir, mreal *m, mreal *w, mreal *s, mreal *k,int)
1458 { mreal mm=0,ww=0,ss=0,kk=0,aa=0;
1459 aa = mgl_data_momentum_val(_DT_,*dir,&mm,&ww,&ss,&kk);
1460 *m=mm; *w=ww; *s=ss; *k=kk; return aa; }
1461 //-----------------------------------------------------------------------------
mgl_data_norm_slice(HMDT d,mreal v1,mreal v2,char dir,long keep_en,long sym)1462 void MGL_EXPORT mgl_data_norm_slice(HMDT d, mreal v1,mreal v2,char dir,long keep_en,long sym)
1463 {
1464 long nx=d->nx, ny=d->ny, nz=d->nz;
1465 mreal *a=d->a;
1466 mglData b(*d);
1467 mreal e0=1, e, m1, m2, aa;
1468 if(sym) { v2 = -v1>v2 ? -v1:v2; v1 = -v2; }
1469 if(dir=='z' && nz>1)
1470 {
1471 //#pragma omp parallel for private(m1,m2,aa,e,e0) // TODO add omp comparison here
1472 for(long k=0;k<nz;k++)
1473 {
1474 m1 = INFINITY; m2 = -INFINITY; e=0;
1475 for(long i=0;i<nx*ny;i++)
1476 {
1477 aa = a[i+nx*ny*k];
1478 m1 = m1<aa ? m1 : aa;
1479 m2 = m2>aa ? m2 : aa;
1480 e += aa*aa;
1481 }
1482 if(m1==m2) m2+=1;
1483 if(sym) { m2 = -m1>m2 ? -m1:m2; m1 = -m2; }
1484 if(keep_en && k) e = e0>0?sqrt(e/e0):1;
1485 else { e0 = e; e=1; }
1486 for(long i=0;i<nx*ny;i++)
1487 b.a[i+nx*ny*k] = (v1 + (v2-v1)*(a[i+nx*ny*k]-m1)/(m2-m1))*e;
1488 }
1489 }
1490 else if(dir=='y' && ny>1)
1491 {
1492 //#pragma omp parallel for private(m1,m2,aa,e,e0) // TODO add omp comparison here
1493 for(long j=0;j<ny;j++)
1494 {
1495 m1 = INFINITY; m2 = -INFINITY; e=0;
1496 for(long k=0;k<nz;k++) for(long i=0;i<nx;i++)
1497 {
1498 aa = a[i+nx*(j+ny*k)];
1499 m1 = m1<aa ? m1 : aa;
1500 m2 = m2>aa ? m2 : aa;
1501 e += aa*aa;
1502 }
1503 if(m1==m2) m2+=1;
1504 if(sym) { m2 = -m1>m2 ? -m1:m2; m1 = -m2; }
1505 if(keep_en && j) e = e0>0?sqrt(e/e0):1;
1506 else { e0 = e; e=1; }
1507 #pragma omp parallel for collapse(2)
1508 for(long k=0;k<nz;k++) for(long i=0;i<nx;i++)
1509 b.a[i+nx*(j+ny*k)] = (v1 + (v2-v1)*(a[i+nx*(j+ny*k)]-m1)/(m2-m1))*e;
1510 }
1511 }
1512 else if(dir=='x' && nx>1)
1513 {
1514 //#pragma omp parallel for private(m1,m2,aa,e,e0) // TODO add omp comparison here
1515 for(long i=0;i<nx;i++)
1516 {
1517 m1 = INFINITY; m2 = -INFINITY; e=0;
1518 for(long k=0;k<ny*nz;k++)
1519 {
1520 aa = a[i+nx*k];
1521 m1 = m1<aa ? m1 : aa;
1522 m2 = m2>aa ? m2 : aa;
1523 e += aa*aa;
1524 }
1525 if(m1==m2) m2+=1;
1526 if(sym) { m2 = -m1>m2 ? -m1:m2; m1 = -m2; }
1527 if(keep_en && i) e = e0>0?sqrt(e/e0):1;
1528 else { e0 = e; e=1; }
1529 #pragma omp parallel for
1530 for(long k=0;k<ny*nz;k++)
1531 b.a[i+nx*k] = (v1 + (v2-v1)*(a[i+nx*k]-m1)/(m2-m1))*e;
1532 }
1533 }
1534 memcpy(d->a, b.a, nx*ny*nz*sizeof(mreal));
1535 }
mgl_data_norm_slice_(uintptr_t * d,mreal * v1,mreal * v2,char * dir,int * keep_en,int * sym,int)1536 void MGL_EXPORT mgl_data_norm_slice_(uintptr_t *d, mreal *v1,mreal *v2,char *dir,int *keep_en,int *sym,int )
1537 { mgl_data_norm_slice(_DT_,*v1,*v2,*dir,*keep_en,*sym); }
1538 //-----------------------------------------------------------------------------
mgl_data_info(HCDT d)1539 MGL_EXPORT const char *mgl_data_info(HCDT d) // NOTE: Not thread safe function!
1540 {
1541 static char buf[512];
1542 char s[128]; buf[0]=0;
1543 snprintf(s,128,"nx = %ld\tny = %ld\tnz = %ld\n",d->GetNx(),d->GetNy(),d->GetNz());
1544 s[127]=0; strcat(buf,s);
1545
1546 long i=0,j=0,k=0;
1547 mreal A=0,Wa=0,X=0,Y=0,Z=0,Wx=0,Wy=0,Wz=0, b;
1548 b = mgl_data_max_int(d,&i,&j,&k);
1549 snprintf(s,128,_("Maximum is %g\t at x = %ld\ty = %ld\tz = %ld\n"), b,i,j,k);
1550 s[127]=0; strcat(buf,s);
1551 b = mgl_data_min_int(d,&i,&j,&k);
1552 snprintf(s,128,_("Minimum is %g\t at x = %ld\ty = %ld\tz = %ld\n"), b,i,j,k);
1553 s[127]=0; strcat(buf,s);
1554
1555 mgl_data_momentum_val(d,'a',&A,&Wa,0,0); mgl_data_momentum_val(d,'x',&X,&Wx,0,0);
1556 mgl_data_momentum_val(d,'y',&Y,&Wy,0,0); mgl_data_momentum_val(d,'z',&Z,&Wz,0,0);
1557 snprintf(s,128,_("Averages are:\n<a> = %g\t<x> = %g\t<y> = %g\t<z> = %g\n"), A,X,Y,Z);
1558 s[127]=0; strcat(buf,s);
1559 snprintf(s,128,_("Widths are:\nWa = %g\tWx = %g\tWy = %g\tWz = %g\n"), Wa,Wx,Wy,Wz);
1560 s[127]=0; strcat(buf,s);
1561 return buf;
1562 }
mgl_data_info_(uintptr_t * d,char * out,int len)1563 int MGL_EXPORT mgl_data_info_(uintptr_t *d, char *out, int len)
1564 {
1565 const char *res = mgl_data_info(_DA_(d));
1566 if(out) mgl_strncpy(out,res,len);
1567 return strlen(res);
1568 }
1569 //-----------------------------------------------------------------------------
mgl_data_insert(HMDT d,char dir,long at,long num)1570 void MGL_EXPORT mgl_data_insert(HMDT d, char dir, long at, long num)
1571 {
1572 if(num<1) return;
1573 at = at<0 ? 0:at;
1574 long nn, nx=d->nx, ny=d->ny, nz=d->nz, nxy=nx*ny;
1575 mglData b;
1576 if(dir=='x')
1577 {
1578 if(at>nx) at=nx;
1579 nn=nx+num; b.Create(nn,ny,nz);
1580 #pragma omp parallel for
1581 for(long k=0;k<ny*nz;k++)
1582 {
1583 if(at>0) memcpy(b.a+nn*k, d->a+nx*k,at*sizeof(mreal));
1584 if(at<nx) memcpy(b.a+at+num+nn*k, d->a+at+nx*k,(nx-at)*sizeof(mreal));
1585 if(at<nx) for(long i=0;i<num;i++) b.a[nn*k+at+i]=d->a[nx*k+at]; // copy values
1586 else for(long i=0;i<num;i++) b.a[nn*k+at+i]=d->a[nx*k+nx-1]; // copy values
1587 }
1588 d->Set(b); nx+=num;
1589 }
1590 if(dir=='y')
1591 {
1592 if(at>ny) at=ny;
1593 nn=num+ny; b.Create(nx,nn,nz);
1594 #pragma omp parallel for
1595 for(long k=0;k<nz;k++)
1596 {
1597 if(at>0) memcpy(b.a+nx*nn*k, d->a+nxy*k,at*nx*sizeof(mreal));
1598 if(at<ny) memcpy(b.a+nx*(at+num+nn*k), d->a+nx*(at+ny*k),(ny-at)*nx*sizeof(mreal));
1599 if(at<ny) for(long i=0;i<num;i++)
1600 memcpy(b.a+nx*(nn*k+at+i),d->a+nx*(ny*k+at),nx*sizeof(mreal));
1601 else for(long i=0;i<num;i++)
1602 memcpy(b.a+nx*(nn*k+at+i),d->a+nx*(ny*k+ny-1),nx*sizeof(mreal));
1603 }
1604 d->Set(b); ny+=num;
1605 }
1606 if(dir=='z')
1607 {
1608 if(at>nz) at=nz;
1609 b.Create(nx,ny,nz+num);
1610 if(at>0) memcpy(b.a, d->a,at*nxy*sizeof(mreal));
1611 if(at<nz) memcpy(b.a+nxy*(at+num), d->a+nxy*at,(nz-at)*nxy*sizeof(mreal));
1612 if(at<nz)
1613 #pragma omp parallel for
1614 for(long i=0;i<num;i++) memcpy(b.a+nxy*(at+i),d->a+nxy*at,nxy*sizeof(mreal));
1615 else
1616 #pragma omp parallel for
1617 for(long i=0;i<num;i++) memcpy(b.a+nxy*(at+i),d->a+nxy*(nz-1),nxy*sizeof(mreal));
1618 d->Set(b);
1619 }
1620 }
1621 //-----------------------------------------------------------------------------
mgl_data_delete(HMDT d,char dir,long at,long num)1622 void MGL_EXPORT mgl_data_delete(HMDT d, char dir, long at, long num)
1623 {
1624 if(num<1 || at<0) return;
1625 mglData b;
1626 long nx=d->nx, ny=d->ny, nz=d->nz, nn;
1627 if(dir=='x')
1628 {
1629 if(at+num>=nx) return;
1630 nn=nx-num; b.Create(nn,ny,nz);
1631 #pragma omp parallel for
1632 for(long k=0;k<ny*nz;k++)
1633 {
1634 if(at>0) memcpy(b.a+nn*k, d->a+nx*k,at*sizeof(mreal));
1635 memcpy(b.a+at+nn*k, d->a+at+num+nx*k,(nx-at-num)*sizeof(mreal));
1636 }
1637 d->Set(b); nx-=num;
1638 }
1639 if(dir=='y')
1640 {
1641 if(at+num>=ny) return;
1642 nn=ny-num; b.Create(nx,nn,nz);
1643 #pragma omp parallel for
1644 for(long k=0;k<nz;k++)
1645 {
1646 if(at>0) memcpy(b.a+nx*nn*k, d->a+nx*ny*k,at*nx*sizeof(mreal));
1647 memcpy(b.a+nx*(at+nn*k), d->a+nx*(at+num+ny*k),(ny-at-num)*nx*sizeof(mreal));
1648 }
1649 d->Set(b); ny-=num;
1650 }
1651 if(dir=='z')
1652 {
1653 if(at+num>=nz) return;
1654 b.Create(nx,ny,nz-num);
1655 if(at>0) memcpy(b.a, d->a,at*nx*ny*sizeof(mreal));
1656 memcpy(b.a+nx*ny*at, d->a+nx*ny*(at+num),(nz-at-num)*nx*ny*sizeof(mreal));
1657 d->Set(b);
1658 }
1659 }
1660 //-----------------------------------------------------------------------------
mgl_data_insert_(uintptr_t * d,const char * dir,int * at,int * num,int)1661 void MGL_EXPORT mgl_data_insert_(uintptr_t *d, const char *dir, int *at, int *num, int)
1662 { mgl_data_insert(_DT_,*dir,*at,*num); }
mgl_data_delete_(uintptr_t * d,const char * dir,int * at,int * num,int)1663 void MGL_EXPORT mgl_data_delete_(uintptr_t *d, const char *dir, int *at, int *num, int)
1664 { mgl_data_delete(_DT_,*dir,*at,*num); }
1665 //-----------------------------------------------------------------------------
1666 #define omod(x,y) (y)*((x)>0?int((x)/(y)+0.5):int((x)/(y)-0.5))
mgl_omod(mreal * a,mreal da,int nx,int n)1667 void static mgl_omod(mreal *a, mreal da, int nx, int n)
1668 {
1669 bool qq=true;
1670 for(long i=1;i<nx;i++)
1671 {
1672 long ii = i*n;
1673 if(mgl_isnan(a[ii-n])) { qq=true; continue; }
1674 if(qq)
1675 {
1676 a[ii] += omod(a[ii-n]-a[ii], da);
1677 qq=false;
1678 }
1679 else
1680 {
1681 mreal q = 2*a[ii-n]-a[ii-2*n];
1682 a[ii] += omod(q-a[ii], da);
1683 }
1684 }
1685 }
1686 //-----------------------------------------------------------------------------
mgl_sew_z(void * par)1687 static void *mgl_sew_z(void *par)
1688 {
1689 mglThreadD *t=(mglThreadD *)par;
1690 long nz=t->p[2], nn=t->n;
1691 #if !MGL_HAVE_PTHREAD
1692 #pragma omp parallel for
1693 #endif
1694 for(long i=t->id;i<nn;i+=mglNumThr) mgl_omod(t->a+i, t->b[0], nz, nn);
1695 return 0;
1696 }
mgl_sew_y(void * par)1697 static void *mgl_sew_y(void *par)
1698 {
1699 mglThreadD *t=(mglThreadD *)par;
1700 long nx=t->p[0], ny=t->p[1], nn=t->n;
1701 #if !MGL_HAVE_PTHREAD
1702 #pragma omp parallel for
1703 #endif
1704 for(long i=t->id;i<nn;i+=mglNumThr) mgl_omod(t->a+(i%nx)+nx*ny*(i/nx), t->b[0], ny, nx);
1705 return 0;
1706 }
mgl_sew_x(void * par)1707 static void *mgl_sew_x(void *par)
1708 {
1709 mglThreadD *t=(mglThreadD *)par;
1710 long nx=t->p[0], nn=t->n;
1711 #if !MGL_HAVE_PTHREAD
1712 #pragma omp parallel for
1713 #endif
1714 for(long i=t->id;i<nn;i+=mglNumThr) mgl_omod(t->a+i*nx, t->b[0], nx, 1);
1715 return 0;
1716 }
mgl_data_sew(HMDT d,const char * dirs,mreal delta)1717 void MGL_EXPORT mgl_data_sew(HMDT d, const char *dirs, mreal delta)
1718 {
1719 if(!dirs || *dirs==0) return;
1720 long nx=d->nx, ny=d->ny, nz=d->nz;
1721 long p[3]={nx,ny,nz};
1722 mreal da = delta;
1723 if(strchr(dirs,'x') && nx>1) mglStartThread(mgl_sew_x,0,nz*ny,d->a,&da,0,p);
1724 if(strchr(dirs,'y') && ny>1) mglStartThread(mgl_sew_y,0,nz*nx,d->a,&da,0,p);
1725 if(strchr(dirs,'z') && nz>1) mglStartThread(mgl_sew_z,0,nx*ny,d->a,&da,0,p);
1726 }
mgl_data_sew_(uintptr_t * d,const char * dirs,mreal * da,int l)1727 void MGL_EXPORT mgl_data_sew_(uintptr_t *d, const char *dirs, mreal *da, int l)
1728 { char *s=new char[l+1]; memcpy(s,dirs,l); s[l]=0;
1729 mgl_data_sew(_DT_,s,*da); delete []s; }
1730 //-----------------------------------------------------------------------------
mgl_data_put_val(HMDT d,mreal val,long xx,long yy,long zz)1731 void MGL_EXPORT mgl_data_put_val(HMDT d, mreal val, long xx, long yy, long zz)
1732 {
1733 long nx=d->nx, ny=d->ny, nz=d->nz;
1734 if(xx>=nx || yy>=ny || zz>=nz) return;
1735 mreal *a=d->a;
1736 if(xx>=0 && yy>=0 && zz>=0) a[xx+nx*(yy+zz*ny)] = val;
1737 else if(xx<0 && yy<0 && zz<0)
1738 #pragma omp parallel for
1739 for(long i=0;i<nx*ny*nz;i++) a[i] = val;
1740 else if(xx<0 && yy<0)
1741 #pragma omp parallel for
1742 for(long i=0;i<nx*ny;i++) a[i+zz*nx*ny] = val;
1743 else if(yy<0 && zz<0)
1744 #pragma omp parallel for
1745 for(long i=0;i<nz*ny;i++) a[xx+i*nx] = val;
1746 else if(xx<0 && zz<0)
1747 #pragma omp parallel for collapse(2)
1748 for(long j=0;j<nz;j++) for(long i=0;i<nx;i++) a[i+nx*(yy+j*ny)] = val;
1749 else if(xx<0)
1750 #pragma omp parallel for
1751 for(long i=0;i<nx;i++) a[i+nx*(yy+zz*ny)] = val;
1752 else if(yy<0)
1753 #pragma omp parallel for
1754 for(long i=0;i<ny;i++) a[xx+nx*(i+zz*ny)] = val;
1755 else //if(zz<0)
1756 #pragma omp parallel for
1757 for(long i=0;i<nz;i++) a[xx+nx*(yy+i*ny)] = val;
1758 }
1759 //-----------------------------------------------------------------------------
mgl_data_put_dat(HMDT d,HCDT v,long xx,long yy,long zz)1760 void MGL_EXPORT mgl_data_put_dat(HMDT d, HCDT v, long xx, long yy, long zz)
1761 {
1762 long nx=d->nx, ny=d->ny, nz=d->nz;
1763 if(xx>=nx || yy>=ny || zz>=nz) return;
1764 const mglData *mv = dynamic_cast<const mglData *>(v);
1765 mreal *a=d->a, vv=v->v(0);
1766 const mreal *b = mv?mv->a:0;
1767 long vx=v->GetNx(), vy=v->GetNy(), vz=v->GetNz();
1768 if(xx<0 && yy<0 && zz<0) // whole array
1769 {
1770 long nn = vx>=nx?nx:vx, mm = vy>=ny?ny:vy, ll = vz>=nz?nz:vz;
1771 if(nn>1 && mm>1 && ll>1)
1772 // #pragma omp parallel for
1773 for(long k=0;k<ll;k++) for(long j=0;j<mm;j++) for(long i=0;i<nn;i++)
1774 a[i+nx*(j+k*ny)] = b?b[i+vx*(j+k*vy)]:v->v(i,j,k);
1775 if(nn>1 && mm>1)
1776 // #pragma omp parallel for
1777 for(long k=0;k<nz;k++) for(long j=0;j<mm;j++) for(long i=0;i<nn;i++)
1778 a[i+nx*(j+k*ny)] = b?b[i+vx*j]:v->v(i,j);
1779 else if(nn>1)
1780 // #pragma omp parallel for
1781 for(long k=0;k<nz;k++) for(long j=0;j<ny;j++) for(long i=0;i<nn;i++)
1782 a[i+nx*(j+k*ny)] = b?b[i]:v->v(i);
1783 else
1784 // #pragma omp parallel for
1785 for(long ii=0;ii<nx*ny*nz;ii++) a[ii] = vv;
1786 }
1787 else if(xx<0 && yy<0) // 2d
1788 {
1789 zz*=nx*ny;
1790 long nn = vx>=nx?nx:vx, mm = vy>=ny?ny:vy;
1791 if(nn>1 && mm>1)
1792 // #pragma omp parallel for
1793 for(long j=0;j<mm;j++) for(long i=0;i<nn;i++)
1794 a[i+j*nx+zz] = b?b[i+vx*j]:v->v(i,j);
1795 else if(nn>1)
1796 // #pragma omp parallel for
1797 for(long j=0;j<ny;j++) for(long i=0;i<nn;i++)
1798 a[i+j*nx+xx] = b?b[i]:v->v(i);
1799 else
1800 // #pragma omp parallel for
1801 for(long ii=0;ii<nx*ny;ii++) a[ii+zz] = vv;
1802 }
1803 else if(yy<0 && zz<0) // 2d
1804 {
1805 long nn = vx>=ny?ny:vx, mm = vy>=nz?nz:vy;
1806 if(nn>1 && mm>1)
1807 // #pragma omp parallel for collapse(2)
1808 for(long j=0;j<mm;j++) for(long i=0;i<nn;i++)
1809 a[(i+ny*j)*nx+xx] = b?b[i+vx*j]:v->v(i,j);
1810 else if(nn>1)
1811 // #pragma omp parallel for collapse(2)
1812 for(long j=0;j<nz;j++) for(long i=0;i<nn;i++)
1813 a[(i+ny*j)*nx+xx] = b?b[i]:v->v(i);
1814 else
1815 // #pragma omp parallel for
1816 for(long ii=0;ii<ny*nz;ii++) a[ii*nx+xx] = vv;
1817 }
1818 else if(xx<0 && zz<0) // 2d
1819 {
1820 yy *= nx; zz = nx*ny;
1821 long nn = vx>=nx?nx:vx, mm = vy>=nz?nz:vy;
1822 if(nn>1 && mm>1)
1823 // #pragma omp parallel for collapse(2)
1824 for(long j=0;j<mm;j++) for(long i=0;i<nn;i++)
1825 a[i+yy+j*zz] = b?b[i+vx*j]:v->v(i,j);
1826 else if(nn>1)
1827 // #pragma omp parallel for collapse(2)
1828 for(long j=0;j<nz;j++) for(long i=0;i<nn;i++)
1829 a[i+yy+j*zz] = b?b[i]:v->v(i);
1830 else
1831 // #pragma omp parallel for
1832 for(long ii=0;ii<nx*nz;ii++)
1833 { long i=ii%nx, j=ii/nx; a[i+yy+j*zz] = vv; }
1834 }
1835 else if(xx<0)
1836 {
1837 xx = nx*(yy+zz*ny);
1838 long nn = vx>=nx?nx:vx;
1839 if(nn>1)
1840 // #pragma omp parallel for
1841 for(long i=0;i<nn;i++) a[i+xx] = b?b[i]:v->v(i);
1842 else
1843 // #pragma omp parallel for
1844 for(long i=0;i<nx;i++) a[i+xx] = vv;
1845 }
1846 else if(yy<0)
1847 {
1848 xx += zz*nx*ny;
1849 long nn = vx>=ny?ny:vx;
1850 if(nn>1)
1851 // #pragma omp parallel for
1852 for(long i=0;i<nn;i++) a[xx+nx*i] = b?b[i]:v->v(i);
1853 else
1854 // #pragma omp parallel for
1855 for(long i=0;i<ny;i++) a[xx+nx*i] = vv;
1856 }
1857 else if(zz<0)
1858 {
1859 xx += nx*yy; yy = nx*ny;
1860 long nn = vx>=nz?nz:vx;
1861 if(nn>1)
1862 // #pragma omp parallel for
1863 for(long i=0;i<nn;i++) a[xx+yy*i] = b?b[i]:v->v(i);
1864 else
1865 // #pragma omp parallel for
1866 for(long i=0;i<nz;i++) a[xx+yy*i] = vv;
1867 }
1868 else a[xx+nx*(yy+ny*zz)] = vv;
1869 }
1870 //-----------------------------------------------------------------------------
mgl_data_put_val_(uintptr_t * d,mreal * val,int * i,int * j,int * k)1871 void MGL_EXPORT mgl_data_put_val_(uintptr_t *d, mreal *val, int *i, int *j, int *k)
1872 { mgl_data_put_val(_DT_,*val, *i,*j,*k); }
mgl_data_put_dat_(uintptr_t * d,uintptr_t * val,int * i,int * j,int * k)1873 void MGL_EXPORT mgl_data_put_dat_(uintptr_t *d, uintptr_t *val, int *i, int *j, int *k)
1874 { mgl_data_put_dat(_DT_,_DA_(val), *i,*j,*k); }
1875 //-----------------------------------------------------------------------------
mgl_diff_3(void * par)1876 static void *mgl_diff_3(void *par)
1877 {
1878 mglThreadD *t=(mglThreadD *)par;
1879 long nx=t->p[0], ny=t->p[1], nz=t->p[2], nn=t->n, n2=nx*ny;
1880 mreal *b=t->a,au,av,aw,xu,xv,xw,yu,yv,yw,zu,zv,zw;
1881 HCDT x=(HCDT)(t->c),y=(HCDT)(t->d),z=(HCDT)(t->e);
1882 const mreal *a=t->b;
1883 #if !MGL_HAVE_PTHREAD
1884 #pragma omp parallel for private(au,av,aw,xu,xv,xw,yu,yv,yw,zu,zv,zw)
1885 #endif
1886 for(long i0=t->id;i0<nn;i0+=mglNumThr)
1887 {
1888 long i=i0%nx, j=((i0/nx)%ny), k=i0/(nx*ny);
1889 if(i==0)
1890 {
1891 au = 3*a[i0]-4*a[i0+1]+a[i0+2];
1892 xu = 3*x->vthr(i0)-4*x->vthr(i0+1)+x->vthr(i0+2);
1893 yu = 3*y->vthr(i0)-4*y->vthr(i0+1)+y->vthr(i0+2);
1894 zu = 3*z->vthr(i0)-4*z->vthr(i0+1)+z->vthr(i0+2);
1895 }
1896 else if(i==nx-1)
1897 {
1898 au = 3*a[i0]-4*a[i0-1]+a[i0-2];
1899 xu = 3*x->vthr(i0)-4*x->vthr(i0-1)+x->vthr(i0-2);
1900 yu = 3*y->vthr(i0)-4*y->vthr(i0-1)+y->vthr(i0-2);
1901 zu = 3*z->vthr(i0)-4*z->vthr(i0-1)+z->vthr(i0-2);
1902 }
1903 else
1904 {
1905 au = a[i0+1]-a[i0-1];
1906 xu = x->vthr(i0+1)-x->vthr(i0-1);
1907 yu = y->vthr(i0+1)-y->vthr(i0-1);
1908 zu = z->vthr(i0+1)-z->vthr(i0-1);
1909 }
1910 if(j==0)
1911 {
1912 av = 3*a[i0]-4*a[i0+nx]+a[i0+2*nx];
1913 xv = 3*x->vthr(i0)-4*x->vthr(i0+nx)+x->vthr(i0+2*nx);
1914 yv = 3*y->vthr(i0)-4*y->vthr(i0+nx)+y->vthr(i0+2*nx);
1915 zv = 3*z->vthr(i0)-4*z->vthr(i0+nx)+z->vthr(i0+2*nx);
1916 }
1917 else if(j==ny-1)
1918 {
1919 av = 3*a[i0]-4*a[i0-nx]+a[i0+(ny-3)*nx];
1920 xv = 3*x->vthr(i0)-4*x->vthr(i0-nx)+x->vthr(i0-2*nx);
1921 yv = 3*y->vthr(i0)-4*y->vthr(i0-nx)+y->vthr(i0-2*nx);
1922 zv = 3*z->vthr(i0)-4*z->vthr(i0-nx)+z->vthr(i0-2*nx);
1923 }
1924 else
1925 {
1926 av = a[i0+nx]-a[i0-nx];
1927 xv = x->vthr(i0+nx)-x->vthr(i0-nx);
1928 yv = y->vthr(i0+nx)-y->vthr(i0-nx);
1929 zv = z->vthr(i0+nx)-z->vthr(i0-nx);
1930 }
1931 if(k==0)
1932 {
1933 aw = 3*a[i0]-4*a[i0+n2]+a[i0+2*n2];
1934 xw = 3*x->vthr(i0)-4*x->vthr(i0+n2)+x->vthr(i0+2*n2);
1935 yw = 3*y->vthr(i0)-4*y->vthr(i0+n2)+y->vthr(i0+2*n2);
1936 zw = 3*z->vthr(i0)-4*z->vthr(i0+n2)+z->vthr(i0+2*n2);
1937 }
1938 else if(k==nz-1)
1939 {
1940 aw = 3*a[i0]-4*a[i0-n2]+a[i0-2*n2];
1941 xw = 3*x->vthr(i0)-4*x->vthr(i0-n2)+x->vthr(i0-2*n2);
1942 yw = 3*y->vthr(i0)-4*y->vthr(i0-n2)+y->vthr(i0-2*n2);
1943 zw = 3*z->vthr(i0)-4*z->vthr(i0-n2)+z->vthr(i0-2*n2);
1944 }
1945 else
1946 {
1947 aw = a[i0+n2]-a[i0-n2];
1948 xw = x->vthr(i0+n2)-x->vthr(i0-n2);
1949 yw = y->vthr(i0+n2)-y->vthr(i0-n2);
1950 zw = z->vthr(i0+n2)-z->vthr(i0-n2);
1951 }
1952 b[i0] = (au*yv*zw-av*yu*zw-au*yw*zv+aw*yu*zv+av*yw*zu-aw*yv*zu) / (xu*yv*zw-xv*yu*zw-xu*yw*zv+xw*yu*zv+xv*yw*zu-xw*yv*zu);
1953 }
1954 return 0;
1955 }
mgl_diff_2(void * par)1956 static void *mgl_diff_2(void *par)
1957 {
1958 mglThreadD *t=(mglThreadD *)par;
1959 long nx=t->p[0], ny=t->p[1], nn=t->n, same=t->p[2];
1960 mreal *b=t->a,au,av,xu,xv,yu,yv;
1961 HCDT x=(HCDT)(t->c),y=(HCDT)(t->d);
1962 const mreal *a=t->b;
1963 #if !MGL_HAVE_PTHREAD
1964 #pragma omp parallel for private(au,av,xu,xv,yu,yv)
1965 #endif
1966 for(long i0=t->id;i0<nn;i0+=mglNumThr)
1967 {
1968 long i=i0%nx, j=((i0/nx)%ny), i1 = same ? i0 : i0%(nx*ny);
1969 if(i==0)
1970 {
1971 au = 3*a[i0]-4*a[i0+1]+a[i0+2];
1972 xu = 3*x->vthr(i1)-4*x->vthr(i1+1)+x->vthr(i1+2);
1973 yu = 3*y->vthr(i1)-4*y->vthr(i1+1)+y->vthr(i1+2);
1974 }
1975 else if(i==nx-1)
1976 {
1977 au = 3*a[i0]-4*a[i0-1]+a[i0-2];
1978 xu = 3*x->vthr(i1)-4*x->vthr(i1-1)+x->vthr(i1-2);
1979 yu = 3*y->vthr(i1)-4*y->vthr(i1-1)+y->vthr(i1-2);
1980 }
1981 else
1982 {
1983 au = a[i0+1]-a[i0-1];
1984 xu = x->vthr(i1+1)-x->vthr(i1-1);
1985 yu = y->vthr(i1+1)-y->vthr(i1-1);
1986 }
1987 if(j==0)
1988 {
1989 av = 3*a[i0]-4*a[i0+nx]+a[i0+2*nx];
1990 xv = 3*x->vthr(i1)-4*x->vthr(i1+nx)+x->vthr(i1+2*nx);
1991 yv = 3*y->vthr(i1)-4*y->vthr(i1+nx)+y->vthr(i1+2*nx);
1992 }
1993 else if(j==ny-1)
1994 {
1995 av = 3*a[i0]-4*a[i0-nx]+a[i0-2*nx];
1996 xv = 3*x->vthr(i1)-4*x->vthr(i1-nx)+x->vthr(i1-2*nx);
1997 yv = 3*y->vthr(i1)-4*y->vthr(i1-nx)+y->vthr(i1-2*nx);
1998 }
1999 else
2000 {
2001 av = a[i0+nx]-a[i0-nx];
2002 xv = x->vthr(i1+nx)-x->vthr(i1-nx);
2003 yv = y->vthr(i1+nx)-y->vthr(i1-nx);
2004 }
2005 b[i0] = (av*yu-au*yv)/(xv*yu-xu*yv);
2006 }
2007 return 0;
2008 }
mgl_diff_1(void * par)2009 static void *mgl_diff_1(void *par)
2010 {
2011 mglThreadD *t=(mglThreadD *)par;
2012 long nx=t->p[0], nn=t->n, same=t->p[1];
2013 mreal *b=t->a,au,xu;
2014 HCDT x=(HCDT)(t->c);
2015 const mreal *a=t->b;
2016 #if !MGL_HAVE_PTHREAD
2017 #pragma omp parallel for private(au,xu)
2018 #endif
2019 for(long i0=t->id;i0<nn;i0+=mglNumThr)
2020 {
2021 long i=i0%nx, i1 = same ? i0 : i;
2022 if(i==0)
2023 {
2024 au = 3*a[i0]-4*a[i0+1]+a[i0+2];
2025 xu = 3*x->vthr(i1)-4*x->vthr(i1+1)+x->vthr(i1+2);
2026 }
2027 else if(i==nx-1)
2028 {
2029 au = 3*a[i0]-4*a[i0-1]+a[i0-2];
2030 xu = 3*x->vthr(i1)-4*x->vthr(i1-1)+x->vthr(i1-2);
2031 }
2032 else
2033 {
2034 au = a[i0+1]-a[i0-1];
2035 xu = x->vthr(i1+1)-x->vthr(i1-1);
2036 }
2037 b[i0] = au/xu;
2038 }
2039 return 0;
2040 }
mgl_data_diff_par(HMDT d,HCDT x,HCDT y,HCDT z)2041 void MGL_EXPORT mgl_data_diff_par(HMDT d, HCDT x, HCDT y, HCDT z)
2042 {
2043 long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz(), nn=nx*ny*nz;
2044 if(nx<2 || ny<2) return;
2045 mreal *b = new mreal[nn]; memset(b,0,nn*sizeof(mreal));
2046 long p[3]={nx,ny,nz};
2047
2048 if(x&&y&&z && x->GetNN()==nn && y->GetNN()==nn && z->GetNN()==nn)
2049 mglStartThread(mgl_diff_3,0,nn,b,d->a,(const mreal *)x,p,0,(const mreal *)y,(const mreal *)z);
2050 else if(x&&y && x->GetNx()*x->GetNy()==nx*ny && y->GetNx()*y->GetNy()==nx*ny)
2051 {
2052 p[2]=(x->GetNz()==nz && y->GetNz()==nz);
2053 mglStartThread(mgl_diff_2,0,nn,b,d->a,(const mreal *)x,p,0,(const mreal *)y);
2054 }
2055 else if(x && x->GetNx()==nx)
2056 {
2057 p[1]=(x->GetNy()*x->GetNz()==ny*nz);
2058 mglStartThread(mgl_diff_1,0,nn,b,d->a,(const mreal *)x,p,0,0);
2059 }
2060 memcpy(d->a,b,nn*sizeof(mreal)); delete []b;
2061 }
mgl_data_diff_par_(uintptr_t * d,uintptr_t * v1,uintptr_t * v2,uintptr_t * v3)2062 void MGL_EXPORT mgl_data_diff_par_(uintptr_t *d, uintptr_t *v1, uintptr_t *v2, uintptr_t *v3)
2063 { mgl_data_diff_par(_DT_,_DA_(v1),_DA_(v2),_DA_(v3)); }
2064 //-----------------------------------------------------------------------------
mgl_data_set_value(HMDT dat,mreal v,long i,long j,long k)2065 void MGL_EXPORT mgl_data_set_value(HMDT dat, mreal v, long i, long j, long k)
2066 { if(i>=0 && i<dat->nx && j>=0 && j<dat->ny && k>=0 && k<dat->nz) dat->a[i+dat->nx*(j+dat->ny*k)]=v; }
mgl_data_set_value_(uintptr_t * d,mreal * v,int * i,int * j,int * k)2067 void MGL_EXPORT mgl_data_set_value_(uintptr_t *d, mreal *v, int *i, int *j, int *k)
2068 { mgl_data_set_value(_DT_,*v,*i,*j,*k); }
2069 //-----------------------------------------------------------------------------
mgl_data_get_value(HCDT dat,long i,long j,long k)2070 mreal MGL_EXPORT mgl_data_get_value(HCDT dat, long i, long j, long k)
2071 { long nx = dat->GetNx(), ny = dat->GetNy();
2072 return (i>=0 && i<nx && j>=0 && j<ny && k>=0 && k<dat->GetNz()) ? dat->vthr(i+nx*(j+ny*k)):NAN; }
mgl_data_get_value_(uintptr_t * d,int * i,int * j,int * k)2073 mreal MGL_EXPORT mgl_data_get_value_(uintptr_t *d, int *i, int *j, int *k)
2074 { return mgl_data_get_value(_DA_(d),*i,*j,*k); }
2075 //-----------------------------------------------------------------------------
mgl_data_data(HMDT dat)2076 MGL_EXPORT_PURE mreal *mgl_data_data(HMDT dat) { return dat->a; }
2077 //-----------------------------------------------------------------------------
mgl_data_value(HMDT dat,long i,long j,long k)2078 MGL_EXPORT mreal *mgl_data_value(HMDT dat, long i,long j,long k)
2079 { long ii=i*dat->nx*(j+dat->ny*k);
2080 return ii>=0 && ii<dat->GetNN() ? dat->a+ii : 0; }
2081 //-----------------------------------------------------------------------------
mgl_data_get_nx(HCDT dat)2082 long MGL_EXPORT mgl_data_get_nx(HCDT dat) { return dat->GetNx(); }
mgl_data_get_ny(HCDT dat)2083 long MGL_EXPORT mgl_data_get_ny(HCDT dat) { return dat->GetNy(); }
mgl_data_get_nz(HCDT dat)2084 long MGL_EXPORT mgl_data_get_nz(HCDT dat) { return dat->GetNz(); }
mgl_data_get_nx_(uintptr_t * d)2085 long MGL_EXPORT mgl_data_get_nx_(uintptr_t *d) { return _DA_(d)->GetNx(); }
mgl_data_get_ny_(uintptr_t * d)2086 long MGL_EXPORT mgl_data_get_ny_(uintptr_t *d) { return _DA_(d)->GetNy(); }
mgl_data_get_nz_(uintptr_t * d)2087 long MGL_EXPORT mgl_data_get_nz_(uintptr_t *d) { return _DA_(d)->GetNz(); }
2088 //-----------------------------------------------------------------------------
mgl_data_join(HMDT d,HCDT v)2089 void MGL_EXPORT mgl_data_join(HMDT d, HCDT v)
2090 {
2091 if(!d || !v) return;
2092 long nx=d->nx, ny=d->ny, nz=d->nz, k=nx*ny*nz;
2093 const mglData *mv = dynamic_cast<const mglData *>(v);
2094 long vx=v->GetNx(), vy=v->GetNy(), vz=v->GetNz(), m = vx*vy*vz;
2095
2096 if(nx==vx && ny==vy && ny>1) d->nz += vz;
2097 else
2098 {
2099 ny *= nz; vy *= vz;
2100 if(nx==vx && nx>1)
2101 { d->nz = 1; d->ny = ny+vy; }
2102 else
2103 { d->ny = d->nz = 1; d->nx = k+m; }
2104 }
2105 mreal *b = new mreal[k+m];
2106 memcpy(b,d->a,k*sizeof(mreal));
2107 if(mv) memcpy(b+k,mv->a,m*sizeof(mreal));
2108 else
2109 #pragma omp parallel for
2110 for(long i=0;i<m;i++) b[k+i] = v->vthr(i);
2111 if(!d->link) delete []d->a;
2112 d->a = b; d->link=false; d->NewId();
2113 }
2114 //-----------------------------------------------------------------------------
mgl_data_join_(uintptr_t * d,uintptr_t * val)2115 void MGL_EXPORT mgl_data_join_(uintptr_t *d, uintptr_t *val)
2116 { mgl_data_join(_DT_,_DA_(val)); }
2117 //-----------------------------------------------------------------------------
mgl_data_refill_gs(HMDT dat,HCDT xdat,HCDT vdat,mreal x1,mreal x2,long sl)2118 void MGL_EXPORT mgl_data_refill_gs(HMDT dat, HCDT xdat, HCDT vdat, mreal x1, mreal x2, long sl)
2119 {
2120 HMDT coef = mgl_gspline_init(xdat, vdat);
2121 if(!coef) return; // incompatible dimensions
2122 const long nx = dat->nx, nn=dat->ny*dat->nz;
2123 mreal x0 = x1-xdat->v(0), dx = (x2-x1)/(nx-1);
2124 #pragma omp parallel for
2125 for(long i=0;i<nx;i++)
2126 {
2127 mreal d = mgl_gspline(coef,x0+dx*i,0,0);
2128 if(sl<0) for(long j=0;j<nn;j++) dat->a[i+j*nx] = d;
2129 else dat->a[i+sl*nx] = d;
2130 }
2131 mgl_delete_data(coef);
2132 }
2133 //-----------------------------------------------------------------------------
mgl_index_1(mreal v,HCDT dat)2134 mreal MGL_NO_EXPORT mgl_index_1(mreal v, HCDT dat)
2135 {
2136 long mx=dat->GetNx();
2137 mreal d,d1=0,d2=mx-1,v1,v2;
2138 v1 = dat->value(d1,0,0);
2139 v2 = dat->value(d2,0,0);
2140 long count=0;
2141
2142 const mreal eps = MGL_EPSILON-1.;
2143 if(fabs(v-v1)<eps) return d1;
2144 if(fabs(v-v2)<eps) return d2;
2145 if((v1-v)*(v2-v)>0) return NAN;
2146 do
2147 {
2148 d = count<10?(d2-d1)*(v-v1)/(v2-v1)+d1:(d1+d2)/2; count++;
2149 mreal val = dat->value(d,0,0);
2150 // if(fabs(val-v)<acx) break;
2151 if(val==v || d2-d<eps) break;
2152 if((v1-v)*(val-v)<0) { v2=val; d2=d; } else { v1=val; d1=d; }
2153 } while(fabs(d2-d1)>1e-5);
2154 return d;
2155 }
2156 //-----------------------------------------------------------------------------
mgl_data_refill_x(HMDT dat,HCDT xdat,HCDT vdat,mreal x1,mreal x2,long sl)2157 void MGL_EXPORT mgl_data_refill_x(HMDT dat, HCDT xdat, HCDT vdat, mreal x1, mreal x2, long sl)
2158 {
2159 long nx=dat->nx,mx=vdat->GetNx(),nn=dat->ny*dat->nz;
2160 if(mx!=xdat->GetNx()) return; // incompatible dimensions
2161 mreal dx = (x2-x1)/(nx-1);
2162 #pragma omp parallel for
2163 for(long i=0;i<nx;i++)
2164 {
2165 mreal u = mgl_index_1(x1+dx*i,xdat);
2166 mreal d = vdat->value(u,0,0);
2167 if(sl<0) for(long j=0;j<nn;j++) dat->a[i+j*nx] = d;
2168 else dat->a[i+sl*nx] = d;
2169 }
2170 }
2171 //-----------------------------------------------------------------------------
mgl_data_refill_xy(HMDT dat,HCDT xdat,HCDT ydat,HCDT vdat,mreal x1,mreal x2,mreal y1,mreal y2,long sl)2172 void MGL_EXPORT mgl_data_refill_xy(HMDT dat, HCDT xdat, HCDT ydat, HCDT vdat, mreal x1, mreal x2, mreal y1, mreal y2, long sl)
2173 {
2174 long nx=dat->nx,ny=dat->ny,nz=dat->nz,mx=vdat->GetNx(),my=vdat->GetNy(),nn=nx*ny;
2175 bool both=(xdat->GetNN()==vdat->GetNN() && ydat->GetNN()==vdat->GetNN());
2176 if(!both && (xdat->GetNx()!=mx || ydat->GetNx()!=my)) return; // incompatible dimensions
2177 mreal dx = (x2-x1)/(nx-1), dy = (y2-y1)/(ny-1);
2178 if(both)
2179 {
2180 #pragma omp parallel for
2181 for(long i=0;i<nn*nz;i++) dat->a[i]=NAN;
2182 #pragma omp parallel for collapse(2)
2183 for(long j=0;j<my-1;j++) for(long i=0;i<mx-1;i++)
2184 {
2185 long i0 = i+mx*j;
2186 mreal vx0 = (xdat->vthr(i0)-x1)/dx, vy0 = (ydat->vthr(i0)-y1)/dy;
2187 mreal vx1 = (xdat->vthr(i0+1)-x1)/dx, vy1 = (ydat->vthr(i0+1)-y1)/dy;
2188 mreal vx2 = (xdat->vthr(i0+mx)-x1)/dx, vy2 = (ydat->vthr(i0+mx)-y1)/dy;
2189 mreal vx3 = (xdat->vthr(i0+mx+1)-x1)/dx, vy3 = (ydat->vthr(i0+mx+1)-y1)/dy;
2190 long xx1 = long(mgl_min( mgl_min(vx0,vx1), mgl_min(vx2,vx3) )); // bounding box
2191 long yy1 = long(mgl_min( mgl_min(vy0,vy1), mgl_min(vy2,vy3) ));
2192 long xx2 = long(mgl_max( mgl_max(vx0,vx1), mgl_max(vx2,vx3) ));
2193 long yy2 = long(mgl_max( mgl_max(vy0,vy1), mgl_max(vy2,vy3) ));
2194 xx1=mgl_imax(xx1,0); xx2=mgl_imin(xx2,nx-1);
2195 yy1=mgl_imax(yy1,0); yy2=mgl_imin(yy2,ny-1);
2196 if(xx1>xx2 || yy1>yy2) continue;
2197
2198 mreal d1x = vx1-vx0, d1y = vy1-vy0;
2199 mreal d2x = vx2-vx0, d2y = vy2-vy0;
2200 mreal d3x = vx3+vx0-vx1-vx2, d3y = vy3+vy0-vy1-vy2;
2201 mreal dd = d1x*d2y-d1y*d2x;
2202 mreal dsx =-4*(d2y*d3x - d2x*d3y)*d1y;
2203 mreal dsy = 4*(d2y*d3x - d2x*d3y)*d1x;
2204
2205 for(long jj=yy1;jj<=yy2;jj++) for(long ii=xx1;ii<=xx2;ii++)
2206 {
2207 mreal xx = (ii-vx0), yy = (jj-vy0);
2208 mreal s = dsx*xx + dsy*yy + (dd+d3y*xx-d3x*yy)*(dd+d3y*xx-d3x*yy);
2209 if(s>=0)
2210 {
2211 s = sqrt(s);
2212 mreal qu = d3x*yy - d3y*xx + dd + s;
2213 mreal qv = d3y*xx - d3x*yy + dd + s;
2214 mreal u = 2.f*(d2y*xx - d2x*yy)/qu;
2215 mreal v = 2.f*(d1x*yy - d1y*xx)/qv;
2216 if(u*(1.f-u)<0.f || v*(1.f-v)<0.f) // first root bad
2217 {
2218 qu = d3x*yy - d3y*xx + dd - s;
2219 qv = d3y*xx - d3x*yy + dd - s;
2220 u = 2.f*(d2y*xx - d2x*yy)/qu;
2221 v = 2.f*(d1x*yy - d1y*xx)/qv;
2222 if(u*(1.f-u)<0.f || v*(1.f-v)<0.f) continue; // second root bad
2223 }
2224 i0 = ii+nx*jj; s = vdat->value(i+u,j+v,0);
2225 if(sl<0) for(long k=0;k<nz;k++) dat->a[i0+k*nn] = s;
2226 else dat->a[i0+sl*nn] = s;
2227 }
2228 }
2229 }
2230 }
2231 else
2232 {
2233 mglData u(nx), v(ny);
2234 #pragma omp parallel for
2235 for(long i=0;i<nx;i++) u.a[i] = mgl_index_1(x1+dx*i,xdat);
2236 #pragma omp parallel for
2237 for(long i=0;i<ny;i++) v.a[i] = mgl_index_1(y1+dy*i,ydat);
2238 #pragma omp parallel for collapse(2)
2239 for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
2240 {
2241 mreal d = vdat->value(u.a[i],v.a[j],0);
2242 long i0=i+nx*j;
2243 if(sl<0) for(long k=0;k<nz;k++) dat->a[i0+k*nn] = d;
2244 else dat->a[i0+sl*nn] = d;
2245 }
2246 }
2247 }
2248 //-----------------------------------------------------------------------------
mgl_data_refill_xyz(HMDT dat,HCDT xdat,HCDT ydat,HCDT zdat,HCDT vdat,mreal x1,mreal x2,mreal y1,mreal y2,mreal z1,mreal z2)2249 void MGL_EXPORT mgl_data_refill_xyz(HMDT dat, HCDT xdat, HCDT ydat, HCDT zdat, HCDT vdat, mreal x1, mreal x2, mreal y1, mreal y2, mreal z1, mreal z2)
2250 {
2251 if(!dat || !xdat || !ydat || !zdat || !vdat) return;
2252 long nx=dat->nx,ny=dat->ny,nz=dat->nz,mx=vdat->GetNx(),my=vdat->GetNy(),mz=vdat->GetNz();
2253 bool both=(xdat->GetNN()==vdat->GetNN() && ydat->GetNN()==vdat->GetNN() && zdat->GetNN()==vdat->GetNN());
2254 if(!both && (xdat->GetNx()!=mx || ydat->GetNx()!=my || zdat->GetNx()!=mz)) return; // incompatible dimensions
2255 const mreal acx=1e-6*fabs(x2-x1), acy=1e-6*fabs(y2-y1), acz=1e-6*fabs(z2-z1);
2256 if(both)
2257 #pragma omp parallel for collapse(3)
2258 for(long k=0;k<nz;k++) for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
2259 {
2260 mreal xx = x1+(x2-x1)*i/(nx-1.),dxx,dxy,dxz,vx,dx=0,dd;
2261 mreal yy = y1+(y2-y1)*j/(ny-1.),dyx,dyy,dyz,vy,dy=0;
2262 mreal zz = z1+(z2-z1)*k/(nz-1.),dzx,dzy,dzz,vz,dz=0;
2263 vx = xdat->valueD(dx,dy,dz,&dxx,&dxy,&dxz);
2264 vy = ydat->valueD(dx,dy,dz,&dyx,&dyy,&dyz);
2265 vz = zdat->valueD(dx,dy,dz,&dzx,&dzy,&dzz);
2266 long count=0;
2267 do // use Newton method to find root
2268 {
2269 if(count>50) { dx=NAN; break; } count++;
2270 dd = -dxx*dyy*dzz+dxy*dyx*dzz+dxx*dyz*dzy-dxz*dyx*dzy-dxy*dyz*dzx+dxz*dyy*dzx;
2271 dx += ((dyz*dzy-dyy*dzz)*(xx-vx)+(dxy*dzz-dxz*dzy)*(yy-vy)+(dxz*dyy-dxy*dyz)*(zz-vz))/dd;
2272 dy += ((dyx*dzz-dyz*dzx)*(xx-vx)+(dxz*dzx-dxx*dzz)*(yy-vy)+(dxx*dyz-dxz*dyx)*(zz-vz))/dd;
2273 dz += ((dyy*dzx-dyx*dzy)*(xx-vx)+(dxx*dzy-dxy*dzx)*(yy-vy)+(dxy*dyx-dxx*dyy)*(zz-vz))/dd;
2274 vx = xdat->valueD(dx,dy,dz,&dxx,&dxy,&dxz);
2275 vy = ydat->valueD(dx,dy,dz,&dyx,&dyy,&dyz);
2276 vz = zdat->valueD(dx,dy,dz,&dzx,&dzy,&dzz);
2277 } while(fabs(xx-vx)>acx && fabs(yy-vy)>acy && fabs(zz-vz)>acz); // this valid for linear interpolation
2278 dat->a[i+nx*(j+ny*k)] = mgl_isnan(dx)?NAN:vdat->value(dx,dy,dz);
2279 }
2280 else
2281 {
2282 mglData u(nx), v(ny), w(nz);
2283 mreal dx = (x2-x1)/(nx-1), dy = (y2-y1)/(ny-1), dz = (z2-z1)/(nz-1);
2284 #pragma omp parallel for
2285 for(long i=0;i<nx;i++) u.a[i] = mgl_index_1(x1+dx*i,xdat);
2286 #pragma omp parallel for
2287 for(long i=0;i<ny;i++) v.a[i] = mgl_index_1(y1+dy*i,ydat);
2288 #pragma omp parallel for
2289 for(long i=0;i<nz;i++) w.a[i] = mgl_index_1(z1+dz*i,zdat);
2290 #pragma omp parallel for collapse(3)
2291 for(long k=0;k<nz;k++) for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
2292 dat->a[i+nx*(j+ny*k)] = vdat->value(u.a[i],v.a[j],w.a[k]);
2293 }
2294 }
2295 //-----------------------------------------------------------------------------
mgl_data_evaluate(HCDT dat,HCDT idat,HCDT jdat,HCDT kdat,int norm)2296 HMDT MGL_EXPORT mgl_data_evaluate(HCDT dat, HCDT idat, HCDT jdat, HCDT kdat, int norm)
2297 {
2298 if(!idat || (jdat && jdat->GetNN()!=idat->GetNN()) || (kdat && kdat->GetNN()!=idat->GetNN())) return 0;
2299 const mglData *dd=dynamic_cast<const mglData *>(dat);
2300 long nx=dat->GetNx(), ny=dat->GetNy(), nz=dat->GetNz();
2301 mglData *r=new mglData(idat->GetNx(),idat->GetNy(),idat->GetNz());
2302 mreal dx = nx-1, dy = ny-1, dz = nz-1;
2303 if(!norm) dx=dy=dz=1;
2304 if(dd)
2305 #pragma omp parallel for
2306 for(long i=0;i<idat->GetNN();i++)
2307 {
2308 mreal x=dx*idat->vthr(i), y=jdat?dy*jdat->vthr(i):0, z=kdat?dz*kdat->vthr(i):0;
2309 r->a[i] = mgl_isnum(x*y*z)?mglSpline3st<mreal>(dd->a,nx,ny,nz, x,y,z):NAN;
2310 }
2311 else
2312 #pragma omp parallel for
2313 for(long i=0;i<idat->GetNN();i++)
2314 {
2315 mreal x=dx*idat->vthr(i), y=jdat?dy*jdat->vthr(i):0, z=kdat?dz*kdat->vthr(i):0;
2316 r->a[i] = mgl_isnum(x*y*z)?dat->linear(x,y,z):NAN;;
2317 }
2318 return r;
2319 }
mgl_data_evaluate_(uintptr_t * d,uintptr_t * idat,uintptr_t * jdat,uintptr_t * kdat,int * norm)2320 uintptr_t MGL_EXPORT mgl_data_evaluate_(uintptr_t *d, uintptr_t *idat, uintptr_t *jdat, uintptr_t *kdat, int *norm)
2321 { return uintptr_t(mgl_data_evaluate(_DT_,_DA_(idat),_DA_(jdat),_DA_(kdat),*norm)); }
2322 //-----------------------------------------------------------------------------
mgl_gspline_init(HCDT x,HCDT v)2323 HMDT MGL_EXPORT mgl_gspline_init(HCDT x, HCDT v)
2324 {
2325 long n = v->GetNx();
2326 if(!x || x->GetNx()!=n) return 0;
2327 mglData *res = new mglData(5*(n-1));
2328 mreal *xx=0, *vv=0;
2329 const mglData *dx = dynamic_cast<const mglData *>(x);
2330 if(!dx)
2331 {
2332 xx = new mreal[n];
2333 for(long i=0;i<n;i++) xx[i] = x->v(i);
2334 }
2335 const mglData *dv = dynamic_cast<const mglData *>(v);
2336 if(!dv)
2337 {
2338 vv = new mreal[n];
2339 for(long i=0;i<n;i++) vv[i] = v->v(i);
2340 }
2341 mgl_gspline_init(n,dx?dx->a:xx,dv?dv->a:vv,res->a);
2342 if(xx) delete []xx;
2343 if(vv) delete []vv;
2344 return res;
2345 }
mgl_gspline_init_(uintptr_t * x,uintptr_t * v)2346 uintptr_t MGL_EXPORT mgl_gspline_init_(uintptr_t *x, uintptr_t *v)
2347 { return uintptr_t(mgl_gspline_init(_DA_(x),_DA_(v))); }
2348 //-----------------------------------------------------------------------------
mgl_gspline(HCDT c,mreal dx,mreal * d1,mreal * d2)2349 mreal MGL_EXPORT mgl_gspline(HCDT c, mreal dx, mreal *d1, mreal *d2)
2350 {
2351 long i=0, n = c->GetNx();
2352 if(n%5 || dx<0) return NAN; // not the table of coefficients
2353 while(dx>c->v(5*i))
2354 {
2355 dx-=c->v(5*i); i++;
2356 if(5*i>=n) return NAN;
2357 }
2358 if(d1) *d1 = c->v(5*i+2)+dx*(2*c->v(5*i+3)+3*dx*c->v(5*i+4));
2359 if(d2) *d2 = 2*c->v(5*i+3)+6*dx*c->v(5*i+4);
2360 return c->v(5*i+1)+dx*(c->v(5*i+2)+dx*(c->v(5*i+3)+dx*c->v(5*i+4)));
2361 }
mgl_gspline_(uintptr_t * c,mreal * dx,mreal * d1,mreal * d2)2362 mreal MGL_EXPORT mgl_gspline_(uintptr_t *c, mreal *dx, mreal *d1, mreal *d2)
2363 { return mgl_gspline(_DA_(c),*dx,d1,d2); }
2364 //-----------------------------------------------------------------------------
pntpnt2365 struct pnt { mreal x,y; pnt(mreal X, mreal Y):x(X),y(Y){} };
mgl_find_pnt(std::vector<mreal> & mpos,mreal est,mreal dj)2366 mreal static mgl_find_pnt(std::vector<mreal> &mpos, mreal est, mreal dj)
2367 {
2368 mreal mv = dj, val=NAN;
2369 size_t n = mpos.size(), kk=0;
2370 for(size_t k=0;k<n;k++)
2371 {
2372 mreal de = mpos[k]-est;
2373 if(fabs(de)<fabs(mv)) { mv=de; val=de+est; kk=k; }
2374 }
2375 if(mgl_isnum(val)) mpos.erase(mpos.begin()+kk);
2376 return val;
2377 }
mgl_data_detect(HCDT d,mreal lvl,mreal dj,mreal di,mreal min_len)2378 HMDT MGL_EXPORT mgl_data_detect(HCDT d, mreal lvl, mreal dj, mreal di, mreal min_len)
2379 {
2380 long nx=d->GetNx(), ny=d->GetNy();
2381 std::vector<mreal> *max_pos = new std::vector<mreal>[nx];
2382 if(di<=0) di=dj;
2383 for(long i=0;i<nx;i++) // first collect maximums for each i
2384 {
2385 for(long j=1;j<ny-1;j++)
2386 {
2387 mreal v = d->v(i,j), v1 = d->v(i,j-1), v2 = d->v(i,j+1);
2388 if(v>lvl && v1<v && v>=v2) // NOTE only one edge is required
2389 // if(v>lvl && ((v1<=v && v>v2) || (v1<v && v>=v2))) // NOTE only edges are required
2390 {
2391 bool c1=false, c2=false;
2392 for(long j1=j-1;j1>=0;j1--)
2393 {
2394 mreal vv = d->v(i,j1);
2395 if(vv>v) break;
2396 if(vv<v/2) { c1=true; break; }
2397 }
2398 for(long j2=j+1;j2<ny;j2++)
2399 {
2400 mreal vv = d->v(i,j2);
2401 if(vv>v) break;
2402 if(vv<v/2) { c2=true; break; }
2403 }
2404 if(c1 && c2) max_pos[i].push_back(j + (v2-v1)/(2*v-v2-v1)/2);
2405 }
2406 }
2407 }
2408 std::vector<pnt> curv;
2409 for(long ii=0;ii<nx-1;ii++) // now join points into curves
2410 {
2411 while(max_pos[ii].size()) // try to start curve
2412 {
2413 mreal vv = max_pos[ii].back();
2414 max_pos[ii].pop_back();
2415 pnt p1(ii,vv), p2(ii-1,vv);
2416 size_t ini = curv.size();
2417 curv.push_back(p1);
2418 for(long i=ii+1;i<nx;i++) // join points to selected curve
2419 {
2420 bool none=true;
2421 for(long k=0;k<di && i+k<nx;k++) // try next points
2422 {
2423 mreal val = mgl_find_pnt(max_pos[i+k],p1.y,dj);
2424 if(mgl_isnum(val)) // first try closest point (for noise data)
2425 { p2=p1; p1=pnt(i+k,val); curv.push_back(p1);
2426 none=false; i+=k; break; }
2427 else // next try linear approximation
2428 {
2429 mreal est = p1.y + (k+1)*(p1.y-p2.y)/(p1.x-p2.x);
2430 val = mgl_find_pnt(max_pos[i+k],est,dj);
2431 if(mgl_isnum(val))
2432 { p2=p1; p1=pnt(i+k,val); curv.push_back(p1);
2433 none=false; i+=k; break; }
2434 }
2435 }
2436 if(none) // break curve
2437 { curv.push_back(pnt(NAN,NAN)); break; }
2438 }
2439 if(mgl_isnum(curv.back().x)) curv.push_back(pnt(NAN,NAN));
2440 if(curv[curv.size()-2].x-curv[ini].x<min_len)
2441 curv.erase(curv.begin()+ini,curv.end());
2442 }
2443 }
2444 size_t nn = curv.size();
2445 HMDT res = new mglData(2,nn);
2446 for(size_t k=0;k<nn;k++)
2447 { res->a[2*k] = curv[k].x; res->a[2*k+1] = curv[k].y; }
2448 delete []max_pos; return res;
2449 }
mgl_data_detect_(uintptr_t * d,mreal * lvl,mreal * dj,mreal * di,mreal * min_len)2450 uintptr_t MGL_EXPORT mgl_data_detect_(uintptr_t *d, mreal *lvl, mreal *dj, mreal *di, mreal *min_len)
2451 { return uintptr_t(mgl_data_detect(_DT_,*lvl,*dj, *di, *min_len)); }
2452 //-----------------------------------------------------------------------------
mgl_data_dilate(HMDT d,mreal val,long step)2453 void MGL_EXPORT mgl_data_dilate(HMDT d, mreal val, long step)
2454 {
2455 long nx = d->GetNx(), ny=d->GetNy(), nz=d->GetNz();
2456 if(step<1 || nx<2) return;
2457 long *dist = new long[nx*ny*nz], nn = nx*ny;
2458 bool done=false;
2459 if(nz>1 && ny>1)
2460 {
2461 done = true;
2462 for(long k=0;k<nz;k++) for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
2463 {
2464 long i0 = i+nx*(j+ny*k);
2465 if(d->vthr(i0)>=val) dist[i0]=0;
2466 else
2467 {
2468 dist[i0] = nx+ny;
2469 if(i>0 && dist[i0-1]+1<dist[i0]) dist[i0]=dist[i0-1]+1;
2470 if(j>0 && dist[i0-nx]+1<dist[i0]) dist[i0]=dist[i0-nx]+1;
2471 if(k>0 && dist[i0-nn]+1<dist[i0]) dist[i0]=dist[i0-nn]+1;
2472 }
2473 }
2474 for(long k=nz-1;k>=0;k--) for(long j=ny-1;j>=0;j--) for(long i=nx-1;i>=0;i--)
2475 {
2476 long i0 = i+nx*(j+ny*k);
2477 if(i<nx-1 && dist[i0+1]+1<dist[i0]) dist[i0]=dist[i0+1]+1;
2478 if(j<ny-1 && dist[i0+nx]+1<dist[i0]) dist[i0]=dist[i0+nx]+1;
2479 if(k<nz-1 && dist[i0+nn]+1<dist[i0]) dist[i0]=dist[i0+nn]+1;
2480 }
2481 }
2482 else if(ny>1)
2483 {
2484 done = true;
2485 for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
2486 {
2487 long i0 = i+nx*j;
2488 if(d->vthr(i0)>=val) dist[i0]=0;
2489 else
2490 {
2491 dist[i0] = nx+ny;
2492 if(i>0 && dist[i0-1]+1<dist[i0]) dist[i0]=dist[i0-1]+1;
2493 if(j>0 && dist[i0-nx]+1<dist[i0]) dist[i0]=dist[i0-nx]+1;
2494 }
2495 }
2496 for(long j=ny-1;j>=0;j--) for(long i=nx-1;i>=0;i--)
2497 {
2498 long i0 = i+nx*j;
2499 if(i<nx-1 && dist[i0+1]+1<dist[i0]) dist[i0]=dist[i0+1]+1;
2500 if(j<ny-1 && dist[i0+nx]+1<dist[i0]) dist[i0]=dist[i0+nx]+1;
2501 }
2502 }
2503 else
2504 {
2505 done = true;
2506 for(long i=0;i<nx;i++)
2507 if(d->v(i)>=val) dist[i]=0;
2508 else
2509 {
2510 dist[i] = nx;
2511 if(i>0 && dist[i-1]+1<dist[i]) dist[i]=dist[i-1]+1;
2512 }
2513 for(long i=nx-2;i>=0;i--)
2514 if(dist[i+1]+1<dist[i]) dist[i]=dist[i+1]+1;
2515 }
2516 if(done) for(long i=0;i<nx*ny*nz;i++) d->a[i] = dist[i]<=step?1:0;
2517 delete []dist;
2518 }
mgl_data_dilate_(uintptr_t * d,mreal * val,int * step)2519 void MGL_EXPORT mgl_data_dilate_(uintptr_t *d, mreal *val, int *step)
2520 { mgl_data_dilate(_DT_,*val,*step); }
2521 //-----------------------------------------------------------------------------
mgl_data_erode(HMDT d,mreal val,long step)2522 void MGL_EXPORT mgl_data_erode(HMDT d, mreal val, long step)
2523 {
2524 long nx = d->GetNx(), ny=d->GetNy(), nz=d->GetNz();
2525 if(step<1 || nx<2) return;
2526 long *dist = new long[nx*ny*nz], nn = nx*ny;
2527 bool done=false;
2528 if(nz>1 && ny>1)
2529 {
2530 done = true;
2531 for(long k=0;k<nz;k++) for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
2532 {
2533 long i0 = i+nx*(j+ny*k);
2534 if(d->vthr(i0)<val) dist[i0]=0;
2535 else
2536 {
2537 dist[i0] = nx+ny;
2538 if(i>0 && dist[i0-1]+1<dist[i0]) dist[i0]=dist[i0-1]+1;
2539 if(j>0 && dist[i0-nx]+1<dist[i0]) dist[i0]=dist[i0-nx]+1;
2540 if(k>0 && dist[i0-nn]+1<dist[i0]) dist[i0]=dist[i0-nn]+1;
2541 }
2542 }
2543 for(long k=nz-1;k>=0;k--) for(long j=ny-1;j>=0;j--) for(long i=nx-1;i>=0;i--)
2544 {
2545 long i0 = i+nx*(j+ny*k);
2546 if(i<nx-1 && dist[i0+1]+1<dist[i0]) dist[i0]=dist[i0+1]+1;
2547 if(j<ny-1 && dist[i0+nx]+1<dist[i0]) dist[i0]=dist[i0+nx]+1;
2548 if(k<nz-1 && dist[i0+nn]+1<dist[i0]) dist[i0]=dist[i0+nn]+1;
2549 }
2550 }
2551 else if(ny>1)
2552 {
2553 done = true;
2554 for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
2555 {
2556 long i0 = i+nx*j;
2557 if(d->vthr(i0)<val) dist[i0]=0;
2558 else
2559 {
2560 dist[i0] = nx+ny;
2561 if(i>0 && dist[i0-1]+1<dist[i0]) dist[i0]=dist[i0-1]+1;
2562 if(j>0 && dist[i0-nx]+1<dist[i0]) dist[i0]=dist[i0-nx]+1;
2563 }
2564 }
2565 for(long j=ny-1;j>=0;j--) for(long i=nx-1;i>=0;i--)
2566 {
2567 long i0 = i+nx*j;
2568 if(i<nx-1 && dist[i0+1]+1<dist[i0]) dist[i0]=dist[i0+1]+1;
2569 if(j<ny-1 && dist[i0+nx]+1<dist[i0]) dist[i0]=dist[i0+nx]+1;
2570 }
2571 }
2572 else
2573 {
2574 done = true;
2575 for(long i=0;i<nx;i++)
2576 if(d->v(i)<val) dist[i]=0;
2577 else
2578 {
2579 dist[i] = nx;
2580 if(i>0 && dist[i-1]+1<dist[i]) dist[i]=dist[i-1]+1;
2581 }
2582 for(long i=nx-2;i>=0;i--)
2583 if(dist[i+1]+1<dist[i]) dist[i]=dist[i+1]+1;
2584 }
2585 if(done) for(long i=0;i<nx*ny*nz;i++) d->a[i] = dist[i]>step?1:0;
2586 delete []dist;
2587 }
mgl_data_erode_(uintptr_t * d,mreal * val,int * step)2588 void MGL_EXPORT mgl_data_erode_(uintptr_t *d, mreal *val, int *step)
2589 { mgl_data_erode(_DT_,*val,*step); }
2590 //-----------------------------------------------------------------------------
2591