Impl ordinary diff equation solver

.vscode/settings.json

@@ -2,9 +2,11 @@
   "cSpell.words": [
     "aiᵢⱼ",
     "archlinux",
+    "Arenstorf",
     "Arioli",
     "BLACS",
     "blasint",
+    "bweuler",
     "Demmel",
     "devel",
     "dgeev",
@@ -17,14 +19,23 @@
     "dgetri",
     "dgssvx",
     "Dᵢⱼₖₗ",
+    "dopri",
+    "Dormand",
     "Dorozhinskii",
     "dpotrf",
     "dscal",
     "dsyev",
     "dsyrk",
     "dtype",
+    "dznrm2",
+    "FACCON",
+    "Fehlberg",
     "Flannery",
+    "FSAL",
+    "fweuler",
     "gfortran",
+    "Gustafsson",
+    "Heun",
     "ifort",
     "IIIₛ",
     "ᵢⱼₖₗ",
@@ -43,27 +54,36 @@
     "lredrhs",
     "lrhs",
     "lsol",
+    "Lubich",
     "lwork",
+    "mdeuler",
     "memcheck",
+    "Merson",
     "msgpass",
     "nelt",
+    "Nørsett",
+    "nstage",
+    "nstep",
     "odyad",
     "oneapi",
     "PERMUT",
     "PREORDER",
     "pthread",
+    "Radau",
     "RINFOG",
     "rowcol",
     "rowcoliter",
     "rowcolrig",
     "rustup",
     "Schur",
+    "substeps",
     "Teukolsky",
     "tgamma",
     "tocsc",
     "tocsr",
     "udyad",
     "unsym",
+    "Verner",
     "Vetterling",
     "zcopy",
     "zgemm",
@@ -74,6 +94,8 @@
     "zgetri",
     "zherk",
     "zlange",
+    "ZMUMPS",
+    "Zonneveld",
     "zpotrf",
     "zscal",
     "zsyrk"

Cargo.toml

@@ -1,10 +1,11 @@
 [workspace]
 
 members = [
-    "russell_lab",
-    "russell_sparse",
-    "russell_stat",
-    "russell_tensor"
+  "russell_lab",
+  "russell_ode",
+  "russell_sparse",
+  "russell_stat",
+  "russell_tensor",
 ]
 
 resolver = "2"

README.md

@@ -7,7 +7,7 @@
 
 ![Bertrand Russell](Bertrand_Russell_1957.jpg)
 
-([CC0](http://creativecommons.org/publicdomain/zero/1.0/deed.en). Photo: [Bertrand Russell](https://en.wikipedia.org/wiki/Bertrand_Russell))
+([CC0](https://creativecommons.org/publicdomain/zero/1.0/deed.en). Photo: [Bertrand Russell](https://en.wikipedia.org/wiki/Bertrand_Russell))
 
 ## Contents
 
@@ -298,7 +298,7 @@ fn main() -> Result<(), StrError> {
     //  . -1 -3  2  .
     //  .  .  1  .  .
     //  . 4  2  .  1
-    let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None, false)?;
+    let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?;
     coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate
     coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate
     coo.put(1, 0, 3.0)?;
@@ -336,7 +336,7 @@ fn main() -> Result<(), StrError> {
     // analysis
     let mut stats = StatsLinSol::new();
     umfpack.update_stats(&mut stats);
-    let (mx, ex) = (stats.determinant.mantissa, stats.determinant.exponent);
+    let (mx, ex) = (stats.determinant.mantissa_real, stats.determinant.exponent);
     println!("det(a) = {:?}", mx * f64::powf(10.0, ex));
     println!("rcond  = {:?}", stats.output.umfpack_rcond_estimate);
     Ok(())

russell_lab/Cargo.toml

@@ -21,7 +21,6 @@ num-traits = "0.2"
 serde = { version = "1.0", features = ["derive"] }
 
 [dev-dependencies]
-rmp-serde = "1.1"
 serde_json = "1.0"
 
 [build-dependencies]

russell_lab/c_code/interface_blas.c

@@ -14,6 +14,8 @@
 #define FN_DGESVD dgesvd_
 #define FN_DGETRF dgetrf_
 #define FN_DGETRI dgetri_
+#define FN_ZGETRF zgetrf_
+#define FN_ZGETRI zgetri_
 #else
 #include "cblas.h"
 #include "lapack.h"
@@ -28,6 +30,8 @@
 #define FN_DGESVD LAPACK_dgesvd
 #define FN_DGETRF LAPACK_dgetrf
 #define FN_DGETRI LAPACK_dgetri
+#define FN_ZGETRF LAPACK_zgetrf
+#define FN_ZGETRI LAPACK_zgetri
 #endif
 
 #include "constants.h"
@@ -78,7 +82,7 @@ int32_t c_get_num_threads() {
 }
 
 // Computes the solution to a system of linear equations
-// <http://www.netlib.org/lapack/explore-html/d8/d72/dgesv_8f.html>
+// <https://www.netlib.org/lapack/explore-html/d8/d72/dgesv_8f.html>
 void c_dgesv(const int32_t *n,
              const int32_t *nrhs,
              double *a,
@@ -91,7 +95,7 @@ void c_dgesv(const int32_t *n,
 }
 
 // Computes the solution to a real system of linear equations (complex version)
-// <http://www.netlib.org/lapack/explore-html/d1/ddc/zgesv_8f.html>
+// <https://www.netlib.org/lapack/explore-html/d1/ddc/zgesv_8f.html>
 void c_zgesv(const int32_t *n,
              const int32_t *nrhs,
              COMPLEX64 *a,
@@ -104,7 +108,7 @@ void c_zgesv(const int32_t *n,
 }
 
 // Computes the matrix norm
-// <http://www.netlib.org/lapack/explore-html/dc/d09/dlange_8f.html>
+// <https://www.netlib.org/lapack/explore-html/dc/d09/dlange_8f.html>
 double c_dlange(int32_t norm_code,
                 const int32_t *m,
                 const int32_t *n,
@@ -119,7 +123,7 @@ double c_dlange(int32_t norm_code,
 }
 
 // Computes the matrix norm (complex version)
-// <http://www.netlib.org/lapack/explore-html/d5/d8f/zlange_8f.html>
+// <https://www.netlib.org/lapack/explore-html/d5/d8f/zlange_8f.html>
 double c_zlange(int32_t norm_code,
                 const int32_t *m,
                 const int32_t *n,
@@ -134,7 +138,7 @@ double c_zlange(int32_t norm_code,
 }
 
 // Computes the Cholesky factorization of a real symmetric positive definite matrix
-// <http://www.netlib.org/lapack/explore-html/d0/d8a/dpotrf_8f.html>
+// <https://www.netlib.org/lapack/explore-html/d0/d8a/dpotrf_8f.html>
 void c_dpotrf(C_BOOL upper,
               const int32_t *n,
               double *a,
@@ -161,7 +165,7 @@ void c_dsyev(C_BOOL calc_v,
 }
 
 // Computes the eigenvalues and eigenvectors of a general matrix
-// <http://www.netlib.org/lapack/explore-html/d9/d28/dgeev_8f.html>
+// <https://www.netlib.org/lapack/explore-html/d9/d28/dgeev_8f.html>
 void c_dgeev(C_BOOL calc_vl,
              C_BOOL calc_vr,
              const int32_t *n,
@@ -182,7 +186,7 @@ void c_dgeev(C_BOOL calc_vl,
 }
 
 // Computes the singular value decomposition (SVD)
-// <http://www.netlib.org/lapack/explore-html/d8/d2d/dgesvd_8f.html>
+// <https://www.netlib.org/lapack/explore-html/d8/d2d/dgesvd_8f.html>
 void c_dgesvd(int32_t jobu_code,
               int32_t jobvt_code,
               const int32_t *m,
@@ -209,7 +213,7 @@ void c_dgesvd(int32_t jobu_code,
 }
 
 // Computes the LU factorization of a general (m,n) matrix
-// <http://www.netlib.org/lapack/explore-html/d3/d6a/dgetrf_8f.html>
+// <https://www.netlib.org/lapack/explore-html/d3/d6a/dgetrf_8f.html>
 void c_dgetrf(const int32_t *m,
               const int32_t *n,
               double *a,
@@ -220,7 +224,7 @@ void c_dgetrf(const int32_t *m,
 }
 
 // Computes the inverse of a matrix using the LU factorization computed by dgetrf
-// <http://www.netlib.org/lapack/explore-html/df/da4/dgetri_8f.html>
+// <https://www.netlib.org/lapack/explore-html/df/da4/dgetri_8f.html>
 void c_dgetri(const int32_t *n,
               double *a,
               const int32_t *lda,
@@ -230,3 +234,26 @@ void c_dgetri(const int32_t *n,
               int32_t *info) {
     FN_DGETRI(n, a, lda, ipiv, work, lwork, info);
 }
+
+// Computes the LU factorization of a general (m,n) matrix
+// <https://www.netlib.org/lapack/explore-html/dd/dd1/zgetrf_8f.html>
+void c_zgetrf(const int32_t *m,
+              const int32_t *n,
+              COMPLEX64 *a,
+              const int32_t *lda,
+              int32_t *ipiv,
+              int32_t *info) {
+    FN_ZGETRF(m, n, a, lda, ipiv, info);
+}
+
+// Computes the inverse of a matrix using the LU factorization computed by zgetrf
+// <https://www.netlib.org/lapack/explore-html/d0/db3/zgetri_8f.html>
+void c_zgetri(const int32_t *n,
+              COMPLEX64 *a,
+              const int32_t *lda,
+              const int32_t *ipiv,
+              COMPLEX64 *work,
+              const int32_t *lwork,
+              int32_t *info) {
+    FN_ZGETRI(n, a, lda, ipiv, work, lwork, info);
+}

russell_lab/examples/ex_matrix_eigenvalues.rs

@@ -58,8 +58,10 @@ fn main() -> Result<(), StrError> {
     // err := a⋅v - v⋅λ
     // ```
     let a = ComplexMatrix::from(&data);
-    let v = complex_mat_zip(&v_real, &v_imag)?;
-    let d = complex_vec_zip(&l_real, &l_imag)?;
+    let mut v = ComplexMatrix::new(m, m);
+    let mut d = ComplexVector::new(m);
+    complex_mat_zip(&mut v, &v_real, &v_imag)?;
+    complex_vec_zip(&mut d, &l_real, &l_imag)?;
     let lam = ComplexMatrix::diagonal(d.as_data());
     let mut a_v = ComplexMatrix::new(m, m);
     let mut v_l = ComplexMatrix::new(m, m);

russell_lab/src/base/auxiliary_blas.rs

@@ -6,7 +6,7 @@ extern "C" {
     fn c_get_num_threads() -> i32;
 
     // Finds the index of the maximum absolute value
-    // <http://www.netlib.org/lapack/explore-html/dd/de0/idamax_8f.html>
+    // <https://www.netlib.org/lapack/explore-html/dd/de0/idamax_8f.html>
     fn cblas_idamax(n: i32, x: *const f64, incx: i32) -> i32;
 }
 

russell_lab/src/base/formatters.rs

@@ -66,11 +66,31 @@ pub fn format_nanoseconds(nanoseconds: u128) -> String {
     buf
 }
 
+/// Formats a number using the scientific notation
+pub fn format_scientific(num: f64, width: usize, precision: usize) -> String {
+    // based on <https://stackoverflow.com/questions/65264069/alignment-of-floating-point-numbers-printed-in-scientific-notation>
+    const EXP_PAD: usize = 2;
+    let mut result = format!("{:.precision$e}", num, precision = precision);
+    let exp = result.split_off(result.find('e').unwrap());
+    let (sign, exp) = if exp.starts_with("e-") {
+        ('-', &exp[2..])
+    } else {
+        ('+', &exp[1..])
+    };
+    result.push_str(&format!("E{}{:0>pad$}", sign, exp, pad = EXP_PAD));
+    format!("{:>width$}", result, width = width)
+}
+
+/// Formats a number using the scientific notation as in Fortran with the ES23.15 format
+pub fn format_fortran(num: f64) -> String {
+    format_scientific(num, 23, 15)
+}
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
 #[cfg(test)]
 mod tests {
-    use super::format_nanoseconds;
+    use super::{format_fortran, format_nanoseconds, format_scientific};
 
     #[test]
     fn format_nanoseconds_works() {
@@ -170,4 +190,32 @@ mod tests {
         res = format_nanoseconds(3_601_100_000_000);
         assert_eq!(res, "1h1.1s");
     }
+
+    #[test]
+    fn format_scientific_works() {
+        assert_eq!(format_scientific(0.1111, 9, 2), " 1.11E-01");
+        assert_eq!(format_scientific(0.02222, 11, 4), " 2.2220E-02");
+        assert_eq!(format_scientific(3333.0, 10, 3), " 3.333E+03");
+        assert_eq!(format_scientific(-44444.0, 9, 1), " -4.4E+04");
+        assert_eq!(format_scientific(0.0, 8, 1), " 0.0E+00");
+        assert_eq!(format_scientific(1.0, 23, 15), "  1.000000000000000E+00");
+        assert_eq!(format_scientific(42.0, 23, 15), "  4.200000000000000E+01");
+        assert_eq!(format_scientific(9999999999.00, 8, 1), " 1.0E+10");
+        assert_eq!(format_scientific(999999999999.00, 23, 15), "  9.999999999990000E+11");
+        assert_eq!(format_scientific(123456789.1011, 11, 4), " 1.2346E+08");
+    }
+
+    #[test]
+    fn format_fortran_works() {
+        assert_eq!(format_fortran(0.1111), "  1.111000000000000E-01");
+        assert_eq!(format_fortran(0.02222), "  2.222000000000000E-02");
+        assert_eq!(format_fortran(3333.0), "  3.333000000000000E+03");
+        assert_eq!(format_fortran(-44444.0), " -4.444400000000000E+04");
+        assert_eq!(format_fortran(0.0), "  0.000000000000000E+00");
+        assert_eq!(format_fortran(1.0), "  1.000000000000000E+00");
+        assert_eq!(format_fortran(42.0), "  4.200000000000000E+01");
+        assert_eq!(format_fortran(9999999999.00), "  9.999999999000000E+09");
+        assert_eq!(format_fortran(999999999999.00), "  9.999999999990000E+11");
+        assert_eq!(format_fortran(123456789.1011), "  1.234567891011000E+08");
+    }
 }