refactor: Update conversions.v, dgetf2.v, dsyev.v, and lapack_notd_vs…

…l_lapack_lapacke.v

- Add functions uplo_from_bool and uplo_to_bool to conversions.v
- Remove unused dlamch_s function from dgetf2.v
- Update dsyev.v to include additional error handling and scaling of matrix
- Update dpotrf function in lapack_notd_vsl_lapack_lapacke.v to use uplo_from_bool for uplo parameter

blas/conversions.v

@@ -21,6 +21,16 @@ pub type Diagonal = blas64.Diagonal
 // Side is used to specify whether a matrix is on the left or right side in a matrix-matrix multiplication.
 pub type Side = blas64.Side
 
+// uplo_from_bool converts a boolean to Uplo.
+pub fn uplo_from_bool(uplo bool) Uplo {
+	return if uplo { .upper } else { .lower }
+}
+
+// uplo_to_bool converts Uplo to a boolean.
+pub fn uplo_to_bool(uplo Uplo) bool {
+	return uplo == .upper
+}
+
 // slice_to_col_major converts nested slice into an array representing a col-major matrix
 //
 // _**NOTE**: make sure to have at least 1x1 item_

lapack/lapack64/dgetf2.v

@@ -24,7 +24,7 @@ pub fn dgetf2(m int, n int, mut a []f64, lda int, mut ipiv []int) {
 		panic(bad_len_ipiv)
 	}
 
-	sfmin := dlamch_s()
+	sfmin := dlamch_s
 
 	for j := 0; j < mn; j++ {
 		// Find a pivot and test for singularity.
@@ -58,9 +58,3 @@ pub fn dgetf2(m int, n int, mut a []f64, lda int, mut ipiv []int) {
 		}
 	}
 }
-
-fn dlamch_s() f64 {
-	// Returns the safe minimum value (sfmin).
-	// This value is used as a threshold for detecting small values in the matrix.
-	return math.ldexp(1.0, -1022) // Smallest positive normal number.
-}

lapack/lapack64/dlansy.v

@@ -0,0 +1,122 @@
+module lapack64
+
+import math
+import vsl.blas
+
+// dlansy returns the value of the specified norm of an n×n symmetric matrix.
+// If norm == MatrixNorm.max_column_sum or norm == MatrixNorm.max_row_sum, work must have length
+// at least n, otherwise work is unused.
+pub fn dlansy(norm MatrixNorm, uplo blas.Uplo, n int, a []f64, lda int, mut work []f64) f64 {
+	if norm != .max_row_sum && norm != .max_column_sum && norm != .frobenius && norm != .max_abs {
+		panic(lapack64.bad_norm)
+	}
+	if uplo != .upper && uplo != .lower {
+		panic(lapack64.bad_uplo)
+	}
+	if n < 0 {
+		panic('lapack: n < 0')
+	}
+	if lda < math.max(1, n) {
+		panic(lapack64.bad_ld_a)
+	}
+
+	// Quick return if possible.
+	if n == 0 {
+		return 0.0
+	}
+
+	if a.len < (n-1) * lda + n {
+		panic(lapack64.short_a)
+	}
+	if (norm == .max_column_sum || norm == .max_row_sum) && work.len < n {
+		panic(lapack64.short_work)
+	}
+
+	match norm {
+		.max_abs {
+			if uplo == .upper {
+				mut max := 0.0
+				for i in 0 .. n {
+					for j in i .. n {
+						v := math.abs(a[i * lda + j])
+						if math.is_nan(v) {
+							return math.nan()
+						}
+						if v > max {
+							max = v
+						}
+					}
+				}
+				return max
+			}
+			mut max := 0.0
+			for i in 0 .. n {
+				for j in 0 .. i + 1 {
+					v := math.abs(a[i * lda + j])
+					if math.is_nan(v) {
+						return math.nan()
+					}
+					if v > max {
+						max = v
+					}
+				}
+			}
+			return max
+		}
+		.max_row_sum, .max_column_sum {
+			// A symmetric matrix has the same 1-norm and ∞-norm.
+			for i in 0 .. n {
+				work[i] = 0.0
+			}
+			if uplo == .upper {
+				for i in 0 .. n {
+					work[i] += math.abs(a[i * lda + i])
+					for j in i + 1 .. n {
+						v := math.abs(a[i * lda + j])
+						work[i] += v
+						work[j] += v
+					}
+				}
+			} else {
+				for i in 0 .. n {
+					for j in 0 .. i {
+						v := math.abs(a[i * lda + j])
+						work[i] += v
+						work[j] += v
+					}
+					work[i] += math.abs(a[i * lda + i])
+				}
+			}
+			mut max := 0.0
+			for i in 0 .. n {
+				v := work[i]
+				if math.is_nan(v) {
+					return math.nan()
+				}
+				if v > max {
+					max = v
+				}
+			}
+			return max
+		}
+		else {
+			// blas.frobenius:
+			mut scale := 0.0
+			mut sum := 1.0
+			// Sum off-diagonals.
+			if uplo == .upper {
+				for i in 0 .. n - 1 {
+					scale, sum = dlassq(n - i - 1, a[i * lda + i + 1..], 1, scale, sum)
+				}
+			} else {
+				for i in 1 .. n {
+					scale, sum = dlassq(i, a[i * lda..], 1, scale, sum)
+				}
+			}
+			sum *= 2.0
+			// Sum diagonal.
+			scale, sum = dlassq(n, a, lda + 1, scale, sum)
+			return scale * math.sqrt(sum)
+		}
+	}
+}

