1 /***************************************************************************
2 * data_new.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 <ctype.h>
21 #include "mgl2/datac.h"
22 #include "mgl2/evalc.h"
23 #include "mgl2/thread.h"
24 #include "interp.hpp"
25 void MGL_NO_EXPORT mgl_txt_funcC(const mreal *x, mreal *dx, void *par);
26 //-----------------------------------------------------------------------------
mgl_datac_trace(HCDT d)27 HADT MGL_EXPORT mgl_datac_trace(HCDT d)
28 {
29 long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
30 const mglDataC *dc = dynamic_cast<const mglDataC *>(d);
31 mglDataC *r=new mglDataC(nx);
32 if(dc)
33 {
34 if(ny>=nx && nz>=nx)
35 #pragma omp parallel for
36 for(long i=0;i<nx;i++) r->a[i] = dc->a[i+nx*(i+ny*i)];
37 else if(ny>=nx)
38 #pragma omp parallel for
39 for(long i=0;i<nx;i++) r->a[i] = dc->a[i+nx*i];
40 else
41 #pragma omp parallel for
42 for(long i=0;i<nx;i++) r->a[i] = dc->a[i];
43 }
44 else if(ny>=nx && nz>=nx)
45 #pragma omp parallel for
46 for(long i=0;i<nx;i++) r->a[i] = d->v(i,i,i);
47 else if(ny>=nx)
48 #pragma omp parallel for
49 for(long i=0;i<nx;i++) r->a[i] = d->v(i,i);
50 else
51 #pragma omp parallel for
52 for(long i=0;i<nx;i++) r->a[i] = d->v(i);
53 return r;
54 }
mgl_datac_trace_(uintptr_t * d)55 uintptr_t MGL_EXPORT mgl_datac_trace_(uintptr_t *d)
56 { return uintptr_t(mgl_datac_trace(_DC_)); }
57 //-----------------------------------------------------------------------------
mgl_datac_subdata_ext(HCDT d,HCDT xx,HCDT yy,HCDT zz)58 HADT MGL_EXPORT mgl_datac_subdata_ext(HCDT d, HCDT xx, HCDT yy, HCDT zz)
59 {
60 if(!xx || !yy || !zz)
61 {
62 mglData tmp; tmp.a[0]=-1;
63 return mgl_datac_subdata_ext(d,xx?xx:&tmp,yy?yy:&tmp,zz?zz:&tmp);
64 }
65
66 long n=0,m=0,l=0,j,k;
67 bool ix=false, iy=false, iz=false;
68 if(xx->GetNz()>1) // 3d data
69 {
70 n = xx->GetNx(); m = xx->GetNy(); l = xx->GetNz();
71 j = yy->GetNN(); if(j>1 && j!=n*m*l) return 0; // wrong sizes
72 k = zz->GetNN(); if(k>1 && k!=n*m*l) return 0; // wrong sizes
73 ix = true; iy = j>1; iz = k>1;
74 }
75 else if(yy->GetNz()>1)
76 {
77 n = yy->GetNx(); m = yy->GetNy(); l = yy->GetNz();
78 j = xx->GetNN(); if(j>1 && j!=n*m*l) return 0; // wrong sizes
79 k = zz->GetNN(); if(k>1 && k!=n*m*l) return 0; // wrong sizes
80 iy = true; ix = j>1; iz = k>1;
81 }
82 else if(zz->GetNz()>1)
83 {
84 n = zz->GetNx(); m = zz->GetNy(); l = zz->GetNz();
85 j = yy->GetNN(); if(j>1 && j!=n*m*l) return 0; // wrong sizes
86 k = xx->GetNN(); if(k>1 && k!=n*m*l) return 0; // wrong sizes
87 iz = true; iy = j>1; ix = k>1;
88 }
89 else if(xx->GetNy()>1) // 2d data
90 {
91 n = xx->GetNx(); m = xx->GetNy(); l = 1;
92 j = yy->GetNx()*yy->GetNy(); if(j>1 && j!=n*m) return 0; // wrong sizes
93 k = zz->GetNx()*zz->GetNy(); if(k>1 && k!=n*m) return 0; // wrong sizes
94 ix = true; iy = j>1; iz = k>1;
95 }
96 else if(yy->GetNy()>1)
97 {
98 n = yy->GetNx(); m = yy->GetNy(); l = 1;
99 j = xx->GetNx()*xx->GetNy(); if(j>1 && j!=n*m) return 0; // wrong sizes
100 k = zz->GetNx()*zz->GetNy(); if(k>1 && k!=n*m) return 0; // wrong sizes
101 iy = true; ix = j>1; iz = k>1;
102 }
103 else if(zz->GetNy()>1)
104 {
105 n = zz->GetNx(); m = zz->GetNy(); l = 1;
106 j = yy->GetNx()*yy->GetNy(); if(j>1 && j!=n*m) return 0; // wrong sizes
107 k = xx->GetNx()*xx->GetNy(); if(k>1 && k!=n*m) return 0; // wrong sizes
108 iz = true; iy = j>1; ix = k>1;
109 }
110 long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
111 long vx=long(xx->v(0)), vy=long(yy->v(0)), vz=long(zz->v(0));
112 const mglDataC *dd = dynamic_cast<const mglDataC *>(d);
113 mglDataC *r;
114 if(n*m*l>1) // this is 2d or 3d data
115 {
116 mglDataV tx(n,m,l),ty(n,m,l),tz(n,m,l);
117 if(!ix) { xx = &tx; if(vx>=0) tx.Fill(vx); else tx.All(); }
118 if(!iy) { yy = &ty; if(vy>=0) ty.Fill(vy); else ty.All(); }
119 if(!iz) { zz = &tz; if(vz>=0) tz.Fill(vz); else tz.All(); }
120 r=new mglDataC(n,m,l);
121 if(dd)
122 #pragma omp parallel for
123 for(long i0=0;i0<n*m*l;i0++)
124 {
125 long x=long(floor(0.5+xx->vthr(i0))), y=long(floor(0.5+yy->vthr(i0))), z=long(floor(0.5+zz->vthr(i0)));
126 r->a[i0] = (x>=0 && x<nx && y>=0 && y<ny && z>=0 && z<nz)?dd->a[x+nx*(y+ny*z)]:NAN;
127 }
128 else
129 #pragma omp parallel for
130 for(long i0=0;i0<n*m*l;i0++)
131 {
132 long x=long(floor(0.5+xx->vthr(i0))), y=long(floor(0.5+yy->vthr(i0))), z=long(floor(0.5+zz->vthr(i0)));
133 r->a[i0] = (x>=0 && x<nx && y>=0 && y<ny && z>=0 && z<nz)?d->v(x,y,z):NAN;
134 }
135 }
136 else // this is 1d data -> try as normal SubData()
137 {
138 mglDataV tx(nx),ty(ny),tz(nz); tx.Fill(0,nx-1); ty.Fill(0,ny-1); tz.Fill(0,nz-1);
139 if(xx->GetNx()>1 || vx>=0) n=xx->GetNx(); else { n=nx; xx = &tx; }
140 if(yy->GetNx()>1 || vy>=0) m=yy->GetNx(); else { m=ny; yy = &ty; }
141 if(zz->GetNx()>1 || vz>=0) l=zz->GetNx(); else { l=nz; zz = &tz; }
142 r=new mglDataC(n,m,l);
143 if(dd)
144 #pragma omp parallel for collapse(3)
145 for(long k=0;k<l;k++) for(long j=0;j<m;j++) for(long i=0;i<n;i++)
146 {
147 long x=long(floor(0.5+xx->v(i))), y=long(floor(0.5+yy->v(j))), z=long(floor(0.5+zz->v(k)));
148 r->a[i+n*(j+m*k)] = (x>=0 && x<nx && y>=0 && y<ny && z>=0 && z<nz)?dd->a[x+nx*(y+ny*z)]:NAN;
149 }
150 else
151 #pragma omp parallel for collapse(3)
152 for(long k=0;k<l;k++) for(long j=0;j<m;j++) for(long i=0;i<n;i++)
153 {
154 long x=long(floor(0.5+xx->v(i))), y=long(floor(0.5+yy->v(j))), z=long(floor(0.5+zz->v(k)));
155 r->a[i+n*(j+m*k)] = (x>=0 && x<nx && y>=0 && y<ny && z>=0 && z<nz)?d->v(x,y,z):NAN;
156 }
157 if(m==1) { r->ny=r->nz; r->nz=1; }// "squeeze" dimensions
158 if(n==1) { r->nx=r->ny; r->ny=r->nz; r->nz=1; r->NewId();}
159 }
160 return r;
161 }
162 //-----------------------------------------------------------------------------
mgl_datac_subdata(HCDT d,long xx,long yy,long zz)163 HADT MGL_EXPORT mgl_datac_subdata(HCDT d, long xx,long yy,long zz)
164 {
165 long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz(), n=1,m=1,l=1;
166 int dx=0,dy=0,dz=0;
167 if(xx<0) { xx=0; dx=1; n=nx; }
168 if(yy<0) { yy=0; dy=1; m=ny; }
169 if(zz<0) { zz=0; dz=1; l=nz; }
170 const mglDataC *dd = dynamic_cast<const mglDataC *>(d);
171 mglDataC *r=new mglDataC(n,m,l);
172 if(xx>=nx || yy>=ny || zz>=nz)
173 #pragma omp parallel for
174 for(long i=0;i<n*m*l;i++) r->a[i] = NAN;
175 else if(dd)
176 #pragma omp parallel for collapse(3)
177 for(long k=0;k<l;k++) for(long j=0;j<m;j++) for(long i=0;i<n;i++)
178 r->a[i+n*(j+m*k)] = dd->a[xx+dx*i + nx*(yy+dy*j + ny*(zz+dz*k))];
179 else
180 #pragma omp parallel for collapse(3)
181 for(long k=0;k<l;k++) for(long j=0;j<m;j++) for(long i=0;i<n;i++)
182 r->a[i+n*(j+m*k)] = d->v(xx+dx*i, yy+dy*j, zz+dz*k);
183 if(m==1) { r->ny=r->nz; r->nz=1; }// "squeeze" dimensions
184 if(n==1) { r->nx=r->ny; r->ny=r->nz; r->nz=1; r->NewId();}
185 return r;
186 }
187 //-----------------------------------------------------------------------------
mgl_datac_subdata_(uintptr_t * d,int * xx,int * yy,int * zz)188 uintptr_t MGL_EXPORT mgl_datac_subdata_(uintptr_t *d, int *xx,int *yy,int *zz)
189 { return uintptr_t(mgl_datac_subdata(_DC_,*xx,*yy,*zz)); }
mgl_datac_subdata_ext_(uintptr_t * d,uintptr_t * xx,uintptr_t * yy,uintptr_t * zz)190 uintptr_t MGL_EXPORT mgl_datac_subdata_ext_(uintptr_t *d, uintptr_t *xx, uintptr_t *yy, uintptr_t *zz)
191 { return uintptr_t(mgl_datac_subdata_ext(_DC_,_DA_(xx),_DA_(yy),_DA_(zz))); }
192 //-----------------------------------------------------------------------------
mgl_cresize(void * par)193 static void *mgl_cresize(void *par)
194 {
195 mglThreadC *t=(mglThreadC *)par;
196 long nx=t->p[0]+0.1, ny=t->p[1]+0.1;
197 long n1=t->p[3]+0.1,n2=t->p[4]+0.1,n3=t->p[5]+0.1;
198 dual *b=t->a;
199 const dual *a=t->b;
200 const mreal *c=(const mreal *)t->v;
201 #if !MGL_HAVE_PTHREAD
202 #pragma omp parallel for
203 #endif
204 for(long i0=t->id;i0<t->n;i0+=mglNumThr)
205 {
206 mreal i=(i0%nx), j=((i0/nx)%ny), k=(i0/(nx*ny));
207 b[i0] = mglSpline3Cs(a,n1,n2,n3, c[0]+i*c[1], c[2]+j*c[3], c[4]+k*c[5]);
208 }
209 return 0;
210 }
mgl_datac_resize_box(HCDT dat,long mx,long my,long mz,mreal x1,mreal x2,mreal y1,mreal y2,mreal z1,mreal z2)211 HADT MGL_EXPORT mgl_datac_resize_box(HCDT dat, long mx,long my,long mz, mreal x1,mreal x2, mreal y1,mreal y2, mreal z1,mreal z2)
212 { // NOTE: only for mglDataC
213 const mglDataC *d=dynamic_cast<const mglDataC *>(dat);
214 if(!d) return 0;
215 long nx = d->nx-1, ny = d->ny-1, nz = d->nz-1;
216 mx = mx<1 ? nx+1:mx; my = my<1 ? ny+1:my; mz = mz<1 ? nz+1:mz;
217 mglDataC *r=new mglDataC(mx,my,mz);
218
219 mreal par[6]={nx*x1,0,ny*y1,0,nz*z1,0};
220 long nn[6]={mx,my,mz,nx+1,ny+1,nz+1};
221 if(mx>1) par[1] = nx*(x2-x1)/(mx-1);
222 if(my>1) par[3] = ny*(y2-y1)/(my-1);
223 if(mz>1) par[5] = nz*(z2-z1)/(mz-1);
224 mglStartThreadC(mgl_cresize,0,mx*my*mz,r->a,d->a,0,nn,par);
225 return r;
226 }
mgl_datac_resize(HCDT d,long mx,long my,long mz)227 HADT MGL_EXPORT mgl_datac_resize(HCDT d, long mx,long my,long mz)
228 { return mgl_datac_resize_box(d, mx,my,mz,0,1,0,1,0,1); }
mgl_datac_resize_(uintptr_t * d,int * mx,int * my,int * mz)229 uintptr_t MGL_EXPORT mgl_datac_resize_(uintptr_t *d, int *mx,int *my,int *mz)
230 { return uintptr_t(mgl_datac_resize(_DC_,*mx,*my,*mz)); }
mgl_datac_resize_box_(uintptr_t * d,int * mx,int * my,int * mz,mreal * x1,mreal * x2,mreal * y1,mreal * y2,mreal * z1,mreal * z2)231 uintptr_t MGL_EXPORT mgl_datac_resize_box_(uintptr_t *d, int *mx,int *my,int *mz, mreal *x1,mreal *x2, mreal *y1,mreal *y2, mreal *z1,mreal *z2)
232 { return uintptr_t(mgl_datac_resize_box(_DC_,*mx,*my,*mz,*x1,*x2,*y1,*y2,*z1,*z2)); }
233 //-----------------------------------------------------------------------------
mgl_datac_combine(HCDT d1,HCDT d2)234 HADT MGL_EXPORT mgl_datac_combine(HCDT d1, HCDT d2)
235 {
236 long n1=d1->GetNy(),n2=d2->GetNx(),nx=d1->GetNx();
237 if(d1->GetNz()>1 || (n1>1 && d2->GetNy()>1) || d2->GetNz()>1) return 0; // wrong dimensions
238 mglDataC *r=new mglDataC;
239 bool dim2=true;
240 if(n1==1) { n1=n2; n2=d2->GetNy(); dim2 = false; }
241 r->Create(nx,n1,n2);
242 if(dim2) n1*=nx; else { n2*=n1; n1=nx; }
243
244 const mglDataC *c1=dynamic_cast<const mglDataC *>(d1);
245 const mglDataC *c2=dynamic_cast<const mglDataC *>(d2);
246 if(c1 && c2)
247 #pragma omp parallel for collapse(2)
248 for(long j=0;j<n2;j++) for(long i=0;i<n1;i++)
249 r->a[i+n1*j] = c1->a[i]*c2->a[j];
250 else if(c1)
251 #pragma omp parallel for collapse(2)
252 for(long j=0;j<n2;j++) for(long i=0;i<n1;i++)
253 r->a[i+n1*j] = c1->a[i]*d2->vthr(j);
254 else if(c2)
255 #pragma omp parallel for collapse(2)
256 for(long j=0;j<n2;j++) for(long i=0;i<n1;i++)
257 r->a[i+n1*j] = d1->vthr(i)*c2->a[j];
258 else
259 #pragma omp parallel for collapse(2)
260 for(long j=0;j<n2;j++) for(long i=0;i<n1;i++)
261 r->a[i+n1*j] = d1->vthr(i)*d2->vthr(j);
262 return r;
263 }
mgl_datac_combine_(uintptr_t * a,uintptr_t * b)264 uintptr_t MGL_EXPORT mgl_datac_combine_(uintptr_t *a, uintptr_t *b)
265 { return uintptr_t(mgl_datac_combine(_DA_(a),_DA_(b))); }
266 //-----------------------------------------------------------------------------
mgl_sumc_z(void * par)267 static void *mgl_sumc_z(void *par)
268 {
269 mglThreadC *t=(mglThreadC *)par;
270 long nz=t->p[2], nn=t->n;
271 dual *b=t->a;
272 const dual *a=t->b;
273 #if !MGL_HAVE_PTHREAD
274 #pragma omp parallel for
275 #endif
276 for(long i=t->id;i<nn;i+=mglNumThr)
277 {
278 b[i]=0;
279 for(long j=0;j<nz;j++) b[i] += a[i+nn*j];
280 b[i] /= nz;
281 }
282 return 0;
283 }
mgl_sumc_y(void * par)284 static void *mgl_sumc_y(void *par)
285 {
286 mglThreadC *t=(mglThreadC *)par;
287 long nx=t->p[0], ny=t->p[1], nn=t->n;
288 dual *b=t->a;
289 const dual *a=t->b;
290 #if !MGL_HAVE_PTHREAD
291 #pragma omp parallel for
292 #endif
293 for(long i=t->id;i<nn;i+=mglNumThr)
294 {
295 long k = (i%nx)+nx*ny*(i/nx); b[i]=0;
296 for(long j=0;j<ny;j++) b[i] += a[k+nx*j];
297 b[i] /= ny;
298 }
299 return 0;
300 }
mgl_sumc_x(void * par)301 static void *mgl_sumc_x(void *par)
302 {
303 mglThreadC *t=(mglThreadC *)par;
304 long nx=t->p[0], nn=t->n;
305 dual *b=t->a;
306 const dual *a=t->b;
307 #if !MGL_HAVE_PTHREAD
308 #pragma omp parallel for
309 #endif
310 for(long i=t->id;i<nn;i+=mglNumThr)
311 {
312 long k = i*nx; b[i]=0;
313 for(long j=0;j<nx;j++) b[i] += a[j+k];
314 b[i] /= nx;
315 }
316 return 0;
317 }
mgl_datac_sum(HCDT dat,const char * dir)318 HADT MGL_EXPORT mgl_datac_sum(HCDT dat, const char *dir)
319 {
320 if(!dir || *dir==0) return 0;
321 long nx=dat->GetNx(),ny=dat->GetNy(),nz=dat->GetNz();
322 long p[3]={nx,ny,nz};
323 dual *b = new dual[nx*ny*nz];
324 dual *c = new dual[nx*ny*nz];
325
326 const mglDataC *d=dynamic_cast<const mglDataC *>(dat);
327 if(d) memcpy(c,d->a,nx*ny*nz*sizeof(dual));
328 else
329 #pragma omp parallel for
330 for(long i=0;i<nx*ny*nz;i++) c[i]=dat->vthr(i);
331
332 if(strchr(dir,'z') && nz>1)
333 {
334 mglStartThreadC(mgl_sumc_z,0,nx*ny,b,c,0,p);
335 memcpy(c,b,nx*ny*sizeof(dual)); p[2] = 1;
336 }
337 if(strchr(dir,'y') && ny>1)
338 {
339 mglStartThreadC(mgl_sumc_y,0,nx*p[2],b,c,0,p);
340 memcpy(c,b,nx*p[2]*sizeof(dual)); p[1] = p[2]; p[2] = 1;
341 }
342 if(strchr(dir,'x') && nx>1)
343 {
344 mglStartThreadC(mgl_sumc_x,0,p[1]*p[2],b,c,0,p);
345 p[0] = p[1]; p[1] = p[2]; p[2] = 1;
346 memcpy(c,b,p[0]*p[1]*sizeof(dual));
347 }
348 mglDataC *r=new mglDataC(p[0],p[1],p[2]);
349 memcpy(r->a,c,p[0]*p[1]*p[2]*sizeof(dual));
350 delete []b; delete []c; return r;
351 }
mgl_datac_sum_(uintptr_t * d,const char * dir,int l)352 uintptr_t MGL_EXPORT mgl_datac_sum_(uintptr_t *d, const char *dir,int l)
353 { char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
354 uintptr_t r=uintptr_t(mgl_datac_sum(_DC_,s)); delete []s; return r; }
355 //-----------------------------------------------------------------------------
mgl_datac_momentum(HCDT dat,char dir,const char * how)356 HADT MGL_EXPORT mgl_datac_momentum(HCDT dat, char dir, const char *how)
357 {
358 if(!how || !(*how) || !strchr("xyz",dir)) return 0;
359 long nx=dat->GetNx(),ny=dat->GetNy(),nz=dat->GetNz();
360 mglDataV x(nx,ny,nz, 0,1,'x'); x.Name(L"x");
361 mglDataV y(nx,ny,nz, 0,1,'y'); y.Name(L"y");
362 mglDataV z(nx,ny,nz, 0,1,'z'); z.Name(L"z");
363 mglDataC u(dat); u.Name(L"u"); // NOTE slow !!!
364 std::vector<mglDataA*> list;
365 list.push_back(&x); list.push_back(&y); list.push_back(&z); list.push_back(&u);
366 HADT res=mglFormulaCalcC(how,list), b=0;
367
368 if(dir=='x')
369 {
370 b=new mglDataC(nx);
371 #pragma omp parallel for
372 for(long i=0;i<nx;i++)
373 {
374 dual i1=0,i0=0;
375 for(long j=0;j<ny*nz;j++)
376 {
377 dual u=dat->vthr(i+nx*j);
378 i0 += u; i1 += u*res->a[i+nx*j];
379 }
380 b->a[i] = i0!=mreal(0) ? i1/i0 : 0;
381 }
382 }
383 if(dir=='y')
384 {
385 b=new mglDataC(ny);
386 #pragma omp parallel for
387 for(long i=0;i<ny;i++)
388 {
389 dual i1=0,i0=0;
390 for(long k=0;k<nz;k++) for(long j=0;j<nx;j++)
391 {
392 dual u=dat->v(j,i,k);
393 i0 += u; i1 += u*res->a[j+nx*(i+ny*k)];
394 }
395 b->a[i] = i0!=mreal(0) ? i1/i0 : 0;
396 }
397 }
398 if(dir=='z')
399 {
400 long nn=nx*ny;
401 b=new mglDataC(nz);
402 #pragma omp parallel for
403 for(long i=0;i<nz;i++)
404 {
405 dual i1=0,i0=0;
406 for(long j=0;j<nn;j++)
407 {
408 dual u=dat->vthr(j+nn*i);
409 i0 += u; i1 += u*res->a[j+nn*i];
410 }
411 b->a[i] = i0!=mreal(0) ? i1/i0 : 0;
412 }
413 }
414 mgl_delete_datac(res); return b;
415 }
mgl_datac_momentum_(uintptr_t * d,char * dir,const char * how,int,int l)416 uintptr_t MGL_EXPORT mgl_datac_momentum_(uintptr_t *d, char *dir, const char *how, int,int l)
417 { char *s=new char[l+1]; memcpy(s,how,l); s[l]=0;
418 uintptr_t r=uintptr_t(mgl_datac_momentum(_DC_,*dir, s)); delete []s; return r; }
419 //-----------------------------------------------------------------------------
mgl_datac_evaluate(HCDT dat,HCDT idat,HCDT jdat,HCDT kdat,int norm)420 HADT MGL_EXPORT mgl_datac_evaluate(HCDT dat, HCDT idat, HCDT jdat, HCDT kdat, int norm)
421 {
422 if(!idat || (jdat && jdat->GetNN()!=idat->GetNN()) || (kdat && kdat->GetNN()!=idat->GetNN())) return 0;
423 const mglData *dd=dynamic_cast<const mglData *>(dat);
424 const mglDataC *dc=dynamic_cast<const mglDataC *>(dat);
425 long nx=dat->GetNx(), ny=dat->GetNy(), nz=dat->GetNz();
426 mglDataC *r=new mglDataC(idat->GetNx(),idat->GetNy(),idat->GetNz());
427 mreal dx = nx-1, dy = ny-1, dz = nz-1;
428 if(!norm) dx=dy=dz=1;
429 if(dd)
430 #pragma omp parallel for
431 for(long i=0;i<idat->GetNN();i++)
432 {
433 mreal x=dx*idat->vthr(i), y=jdat?dy*jdat->vthr(i):0, z=kdat?dz*kdat->vthr(i):0;
434 r->a[i] = mgl_isnum(x*y*z)?mglSpline3st<mreal>(dd->a,nx,ny,nz, x,y,z):NAN;
435 }
436 else if(dc)
437 #pragma omp parallel for
438 for(long i=0;i<idat->GetNN();i++)
439 {
440 mreal x=dx*idat->vthr(i), y=jdat?dy*jdat->vthr(i):0, z=kdat?dz*kdat->vthr(i):0;
441 r->a[i] = mgl_isnum(x*y*z)?mglSpline3st<dual>(dc->a,nx,ny,nz, x,y,z):NAN;
442 }
443 else
444 #pragma omp parallel for
445 for(long i=0;i<idat->GetNN();i++)
446 {
447 mreal x=dx*idat->vthr(i), y=jdat?dy*jdat->vthr(i):0, z=kdat?dz*kdat->vthr(i):0;
448 r->a[i] = mgl_isnum(x*y*z)?dat->linear(x,y,z):NAN;;
449 }
450 return r;
451 }
mgl_datac_evaluate_(uintptr_t * d,uintptr_t * idat,uintptr_t * jdat,uintptr_t * kdat,int * norm)452 uintptr_t MGL_EXPORT mgl_datac_evaluate_(uintptr_t *d, uintptr_t *idat, uintptr_t *jdat, uintptr_t *kdat, int *norm)
453 { return uintptr_t(mgl_datac_evaluate(_DC_,_DA_(idat),_DA_(jdat),_DA_(kdat),*norm)); }
454 //-----------------------------------------------------------------------------
mgl_datac_column(HCDT dat,const char * eq)455 HADT MGL_EXPORT mgl_datac_column(HCDT dat, const char *eq)
456 {
457 std::vector<mglDataA*> list;
458 const char *id = dat->GetColumnId();
459 size_t len = strlen(id);
460 for(size_t i=0;i<len;i++)
461 {
462 mglDataT *col = new mglDataT(*dat);
463 col->SetInd(i,id[i]);
464 list.push_back(col);
465 }
466 if(list.size()==0) return 0; // no named columns
467 mglDataV *t = new mglDataV(dat->GetNy(),dat->GetNz());
468 t->Name(L"#$mgl"); list.push_back(t);
469 HADT r = mglFormulaCalcC(eq,list);
470 for(size_t i=0;i<list.size();i++) delete list[i];
471 return r;
472 }
mgl_datac_column_(uintptr_t * d,const char * eq,int l)473 uintptr_t MGL_EXPORT mgl_datac_column_(uintptr_t *d, const char *eq,int l)
474 { char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
475 uintptr_t r = uintptr_t(mgl_datac_column(_DC_,s));
476 delete []s; return r; }
477 //-----------------------------------------------------------------------------
mgl_datac_mul_dat(HADT d,HCDT a)478 void MGL_EXPORT mgl_datac_mul_dat(HADT d, HCDT a)
479 {
480 long nx=d->nx, ny=d->ny, nz=d->nz;
481 long mx=a->GetNx(), my=a->GetNy(), mz=a->GetNz();
482 const mglDataC *c = dynamic_cast<const mglDataC*>(a);
483
484 if(mz*my*mx==1)
485 {
486 dual v=c?c->a[0]:a->v(0);
487 #pragma omp parallel for
488 for(long i=0;i<nx*ny*nz;i++) d->a[i] += v;
489 }
490 else
491 {
492 long n=0, m=0;
493 if(nz*ny*nx==mz*my*mx) { n=nx*ny*nz; m=1; }
494 else if(ny*nx==my*mx) { n=nx*ny; m=nz; }
495 else if(nx==mx) { n=nx; m=ny*nz; }
496 if(c)
497 #pragma omp parallel for collapse(2)
498 for(long k=0;k<m;k++) for(long i=0;i<n;i++) d->a[i+n*k] *= c->a[i];
499 else
500 #pragma omp parallel for collapse(2)
501 for(long k=0;k<m;k++) for(long i=0;i<n;i++) d->a[i+n*k] *= a->vthr(i);
502 }
503 }
504 //-----------------------------------------------------------------------------
mgl_datac_mul_num(HADT d,mdual a)505 void MGL_EXPORT mgl_datac_mul_num(HADT d, mdual a)
506 {
507 long n=d->GetNN();
508 #pragma omp parallel for
509 for(long i=0;i<n;i++) d->a[i] *= dual(a);
510 }
511 //-----------------------------------------------------------------------------
mgl_datac_div_dat(HADT d,HCDT a)512 void MGL_EXPORT mgl_datac_div_dat(HADT d, HCDT a)
513 {
514 long nx=d->nx, ny=d->ny, nz=d->nz;
515 long mx=a->GetNx(), my=a->GetNy(), mz=a->GetNz();
516 const mglDataC *c = dynamic_cast<const mglDataC*>(a);
517
518 if(mz*my*mx==1)
519 {
520 dual v=c?c->a[0]:a->v(0);
521 #pragma omp parallel for
522 for(long i=0;i<nx*ny*nz;i++) d->a[i] /= v;
523 }
524 else
525 {
526 long n=0, m=0;
527 if(nz*ny*nx==mz*my*mx) { n=nx*ny*nz; m=1; }
528 else if(ny*nx==my*mx) { n=nx*ny; m=nz; }
529 else if(nx==mx) { n=nx; m=ny*nz; }
530 if(c)
531 #pragma omp parallel for collapse(2)
532 for(long k=0;k<m;k++) for(long i=0;i<n;i++) d->a[i+n*k] /= c->a[i];
533 else
534 #pragma omp parallel for collapse(2)
535 for(long k=0;k<m;k++) for(long i=0;i<n;i++) d->a[i+n*k] /= a->vthr(i);
536 }
537 }
538 //-----------------------------------------------------------------------------
mgl_datac_div_num(HADT d,mdual a)539 void MGL_EXPORT mgl_datac_div_num(HADT d, mdual a)
540 {
541 long n=d->GetNN();
542 #pragma omp parallel for
543 for(long i=0;i<n;i++) d->a[i] /= dual(a);
544 }
545 //-----------------------------------------------------------------------------
mgl_datac_add_dat(HADT d,HCDT a)546 void MGL_EXPORT mgl_datac_add_dat(HADT d, HCDT a)
547 {
548 long nx=d->nx, ny=d->ny, nz=d->nz;
549 long mx=a->GetNx(), my=a->GetNy(), mz=a->GetNz();
550 const mglDataC *c = dynamic_cast<const mglDataC*>(a);
551
552 if(mz*my*mx==1)
553 {
554 dual v=c?c->a[0]:a->v(0);
555 #pragma omp parallel for
556 for(long i=0;i<nx*ny*nz;i++) d->a[i] += v;
557 }
558 else
559 {
560 long n=0, m=0;
561 if(nz*ny*nx==mz*my*mx) { n=nx*ny*nz; m=1; }
562 else if(ny*nx==my*mx) { n=nx*ny; m=nz; }
563 else if(nx==mx) { n=nx; m=ny*nz; }
564 if(c)
565 #pragma omp parallel for collapse(2)
566 for(long k=0;k<m;k++) for(long i=0;i<n;i++) d->a[i+n*k] += c->a[i];
567 else
568 #pragma omp parallel for collapse(2)
569 for(long k=0;k<m;k++) for(long i=0;i<n;i++) d->a[i+n*k] += a->vthr(i);
570 }
571 }
572 //-----------------------------------------------------------------------------
mgl_datac_add_num(HADT d,mdual a)573 void MGL_EXPORT mgl_datac_add_num(HADT d, mdual a)
574 {
575 long n=d->GetNN();
576 #pragma omp parallel for
577 for(long i=0;i<n;i++) d->a[i] += dual(a);
578 }
579 //-----------------------------------------------------------------------------
mgl_datac_sub_dat(HADT d,HCDT a)580 void MGL_EXPORT mgl_datac_sub_dat(HADT d, HCDT a)
581 {
582 long nx=d->nx, ny=d->ny, nz=d->nz;
583 long mx=a->GetNx(), my=a->GetNy(), mz=a->GetNz();
584 const mglDataC *c = dynamic_cast<const mglDataC*>(a);
585
586 if(mz*my*mx==1)
587 {
588 dual v=c?c->a[0]:a->v(0);
589 #pragma omp parallel for
590 for(long i=0;i<nx*ny*nz;i++) d->a[i] -= v;
591 }
592 else
593 {
594 long n=0, m=0;
595 if(nz*ny*nx==mz*my*mx) { n=nx*ny*nz; m=1; }
596 else if(ny*nx==my*mx) { n=nx*ny; m=nz; }
597 else if(nx==mx) { n=nx; m=ny*nz; }
598 if(c)
599 #pragma omp parallel for collapse(2)
600 for(long k=0;k<m;k++) for(long i=0;i<n;i++) d->a[i+n*k] -= c->a[i];
601 else
602 #pragma omp parallel for collapse(2)
603 for(long k=0;k<m;k++) for(long i=0;i<n;i++) d->a[i+n*k] -= a->vthr(i);
604 }
605 }
606 //-----------------------------------------------------------------------------
mgl_datac_sub_num(HADT d,mdual a)607 void MGL_EXPORT mgl_datac_sub_num(HADT d, mdual a)
608 {
609 long n=d->GetNN();
610 #pragma omp parallel for
611 for(long i=0;i<n;i++) d->a[i] -= dual(a);
612 }
613 //-----------------------------------------------------------------------------
mgl_datac_mul_dat_(uintptr_t * d,uintptr_t * b)614 void MGL_EXPORT mgl_datac_mul_dat_(uintptr_t *d, uintptr_t *b) { mgl_datac_mul_dat(_DC_, _DA_(b)); }
mgl_datac_div_dat_(uintptr_t * d,uintptr_t * b)615 void MGL_EXPORT mgl_datac_div_dat_(uintptr_t *d, uintptr_t *b) { mgl_datac_div_dat(_DC_, _DA_(b)); }
mgl_datac_add_dat_(uintptr_t * d,uintptr_t * b)616 void MGL_EXPORT mgl_datac_add_dat_(uintptr_t *d, uintptr_t *b) { mgl_datac_add_dat(_DC_, _DA_(b)); }
mgl_datac_sub_dat_(uintptr_t * d,uintptr_t * b)617 void MGL_EXPORT mgl_datac_sub_dat_(uintptr_t *d, uintptr_t *b) { mgl_datac_sub_dat(_DC_, _DA_(b)); }
mgl_datac_mul_num_(uintptr_t * d,mdual * b)618 void MGL_EXPORT mgl_datac_mul_num_(uintptr_t *d, mdual *b) { mgl_datac_mul_num(_DC_, *b); }
mgl_datac_div_num_(uintptr_t * d,mdual * b)619 void MGL_EXPORT mgl_datac_div_num_(uintptr_t *d, mdual *b) { mgl_datac_div_num(_DC_, *b); }
mgl_datac_add_num_(uintptr_t * d,mdual * b)620 void MGL_EXPORT mgl_datac_add_num_(uintptr_t *d, mdual *b) { mgl_datac_add_num(_DC_, *b); }
mgl_datac_sub_num_(uintptr_t * d,mdual * b)621 void MGL_EXPORT mgl_datac_sub_num_(uintptr_t *d, mdual *b) { mgl_datac_sub_num(_DC_, *b); }
622 //-----------------------------------------------------------------------------
mgl_datac_section(HCDT dat,HCDT ids,char dir,mreal val)623 HADT MGL_EXPORT mgl_datac_section(HCDT dat, HCDT ids, char dir, mreal val)
624 {
625 long di = 1, n = dat->GetNx();
626 if(dir=='y') { di = dat->GetNx(); n = dat->GetNy(); }
627 if(dir=='z') { di = dat->GetNx()*dat->GetNy(); n = dat->GetNz(); }
628 // first collect position of key values
629 std::vector<long> pos; pos.push_back(0);
630 if(mgl_isnan(val)) for(long i=1;i<n;i++)
631 {
632 if(mgl_isnan(dat->vthr(i*di))) pos.push_back(i);
633 }
634 else for(long i=0;i<n;i++)
635 {
636 if(dat->vthr(i*di)==val) pos.push_back(i);
637 }
638 pos.push_back(n); // add last point (size of data)
639 // now collect required position from section and its lengths
640 std::vector<long> ls, ps;
641 long np = pos.size()-1, nl=0;
642 if(np<1) return NULL; // nothing to do
643 for(long i=0;i<ids->GetNN();i++)
644 {
645 long j = mgl_int(ids->vthr(i)+0.5); j = j<0?np+j:j;
646 if(j>=0 && j<np)
647 { long l = pos[j+1]-pos[j]; nl += l;
648 ls.push_back(l); ps.push_back(pos[j]); }
649 }
650 if(nl==0) return NULL;
651 mglDataC *r=0;
652 size_t ns = ps.size();
653 if(dir=='y')
654 {
655 long nx=dat->GetNx(), nz=dat->GetNz(), sh=0;
656 r = new mglDataC(nx,nl,nz);
657 for(size_t s=0;s<ns;s++)
658 {
659 long pp = ps[s];
660 #pragma omp parallel for collapse(3)
661 for(long k=0;k<nz;k++) for(long j=0;j<ls[s];j++) for(long i=0;i<nx;i++)
662 r->a[i+nx*(sh+j+nl*k)] = dat->vc(i,pp+j,k);
663 sh += ls[s];
664 }
665 }
666 else if(dir=='x')
667 {
668 long ny=dat->GetNy(), nz=dat->GetNz(), sh=0;
669 r = new mglDataC(nl,ny,nz);
670 for(size_t s=0;s<ns;s++)
671 {
672 long pp = ps[s];
673 #pragma omp parallel for collapse(3)
674 for(long k=0;k<nz;k++) for(long j=0;j<ny;j++) for(long i=0;i<ls[s];i++)
675 r->a[sh+i+nl*(j+ny*k)] = dat->vc(pp+i,j,k);
676 sh += ls[s];
677 }
678 }
679 else if(dir=='z')
680 {
681 long nx=dat->GetNx(), ny=dat->GetNy(), sh=0;
682 r = new mglDataC(nx,ny,nl);
683 for(size_t s=0;s<ns;s++)
684 {
685 long pp = ps[s];
686 #pragma omp parallel for collapse(3)
687 for(long k=0;k<ls[s];k++) for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
688 r->a[i+nx*(j+ny*(sh+k))] = dat->vc(i,j,pp+k);
689 sh += ls[s];
690 }
691 }
692 return r;
693 }
mgl_datac_section_val(HCDT dat,long id,char dir,mreal val)694 HADT MGL_EXPORT mgl_datac_section_val(HCDT dat, long id, char dir, mreal val)
695 { mglData v; v.a[0]=id; return mgl_datac_section(dat,&v,dir,val); }
mgl_datac_section_(uintptr_t * d,uintptr_t * ids,const char * dir,mreal * val,int)696 uintptr_t MGL_EXPORT mgl_datac_section_(uintptr_t *d, uintptr_t *ids, const char *dir, mreal *val,int)
697 { return uintptr_t(mgl_datac_section(_DT_,_DA_(ids),dir[0],*val)); }
mgl_datac_section_val_(uintptr_t * d,int * id,const char * dir,mreal * val,int)698 uintptr_t MGL_EXPORT mgl_datac_section_val_(uintptr_t *d, int *id, const char *dir, mreal *val,int)
699 { return uintptr_t(mgl_datac_section_val(_DT_,*id,dir[0],*val)); }
700 //-----------------------------------------------------------------------------
mgl_find_roots_txt_c(const char * func,const char * vars,HCDT ini)701 HADT MGL_EXPORT mgl_find_roots_txt_c(const char *func, const char *vars, HCDT ini)
702 {
703 if(!vars || !(*vars) || !func || !ini) return 0;
704 mglEqTxT par;
705 par.var=vars; par.FillCmplx(func);
706 size_t n = par.str.size();
707 if(ini->GetNx()!=long(n)) return 0;
708 mreal *xx = new mreal[2*n];
709 mglDataC *res = new mglDataC(ini);
710 for(long j=0;j<ini->GetNy()*ini->GetNz();j++)
711 {
712 for(size_t i=0;i<n;i++)
713 { dual c = ini->vcthr(i+n*j); xx[2*i] = real(c); xx[2*i+1] = imag(c); }
714 bool ok=mgl_find_roots(2*n,mgl_txt_funcC,xx,&par);
715 for(size_t i=0;i<n;i++) res->a[i+n*j] = ok?dual(xx[2*i],xx[2*i+1]):NAN;
716 }
717 delete []xx; return res;
718 }
mgl_find_roots_txt_c_(const char * func,const char * vars,uintptr_t * ini,int l,int m)719 uintptr_t MGL_EXPORT mgl_find_roots_txt_c_(const char *func, const char *vars, uintptr_t *ini,int l,int m)
720 { char *s=new char[l+1]; memcpy(s,func,l); s[l]=0;
721 char *v=new char[m+1]; memcpy(v,vars,m); v[m]=0;
722 uintptr_t r = uintptr_t(mgl_find_roots_txt_c(s,v,_DA_(ini)));
723 delete []s; delete []v; return r; }
724 //-----------------------------------------------------------------------------
mgl_datac_keep(HADT dat,const char * how,long i,long j)725 void MGL_EXPORT mgl_datac_keep(HADT dat, const char *how, long i, long j)
726 {
727 const long nx = dat->GetNx(), ny = dat->GetNy(), nz = dat->GetNz();
728 const bool phase = !mglchr(how, 'a');
729 if(mglchr(how,'z'))
730 {
731 long ix = i, iy = j;
732 if(ix<0 || ix>=nx) ix = 0;
733 if(iy<0 || iy>=ny) iy = 0;
734 const long i0 = ix+nx*iy, nn = nx*ny;
735 const dual v0 = dat->a[i0];
736 for(long k=0;k<nz;k++)
737 {
738 dual v = dat->a[i0+nn*k], f = v0/v;
739 if(phase) f /= abs(f);
740 for(long ii=0;ii<nn;ii++) dat->a[ii+nn*k] *= f;
741 }
742 }
743 else if(mglchr(how,'x'))
744 {
745 long iy = i, iz = j;
746 if(iz<0 || iz>=nz) iz = 0;
747 if(iy<0 || iy>=ny) iy = 0;
748 const long i0 = nx*(iy+ny*iz), nn = ny*nz;
749 const dual v0 = dat->a[i0];
750 for(long k=0;k<nx;k++)
751 {
752 dual v = dat->a[i0+k], f = v0/v;
753 if(phase) f /= abs(f);
754 for(long ii=0;ii<nn;ii++) dat->a[k+nx*ii] *= f;
755 }
756 }
757 else // default is "y"
758 {
759 long ix = i, iz = j;
760 if(ix<0 || ix>=nx) ix = 0;
761 if(iz<0 || iz>=nz) iz = 0;
762 const long i0 = ix+nx*ny*iz;
763 const dual v0 = dat->a[i0];
764 for(long k=0;k<ny;k++)
765 {
766 dual v = dat->a[i0+nx*k], f = v0/v;
767 if(phase) f /= abs(f);
768 for(long ii=0;ii<nz;ii++) for(long jj=0;jj<nx;jj++) dat->a[jj+nx*(k+ny*ii)] *= f;
769 }
770 }
771 }
mgl_datac_keep_(uintptr_t * d,const char * how,long * i,long * j,int l)772 void MGL_EXPORT mgl_datac_keep_(uintptr_t *d, const char *how, long *i, long *j, int l)
773 { char *s=new char[l+1]; memcpy(s,how,l); s[l]=0;
774 mgl_datac_keep(_DC_, s, *i, *j); delete []s; }
775 //-----------------------------------------------------------------------------
776