summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/runarray.in
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/runarray.in')
-rw-r--r--Build/source/utils/asymptote/runarray.in218
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));