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.in171
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;
}