--- /dev/null
+/***************************************************************************\r
+ * mgl_data_io.cpp is part of Math Graphic Library\r
+ * Copyright (C) 2007 Alexey Balakin <balakin@appl.sci-nnov.ru> *\r
+ * *\r
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU Library General Public License as *
+ * published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU Library General Public *
+ * License along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#include <stdlib.h>\r
+#include <ctype.h>\r
+#include <math.h>\r
+#include <string.h>\r
+#include <zlib.h>
+#ifdef HAVE_HDF5\r
+#include <hdf5.h>\r
+#endif\r
+#ifdef HAVE_HDF4\r
+#define intf hdf4_intf
+#include <hdf/mfhdf.h>
+#undef intf
+#endif\r
+\r
+#ifndef WIN32\r
+#include <glob.h>\r
+#endif\r
+\r
+#include "mgl/mgl_eval.h"\r
+#include "mgl/mgl_data.h"\r
+\r
+//#define isn(ch) ((ch)<' ' && (ch)!='\t')\r
+#define isn(ch) ((ch)=='\n')\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(const char *v,int NX,int NY,int NZ)\r
+{\r
+ if(NX<1 || NY <1 || NZ<1) return;\r
+ register int i,j=0,m=NX*NY*NZ;\r
+ mreal *b = new mreal[m]; memset(b,0,m*sizeof(mreal));\r
+ for(i=0;i<m;i++)\r
+ {\r
+ while(isspace(v[j]) && v[j]) j++;\r
+ if(v[j]==0) break;\r
+ b[i] = atof(v+j);\r
+ while(!isspace(v[j])&& v[j]) j++;\r
+ }\r
+ delete []a;\r
+ a=b; nx=NX; ny=NY; nz=NZ; NewId();\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(gsl_vector *v)\r
+{\r
+#ifndef NO_GSL\r
+ if(!v || v->size<1) return;\r
+ Create(v->size);\r
+ for(long i=0;i<nx;i++) a[i] = v->data[i*v->stride];\r
+#endif\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(gsl_matrix *m)\r
+{\r
+#ifndef NO_GSL\r
+ if(!m || m->size1<1 || m->size2<1) return;\r
+ Create(m->size1,m->size2);\r
+ register long i,j;\r
+ for(i=0;i<nx;i++) for(j=0;j<ny;j++)\r
+ a[i+j*nx] = m->data[i * m->tda + j];\r
+#endif\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(const float *A,int NX,int NY,int NZ)\r
+{\r
+ if(NX<=0 || NY<=0 || NZ<=0) return;\r
+ Create(NX,NY,NZ);\r
+#if(MGL_USE_DOUBLE==1)\r
+ for(long i=0;i<NX*NY*NZ;i++) a[i] = A[i];\r
+#else\r
+ memcpy(a,A,NX*NY*NZ*sizeof(float));\r
+#endif\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(const double *A,int NX,int NY,int NZ)\r
+{\r
+ if(NX<=0 || NY<=0 || NZ<=0) return;\r
+ Create(NX,NY,NZ);\r
+#if(MGL_USE_DOUBLE==1)\r
+ memcpy(a,A,NX*NY*NZ*sizeof(double));\r
+#else\r
+ for(long i=0;i<NX*NY*NZ;i++) a[i] = A[i];\r
+#endif\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(const float **A,int N1,int N2)\r
+{\r
+ if(N1<=0 || N2<=0) return;\r
+ Create(N2,N1);\r
+#if(MGL_USE_DOUBLE==1)\r
+ for(long i=0;i<N1;i++) for(long j=0;j<N2;j++) a[j+i*N2] = A[i][j];\r
+#else\r
+ for(long i=0;i<N1;i++)\r
+ memcpy(a+i*N2,A[i],N2*sizeof(float));\r
+#endif\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(const double **A,int N1,int N2)\r
+{\r
+ if(N1<=0 || N2<=0) return;\r
+ Create(N2,N1);\r
+#if(MGL_USE_DOUBLE==1)\r
+ for(long i=0;i<N1;i++)\r
+ memcpy(a+i*N2,A[i],N2*sizeof(double));\r
+#else\r
+ for(long i=0;i<N1;i++) for(long j=0;j<N2;j++) a[j+i*N2] = A[i][j];\r
+#endif\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(const float ***A,int N1,int N2,int N3)\r
+{\r
+ if(N1<=0 || N2<=0 || N3<=0) return;\r
+ Create(N3,N2,N1);\r
+#if(MGL_USE_DOUBLE==1)\r
+ for(long i=0;i<N1;i++) for(long j=0;j<N2;j++) for(long k=0;k<N3;k++)\r
+ a[k+N3*(j+i*N2)] = A[i][j][k];\r
+#else\r
+ for(long i=0;i<N1;i++) for(long j=0;j<N2;j++)\r
+ memcpy(a+N3*(j+i*N2),A[i][j],N3*sizeof(float));\r
+#endif\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(const double ***A,int N1,int N2,int N3)\r
+{\r
+ if(N1<=0 || N2<=0 || N3<=0) return;\r
+ Create(N3,N2,N1);\r
+#if(MGL_USE_DOUBLE==1)\r
+ for(long i=0;i<N1;i++) for(long j=0;j<N2;j++)\r
+ memcpy(a+N3*(j+i*N2),A[i][j],N3*sizeof(double));\r
+#else\r
+ for(long i=0;i<N1;i++) for(long j=0;j<N2;j++) for(long k=0;k<N3;k++)\r
+ a[k+N3*(j+i*N2)] = A[i][j][k];\r
+#endif\r
+}\r
+//-----------------------------------------------------------------------------\r
+mglData mglData::Trace() const\r
+{\r
+ mglData r(nx);\r
+ register long i;\r
+ if(ny>=nx && nz>=nx) for(i=0;i<nx;i++) r.a[i] = a[i+nx*(i+ny*i)];\r
+ else if(ny>=nx) for(i=0;i<nx;i++) r.a[i] = a[i+nx*i];\r
+ else memcpy(r.a,a,nx*sizeof(mreal));\r
+ return r;\r
+}\r
+//-----------------------------------------------------------------------------\r
+mglData mglData::SubData(const mglData &xx, const mglData &yy, const mglData &zz) const\r
+{\r
+ long n=0,m=0,l=0,i,j,k,i0,x,y,z;\r
+ mglData d;\r
+ bool ix=false, iy=false, iz=false;\r
+ if(xx.nz>1) // 3d data\r
+ {\r
+ n = xx.nx; m = xx.ny; l = xx.nz;\r
+ j = yy.nx*yy.ny*yy.nz; if(j>1 && j!=n*m*l) return d; // wrong sizes\r
+ k = zz.nx*zz.ny*zz.nz; if(k>1 && k!=n*m*l) return d; // wrong sizes\r
+ ix = true; iy = j>1; iz = k>1;\r
+ }\r
+ else if(yy.nz>1)\r
+ {\r
+ n = yy.nx; m = yy.ny; l = yy.nz;\r
+ j = xx.nx*xx.ny*xx.nz; if(j>1 && j!=n*m*l) return d; // wrong sizes\r
+ k = zz.nx*zz.ny*zz.nz; if(k>1 && k!=n*m*l) return d; // wrong sizes\r
+ iy = true; ix = j>1; iz = k>1;\r
+ }\r
+ else if(zz.nz>1)\r
+ {\r
+ n = zz.nx; m = zz.ny; l = zz.nz;\r
+ j = yy.nx*yy.ny*yy.nz; if(j>1 && j!=n*m*l) return d; // wrong sizes\r
+ k = xx.nx*xx.ny*xx.nz; if(k>1 && k!=n*m*l) return d; // wrong sizes\r
+ iz = true; iy = j>1; ix = k>1;\r
+ }\r
+ else if(xx.ny>1) // 2d data\r
+ {\r
+ n = xx.nx; m = xx.ny; l = 1;\r
+ j = yy.nx*yy.ny; if(j>1 && j!=n*m) return d; // wrong sizes\r
+ k = zz.nx*zz.ny; if(k>1 && k!=n*m) return d; // wrong sizes\r
+ ix = true; iy = j>1; iz = k>1;\r
+ }\r
+ else if(yy.ny>1)\r
+ {\r
+ n = yy.nx; m = yy.ny; l = 1;\r
+ j = xx.nx*xx.ny; if(j>1 && j!=n*m) return d; // wrong sizes\r
+ k = zz.nx*zz.ny; if(k>1 && k!=n*m) return d; // wrong sizes\r
+ iy = true; ix = j>1; iz = k>1;\r
+ }\r
+ else if(zz.ny>1)\r
+ {\r
+ n = zz.nx; m = zz.ny; l = 1;\r
+ j = yy.nx*yy.ny; if(j>1 && j!=n*m) return d; // wrong sizes\r
+ k = xx.nx*xx.ny; if(k>1 && k!=n*m) return d; // wrong sizes\r
+ iz = true; iy = j>1; ix = k>1;\r
+ }\r
+ if(n*m*l>1) // this is 2d or 3d data\r
+ {\r
+ d.Create(n,m,l);\r
+ for(i0=0;i0<n*m*l;i0++)\r
+ {\r
+ i = int((ix?xx.a[i0]:xx.a[0])+0.5); if(i<0)i=0; if(i>=nx)i=nx-1;\r
+ j = int((iy?yy.a[i0]:yy.a[0])+0.5); if(j<0)j=0; if(j>=ny)j=ny-1;\r
+ k = int((iz?zz.a[i0]:zz.a[0])+0.5); if(k<0)k=0; if(k>=nz)k=nz-1;\r
+ d.a[i0] = a[i+nx*(j+ny*k)];\r
+ }\r
+ return d;\r
+ }\r
+ // this is 1d data -> try as normal SubData()\r
+ if(xx.nx>1 || xx.a[0]>=0) { n=xx.nx; ix=true; }\r
+ else { n=nx; ix=false; }\r
+ if(yy.nx>1 || yy.a[0]>=0) { m=yy.nx; iy=true; }\r
+ else { m=ny; iy=false; }\r
+ if(zz.nx>1 || zz.a[0]>=0) { l=zz.nx; iz=true; }\r
+ else { l=nz; iz=false; }\r
+ d.Create(n,m,l);\r
+ for(i=0;i<n;i++) for(j=0;j<m;j++) for(k=0;k<l;k++)\r
+ {\r
+ x = ix?int(xx.a[i]+0.5):i; if(x<0)x=0; if(x>=nx)x=nx-1;\r
+ y = iy?int(yy.a[j]+0.5):j; if(y<0)y=0; if(y>=ny)y=ny-1;\r
+ z = iz?int(zz.a[k]+0.5):k; if(z<0)z=0; if(z>=nz)z=nz-1;\r
+ d.a[i+n*(j+m*k)] = a[x+nx*(y+ny*z)];\r
+ }\r
+ if(m==1) { d.ny=d.nz; d.nz=1; }// "squeeze" dimensions\r
+ if(n==1) { d.nx=d.ny; d.ny=d.nz; d.nz=1; d.NewId();}\r
+ return d;\r
+}\r
+//-----------------------------------------------------------------------------\r
+mglData mglData::SubData(int xx,int yy,int zz) const\r
+{\r
+ mglData x,y,z;\r
+ x.a[0]=xx; y.a[0]=yy; z.a[0]=zz;\r
+ return SubData(x,y,z);\r
+}\r
+//-----------------------------------------------------------------------------\r
+mglData mglData::Column(const char *eq)\r
+{\r
+ mglFormula f(eq);\r
+ mglData d;\r
+ d.Create(ny,nz);\r
+ mreal var[MGL_VS];\r
+ memset(var,0,('z'-'a')*sizeof(mreal));\r
+ register long i,j;\r
+ for(i=0;i<ny*nz;i++)\r
+ {\r
+ for(j=0;j<nx;j++)\r
+ if(id[j]>='a' && id[j]<='z')\r
+ var[id[j]-'a'] = a[j+nx*i];\r
+ d.a[i] = f.Calc(var);\r
+ }\r
+ return d;\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::SetColumnId(const char *ids)\r
+{\r
+ NewId(); // clearing + be sure about correct length\r
+ if(ids) for(long i=0;i<nx && ids[i]!=0;i++) id[i]=ids[i];\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Save(const char *fname,int ns) const\r
+{\r
+ FILE *fp;\r
+ fp = fopen(fname,"w");\r
+ if(ns<0 || (ns>=nz && nz>1)) for(long k=0;k<nz;k++)\r
+ { // ñîõðàíÿåì âñå äàííûå\r
+ for(long i=0;i<ny;i++)\r
+ {\r
+ for(long j=0;j<nx;j++)\r
+ fprintf(fp,"%g\t",a[j+nx*(i+ny*k)]);\r
+ fprintf(fp,"\n");\r
+ }\r
+ fprintf(fp,"\n");\r
+ }\r
+ else\r
+ { // ñîõðàíÿåì òîëüêî ñðåç\r
+ if(nz>1) // äëÿ 3D -- ïëîñêîñòü\r
+ {\r
+ for(long i=0;i<ny;i++)\r
+ {\r
+ for(long j=0;j<nx;j++)\r
+ fprintf(fp,"%g\t",a[j+nx*(i+ny*ns)]);\r
+ fprintf(fp,"\n");\r
+ }\r
+ }\r
+ else if(ns<ny) // äëÿ 2D -- ëèíèÿ\r
+ {\r
+ for(long j=0;j<nx;j++)\r
+ fprintf(fp,"%g\t",a[j+nx*ns]);\r
+ }\r
+ }\r
+ fclose(fp);\r
+}\r
+//-----------------------------------------------------------------------------
+char *mgl_read_gz(gzFile fp)
+{
+ long size=1024,n=0,m;
+ char *buf=(char*)malloc(size);
+ while((m=gzread(fp,buf+size*n,size))>0)
+ {
+ if(m<size) { buf[size*n+m]=0; break; }
+ n++; buf=(char*)realloc(buf,size*(n+1));
+ }
+ return buf;
+}\r
+//-----------------------------------------------------------------------------\r
+bool mglData::Read(const char *fname)\r
+{\r
+ long l=1,m=1,k=1;\r
+ long nb,i;\r
+ gzFile fp = gzopen(fname,"r");\r
+ if(!fp)\r
+ {
+ if(!a) Create(1,1,1);\r
+ return false;\r
+ }\r
+ char *buf = mgl_read_gz(fp);
+ nb = strlen(buf); gzclose(fp);\r
+\r
+ bool first=false,com=false;\r
+ register char ch;\r
+ for(i=nb-1;i>=0;i--) if(buf[i]>' ') break;\r
+ buf[i+1]=0; nb = i; // remove tailing spaces
+ for(i=0;i<nb-1 && !isn(buf[i]);i++) // determine nx\r
+ {\r
+ if(buf[i]=='#') while(!isn(buf[i]) && i<nb) i++;\r
+ ch = buf[i];\r
+ if(ch>' ' && !first) first=true;\r
+ if(first && (ch==' ' || ch=='\t') && buf[i+1]>' ') k++;
+ }\r
+ first = false;\r
+ for(i=0;i<nb-1;i++) // determine ny\r
+ {\r
+ ch = buf[i];\r
+ if(ch=='#') while(!isn(buf[i]) && i<nb) i++;\r
+ if(isn(ch))\r
+ {\r
+ if(isn(buf[i+1])) {first=true; break; }\r
+ m++;\r
+ }\r
+ if(ch=='\f') break;\r
+ }\r
+ if(first) for(i=0;i<nb-1;i++) // determine nz\r
+ {\r
+ ch = buf[i];\r
+ if(ch=='#') com = true; // comment\r
+ if(isn(ch))\r
+ {\r
+ if(com) { com=false; continue; }\r
+ if(isn(buf[i+1])) l++;\r
+ }\r
+ }\r
+ else for(i=0;i<nb-1;i++) if(buf[i]=='\f') l++;\r
+ free(buf);
+ return Read(fname,k,m,l);\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Create(int mx,int my,int mz)\r
+{\r
+ nx = mx>0 ? mx:1; ny = my>0 ? my:1; nz = mz>0 ? mz:1;\r
+ if(a) { delete []a; delete []id; }\r
+ a = new mreal[nx*ny*nz];\r
+ id = new char[nx];\r
+ memset(a,0,nx*ny*nz*sizeof(mreal));\r
+ memset(id,0,nx*sizeof(char));\r
+}\r
+//-----------------------------------------------------------------------------\r
+bool mglData::Read(const char *fname,int mx,int my,int mz)\r
+{\r
+ if(mx<=0 || my<=0 || mz<=0) return false;\r
+ gzFile fp = gzopen(fname,"r");\r
+ if(!fp) return false;\r
+ Create(mx,my,mz);\r
+ char *buf = mgl_read_gz(fp);
+ long nb = strlen(buf); gzclose(fp);\r
+\r
+ register long i=0, j=0, k=0;\r
+ while(j<nb)\r
+ {\r
+ while(buf[j]<=' ' && j<nb) j++;\r
+ while(buf[j]=='#') // skip comment\r
+ {\r
+ if(i>0 || buf[j+1]!='#') // this is columns id\r
+ while(!isn(buf[j]) && j<nb) j++;\r
+ else\r
+ {\r
+ while(!isn(buf[j]) && j<nb)\r
+ {\r
+ if(buf[j]>='a' && buf[j]<='z')\r
+ id[k++] = buf[j];\r
+ j++;\r
+ }\r
+ }\r
+// while(buf[j]!='\n' && j<nb) j++;\r
+ while(buf[j]<=' ' && j<nb) j++;\r
+ }\r
+ char *s=buf+j;\r
+ while(buf[j]>' ' && j<nb) j++;\r
+ buf[j]=0;\r
+ a[i] = atof(s);
+ i++; if(i>=nx*ny*nz) break;
+ }\r
+ free(buf);\r
+ return true;\r
+}\r
+//-----------------------------------------------------------------------------\r
+bool mglData::ReadMat(const char *fname,int dim)\r
+{\r
+ if(dim<=0 || dim>3) return false;\r
+ gzFile fp = gzopen(fname,"r");\r
+ if(!fp) return false;\r
+ nx = ny = nz = 1; NewId();\r
+ char *buf = mgl_read_gz(fp);
+ long nb = strlen(buf); gzclose(fp);\r
+\r
+ register long i=0,j=0;\r
+ while(j<nb)\r
+ {\r
+ if(buf[j]=='#') while(!isn(buf[j])) j++; // skip comment\r
+ while(buf[j]<=' ' && j<nb) j++;\r
+ break;\r
+ }\r
+ if(dim==1)\r
+ {\r
+ sscanf(buf+j,"%ld",&nx);\r
+ while(buf[j]>' ') j++;\r
+ }\r
+ else if(dim==2)\r
+ {\r
+ sscanf(buf+j,"%ld%ld",&nx,&ny);\r
+ while(buf[j]>' ' && j<nb) j++;\r
+ while(buf[j]<=' ' && j<nb) j++;\r
+ while(buf[j]>' ' && j<nb) j++;\r
+ }\r
+ else if(dim==3)\r
+ {\r
+ sscanf(buf+j,"%ld%ld%ld",&nx,&ny,&nz);\r
+ while(buf[j]>' ' && j<nb) j++;\r
+ while(buf[j]<=' ' && j<nb) j++;\r
+ while(buf[j]>' ' && j<nb) j++;\r
+ while(buf[j]<=' ' && j<nb) j++;\r
+ while(buf[j]>' ' && j<nb) j++;\r
+ }\r
+ Create(nx,ny,nz);\r
+ while(j<nb)\r
+ {\r
+ while(buf[j]<=' ' && j<nb) j++;\r
+ while(buf[j]=='#') // skip comment\r
+ {\r
+ while(!isn(buf[j]) && j<nb) j++;\r
+ while(buf[j]<=' ' && j<nb) j++;\r
+ }\r
+ a[i] = atof(buf+j); i++;\r
+ if(i>=nx*ny*nz) break;\r
+ while(buf[j]>' ' && j<nb) j++;\r
+ }\r
+ free(buf);\r
+ return true;\r
+}\r
+//-----------------------------------------------------------------------------\r
+mreal mglData::v(int i,int j,int k) const\r
+{\r
+ bool not_ok = i<0 || i>=nx || j<0 || j>=ny || k<0 || k>=nz;\r
+ if(not_ok) return 0;\r
+ return a[i+nx*(j+ny*k)];\r
+}\r
+//-----------------------------------------------------------------------------\r
+mglData mglData::Resize(int mx, int my, int mz, mreal x1, mreal x2,\r
+ mreal y1, mreal y2, mreal z1, mreal z2) const\r
+{\r
+ register long i,j,k;\r
+ mglData d;\r
+ mx = mx<1 ? 1:mx; my = my<1 ? 1:my; mz = mz<1 ? 1:mz;\r
+ d.Create(mx,my,mz);\r
+ mreal dx, dy, dz;\r
+ dx = mx>1 ? (x2-x1)/(mx-1):0;\r
+ dy = my>1 ? (y2-y1)/(my-1):0;\r
+ dz = mz>1 ? (z2-z1)/(mz-1):0;\r
+ for(i=0;i<mx;i++) for(j=0;j<my;j++) for(k=0;k<mz;k++)\r
+ d.a[i+mx*(j+my*k)] = Spline1(x1+i*dx, y1+j*dy, z1+k*dz);\r
+ return d;\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::operator*=(const mglData &d)\r
+{\r
+ register long i,j;\r
+ if(d.nz==1 && d.ny==1 && nx==d.nx)\r
+ {\r
+ for(j=0;j<ny*nz;j++) for(i=0;i<nx;i++) a[i+nx*j] *= d.a[i];\r
+ }\r
+ else if(d.nz==1 && ny==d.ny && nx==d.nx)\r
+ {\r
+ for(j=0;j<nz;j++) for(i=0;i<ny*nx;i++) a[i+nx*ny*j] *= d.a[i];\r
+ }\r
+ else if(nz==d.nz && d.ny==ny && nx==d.nx)\r
+ {\r
+ for(i=0;i<ny*nz*nx;i++) a[i] *= d.a[i];\r
+ }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::operator/=(const mglData &d)\r
+{\r
+ register long i,j;\r
+ if(d.nz==1 && d.ny==1 && nx==d.nx)\r
+ {\r
+ for(j=0;j<ny*nz;j++) for(i=0;i<nx;i++)\r
+ a[i+nx*j] = d.a[i] ? a[i+nx*j]/d.a[i] : 0;\r
+ }\r
+ else if(d.nz==1 && ny==d.ny && nx==d.nx)\r
+ {\r
+ for(j=0;j<nz;j++) for(i=0;i<ny*nx;i++)\r
+ a[i+nx*ny*j] = d.a[i] ? a[i+nx*ny*j]/d.a[i] : 0;\r
+ }\r
+ else if(nz==d.nz && d.ny==ny && nx==d.nx)\r
+ {\r
+ for(i=0;i<ny*nz*nx;i++)\r
+ a[i] = d.a[i] ? a[i]/d.a[i] : 0;\r
+ }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::operator+=(const mglData &d)\r
+{\r
+ register long i,j;\r
+ if(d.nz==1 && d.ny==1 && nx==d.nx)\r
+ {\r
+ for(j=0;j<ny*nz;j++) for(i=0;i<nx;i++)\r
+ a[i+nx*j] += d.a[i];\r
+ }\r
+ else if(d.nz==1 && ny==d.ny && nx==d.nx)\r
+ {\r
+ for(j=0;j<nz;j++) for(i=0;i<ny*nx;i++)\r
+ a[i+nx*ny*j] += d.a[i];\r
+ }\r
+ else if(nz==d.nz && d.ny==ny && nx==d.nx)\r
+ {\r
+ for(i=0;i<ny*nz*nx;i++)\r
+ a[i] += d.a[i];\r
+ }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::operator-=(const mglData &d)\r
+{\r
+ register long i,j;\r
+ if(d.nz==1 && d.ny==1 && nx==d.nx)\r
+ {\r
+ for(j=0;j<ny*nz;j++) for(i=0;i<nx;i++)\r
+ a[i+nx*j] -= d.a[i];\r
+ }\r
+ else if(d.nz==1 && ny==d.ny && nx==d.nx)\r
+ {\r
+ for(j=0;j<nz;j++) for(i=0;i<ny*nx;i++)\r
+ a[i+nx*ny*j] -= d.a[i];\r
+ }\r
+ else if(nz==d.nz && d.ny==ny && nx==d.nx)\r
+ {\r
+ for(i=0;i<ny*nz*nx;i++)\r
+ a[i] -= d.a[i];\r
+ }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::operator-=(mreal d)\r
+{\r
+ for(long i=0;i<ny*nz*nx;i++) a[i] -= d;\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::operator+=(mreal d)\r
+{\r
+ for(long i=0;i<ny*nz*nx;i++) a[i] += d;\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::operator*=(mreal d)\r
+{\r
+ for(long i=0;i<ny*nz*nx;i++) a[i] *= d;\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::operator/=(mreal d)\r
+{\r
+ for(long i=0;i<ny*nz*nx;i++) a[i] = d ? a[i]/d : 0;\r
+}\r
+//-----------------------------------------------------------------------------\r
+mreal mglData::Maximal() const\r
+{\r
+ register mreal m=-1e10;\r
+ for(long i=0;i<nx*ny*nz;i++)\r
+ m = m>a[i] ? m : a[i];\r
+ return m;\r
+}\r
+//-----------------------------------------------------------------------------\r
+mreal mglData::Minimal() const\r
+{\r
+ register mreal m=1e10;\r
+ for(long i=0;i<nx*ny*nz;i++)\r
+ m = m<a[i] ? m : a[i];\r
+ return m;\r
+}\r
+//-----------------------------------------------------------------------------\r
+mreal mglData::Maximal(int &im,int &jm,int &km) const\r
+{\r
+ register mreal m=-1e10;\r
+ for(long i=0;i<nx*ny*nz;i++)\r
+ if(!isnan(a[i]) && m<a[i])\r
+ { m=a[i]; im=i%nx; jm=(i/nx)%ny; km=i/(nx*ny); }\r
+ return m;\r
+}\r
+//-----------------------------------------------------------------------------\r
+mreal mglData::Minimal(int &im,int &jm,int &km) const\r
+{\r
+ register mreal m=1e10;\r
+ for(long i=0;i<nx*ny*nz;i++)\r
+ if(!isnan(a[i]) && m>a[i])\r
+ { m=a[i]; im=i%nx; jm=(i/nx)%ny; km=i/(nx*ny); }\r
+ return m;\r
+}\r
+//-----------------------------------------------------------------------------\r
+mreal mglData::Maximal(mreal &x,mreal &y,mreal &z) const\r
+{\r
+ int im=-1,jm=-1,km=-1;\r
+ register long tm,i;\r
+ mreal m=Maximal(im,jm,km);\r
+ x=im; y=jm; z=km;\r
+\r
+ if(nx>2)\r
+ {\r
+ if(im==0) im=1;\r
+ if(im==nx-1)im=nx-2;\r
+ x = (a[im+1]+a[im-1]-2*a[im])==0 ? im : im+(a[im+1]-a[im-1])/(a[im+1]+a[im-1]-2*a[im])/2;\r
+ }\r
+ if(ny>2)\r
+ {\r
+ if(jm==0) jm=1;\r
+ if(jm==ny-1)jm=ny-2;\r
+ i=nx; tm = jm*nx;\r
+ y = (a[tm+i]+a[tm-i]-2*a[tm])==0? jm : jm+(a[tm+i]-a[tm-i])/(a[tm+i]+a[tm-i]-2*a[tm])/2;\r
+ }\r
+ if(nz>2)\r
+ {\r
+ if(km==0) km=1;\r
+ if(km==nz-1)km=nz-2;\r
+ i=nx*ny; tm = km*i;\r
+ z = (a[tm+i]+a[tm-i]-2*a[tm])==0? km : km+(a[tm+i]-a[tm-i])/(a[tm+i]+a[tm-i]-2*a[tm])/2;\r
+ }\r
+ return m;\r
+}\r
+//-----------------------------------------------------------------------------\r
+mreal mglData::Minimal(mreal &x,mreal &y,mreal &z) const\r
+{\r
+ int im=-1,jm=-1,km=-1;\r
+ register long tm,i;\r
+ mreal m=Minimal(im,jm,km);\r
+ x=im; y=jm; z=km;\r
+ if(nx>2)\r
+ {\r
+ if(im==0) im=1;\r
+ if(im==nx-1)im=nx-2;\r
+ x = im+(a[im+1]-a[im-1])/(a[im+1]+a[im-1]-2*a[im])/2;\r
+ }\r
+ if(ny>2)\r
+ {\r
+ if(jm==0) jm=1;\r
+ if(jm==ny-1)jm=ny-2;\r
+ i=nx; tm = jm*nx;\r
+ y = jm+(a[tm+i]-a[tm-i])/(a[tm+i]+a[tm-i]-2*a[tm])/2;\r
+ }\r
+ if(nz>2)\r
+ {\r
+ if(km==0) km=1;\r
+ if(km==nz-1)km=nz-2;\r
+ i=nx*ny; tm = km*i;\r
+ z = km+(a[tm+i]-a[tm-i])/(a[tm+i]+a[tm-i]-2*a[tm])/2;\r
+ }\r
+ return m;\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Modify(const char *eq,int dim)\r
+{\r
+ long i,j,k;\r
+ mreal y,z,dx=nx>1?1/(nx-1.):0,dy=ny>1?1/(ny-1.):0, *aa;\r
+ mglFormula eqs(eq);\r
+ if(dim<0) dim=0;\r
+ if(nz>1) // 3D array\r
+ {\r
+ for(k=dim;k<nz;k++)\r
+ {\r
+ z = (nz>dim+1) ? (k-dim)/(nz-dim-1.) : 0;
+ aa = a+nx*ny*k;
+//#pragma omp parallel for
+ for(i=0;i<nx*ny;i++)
+ aa[i] = eqs.Calc(dx*(i%nx),dy*(i/nx),z,aa[i]);\r
+ }\r
+ }\r
+ else // 2D or 1D array\r
+ {\r
+ if(ny==1) dim = 0;\r
+ dy = ny>dim+1 ? 1/(ny-dim-1.) : 0;\r
+ for(j=dim;j<ny;j++)
+ {
+ y = dy*(j-dim); aa = a+nx*j;
+//#pragma omp parallel for
+ for(i=0;i<nx;i++) aa[i] = eqs.Calc(dx*i,y,0,aa[i]);\r
+ }\r
+ }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Fill(mreal x1,mreal x2,char dir)\r
+{\r
+ long i,j,k;\r
+ register mreal x;\r
+ if(isnan(x2)) x2=x1;\r
+ if(dir<'x' || dir>'z') dir='x';\r
+ for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++)\r
+ {\r
+ x = nx>1 ? i/(nx-1.):0;\r
+ if(dir=='y') x = ny>1 ? j/(ny-1.):0;\r
+ if(dir=='z') x = nz>1 ? k/(nz-1.):0;\r
+ a[i+nx*(j+ny*k)] = x1+(x2-x1)*x;\r
+ }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Norm(mreal v1,mreal v2,bool sym,int dim)\r
+{\r
+ long i,s,nn=nx*ny*nz;\r
+ mreal a1=1e20,a2=-1e20,v;\r
+ if(nz>1) s = dim*nx*ny;\r
+ else s = dim*ny;\r
+ for(i=s;i<nn;i++) // determines borders of existing data\r
+ {\r
+ if(isnan(a[i])) continue;\r
+ a1 = (a1<a[i] ? a1 : a[i]); a2 = (a2>a[i] ? a2 : a[i]);\r
+ }\r
+ if(a1==a2) { if(a1!=0) a1=0.; else a2=1; }\r
+ if(v1>v2) { v=v1; v1=v2; v2=v; } // swap if uncorrect\r
+ if(sym) // use symmetric\r
+ {\r
+ v2 = -v1>v2 ? -v1:v2; v1 = -v2;\r
+ a2 = -a1>a2 ? -a1:a2; a1 = -a2;\r
+ }\r
+ for(i=s;i<nn;i++) // normalize\r
+ {\r
+ a[i] = v1 + (v2-v1)*(a[i]-a1)/(a2-a1);\r
+ }\r
+\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Squeeze(int rx,int ry,int rz,bool smooth)\r
+{\r
+ long kx,ky,kz,i,j,k;\r
+ mreal *b;\r
+\r
+ // simple checking\r
+ if(rx>=nx) rx=nx-1; if(rx<1) rx=1;\r
+ if(ry>=ny) ry=ny-1; if(ry<1) ry=1;\r
+ if(rz>=nz) rz=nz-1; if(rz<1) rz=1;\r
+ // new sizes\r
+ kx = 1+(nx-1)/rx; ky = 1+(ny-1)/ry; kz = 1+(nz-1)/rz;\r
+ b = new mreal[kx*ky*kz];\r
+ if(!smooth) for(i=0;i<kx;i++) for(j=0;j<ky;j++) for(k=0;k<kz;k++)\r
+ b[i+kx*(j+ky*k)] = a[i*rx+nx*(j*ry+ny*rz*k)];\r
+ else for(i=0;i<kx;i++) for(j=0;j<ky;j++) for(k=0;k<kz;k++)\r
+ {\r
+ long dx,dy,dz,i1,j1,k1;\r
+ dx = (i+1)*rx<=nx ? rx : nx-i*rx;\r
+ dy = (j+1)*ry<=ny ? ry : ny-j*ry;\r
+ dz = (k+1)*rz<=nz ? rz : nz-k*rz;
+ mreal s = 0;\r
+ for(i1=i*rx;i1<i*rx+dx;i1++) for(j1=j*ry;j1<j*ry+dz;j1++) for(k1=k*rz;k1<k*rz+dz;k1++)\r
+ s += a[i1+nx*(j1+ny*k1)];
+ b[i+kx*(j+ky*k)] = s/dx*dy*dz;\r
+ }\r
+ delete []a; a=b;\r
+ nx = kx; ny = ky; nz = kz; NewId();\r
+}\r
+//-----------------------------------------------------------------------------\r
+mglData mglData::Combine(const mglData &b) const\r
+{\r
+ mglData d;\r
+ d.Create(1,1,1);\r
+ if(nz>1 || (ny>1 && b.ny>1) || b.nz>1) return d;\r
+ long n1=ny,n2=b.nx;\r
+ bool dim2=true;\r
+ if(ny==1) { n1 = b.nx; n2 = b.ny; dim2 = false; }\r
+ d.Create(nx,n1,n2);\r
+ register long i,j;\r
+ if(dim2) n1=nx*ny; else { n1=nx; n2=b.nx*b.ny; }\r
+ for(i=0;i<n1;i++) for(j=0;j<n2;j++) d.a[i+n1*j] = a[i]*b.a[j];\r
+ return d;\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Extend(int n1, int n2)\r
+{\r
+ if(nz>2 || n1==0) return;\r
+ long mx,my,mz;\r
+ mreal *b=0;\r
+ register long i,j;\r
+ if(n1>0) // extend to higher dimension(s)\r
+ {\r
+ n2 = n2>0 ? n2:1;\r
+ mx = nx; my = ny>1?ny:n1; mz = ny>1 ? n1 : n2;\r
+ b = new mreal[mx*my*mz];\r
+ if(ny>1) for(i=0;i<n1;i++)\r
+ memcpy(b+i*nx*ny, a, nx*ny*sizeof(mreal));\r
+ else for(i=0;i<n1*n2;i++)\r
+ memcpy(b+i*nx, a, nx*sizeof(mreal));\r
+ }\r
+ else\r
+ {\r
+ mx = -n1; my = n2<0 ? -n2 : nx; mz = n2<0 ? nx : ny;\r
+ if(n2>0 && ny==1) mz = n2;\r
+ b = new mreal[mx*my*mz];\r
+ if(n2<0) for(i=0;i<mx*my;i++) for(j=0;j<nx;j++)\r
+ b[i+mx*my*j] = a[j];\r
+ else for(i=0;i<mx;i++) for(j=0;j<nx*ny;j++)\r
+ b[i+mx*j] = a[j];\r
+ if(n2>0 && ny==1) for(i=0;i<n2;i++)\r
+ memcpy(b+i*mx*my, a, mx*my*sizeof(mreal));\r
+ }\r
+ if(b) { delete []a; a=b; nx=mx; ny=my; nz=mz; NewId(); }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Transpose(const char *dim)\r
+{\r
+ mreal *b=new mreal[nx*ny*nz];\r
+ register long i,j,k,n;\r
+ if(!strcmp(dim,"xyz")) memcpy(b,a,nx*ny*nz*sizeof(mreal));\r
+ else if(!strcmp(dim,"xzy") || !strcmp(dim,"zy"))\r
+ {\r
+ for(i=0;i<nx;i++) for(j=0;j<ny;j++) for(k=0;k<nz;k++)\r
+ b[i+nx*(k+nz*j)] = a[i+nx*(j+ny*k)];\r
+ n=nz; nz=ny; ny=n;\r
+ }\r
+ else if(!strcmp(dim,"yxz") || !strcmp(dim,"yx"))\r
+ {\r
+ for(i=0;i<nx;i++) for(j=0;j<ny;j++) for(k=0;k<nz;k++)\r
+ b[j+ny*(i+nx*k)] = a[i+nx*(j+ny*k)];\r
+ n=nx; nx=ny; ny=n; NewId();\r
+ }\r
+ else if(!strcmp(dim,"yzx"))\r
+ {\r
+ for(i=0;i<nx;i++) for(j=0;j<ny;j++) for(k=0;k<nz;k++)\r
+ b[j+ny*(k+nz*i)] = a[i+nx*(j+ny*k)];\r
+ n=nx; nx=ny; ny=nz; nz=n; NewId();\r
+ }\r
+ else if(!strcmp(dim,"zxy"))\r
+ {\r
+ for(i=0;i<nx;i++) for(j=0;j<ny;j++) for(k=0;k<nz;k++)\r
+ b[k+nz*(i+nx*j)] = a[i+nx*(j+ny*k)];\r
+ n=nx; nx=nz; nz=ny; ny=n; NewId();\r
+ }\r
+ else if(!strcmp(dim,"zyx") || !strcmp(dim,"zx"))\r
+ {\r
+ for(i=0;i<nx;i++) for(j=0;j<ny;j++) for(k=0;k<nz;k++)\r
+ b[k+nz*(j+ny*i)] = a[i+nx*(j+ny*k)];\r
+ n=nz; nz=nx; nx=n; NewId();\r
+ }\r
+ delete []a; a = b;\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Modify(const char *eq, const mglData &v, const mglData &w)\r
+{\r
+ if(v.nx*v.ny*v.nz!=nx*ny*nz || w.nx*w.ny*w.nz!=nx*ny*nz)\r
+ return;\r
+ long i,j,k,i0;\r
+ 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;\r
+ mglFormula eqs(eq);\r
+ for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++)\r
+ {\r
+ x = dx*i; y = dy*j; z = dz*k;\r
+ i0 = i+nx*(j+ny*k);\r
+ a[i0] = eqs.Calc(x,y,z,a[i0],v.a[i0],w.a[i0]);\r
+ }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Modify(const char *eq, const mglData &v)\r
+{\r
+ if(v.nx*v.ny*v.nz!=nx*ny*nz) return;\r
+ long i,j,k,i0;\r
+ 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;\r
+ mglFormula eqs(eq);\r
+ for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++)\r
+ {\r
+ x = dx*i; y = dy*j; z = dz*k;\r
+ i0 = i+nx*(j+ny*k);\r
+ a[i0] = eqs.Calc(x,y,z,a[i0],v.a[i0],0);\r
+ }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::Fill(const char *eq, mglPoint r1, mglPoint r2, const mglData *v, const mglData *w)\r
+{\r
+ if(v && v->nx*v->ny*v->nz!=nx*ny*nz) return;\r
+ if(w && w->nx*w->ny*w->nz!=nx*ny*nz) return;\r
+ long i,j,k,i0;\r
+ mreal x,y,z,dx=nx>1?(r2.x-r1.x)/(nx-1.):0;\r
+ mreal dy=ny>1?(r2.y-r1.y)/(ny-1.):0;\r
+ mreal dz=nz>1?(r2.z-r1.z)/(nz-1.):0;\r
+ mglFormula eqs(eq);\r
+ for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++)\r
+ {\r
+ x = r1.x+dx*i; y = r1.y+dy*j; z = r1.z+dz*k;\r
+ i0 = i+nx*(j+ny*k);\r
+ a[i0] = eqs.Calc(x,y,z,a[i0], v?v->a[i0]:0, w?w->a[i0]:0);\r
+ }\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::ReadHDF4(const char *fname,const char *data)\r
+{
+#ifdef HAVE_HDF4
+ int sd = SDstart(fname,DFACC_READ), nn, i;
+ if(sd==-1) return; // is not a HDF4 file
+ char name[64];
+ SDfileinfo(sd,&nn,&i);
+ for(i=0;i<nn;i++)
+ {
+ int sds, rank, dims[32], type, attr, in[2]={0,0};
+ sds = SDselect(sd,i);
+ SDgetinfo(sds,name,&rank,dims,&type,&attr);
+ if(!strcmp(name,data)) // as I understand there are possible many datas with the same name
+ {
+ if(rank==1) Create(dims[0]);
+ else if(rank==2) Create(dims[1],dims[0]);
+ else if(rank==3) Create(dims[3],dims[1],dims[0]);
+ else continue;
+ if(type==DFNT_FLOAT32)
+ {
+ float *b = new float[nx*ny*nz];
+ SDreaddata(sds,in,0,dims,b);
+ for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];
+ delete []b;
+ }
+ if(type==DFNT_FLOAT64)
+ {
+ double *b = new double[nx*ny*nz];
+ SDreaddata(sds,in,0,dims,b);
+ for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];
+ delete []b;
+ }
+ }
+ SDendaccess(sds);
+ }
+ SDend(sd);
+#endif
+}\r
+//-----------------------------------------------------------------------------\r
+#ifdef HAVE_HDF5\r
+void mglData::SaveHDF(const char *fname,const char *data,bool rewrite) const\r
+{\r
+ hid_t hf,hd,hs;\r
+ hsize_t dims[3];\r
+ long rank = 3, res;\r
+#ifndef H5_USE_16_API\r
+ H5Eset_auto(H5E_DEFAULT,0,0);\r
+#else\r
+ H5Eset_auto(0,0);\r
+#endif\r
+ res=H5Fis_hdf5(fname);\r
+ if(res>0 && !rewrite) hf = H5Fopen(fname, H5F_ACC_RDWR, H5P_DEFAULT);\r
+ else hf = H5Fcreate(fname, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);\r
+ if(hf<0) return;\r
+ if(nz==1 && ny == 1) { rank = 1; dims[0] = nx; }\r
+ else if(nz==1) { rank = 2; dims[0] = ny; dims[1] = nx; }\r
+ else { rank = 3; dims[0] = nz; dims[1] = ny; dims[2] = nx; }\r
+ hs = H5Screate_simple(rank, dims, 0);\r
+#if(MGL_USE_DOUBLE==1)\r
+#ifndef H5_USE_16_API\r
+ hd = H5Dcreate(hf, data, H5T_NATIVE_DOUBLE, hs, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);\r
+#else /* ! HAVE_HDF5_18 */\r
+ hd = H5Dcreate(hf, data, H5T_NATIVE_DOUBLE, hs, H5P_DEFAULT);\r
+#endif /* HAVE_HDF5_18 */\r
+ H5Dwrite(hd, H5T_NATIVE_DOUBLE, hs, hs, H5P_DEFAULT, a);\r
+#else\r
+#ifndef H5_USE_16_API\r
+ hd = H5Dcreate(hf, data, H5T_NATIVE_FLOAT, hs, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);\r
+#else /* ! HAVE_HDF5_18 */\r
+ hd = H5Dcreate(hf, data, H5T_NATIVE_FLOAT, hs, H5P_DEFAULT);\r
+#endif /* HAVE_HDF5_18 */\r
+ H5Dwrite(hd, H5T_NATIVE_FLOAT, hs, hs, H5P_DEFAULT, a);\r
+#endif\r
+ H5Dclose(hd); H5Sclose(hs); H5Fclose(hf);\r
+}\r
+//-----------------------------------------------------------------------------\r
+void mglData::ReadHDF(const char *fname,const char *data)\r
+{\r
+ hid_t hf,hd,hs;\r
+ hsize_t dims[3];\r
+ long rank, res = H5Fis_hdf5(fname);
+ if(res<=0) { ReadHDF4(fname,data); return; }\r
+ hf = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT); if(hf<0) return;\r
+#ifndef H5_USE_16_API\r
+ hd = H5Dopen(hf,data,H5P_DEFAULT);\r
+#else\r
+ hd = H5Dopen(hf,data);\r
+#endif
+ if(hd<0) return;\r
+ hs = H5Dget_space(hd);\r
+ rank = H5Sget_simple_extent_ndims(hs);\r
+ if(rank>0 && rank<=3)\r
+ {\r
+ H5Sget_simple_extent_dims(hs,dims,0);\r
+ nx = ny = nz = 1;\r
+ switch(rank)\r
+ {\r
+ case 1: nx = dims[0]; break;\r
+ case 2: nx = dims[1]; ny = dims[0]; break;\r
+ case 3: nx = dims[2]; ny = dims[1]; nz = dims[0]; break;\r
+ }\r
+ delete []a; a = new mreal[nx*ny*nz]; NewId();\r
+#if(MGL_USE_DOUBLE==1)\r
+ H5Dread(hd, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, a);\r
+#else\r
+ H5Dread(hd, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, a);\r
+#endif\r
+ }\r
+ H5Dclose(hd); H5Sclose(hs); H5Fclose(hf);\r
+}\r
+#else\r
+void mglData::SaveHDF(const char *,const char *,bool ) const {}\r
+void mglData::ReadHDF(const char *,const char *) {}\r
+#endif\r
+//-----------------------------------------------------------------------------\r
+bool mgl_add_file(long &kx,long &ky, long &kz, mreal *&b, mglData &d,bool as_slice)\r
+{\r
+ if(as_slice && d.nz==1)\r
+ {\r
+ if(kx==d.nx && d.ny==1)\r
+ {\r
+ b = (mreal *)realloc(b,kx*(ky+1)*sizeof(mreal));\r
+ memcpy(b+kx*ky,d.a,kx*sizeof(mreal)); ky++;\r
+ }\r
+ else if(kx==d.nx && ky==d.ny)\r
+ {\r
+ b = (mreal *)realloc(b,kx*ky*(kz+1)*sizeof(mreal));\r
+ memcpy(b+kx*ky*kz,d.a,kx*ky*sizeof(mreal)); kz++;\r
+ }\r
+ else return false;\r
+ }\r
+ else\r
+ {\r
+ if(d.ny*d.nz==1 && ky*kz==1)\r
+ {\r
+ b = (mreal *)realloc(b,(kx+d.nx)*sizeof(mreal));\r
+ memcpy(b+kx,d.a,d.nx*sizeof(mreal)); kx+=d.nx;\r
+ }\r
+ else if(kx==d.nx && kz==1 && d.nz==1)\r
+ {\r
+ b = (mreal *)realloc(b,kx*(ky+d.ny)*sizeof(mreal));\r
+ memcpy(b+kx*ky,d.a,kx*d.ny*sizeof(mreal)); ky+=d.ny;\r
+ }\r
+ else if(kx==d.nx && ky==d.ny)\r
+ {\r
+ b = (mreal *)realloc(b,kx*kx*(kz+d.nz)*sizeof(mreal));\r
+ memcpy(b+kx*ky*kz,d.a,kx*ky*d.nz*sizeof(mreal)); kz+=d.nz;\r
+ }\r
+ else return false;\r
+ }\r
+ return true;\r
+}\r
+//-----------------------------------------------------------------------------\r
+bool mglData::ReadRange(const char *templ, mreal from, mreal to, mreal step, bool as_slice)\r
+{\r
+ mglData d;\r
+ mreal t = from, *b;\r
+ long kx,ky,kz;\r
+ char *fname = new char[strlen(templ)+20];\r
+\r
+ //read first file\r
+ do{ sprintf(fname,templ,t); t+= step; } while(!d.Read(fname) && t<=to);\r
+\r
+ if(t>to) return false;\r
+ kx = d.nx; ky = d.ny; kz = d.nz;\r
+ b = (mreal *)malloc(kx*ky*kz*sizeof(mreal));\r
+ memcpy(b,d.a,kx*ky*kz*sizeof(mreal));\r
+\r
+ // read other files\r
+ for(;t<=to;t+=step)\r
+ {\r
+ sprintf(fname,templ,t);\r
+ if(d.Read(fname))\r
+ if(!mgl_add_file(kx,ky,kz,b,d,as_slice))\r
+ return false;\r
+ }\r
+ Set(b,kx,ky,kz);\r
+ delete []fname; free(b);\r
+ return true;\r
+}\r
+//-----------------------------------------------------------------------------\r
+bool mglData::ReadAll(const char *templ, bool as_slice)\r
+{\r
+#ifndef WIN32\r
+ mglData d;\r
+ glob_t res;\r
+ unsigned long i;\r
+ mreal *b;\r
+ long kx,ky,kz;\r
+ char *fname = new char[256];\r
+ glob (templ, GLOB_TILDE, NULL, &res);\r
+\r
+ //read first file\r
+ for(i=0;i<res.gl_pathc;i++)\r
+ if(d.Read(res.gl_pathv[i])) break;\r
+\r
+ if(i>=res.gl_pathc) { delete []fname; return false; }\r
+ kx = d.nx; ky = d.ny; kz = d.nz;\r
+ b = (mreal *)malloc(kx*ky*kz*sizeof(mreal));\r
+ memcpy(b,d.a,kx*ky*kz*sizeof(mreal));\r
+\r
+ for(;i<res.gl_pathc;i++)\r
+ {\r
+ if(d.Read(res.gl_pathv[i]))\r
+ if(!mgl_add_file(kx,ky,kz,b,d,as_slice))
+ { delete []fname; return false; }\r
+ }\r
+ Set(b,kx,ky,kz);\r
+\r
+ globfree (&res);\r
+ delete []fname; free(b);\r
+ return true;\r
+#else\r
+ return false;\r
+#endif\r
+}\r
+//-----------------------------------------------------------------------------\r
+mglData operator*(const mglData &b, const mglData &d)\r
+{ mglData a(b); a*=d; return a; }\r
+mglData operator*(mreal b, const mglData &d)\r
+{ mglData a(d); a*=b; return a; }\r
+mglData operator*(const mglData &d, mreal b)\r
+{ mglData a(d); a*=b; return a; }\r
+mglData operator-(const mglData &b, const mglData &d)\r
+{ mglData a(b); a-=d; return a; }\r
+mglData operator-(mreal b, const mglData &d)\r
+{ mglData a(d); a-=b; return a; }\r
+mglData operator-(const mglData &d, mreal b)\r
+{ mglData a(d); a-=b; return a; }\r
+mglData operator+(const mglData &b, const mglData &d)\r
+{ mglData a(b); a+=d; return a; }\r
+mglData operator+(mreal b, const mglData &d)\r
+{ mglData a(d); a+=b; return a; }\r
+mglData operator+(const mglData &d, mreal b)\r
+{ mglData a(d); a+=b; return a; }\r
+mglData operator/(const mglData &b, const mglData &d)\r
+{ mglData a(b); a/=d; return a; }\r
+mglData operator/(const mglData &d, mreal b)\r
+{ mglData a(d); a/=b; return a; }\r
+void mglData::operator=(mreal v)\r
+{ for(long i=0;i<nx*ny*nz;i++) a[i]=v; }\r
+//-----------------------------------------------------------------------------\r
+void mglData::Set(const std::vector<int> &d)\r
+{ Create(d.size()); for(long i=0;i<nx;i++) a[i] = d[i]; }\r
+void mglData::Set(const std::vector<float> &d)\r
+{ Create(d.size()); for(long i=0;i<nx;i++) a[i] = d[i]; }\r
+void mglData::Set(const std::vector<double> &d)\r
+{ Create(d.size()); for(long i=0;i<nx;i++) a[i] = d[i]; }\r
+//-----------------------------------------------------------------------------\r
+void mglData::NewId()\r
+{ delete []id; id=new char[nx]; memset(id,0,nx*sizeof(char)); }\r
+//-----------------------------------------------------------------------------\r
--- /dev/null
+Upstream patch for FTBFS on i386
+--- a/mgl/mgl_data_io.cpp
++++ b/mgl/mgl_data_io.cpp
+@@ -2,33 +2,33 @@
+ * mgl_data_io.cpp is part of Math Graphic Library\r
+ * Copyright (C) 2007 Alexey Balakin <balakin@appl.sci-nnov.ru> *\r
+ * *\r
+- * This program is free software; you can redistribute it and/or modify *
+- * it under the terms of the GNU Library General Public License as *
+- * published by the Free Software Foundation; either version 3 of the *
+- * License, or (at your option) any later version. *
+- * *
+- * This program is distributed in the hope that it will be useful, *
+- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+- * GNU General Public License for more details. *
+- * *
+- * You should have received a copy of the GNU Library General Public *
+- * License along with this program; if not, write to the *
+- * Free Software Foundation, Inc., *
+- * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+- ***************************************************************************/
++ * This program is free software; you can redistribute it and/or modify *\r
++ * it under the terms of the GNU Library General Public License as *\r
++ * published by the Free Software Foundation; either version 3 of the *\r
++ * License, or (at your option) any later version. *\r
++ * *\r
++ * This program is distributed in the hope that it will be useful, *\r
++ * but WITHOUT ANY WARRANTY; without even the implied warranty of *\r
++ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *\r
++ * GNU General Public License for more details. *\r
++ * *\r
++ * You should have received a copy of the GNU Library General Public *\r
++ * License along with this program; if not, write to the *\r
++ * Free Software Foundation, Inc., *\r
++ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *\r
++ ***************************************************************************/\r
+ #include <stdlib.h>\r
+ #include <ctype.h>\r
+ #include <math.h>\r
+ #include <string.h>\r
+-#include <zlib.h>
++#include <zlib.h>\r
+ #ifdef HAVE_HDF5\r
+ #include <hdf5.h>\r
+ #endif\r
+ #ifdef HAVE_HDF4\r
+-#define intf hdf4_intf
+-#include <hdf/mfhdf.h>
+-#undef intf
++#define intf hdf4_intf\r
++#include <hdf/mfhdf.h>\r
++#undef intf\r
+ #endif\r
+ \r
+ #ifndef WIN32\r
+@@ -302,17 +302,17 @@
+ }\r
+ fclose(fp);\r
+ }\r
+-//-----------------------------------------------------------------------------
+-char *mgl_read_gz(gzFile fp)
+-{
+- long size=1024,n=0,m;
+- char *buf=(char*)malloc(size);
+- while((m=gzread(fp,buf+size*n,size))>0)
+- {
+- if(m<size) { buf[size*n+m]=0; break; }
+- n++; buf=(char*)realloc(buf,size*(n+1));
+- }
+- return buf;
++//-----------------------------------------------------------------------------\r
++char *mgl_read_gz(gzFile fp)\r
++{\r
++ long size=1024,n=0,m;\r
++ char *buf=(char*)malloc(size);\r
++ while((m=gzread(fp,buf+size*n,size))>0)\r
++ {\r
++ if(m<size) { buf[size*n+m]=0; break; }\r
++ n++; buf=(char*)realloc(buf,size*(n+1));\r
++ }\r
++ return buf;\r
+ }\r
+ //-----------------------------------------------------------------------------\r
+ bool mglData::Read(const char *fname)\r
+@@ -321,23 +321,23 @@
+ long nb,i;\r
+ gzFile fp = gzopen(fname,"r");\r
+ if(!fp)\r
+- {
++ {\r
+ if(!a) Create(1,1,1);\r
+ return false;\r
+ }\r
+- char *buf = mgl_read_gz(fp);
++ char *buf = mgl_read_gz(fp);\r
+ nb = strlen(buf); gzclose(fp);\r
+ \r
+ bool first=false,com=false;\r
+ register char ch;\r
+ for(i=nb-1;i>=0;i--) if(buf[i]>' ') break;\r
+- buf[i+1]=0; nb = i; // remove tailing spaces
++ buf[i+1]=0; nb = i; // remove tailing spaces\r
+ for(i=0;i<nb-1 && !isn(buf[i]);i++) // determine nx\r
+ {\r
+ if(buf[i]=='#') while(!isn(buf[i]) && i<nb) i++;\r
+ ch = buf[i];\r
+ if(ch>' ' && !first) first=true;\r
+- if(first && (ch==' ' || ch=='\t') && buf[i+1]>' ') k++;
++ if(first && (ch==' ' || ch=='\t') && buf[i+1]>' ') k++;\r
+ }\r
+ first = false;\r
+ for(i=0;i<nb-1;i++) // determine ny\r
+@@ -362,7 +362,7 @@
+ }\r
+ }\r
+ else for(i=0;i<nb-1;i++) if(buf[i]=='\f') l++;\r
+- free(buf);
++ free(buf);\r
+ return Read(fname,k,m,l);\r
+ }\r
+ //-----------------------------------------------------------------------------\r
+@@ -382,7 +382,7 @@
+ gzFile fp = gzopen(fname,"r");\r
+ if(!fp) return false;\r
+ Create(mx,my,mz);\r
+- char *buf = mgl_read_gz(fp);
++ char *buf = mgl_read_gz(fp);\r
+ long nb = strlen(buf); gzclose(fp);\r
+ \r
+ register long i=0, j=0, k=0;\r
+@@ -408,8 +408,8 @@
+ char *s=buf+j;\r
+ while(buf[j]>' ' && j<nb) j++;\r
+ buf[j]=0;\r
+- a[i] = atof(s);
+- i++; if(i>=nx*ny*nz) break;
++ a[i] = atof(s);\r
++ i++; if(i>=nx*ny*nz) break;\r
+ }\r
+ free(buf);\r
+ return true;\r
+@@ -421,7 +421,7 @@
+ gzFile fp = gzopen(fname,"r");\r
+ if(!fp) return false;\r
+ nx = ny = nz = 1; NewId();\r
+- char *buf = mgl_read_gz(fp);
++ char *buf = mgl_read_gz(fp);\r
+ long nb = strlen(buf); gzclose(fp);\r
+ \r
+ register long i=0,j=0;\r
+@@ -692,10 +692,10 @@
+ {\r
+ for(k=dim;k<nz;k++)\r
+ {\r
+- z = (nz>dim+1) ? (k-dim)/(nz-dim-1.) : 0;
+- aa = a+nx*ny*k;
+-//#pragma omp parallel for
+- for(i=0;i<nx*ny;i++)
++ z = (nz>dim+1) ? (k-dim)/(nz-dim-1.) : 0;\r
++ aa = a+nx*ny*k;\r
++//#pragma omp parallel for\r
++ for(i=0;i<nx*ny;i++)\r
+ aa[i] = eqs.Calc(dx*(i%nx),dy*(i/nx),z,aa[i]);\r
+ }\r
+ }\r
+@@ -703,10 +703,10 @@
+ {\r
+ if(ny==1) dim = 0;\r
+ dy = ny>dim+1 ? 1/(ny-dim-1.) : 0;\r
+- for(j=dim;j<ny;j++)
+- {
+- y = dy*(j-dim); aa = a+nx*j;
+-//#pragma omp parallel for
++ for(j=dim;j<ny;j++)\r
++ {\r
++ y = dy*(j-dim); aa = a+nx*j;\r
++//#pragma omp parallel for\r
+ for(i=0;i<nx;i++) aa[i] = eqs.Calc(dx*i,y,0,aa[i]);\r
+ }\r
+ }\r
+@@ -771,10 +771,10 @@
+ long dx,dy,dz,i1,j1,k1;\r
+ dx = (i+1)*rx<=nx ? rx : nx-i*rx;\r
+ dy = (j+1)*ry<=ny ? ry : ny-j*ry;\r
+- dz = (k+1)*rz<=nz ? rz : nz-k*rz;
++ dz = (k+1)*rz<=nz ? rz : nz-k*rz;\r
+ mreal s = 0;\r
+ for(i1=i*rx;i1<i*rx+dx;i1++) for(j1=j*ry;j1<j*ry+dz;j1++) for(k1=k*rz;k1<k*rz+dz;k1++)\r
+- s += a[i1+nx*(j1+ny*k1)];
++ s += a[i1+nx*(j1+ny*k1)];\r
+ b[i+kx*(j+ky*k)] = s/dx*dy*dz;\r
+ }\r
+ delete []a; a=b;\r
+@@ -912,42 +912,42 @@
+ }\r
+ //-----------------------------------------------------------------------------\r
+ void mglData::ReadHDF4(const char *fname,const char *data)\r
+-{
+-#ifdef HAVE_HDF4
+- int sd = SDstart(fname,DFACC_READ), nn, i;
+- if(sd==-1) return; // is not a HDF4 file
+- char name[64];
+- SDfileinfo(sd,&nn,&i);
+- for(i=0;i<nn;i++)
+- {
+- int sds, rank, dims[32], type, attr, in[2]={0,0};
+- sds = SDselect(sd,i);
+- SDgetinfo(sds,name,&rank,dims,&type,&attr);
+- if(!strcmp(name,data)) // as I understand there are possible many datas with the same name
+- {
+- if(rank==1) Create(dims[0]);
+- else if(rank==2) Create(dims[1],dims[0]);
+- else if(rank==3) Create(dims[3],dims[1],dims[0]);
+- else continue;
+- if(type==DFNT_FLOAT32)
+- {
+- float *b = new float[nx*ny*nz];
+- SDreaddata(sds,in,0,dims,b);
+- for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];
+- delete []b;
+- }
+- if(type==DFNT_FLOAT64)
+- {
+- double *b = new double[nx*ny*nz];
+- SDreaddata(sds,in,0,dims,b);
+- for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];
+- delete []b;
+- }
+- }
+- SDendaccess(sds);
+- }
+- SDend(sd);
+-#endif
++{\r
++#ifdef HAVE_HDF4\r
++ long sd = SDstart(fname,DFACC_READ), nn, i;\r
++ if(sd==-1) return; // is not a HDF4 file\r
++ char name[64];\r
++ SDfileinfo(sd,&nn,&i);\r
++ for(i=0;i<nn;i++)\r
++ {\r
++ long sds, rank, dims[32], type, attr, in[2]={0,0};\r
++ sds = SDselect(sd,i);\r
++ SDgetinfo(sds,name,&rank,dims,&type,&attr);\r
++ if(!strcmp(name,data)) // as I understand there are possible many datas with the same name\r
++ {\r
++ if(rank==1) Create(dims[0]);\r
++ else if(rank==2) Create(dims[1],dims[0]);\r
++ else if(rank==3) Create(dims[3],dims[1],dims[0]);\r
++ else continue;\r
++ if(type==DFNT_FLOAT32)\r
++ {\r
++ float *b = new float[nx*ny*nz];\r
++ SDreaddata(sds,in,0,dims,b);\r
++ for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];\r
++ delete []b;\r
++ }\r
++ if(type==DFNT_FLOAT64)\r
++ {\r
++ double *b = new double[nx*ny*nz];\r
++ SDreaddata(sds,in,0,dims,b);\r
++ for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];\r
++ delete []b;\r
++ }\r
++ }\r
++ SDendaccess(sds);\r
++ }\r
++ SDend(sd);\r
++#endif\r
+ }\r
+ //-----------------------------------------------------------------------------\r
+ #ifdef HAVE_HDF5\r
+@@ -991,14 +991,14 @@
+ {\r
+ hid_t hf,hd,hs;\r
+ hsize_t dims[3];\r
+- long rank, res = H5Fis_hdf5(fname);
++ long rank, res = H5Fis_hdf5(fname);\r
+ if(res<=0) { ReadHDF4(fname,data); return; }\r
+ hf = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT); if(hf<0) return;\r
+ #ifndef H5_USE_16_API\r
+ hd = H5Dopen(hf,data,H5P_DEFAULT);\r
+ #else\r
+ hd = H5Dopen(hf,data);\r
+-#endif
++#endif\r
+ if(hd<0) return;\r
+ hs = H5Dget_space(hd);\r
+ rank = H5Sget_simple_extent_ndims(hs);\r
+@@ -1115,7 +1115,7 @@
+ for(;i<res.gl_pathc;i++)\r
+ {\r
+ if(d.Read(res.gl_pathv[i]))\r
+- if(!mgl_add_file(kx,ky,kz,b,d,as_slice))
++ if(!mgl_add_file(kx,ky,kz,b,d,as_slice))\r
+ { delete []fname; return false; }\r
+ }\r
+ Set(b,kx,ky,kz);\r
* mgl_data_io.cpp is part of Math Graphic Library\r
* Copyright (C) 2007 Alexey Balakin <balakin@appl.sci-nnov.ru> *\r
* *\r
- * This program is free software; you can redistribute it and/or modify *
- * it under the terms of the GNU Library General Public License as *
- * published by the Free Software Foundation; either version 3 of the *
- * License, or (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU Library General Public *
- * License along with this program; if not, write to the *
- * Free Software Foundation, Inc., *
- * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
- ***************************************************************************/
+ * This program is free software; you can redistribute it and/or modify *\r
+ * it under the terms of the GNU Library General Public License as *\r
+ * published by the Free Software Foundation; either version 3 of the *\r
+ * License, or (at your option) any later version. *\r
+ * *\r
+ * This program is distributed in the hope that it will be useful, *\r
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *\r
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *\r
+ * GNU General Public License for more details. *\r
+ * *\r
+ * You should have received a copy of the GNU Library General Public *\r
+ * License along with this program; if not, write to the *\r
+ * Free Software Foundation, Inc., *\r
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *\r
+ ***************************************************************************/\r
#include <stdlib.h>\r
#include <ctype.h>\r
#include <math.h>\r
#include <string.h>\r
-#include <zlib.h>
+#include <zlib.h>\r
#ifdef HAVE_HDF5\r
#include <hdf5.h>\r
#endif\r
#ifdef HAVE_HDF4\r
-#define intf hdf4_intf
-#include <hdf/mfhdf.h>
-#undef intf
+#define intf hdf4_intf\r
+#include <hdf/mfhdf.h>\r
+#undef intf\r
#endif\r
\r
#ifndef WIN32\r
}\r
fclose(fp);\r
}\r
-//-----------------------------------------------------------------------------
-char *mgl_read_gz(gzFile fp)
-{
- long size=1024,n=0,m;
- char *buf=(char*)malloc(size);
- while((m=gzread(fp,buf+size*n,size))>0)
- {
- if(m<size) { buf[size*n+m]=0; break; }
- n++; buf=(char*)realloc(buf,size*(n+1));
- }
- return buf;
+//-----------------------------------------------------------------------------\r
+char *mgl_read_gz(gzFile fp)\r
+{\r
+ long size=1024,n=0,m;\r
+ char *buf=(char*)malloc(size);\r
+ while((m=gzread(fp,buf+size*n,size))>0)\r
+ {\r
+ if(m<size) { buf[size*n+m]=0; break; }\r
+ n++; buf=(char*)realloc(buf,size*(n+1));\r
+ }\r
+ return buf;\r
}\r
//-----------------------------------------------------------------------------\r
bool mglData::Read(const char *fname)\r
long nb,i;\r
gzFile fp = gzopen(fname,"r");\r
if(!fp)\r
- {
+ {\r
if(!a) Create(1,1,1);\r
return false;\r
}\r
- char *buf = mgl_read_gz(fp);
+ char *buf = mgl_read_gz(fp);\r
nb = strlen(buf); gzclose(fp);\r
\r
bool first=false,com=false;\r
register char ch;\r
for(i=nb-1;i>=0;i--) if(buf[i]>' ') break;\r
- buf[i+1]=0; nb = i; // remove tailing spaces
+ buf[i+1]=0; nb = i; // remove tailing spaces\r
for(i=0;i<nb-1 && !isn(buf[i]);i++) // determine nx\r
{\r
if(buf[i]=='#') while(!isn(buf[i]) && i<nb) i++;\r
ch = buf[i];\r
if(ch>' ' && !first) first=true;\r
- if(first && (ch==' ' || ch=='\t') && buf[i+1]>' ') k++;
+ if(first && (ch==' ' || ch=='\t') && buf[i+1]>' ') k++;\r
}\r
first = false;\r
for(i=0;i<nb-1;i++) // determine ny\r
}\r
}\r
else for(i=0;i<nb-1;i++) if(buf[i]=='\f') l++;\r
- free(buf);
+ free(buf);\r
return Read(fname,k,m,l);\r
}\r
//-----------------------------------------------------------------------------\r
gzFile fp = gzopen(fname,"r");\r
if(!fp) return false;\r
Create(mx,my,mz);\r
- char *buf = mgl_read_gz(fp);
+ char *buf = mgl_read_gz(fp);\r
long nb = strlen(buf); gzclose(fp);\r
\r
register long i=0, j=0, k=0;\r
char *s=buf+j;\r
while(buf[j]>' ' && j<nb) j++;\r
buf[j]=0;\r
- a[i] = atof(s);
- i++; if(i>=nx*ny*nz) break;
+ a[i] = atof(s);\r
+ i++; if(i>=nx*ny*nz) break;\r
}\r
free(buf);\r
return true;\r
gzFile fp = gzopen(fname,"r");\r
if(!fp) return false;\r
nx = ny = nz = 1; NewId();\r
- char *buf = mgl_read_gz(fp);
+ char *buf = mgl_read_gz(fp);\r
long nb = strlen(buf); gzclose(fp);\r
\r
register long i=0,j=0;\r
{\r
for(k=dim;k<nz;k++)\r
{\r
- z = (nz>dim+1) ? (k-dim)/(nz-dim-1.) : 0;
- aa = a+nx*ny*k;
-//#pragma omp parallel for
- for(i=0;i<nx*ny;i++)
+ z = (nz>dim+1) ? (k-dim)/(nz-dim-1.) : 0;\r
+ aa = a+nx*ny*k;\r
+//#pragma omp parallel for\r
+ for(i=0;i<nx*ny;i++)\r
aa[i] = eqs.Calc(dx*(i%nx),dy*(i/nx),z,aa[i]);\r
}\r
}\r
{\r
if(ny==1) dim = 0;\r
dy = ny>dim+1 ? 1/(ny-dim-1.) : 0;\r
- for(j=dim;j<ny;j++)
- {
- y = dy*(j-dim); aa = a+nx*j;
-//#pragma omp parallel for
+ for(j=dim;j<ny;j++)\r
+ {\r
+ y = dy*(j-dim); aa = a+nx*j;\r
+//#pragma omp parallel for\r
for(i=0;i<nx;i++) aa[i] = eqs.Calc(dx*i,y,0,aa[i]);\r
}\r
}\r
long dx,dy,dz,i1,j1,k1;\r
dx = (i+1)*rx<=nx ? rx : nx-i*rx;\r
dy = (j+1)*ry<=ny ? ry : ny-j*ry;\r
- dz = (k+1)*rz<=nz ? rz : nz-k*rz;
+ dz = (k+1)*rz<=nz ? rz : nz-k*rz;\r
mreal s = 0;\r
for(i1=i*rx;i1<i*rx+dx;i1++) for(j1=j*ry;j1<j*ry+dz;j1++) for(k1=k*rz;k1<k*rz+dz;k1++)\r
- s += a[i1+nx*(j1+ny*k1)];
+ s += a[i1+nx*(j1+ny*k1)];\r
b[i+kx*(j+ky*k)] = s/dx*dy*dz;\r
}\r
delete []a; a=b;\r
}\r
//-----------------------------------------------------------------------------\r
void mglData::ReadHDF4(const char *fname,const char *data)\r
-{
-#ifdef HAVE_HDF4
- int sd = SDstart(fname,DFACC_READ), nn, i;
- if(sd==-1) return; // is not a HDF4 file
- char name[64];
- SDfileinfo(sd,&nn,&i);
- for(i=0;i<nn;i++)
- {
- int sds, rank, dims[32], type, attr, in[2]={0,0};
- sds = SDselect(sd,i);
- SDgetinfo(sds,name,&rank,dims,&type,&attr);
- if(!strcmp(name,data)) // as I understand there are possible many datas with the same name
- {
- if(rank==1) Create(dims[0]);
- else if(rank==2) Create(dims[1],dims[0]);
- else if(rank==3) Create(dims[3],dims[1],dims[0]);
- else continue;
- if(type==DFNT_FLOAT32)
- {
- float *b = new float[nx*ny*nz];
- SDreaddata(sds,in,0,dims,b);
- for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];
- delete []b;
- }
- if(type==DFNT_FLOAT64)
- {
- double *b = new double[nx*ny*nz];
- SDreaddata(sds,in,0,dims,b);
- for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];
- delete []b;
- }
- }
- SDendaccess(sds);
- }
- SDend(sd);
-#endif
+{\r
+#ifdef HAVE_HDF4\r
+ long sd = SDstart(fname,DFACC_READ), nn, i;\r
+ if(sd==-1) return; // is not a HDF4 file\r
+ char name[64];\r
+ SDfileinfo(sd,&nn,&i);\r
+ for(i=0;i<nn;i++)\r
+ {\r
+ long sds, rank, dims[32], type, attr, in[2]={0,0};\r
+ sds = SDselect(sd,i);\r
+ SDgetinfo(sds,name,&rank,dims,&type,&attr);\r
+ if(!strcmp(name,data)) // as I understand there are possible many datas with the same name\r
+ {\r
+ if(rank==1) Create(dims[0]);\r
+ else if(rank==2) Create(dims[1],dims[0]);\r
+ else if(rank==3) Create(dims[3],dims[1],dims[0]);\r
+ else continue;\r
+ if(type==DFNT_FLOAT32)\r
+ {\r
+ float *b = new float[nx*ny*nz];\r
+ SDreaddata(sds,in,0,dims,b);\r
+ for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];\r
+ delete []b;\r
+ }\r
+ if(type==DFNT_FLOAT64)\r
+ {\r
+ double *b = new double[nx*ny*nz];\r
+ SDreaddata(sds,in,0,dims,b);\r
+ for(long j=0;j<nx*ny*nz;j++) a[j]=b[j];\r
+ delete []b;\r
+ }\r
+ }\r
+ SDendaccess(sds);\r
+ }\r
+ SDend(sd);\r
+#endif\r
}\r
//-----------------------------------------------------------------------------\r
#ifdef HAVE_HDF5\r
{\r
hid_t hf,hd,hs;\r
hsize_t dims[3];\r
- long rank, res = H5Fis_hdf5(fname);
+ long rank, res = H5Fis_hdf5(fname);\r
if(res<=0) { ReadHDF4(fname,data); return; }\r
hf = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT); if(hf<0) return;\r
#ifndef H5_USE_16_API\r
hd = H5Dopen(hf,data,H5P_DEFAULT);\r
#else\r
hd = H5Dopen(hf,data);\r
-#endif
+#endif\r
if(hd<0) return;\r
hs = H5Dget_space(hd);\r
rank = H5Sget_simple_extent_ndims(hs);\r
for(;i<res.gl_pathc;i++)\r
{\r
if(d.Read(res.gl_pathv[i]))\r
- if(!mgl_add_file(kx,ky,kz,b,d,as_slice))
+ if(!mgl_add_file(kx,ky,kz,b,d,as_slice))\r
{ delete []fname; return false; }\r
}\r
Set(b,kx,ky,kz);\r
projects / mathgl.git / commitdiff