lapack/lapack64/dlassq.v

@@ -0,0 +1,120 @@
+module lapack64
+
+import math
+
+// dlassq updates a sum of squares represented in scaled form. It returns
+// the values scl and smsq such that
+//
+// 	scl^2*smsq = X[0]^2 + ... + X[n-1]^2 + scale^2*sumsq
+//
+// The value of sumsq is assumed to be non-negative.
+pub fn dlassq(n int, x []f64, incx int, scale f64, sumsq f64) (f64, f64) {
+	if n < 0 {
+		panic('lapack: n < 0')
+	}
+	if incx <= 0 {
+		panic('lapack: increment not one or negative one')
+	}
+	if x.len < 1 + (n - 1) * incx {
+		panic('lapack: insufficient length of x')
+	}
+
+	if math.is_nan(scale) || math.is_nan(sumsq) {
+		return scale, sumsq
+	}
+
+	mut scl := scale
+	mut smsq := sumsq
+
+	if smsq == 0.0 {
+		scl = 1.0
+	}
+	if scl == 0.0 {
+		scl = 1.0
+		smsq = 0.0
+	}
+
+	if n == 0 {
+		return scl, smsq
+	}
+
+	// Compute the sum of squares in 3 accumulators:
+	//  - abig: sum of squares scaled down to avoid overflow
+	//  - asml: sum of squares scaled up to avoid underflow
+	//  - amed: sum of squares that do not require scaling
+	// The thresholds and multipliers are:
+	//  - values bigger than dtbig are scaled down by dsbig
+	//  - values smaller than dtsml are scaled up by dssml
+	mut is_big := false
+	mut asml, mut amed, mut abig := 0.0, 0.0, 0.0
+	mut ix := 0
+	for _ in 0 .. n {
+		mut ax := math.abs(x[ix])
+		if ax > dtbig {
+			ax *= dsbig
+			abig += ax * ax
+			is_big = true
+		} else if ax < dtsml {
+			if !is_big {
+				ax *= dssml
+				asml += ax * ax
+			}
+		} else {
+			amed += ax * ax
+		}
+		ix += incx
+	}
+	// Put the existing sum of squares into one of the accumulators.
+	if smsq > 0.0 {
+		ax := scl * math.sqrt(smsq)
+		if ax > dtbig {
+			if scl > 1.0 {
+				scl *= dsbig
+				abig += scl * scl * smsq
+			} else {
+				// sumsq > dtbig^2 => (dsbig * (dsbig * sumsq)) is representable.
+				abig += scl * scl * dsbig * dsbig * smsq
+			}
+		} else if ax < dtsml {
+			if !is_big {
+				if scl < 1.0 {
+					scl *= dssml
+					asml += scl * scl * smsq
+				} else {
+					// sumsq < dtsml^2 => (dssml * (dssml * sumsq)) is representable.
+					asml += scl * scl * dssml * dssml * smsq
+				}
+			}
+		} else {
+			amed += scl * scl * smsq
+		}
+	}
+	// Combine abig and amed or amed and asml if more than one accumulator was used.
+	if abig > 0.0 {
+		// Combine abig and amed:
+		if amed > 0.0 || math.is_nan(amed) {
+			abig += amed * dsbig * dsbig
+		}
+		scl = 1.0 / dsbig
+		smsq = abig
+	} else if asml > 0.0 {
+		// Combine amed and asml:
+		if amed > 0.0 || math.is_nan(amed) {
+			amed = math.sqrt(amed)
+			asml = math.sqrt(asml) / dssml
+			mut ymin, mut ymax := asml, amed
+			if asml > amed {
+				ymin, ymax = amed, asml
+			}
+			scl = 1.0
+			smsq = ymax * ymax * (1.0 + (ymin / ymax) * (ymin / ymax))
+		} else {
+			scl = 1.0 / dssml
+			smsq = asml
+		}
+	} else {
+		scl = 1.0
+		smsq = amed
+	}
+	return scl, smsq
+}