[Octave-bug-tracker] [bug #60818] delaunayn - 2D code path vectorization doesn't match nD algorithm

[anonymous]

Follow-up Comment #1, bug #60818 (project octave):

I have worked out a way to speed up the check for sliver simplexes without
using loops for any number of dimensions. It appears to be faster than the
existing method for all but the 2d case and can remove the same simplexes. The
check is as follows:


edgvec=pts(T(:,2:end).'(:),:)-kron(pts(T(:,1),:),ones(nd-1,1));
eqs=sparse((nd-1)*nt,(nd-1)*nt);
eqs(logical(kron(speye(nt,nt),true(nd-1))))=edgvec.'(:);
[l u p q]=lu(eqs,"vector");
R=abs(diag(u));
reorderdtriidx=kron(1:nt,ones(1,nd-1))(p)(q);
idx=unique(reorderdtriidx(R<100*(nd-1)*nt*eps(max(R))));


The check is as follows:
* Calculate the edge vectors of the simplex

edgvec=pts(T(:,2:end).'(:),:)-kron(pts(T(:,1),:),ones(nd-1,1));

* Place the edge vectors into a block diagonal matrix where each block is for
each simplex

eqs=sparse((nd-1)*nt,(nd-1)*nt);
eqs(logical(kron(speye(nt,nt),true(nd-1))))=edgvec.'(:);

* Find the lu factorization and extract the diagonal, the operations used in
the lu factorization do not change the determinant (other than the sign and we
only care about the absolute value), and the product of u in the diagonal of
the lu factorization equals the determinant.

[l u p q]=lu(eqs,"vector");
R=abs(diag(u));

* The lu factorization reorders the matrix so we can reorder the indexes to
the simplex blocks to match

reorderdtriidx=kron(1:nt,ones(1,nd-1))(p)(q);

* Perform the check, the 100 is arbitrary and eps(max(R)) could be replaced by
a matrix norm. With 

idx=unique(reorderdtriidx(R<100*(nd-1)*nt*eps(max(R))));

using

x=rand(10,3);
x=[x;x-1000*eps];
delaunayn(x);

I get the same indexes as the existing method with the 100*, without it I get
a subset of the indexes.

This method is similar to a method I use to find which rows of a matrix are
causing it to be singular. I normally use a QR factorization but the QR
factorization was significantly slower than the existing method.

If you want to use the existing method the determinant of a particular
triangle (the first on in this case) is 

prod(R(reorderdtriidx==1))

or for all determinants in a vector

absdet=prod(reshape(kron(R,ones(1,nt))(reorderdtriidx.'==(1:nt)),nd-1,nt),1)

where R=abs(diag(u)) from the lu factorization.

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60818>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/