diff options
Diffstat (limited to 'Build/source/utils/asymptote/runarray.in')
-rw-r--r-- | Build/source/utils/asymptote/runarray.in | 171 |
1 files changed, 117 insertions, 54 deletions
diff --git a/Build/source/utils/asymptote/runarray.in b/Build/source/utils/asymptote/runarray.in index ca2f0bc6d39..24859127bfa 100644 --- a/Build/source/utils/asymptote/runarray.in +++ b/Build/source/utils/asymptote/runarray.in @@ -110,8 +110,29 @@ array *Identity(Int n) static const char *incommensurate="Incommensurate matrices"; static const char *singular="Singular matrix"; +static const char *invalidarraylength="Invalid array length: "; static size_t *pivot,*Row,*Col; +bound_double *bounddouble(int N) +{ + if(N == 16) return bound; + if(N == 10) return boundtri; + ostringstream buf; + buf << invalidarraylength << " " << N; + error(buf); + return NULL; +} + +bound_triple *boundtriple(int N) +{ + if(N == 16) return bound; + if(N == 10) return boundtri; + ostringstream buf; + buf << invalidarraylength << " " << N; + error(buf); + return NULL; +} + static inline void inverseAllocate(size_t n) { pivot=new size_t[n]; @@ -152,39 +173,49 @@ array *copyArray2(array *a) return c; } -double *copyTripleArray2Components(array *a, bool square, size_t dim2, - GCPlacement placement) +double *copyTripleArray2Components(array *a, size_t &N, GCPlacement placement) { size_t n=checkArray(a); - size_t m=(square || n == 0) ? n : checkArray(read<array*>(a,0)); - if(n > 0 && dim2 && m != dim2) { - ostringstream buf; - buf << "second matrix dimension must be " << dim2; - error(buf); + N=0; + for(size_t i=0; i < n; i++) + N += checkArray(read<array*>(a,i)); + + double *A=(placement == NoGC) ? new double [3*N] : + new(placement) double[3*N]; + double *p=A; + + for(size_t i=0; i < n; i++) { + array *ai=read<array*>(a,i); + size_t m=checkArray(ai); + for(size_t j=0; j < m; j++) { + triple v=read<triple>(ai,j); + *p=v.getx(); + *(p+N)=v.gety(); + *(p+2*N)=v.getz(); + ++p; + } } + return A; +} + +triple *copyTripleArray2C(array *a, size_t &N, GCPlacement placement) +{ + size_t n=checkArray(a); + N=0; + for(size_t i=0; i < n; i++) + N += checkArray(read<array*>(a,i)); + + triple *A=(placement == NoGC) ? new triple [N] : + new(placement) triple[N]; + triple *p=A; - size_t nm=n*m; - double *cx=(placement == NoGC) ? new double [3*nm] : - new(placement) double[3*nm]; - double *cy=cx+nm; - double *cz=cx+2*nm; for(size_t i=0; i < n; i++) { array *ai=read<array*>(a,i); - size_t aisize=checkArray(ai); - if(aisize == m) { - double *xi=cx+i*m; - double *yi=cy+i*m; - double *zi=cz+i*m; - for(size_t j=0; j < m; j++) { - triple v=read<triple>(ai,j); - xi[j]=v.getx(); - yi[j]=v.gety(); - zi[j]=v.getz(); - } - } else - error(square ? "matrix must be square" : "matrix must be rectangular"); + size_t m=checkArray(ai); + for(size_t j=0; j < m; j++) + *(p++)=read<triple>(ai,j); } - return cx; + return A; } triple operator *(const array& t, const triple& v) @@ -254,6 +285,35 @@ array *mult(array *a, array *b) return c; } +// Compute transpose(A)*A where A is an n x m matrix. +template<class T> +array *AtA(array *a) +{ + size_t n=checkArray(a); + size_t m=n == 0 ? 0 : checkArray(read<array*>(a,0)); + + array *c=new array(m); + + T *A; + copyArray2C(A,a,false); + + for(size_t i=0; i < m; ++i) { + array *ci=new array(m); + (*c)[i]=ci; + for(size_t j=0; j < m; ++j) { + T sum=T(); + size_t kj=j; + size_t ki=i; + for(size_t k=0; k < n; ++k, kj += m, ki += m) + sum += A[ki]*A[kj]; + (*ci)[j]=sum; + } + } + + delete[] A; + return c; +} + double norm(double *a, size_t n) { if(n == 0) return 0.0; @@ -392,14 +452,6 @@ bool compareFunction(const vm::item& i, const vm::item& j) return pop<bool>(FuncStack); } -void checkSquare(array *a) -{ - size_t n=checkArray(a); - for(size_t i=0; i < n; i++) - if(checkArray(read<array*>(a,i)) != n) - error("matrix a must be square"); -} - // Crout's algorithm for computing the LU decomposition of a square matrix. // cf. routine ludcmp (Press et al., Numerical Recipes, 1991). Int LUdecompose(double *a, size_t n, size_t* index, bool warn=true) @@ -1316,6 +1368,11 @@ triple Operator *(realarray2 *t, triple v) return *t*v; } +realarray2 *AtA(realarray2 *a) +{ + return AtA<real>(a); +} + pair project(triple v, realarray2 *t) { size_t n=checkArray(t); @@ -1615,7 +1672,7 @@ pairarray* fft(pairarray *a, Int sign=1) unsigned n=(unsigned) checkArray(a); array *c=new array(n); if(n) { - Complex *f=fftwpp::ComplexAlign(n); + Complex *f=utils::ComplexAlign(n); fftwpp::fft1d Forward(n,intcast(sign),f); for(size_t i=0; i < n; i++) { @@ -1628,7 +1685,7 @@ pairarray* fft(pairarray *a, Int sign=1) Complex z=f[i]; (*c)[i]=pair(z.real(),z.imag()); } - fftwpp::deleteAlign(f); + utils::deleteAlign(f); } #else unused(&sign); @@ -1738,42 +1795,48 @@ real change2(triplearray2 *a) triple minbezier(triplearray2 *P, triple b) { - real *A=copyTripleArray2Components(P,true,4); - b=triple(bound(A,::min,b.getx(),sqrtFuzz*norm(A,16)), - bound(A+16,::min,b.gety(),sqrtFuzz*norm(A+16,16)), - bound(A+32,::min,b.getz(),sqrtFuzz*norm(A+32,16))); + size_t N; + real *A=copyTripleArray2Components(P,N); + bound_double *B=bounddouble(N); + b=triple(B(A,::min,b.getx(),sqrtFuzz*norm(A,N),maxdepth), + B(A+N,::min,b.gety(),sqrtFuzz*norm(A+N,N),maxdepth), + B(A+2*N,::min,b.getz(),sqrtFuzz*norm(A+2*N,N),maxdepth)); delete[] A; return b; } triple maxbezier(triplearray2 *P, triple b) { - real *A=copyTripleArray2Components(P,true,4); - b=triple(bound(A,::max,b.getx(),sqrtFuzz*norm(A,16)), - bound(A+16,::max,b.gety(),sqrtFuzz*norm(A+16,16)), - bound(A+32,::max,b.getz(),sqrtFuzz*norm(A+32,16))); + size_t N; + real *A=copyTripleArray2Components(P,N); + bound_double *B=bounddouble(N); + b=triple(B(A,::max,b.getx(),sqrtFuzz*norm(A,N),maxdepth), + B(A+N,::max,b.gety(),sqrtFuzz*norm(A+N,N),maxdepth), + B(A+2*N,::max,b.getz(),sqrtFuzz*norm(A+2*N,N),maxdepth)); delete[] A; return b; } pair minratio(triplearray2 *P, pair b) { - triple *A; - copyArray2C(A,P,true,4); - real fuzz=sqrtFuzz*norm(A,16); - b=pair(bound(A,::min,xratio,b.getx(),fuzz), - bound(A,::min,yratio,b.gety(),fuzz)); + size_t N; + triple *A=copyTripleArray2C(P,N); + real fuzz=sqrtFuzz*norm(A,N); + bound_triple *B=boundtriple(N); + b=pair(B(A,::min,xratio,b.getx(),fuzz,maxdepth), + B(A,::min,yratio,b.gety(),fuzz,maxdepth)); delete[] A; return b; } pair maxratio(triplearray2 *P, pair b) { - triple *A; - copyArray2C(A,P,true,4); - real fuzz=sqrtFuzz*norm(A,16); - b=pair(bound(A,::max,xratio,b.getx(),fuzz), - bound(A,::max,yratio,b.gety(),fuzz)); + size_t N; + triple *A=copyTripleArray2C(P,N); + bound_triple *B=boundtriple(N); + real fuzz=sqrtFuzz*norm(A,N); + b=pair(B(A,::max,xratio,b.getx(),fuzz,maxdepth), + B(A,::max,yratio,b.gety(),fuzz,maxdepth)); delete[] A; return b; } |