mat64: implement triangle shadow detection

mat64/shadow.go

@@ -4,7 +4,10 @@
 
 package mat64
 
-import "github.com/gonum/blas/blas64"
+import (
+	"github.com/gonum/blas"
+	"github.com/gonum/blas/blas64"
+)
 
 const (
 	// regionOverlap is the panic string used for the general case
@@ -66,10 +69,6 @@ func (m *Dense) checkOverlap(a blas64.General) bool {
 	return false
 }
 
-// BUG(kortschak): Overlap detection for symmetric and triangular matrices is not
-// precise; currently overlap is detected if the bounding rectangles overlap rather
-// than exact overlap between visible elements.
-
 func (s *SymDense) checkOverlap(a blas64.Symmetric) bool {
 	mat := s.RawSymmetric()
 	if cap(mat.Data) == 0 || cap(a.Data) == 0 {
@@ -100,12 +99,11 @@ func (s *SymDense) checkOverlap(a blas64.Symmetric) bool {
 		panic(mismatchedStrides)
 	}
 
-	// TODO(kortschak) Make this analysis more precise.
 	if off > 0 {
 		off = -off
 		mat.N, a.N = a.N, mat.N
 	}
-	if rectanglesOverlap(off, mat.N, a.N, mat.Stride) {
+	if trianglesOverlap(off, mat.N, a.N, mat.Stride, mat.Uplo == blas.Upper, a.Uplo == blas.Upper) {
 		panic(regionOverlap)
 	}
 	return false
@@ -141,12 +139,11 @@ func (t *TriDense) checkOverlap(a blas64.Triangular) bool {
 		panic(mismatchedStrides)
 	}
 
-	// TODO(kortschak) Make this analysis more precise.
 	if off > 0 {
 		off = -off
 		mat.N, a.N = a.N, mat.N
 	}
-	if rectanglesOverlap(off, mat.N, a.N, mat.Stride) {
+	if trianglesOverlap(off, mat.N, a.N, mat.Stride, mat.Uplo == blas.Upper, a.Uplo == blas.Upper) {
 		panic(regionOverlap)
 	}
 	return false
@@ -223,3 +220,55 @@ func rectanglesOverlap(off, aCols, bCols, stride int) bool {
 	// b strictly wraps and so must overlap with a.
 	return true
 }
+
+// trianglesOverlap returns whether the strided triangles a and b overlap
+// when b is offset by off elements after a but has at least one element before
+// the end of a. a and b have aSize and bSize respectively.
+func trianglesOverlap(off, aSize, bSize, stride int, aUpper, bUpper bool) bool {
+	if !rectanglesOverlap(off, aSize, bSize, stride) {
+		// Fast return if bounding rectangles do not overlap.
+		return false
+	}
+
+	if aUpper {
+		// Test overlap by projecting diagonals from the b upper
+		// left corner, and upper right or lower left corner onto
+		// the top row of a.
+		rows := off / stride
+		proj := off - rows*(stride+1)
+		if proj >= 0 && proj < aSize {
+			// Top left corner of b projects to element in top
+			// row of a.
+			return true
+		}
+		if bUpper {
+			// Project from upper right corner of b.
+			proj += bSize
+		} else {
+			// Project from lower left corner of b.
+			proj -= bSize
+		}
+		// Return whether corner projects to element in top row of a.
+		return proj >= 0 && proj < aSize
+	}
+
+	// Test overlap by projecting diagonals from the b upper
+	// left corner, and upper right or lower left corner onto
+	// the left column of a.
+	cols := off % stride
+	proj := off - cols*(stride+1)
+	if proj >= 0 && proj < aSize*stride {
+		// Top left corner of b projects to element in left
+		// column of a.
+		return true
+	}
+	if bUpper {
+		// Project from upper right corner of b.
+		proj -= bSize * stride
+	} else {
+		// Project from lower left corner of b.
+		proj += bSize * stride
+	}
+	// Return whether corner projects to element in left column of a.
+	return proj >= 0 && proj < aSize*stride
+}