Use CBLAS for double precision vector norm

blas/blas_test.v

@@ -295,6 +295,15 @@ fn test_dnrm2() {
 	}
 }
 
+fn test_dnrm2_handles_stride_and_extreme_magnitudes() {
+	strided := dnrm2(2, [3.0, 99.0, 4.0], 2)
+	assert float64.tolerance(strided, 5.0, test_tol)
+	large := dnrm2(2, [1e308, 1e308], 1)
+	assert float64.tolerance(large / 1e308, math.sqrt(2.0), test_tol)
+	small := dnrm2(2, [1e-308, 1e-308], 1)
+	assert float64.tolerance(small / 1e-308, math.sqrt(2.0), test_tol)
+}
+
 fn test_idamax() {
 	for case in level1_test_cases {
 		result := idamax(case.n, case.x, case.incx)

blas/oblas_d_vsl_blas_cblas.v

@@ -292,7 +292,10 @@ pub fn dnrm2(n int, x []f64, incx int) f64 {
 	if n <= 0 || x.len == 0 {
 		return 0.0
 	}
-	return blas64.dnrm2(n, x, incx)
+	if incx < 1 || x.len <= (n - 1) * incx {
+		return blas64.dnrm2(n, x, incx)
+	}
+	return C.cblas_dnrm2(n, unsafe { &x[0] }, incx)
 }
 
 // scnrm2 exposes this operation as part of the public API.