From f1e79c3669f68804c471699badc56a4a33a23246 Mon Sep 17 00:00:00 2001 From: Norbert Preining Date: Mon, 8 Aug 2022 03:05:28 +0000 Subject: CTAN sync 202208080305 --- graphics/asymptote/runarray.in | 116 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) (limited to 'graphics/asymptote/runarray.in') diff --git a/graphics/asymptote/runarray.in b/graphics/asymptote/runarray.in index 5e45cfd4bf..fdc28f2d3c 100644 --- a/graphics/asymptote/runarray.in +++ b/graphics/asymptote/runarray.in @@ -12,6 +12,7 @@ Intarray* => IntArray() Intarray2* => IntArray2() realarray* => realArray() realarray2* => realArray2() +realarray3* => realArray3() pairarray* => pairArray() pairarray2* => pairArray2() pairarray3* => pairArray3() @@ -34,6 +35,19 @@ static const char *installFFTW= "Please install fftw3, then ./configure; make"; #endif +#ifdef HAVE_EIGEN_DENSE +#include +typedef std::complex Complex; +static const char *square="matrix must be square"; +using Eigen::MatrixXd; +using Eigen::MatrixXcd; +using Eigen::RealSchur; +using Eigen::ComplexSchur; +#else +static const char *installEIGEN= + "Please install eigen3, then ./configure; make"; +#endif + using namespace camp; using namespace vm; @@ -46,6 +60,7 @@ typedef array Intarray; typedef array Intarray2; typedef array realarray; typedef array realarray2; +typedef array realarray3; typedef array pairarray; typedef array pairarray2; typedef array pairarray3; @@ -56,6 +71,7 @@ using types::IntArray; using types::IntArray2; using types::realArray; using types::realArray2; +using types::realArray3; using types::pairArray; using types::pairArray2; using types::pairArray3; @@ -1964,6 +1980,106 @@ pairarray3* fft(pairarray3 *a, Int sign=1) return c; } +// Compute the real Schur decomposition of a 2D pair array +realarray3* _schur(realarray2 *a) +{ +#ifdef HAVE_EIGEN_DENSE + size_t n=checkArray(a); + + MatrixXd A=MatrixXd::Random(n,n); + RealSchur schur(n); + + array *S=new array(2); + + if(n) { + for(size_t i=0; i < n; ++i) { + array *ai=read(a,i); + size_t aisize=checkArray(ai); + if(aisize != n) error(square); + for(size_t j=0; j < n; ++j) + A(i,j)=read(ai,j); + } + + schur.compute(A); + MatrixXd U=schur.matrixU(); + MatrixXd T=schur.matrixT(); + + array *u=new array(n); + array *t=new array(n); + (*S)[0]=u; + (*S)[1]=t; + + for(size_t i=0; i < n; ++i) { + array *ui=new array(n); + array *ti=new array(n); + (*u)[i]=ui; + (*t)[i]=ti; + for(size_t j=0; j < n; ++j) { + (*ui)[j]=U(i,j); + (*ti)[j]=T(i,j); + } + } + } +#else + unused(a); + array *S=new array(0); + error(installEIGEN); +#endif // HAVE_EIGEN_DENSE + return S; +} + +// Compute the Schur decomposition of a 2D pair array +pairarray3* _schur(pairarray2 *a) +{ +#ifdef HAVE_EIGEN_DENSE + size_t n=checkArray(a); + + MatrixXcd A=MatrixXcd::Random(n,n); + ComplexSchur schur(n); + + array *S=new array(2); + + if(n) { + for(size_t i=0; i < n; ++i) { + array *ai=read(a,i); + size_t aisize=checkArray(ai); + if(aisize != n) error(square); + for(size_t j=0; j < n; ++j) { + pair z=read(ai,j); + A(i,j)=Complex(z.getx(),z.gety()); + } + } + + schur.compute(A); + MatrixXcd U=schur.matrixU(); + MatrixXcd T=schur.matrixT(); + + array *u=new array(n); + array *t=new array(n); + (*S)[0]=u; + (*S)[1]=t; + + for(size_t i=0; i < n; ++i) { + array *ui=new array(n); + array *ti=new array(n); + (*u)[i]=ui; + (*t)[i]=ti; + for(size_t j=0; j < n; ++j) { + Complex z=U(i,j); + Complex w=T(i,j); + (*ui)[j]=pair(z.real(),z.imag()); + (*ti)[j]=pair(w.real(),w.imag()); + } + } + } +#else + unused(a); + array *S=new array(0); + error(installEIGEN); +#endif // HAVE_EIGEN_DENSE + return S; +} + Intarray2 *triangulate(pairarray *z) { size_t nv=checkArray(z); -- cgit v1.2.3