diff options
Diffstat (limited to 'Build/source/utils/asymptote/runarray.in')
-rw-r--r-- | Build/source/utils/asymptote/runarray.in | 218 |
1 files changed, 82 insertions, 136 deletions
diff --git a/Build/source/utils/asymptote/runarray.in b/Build/source/utils/asymptote/runarray.in index ac7c6ed6781..7e99a875f4b 100644 --- a/Build/source/utils/asymptote/runarray.in +++ b/Build/source/utils/asymptote/runarray.in @@ -13,6 +13,7 @@ Intarray2* => IntArray2() realarray* => realArray() realarray2* => realArray2() pairarray* => pairArray() +pairarray2* => pairArray2() triplearray2* => tripleArray2() callableReal* => realRealFunction() @@ -31,12 +32,17 @@ callableReal* => realRealFunction() using namespace camp; using namespace vm; +namespace run { +extern pair zero; +} + typedef array boolarray; typedef array Intarray; typedef array Intarray2; typedef array realarray; typedef array realarray2; typedef array pairarray; +typedef array pairarray2; typedef array triplearray2; using types::booleanArray; @@ -45,6 +51,7 @@ using types::IntArray2; using types::realArray; using types::realArray2; using types::pairArray; +using types::pairArray2; using types::tripleArray2; typedef callable callableReal; @@ -116,36 +123,6 @@ array *copyArray(array *a) return c; } -inline size_t checkdimension(const array *a, size_t dim) -{ - size_t size=checkArray(a); - if(dim && size != dim) { - ostringstream buf; - buf << "array of length " << dim << " expected"; - error(buf); - } - return size; -} - -double *copyArrayC(const array *a, size_t dim, GCPlacement placement) -{ - size_t size=checkdimension(a,dim); - double *c=(placement == NoGC) ? new double [size] : - new(placement) double[size]; - for(size_t i=0; i < size; i++) - c[i]=read<double>(a,i); - return c; -} - -triple *copyTripleArrayC(const array *a, size_t dim) -{ - size_t size=checkdimension(a,dim); - triple *c=new triple[size]; - for(size_t i=0; i < size; i++) - c[i]=read<triple>(a,i); - return c; -} - array *copyArray2(array *a) { size_t size=checkArray(a); @@ -182,55 +159,6 @@ array *copyArray3(array *a) return c; } -double *copyArray2C(const array *a, bool square, size_t dim2, - 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); - } - - double *c=(placement == NoGC) ? new double [n*m] : new(placement) double[n*m]; - for(size_t i=0; i < n; i++) { - array *ai=read<array*>(a,i); - size_t aisize=checkArray(ai); - if(aisize == m) { - double *ci=c+i*m; - for(size_t j=0; j < m; j++) - ci[j]=read<double>(ai,j); - } else - error(square ? "matrix must be square" : "matrix must be rectangular"); - } - return c; -} - -triple *copyTripleArray2C(const array *a, bool square, size_t dim2) -{ - 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); - } - - triple *c=new triple[n*m]; - for(size_t i=0; i < n; i++) { - array *ai=read<array*>(a,i); - size_t aisize=checkArray(ai); - if(aisize == m) { - triple *ci=c+i*m; - for(size_t j=0; j < m; j++) - ci[j]=read<triple>(ai,j); - } else - error(square ? "matrix must be square" : "matrix must be rectangular"); - } - return c; -} - double *copyTripleArray2Components(array *a, bool square, size_t dim2, GCPlacement placement) { @@ -323,6 +251,43 @@ triple multshiftless(const array& t, const triple& v) (read<real>(t2,0)*x+read<real>(t2,1)*y+read<real>(t2,2)*z)*f); } +template<class T> +array *mult(array *a, array *b) +{ + size_t n=checkArray(a); + + size_t nb=checkArray(b); + size_t na0=n == 0 ? 0 : checkArray(read<array*>(a,0)); + if(na0 != nb) + error(incommensurate); + + size_t nb0=nb == 0 ? 0 : checkArray(read<array*>(b,0)); + + array *c=new array(n); + + T *A,*B; + copyArray2C(A,a,false); + copyArray2C(B,b,false); + + for(size_t i=0; i < n; ++i) { + T *Ai=A+i*nb; + array *ci=new array(nb0); + (*c)[i]=ci; + for(size_t j=0; j < nb0; ++j) { + T sum=0.0; + size_t kj=j; + for(size_t k=0; k < nb; ++k, kj += nb0) + sum += Ai[k]*B[kj]; + (*ci)[j]=sum; + } + } + + delete[] B; + delete[] A; + + return c; +} + double norm(double *a, size_t n) { if(n == 0) return 0.0; @@ -1064,23 +1029,6 @@ realarray2 *identity(Int n) return Identity(n); } -// Return the diagonal matrix with diagonal entries given by a. -realarray2* :diagonal(realarray *a) -{ - size_t n=checkArray(a); - array *c=new array(n); - for(size_t i=0; i < n; ++i) { - array *ci=new array(n); - (*c)[i]=ci; - for(size_t j=0; j < i; ++j) - (*ci)[j]=0.0; - (*ci)[i]=read<real>(a,i); - for(size_t j=i+1; j < n; ++j) - (*ci)[j]=0.0; - } - return c; -} - // Return the inverse of an n x n matrix a using Gauss-Jordan elimination. realarray2 *inverse(realarray2 *a) { @@ -1185,7 +1133,8 @@ realarray *solve(realarray2 *a, realarray *b, bool warn=true) size_t m=checkArray(b); if(m != n) error(incommensurate); - real *A=copyArray2C(a); + real *A; + copyArray2C(A,a); size_t *index=new size_t[n]; if(LUdecompose(A,n,index,warn) == 0) @@ -1193,7 +1142,8 @@ realarray *solve(realarray2 *a, realarray *b, bool warn=true) array *x=new array(n); - real *B=copyArrayC(b); + real *B; + copyArrayC(B,b); for(size_t i=0; i < n; ++i) { size_t ip=index[i]; @@ -1236,8 +1186,9 @@ realarray2 *solve(realarray2 *a, realarray2 *b, bool warn=true) if(checkArray(b) != n) error(incommensurate); size_t m=checkArray(read<array*>(b,0)); - real *A=copyArray2C(a); - real *B=copyArray2C(b,false); + real *A,*B; + copyArray2C(A,a); + copyArray2C(B,b,false); size_t *index=new size_t[n]; @@ -1291,7 +1242,8 @@ realarray2 *solve(realarray2 *a, realarray2 *b, bool warn=true) // Compute the determinant of an n x n matrix. real determinant(realarray2 *a) { - real *A=copyArray2C(a); + real *A; + copyArray2C(A,a); size_t n=checkArray(a); real det=LUdecompose(A,n,NULL,false); @@ -1309,7 +1261,8 @@ realarray *Operator *(realarray2 *a, realarray *b) size_t n=checkArray(a); size_t m=checkArray(b); array *c=new array(n); - real *B=copyArrayC(b); + real *B; + copyArrayC(B,b); for(size_t i=0; i < n; ++i) { array *ai=read<array*>(a,i); if(checkArray(ai) != m) error(incommensurate); @@ -1326,7 +1279,8 @@ realarray *Operator *(realarray *a, realarray2 *b) { size_t n=checkArray(a); if(n != checkArray(b)) error(incommensurate); - real *A=copyArrayC(a); + real *A; + copyArrayC(A,a); array **B=new array*[n]; array *bk=read<array *>(b,0); @@ -1350,39 +1304,19 @@ realarray *Operator *(realarray *a, realarray2 *b) return c; } +Intarray2 *Operator *(Intarray2 *a, Intarray2 *b) +{ + return mult<Int>(a,b); +} + realarray2 *Operator *(realarray2 *a, realarray2 *b) { - size_t n=checkArray(a); - - size_t nb=checkArray(b); - size_t na0=n == 0 ? 0 : checkArray(read<array*>(a,0)); - if(na0 != nb) - error(incommensurate); - - size_t nb0=nb == 0 ? 0 : checkArray(read<array*>(b,0)); - - array *c=new array(n); - - real *A=copyArray2C(a,false); - real *B=copyArray2C(b,false); - - for(size_t i=0; i < n; ++i) { - real *Ai=A+i*nb; - array *ci=new array(nb0); - (*c)[i]=ci; - for(size_t j=0; j < nb0; ++j) { - real sum=0.0; - size_t kj=j; - for(size_t k=0; k < nb; ++k, kj += nb0) - sum += Ai[k]*B[kj]; - (*ci)[j]=sum; - } - } - - delete[] B; - delete[] A; - - return c; + return mult<real>(a,b); +} + +pairarray2 *Operator *(pairarray2 *a, pairarray2 *b) +{ + return mult<pair>(a,b); } triple Operator *(realarray2 *t, triple v) @@ -1425,6 +1359,16 @@ real dot(realarray *a, realarray *b) return sum; } +// Compute the complex dot product of vectors a and b. +pair dot(pairarray *a, pairarray *b) +{ + size_t n=checkArrays(a,b); + pair sum=zero; + for(size_t i=0; i < n; ++i) + sum += read<pair>(a,i)*conj(read<pair>(b,i)); + return sum; +} + // Solve the problem L\inv f, where f is an n vector and L is the n x n matrix // // [ b[0] c[0] a[0] ] @@ -1821,7 +1765,8 @@ triple maxbezier(triplearray2 *P, triple b) pair minratio(triplearray2 *P, pair b) { - triple *A=copyTripleArray2C(P,true,4); + 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)); @@ -1831,7 +1776,8 @@ pair minratio(triplearray2 *P, pair b) pair maxratio(triplearray2 *P, pair b) { - triple *A=copyTripleArray2C(P,true,4); + 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)); |