testlapack: add test for Dtbtrs

lapack/testlapack/general.go

@@ -1766,3 +1766,136 @@ func dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, m, n int, a
 		panic("not implemented")
 	}
 }
+
+// dlantb is a local implementation of Dlantb to keep code paths independent.
+func dlantb(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n, k int, a []float64, lda int, work []float64) float64 {
+	if n == 0 {
+		return 0
+	}
+	var value float64
+	switch norm {
+	case lapack.MaxAbs:
+		if uplo == blas.Upper {
+			var jfirst int
+			if diag == blas.Unit {
+				value = 1
+				jfirst = 1
+			}
+			for i := 0; i < n; i++ {
+				for _, aij := range a[i*lda+jfirst : i*lda+min(n-i, k+1)] {
+					if math.IsNaN(aij) {
+						return aij
+					}
+					aij = math.Abs(aij)
+					if aij > value {
+						value = aij
+					}
+				}
+			}
+		} else {
+			jlast := k + 1
+			if diag == blas.Unit {
+				value = 1
+				jlast = k
+			}
+			for i := 0; i < n; i++ {
+				for _, aij := range a[i*lda+max(0, k-i) : i*lda+jlast] {
+					if math.IsNaN(aij) {
+						return math.NaN()
+					}
+					aij = math.Abs(aij)
+					if aij > value {
+						value = aij
+					}
+				}
+			}
+		}
+	case lapack.MaxRowSum:
+		var sum float64
+		if uplo == blas.Upper {
+			var jfirst int
+			if diag == blas.Unit {
+				jfirst = 1
+			}
+			for i := 0; i < n; i++ {
+				sum = 0
+				if diag == blas.Unit {
+					sum = 1
+				}
+				for _, aij := range a[i*lda+jfirst : i*lda+min(n-i, k+1)] {
+					sum += math.Abs(aij)
+				}
+				if math.IsNaN(sum) {
+					return math.NaN()
+				}
+				if sum > value {
+					value = sum
+				}
+			}
+		} else {
+			jlast := k + 1
+			if diag == blas.Unit {
+				jlast = k
+			}
+			for i := 0; i < n; i++ {
+				sum = 0
+				if diag == blas.Unit {
+					sum = 1
+				}
+				for _, aij := range a[i*lda+max(0, k-i) : i*lda+jlast] {
+					sum += math.Abs(aij)
+				}
+				if math.IsNaN(sum) {
+					return math.NaN()
+				}
+				if sum > value {
+					value = sum
+				}
+			}
+		}
+	case lapack.MaxColumnSum:
+		work = work[:n]
+		if diag == blas.Unit {
+			for i := range work {
+				work[i] = 1
+			}
+		} else {
+			for i := range work {
+				work[i] = 0
+			}
+		}
+		if uplo == blas.Upper {
+			var jfirst int
+			if diag == blas.Unit {
+				jfirst = 1
+			}
+			for i := 0; i < n; i++ {
+				for j, aij := range a[i*lda+jfirst : i*lda+min(n-i, k+1)] {
+					work[i+jfirst+j] += math.Abs(aij)
+				}
+			}
+		} else {
+			jlast := k + 1
+			if diag == blas.Unit {
+				jlast = k
+			}
+			for i := 0; i < n; i++ {
+				off := max(0, k-i)
+				for j, aij := range a[i*lda+off : i*lda+jlast] {
+					work[i+j+off-k] += math.Abs(aij)
+				}
+			}
+		}
+		for _, wi := range work {
+			if math.IsNaN(wi) {
+				return math.NaN()
+			}
+			if wi > value {
+				value = wi
+			}
+		}
+	case lapack.Frobenius:
+		panic("not implemented")
+	}
+	return value
+}