Merge branch 'features/O3' into dev

Author: Axect

src/fuga/mod.rs

@@ -170,6 +170,7 @@ pub use crate::statistics::stat::QType::{
 pub use crate::structure::matrix::{
     Form::{Diagonal, Identity},
     SolveKind::{LU, WAZ},
+    UPLO::{Upper, Lower}
 };
 pub use crate::structure::dataframe::DType::*;
 pub use crate::structure::ad::AD::*;

src/prelude/simpler.rs

@@ -27,13 +27,15 @@ pub trait SimplerLinearAlgebra {
     fn waz_diag(&self) -> Option<matrix::WAZD>;
     fn waz(&self) -> Option<matrix::WAZD>;
     fn qr(&self) -> matrix::QR;
+    fn cholesky(&self) -> Matrix;
     fn rref(&self) -> Matrix;
     fn det(&self) -> f64;
     fn block(&self) -> (Matrix, Matrix, Matrix, Matrix);
     fn inv(&self) -> Matrix;
     fn pseudo_inv(&self) -> Matrix;
     fn solve(&self, b: &Vec<f64>) -> Vec<f64>;
     fn solve_mat(&self, m: &Matrix) -> Matrix;
+    fn is_symmetric(&self) -> bool;
 }
 
 /// Simple Eigenpair
@@ -115,6 +117,14 @@ impl SimplerLinearAlgebra for Matrix {
     fn solve_mat(&self, m: &Matrix) -> Matrix {
         matrix::LinearAlgebra::solve_mat(self, m, matrix::SolveKind::LU)
     }
+
+    fn cholesky(&self) -> Matrix {
+        matrix::LinearAlgebra::cholesky(self, matrix::UPLO::Lower)
+    }
+
+    fn is_symmetric(&self) -> bool {
+        matrix::LinearAlgebra::is_symmetric(self)
+    }
 }
 
 /// Simple solve

src/structure/matrix.rs

@@ -594,7 +594,7 @@ use ::matrixmultiply;
 #[cfg(feature = "O3")]
 use blas::{daxpy, dgemm, dgemv};
 #[cfg(feature = "O3")]
-use lapack::{dgecon, dgeqrf, dgetrf, dgetri, dgetrs, dorgqr, dgesvd};
+use lapack::{dgecon, dgeqrf, dgetrf, dgetri, dgetrs, dorgqr, dgesvd, dpotrf};
 #[cfg(feature = "O3")]
 use std::f64::NAN;
 
@@ -617,7 +617,7 @@ use std::convert;
 pub use std::error::Error;
 use std::fmt;
 use std::ops::{Add, Div, Index, IndexMut, Mul, Neg, Sub};
-use crate::traits::sugar::{Scalable, ScalableMut};
+use crate::traits::sugar::ScalableMut;
 
 pub type Perms = Vec<(usize, usize)>;
 
@@ -2770,13 +2770,15 @@ pub trait LinearAlgebra {
     fn waz(&self, d_form: Form) -> Option<WAZD>;
     fn qr(&self) -> QR;
     fn svd(&self) -> SVD;
+    fn cholesky(&self, uplo: UPLO) -> Matrix;
     fn rref(&self) -> Matrix;
     fn det(&self) -> f64;
     fn block(&self) -> (Matrix, Matrix, Matrix, Matrix);
     fn inv(&self) -> Matrix;
     fn pseudo_inv(&self) -> Matrix;
     fn solve(&self, b: &Vec<f64>, sk: SolveKind) -> Vec<f64>;
     fn solve_mat(&self, m: &Matrix, sk: SolveKind) -> Matrix;
+    fn is_symmetric(&self) -> bool;
 }
 
 pub fn diag(n: usize) -> Matrix {
@@ -3221,6 +3223,55 @@ impl LinearAlgebra for Matrix {
         }
     }
 
+    /// Cholesky Decomposition
+    ///
+    /// # Examples
+    /// ```
+    /// extern crate peroxide;
+    /// use peroxide::fuga::*;
+    ///
+    /// fn main() {
+    ///     let a = ml_matrix("1 2;2 5");
+    ///     #[cfg(feature = "O3")]
+    ///     {
+    ///         let u = a.cholesky(Upper);
+    ///         let l = a.cholesky(Lower);
+    ///
+    ///         assert_eq!(u, ml_matrix("1 2;0 1"));
+    ///         assert_eq!(l, ml_matrix("1 0;2 1"));
+    ///     }
+    ///     a.print();
+    /// }
+    /// ```
+    fn cholesky(&self, uplo: UPLO) -> Matrix {
+        match () {
+            #[cfg(feature = "O3")]
+            () => {
+                if !self.is_symmetric() {
+                    panic!("Cholesky Error: Matrix is not symmetric!");
+                }
+                let dpotrf = lapack_dpotrf(self, uplo);
+                match dpotrf {
+                    None => panic!("Cholesky Error: Not symmetric or not positive definite."),
+                    Some(x) => {
+                        match x.status {
+                            POSITIVE_STATUS::Failed(i) => panic!("Cholesky Error: the leading minor of order {} is not positive definite", i),
+                            POSITIVE_STATUS::Success => {
+                                match uplo {
+                                    UPLO::Upper => x.get_U().unwrap(),
+                                    UPLO::Lower => x.get_L().unwrap()
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+            _ => {
+                unimplemented!()
+            }
+        }
+    }
+
     /// Reduced Row Echelon Form
     ///
     /// Implementation of [RosettaCode](https://rosettacode.org/wiki/Reduced_row_echelon_form)
@@ -3551,6 +3602,21 @@ impl LinearAlgebra for Matrix {
             }
         }
     }
+
+    fn is_symmetric(&self) -> bool {
+        if self.row != self.col {
+            return false;
+        }
+
+        for i in 0 .. self.row {
+            for j in i .. self.col {
+                if !nearly_eq(self[(i,j)], self[(j,i)]) {
+                    return false;
+                }
+            }
+        }
+        true
+    }
 }
 
 #[allow(non_snake_case)]
@@ -4038,6 +4104,20 @@ pub enum SVD_STATUS {
     Diverge(i32),
 }
 
+#[allow(non_camel_case_types)]
+#[derive(Debug, Copy, Clone, Eq, PartialEq)]
+pub enum POSITIVE_STATUS {
+    Success,
+    Failed(i32),
+}
+
+#[allow(non_camel_case_types)]
+#[derive(Debug, Copy, Clone, Eq, PartialEq)]
+pub enum UPLO {
+    Upper,
+    Lower
+}
+
 /// Temporary data structure from `dgetrf`
 #[derive(Debug, Clone)]
 pub struct DGETRF {
@@ -4062,6 +4142,13 @@ pub struct DGESVD {
     pub status: SVD_STATUS,
 }
 
+#[derive(Debug, Clone)]
+pub struct DPOTRF {
+    pub fact_mat: Matrix,
+    pub uplo: UPLO,
+    pub status: POSITIVE_STATUS
+}
+
 ///// Temporary data structure from `dgeev`
 //#[derive(Debug, Clone)]
 //pub struct DGEEV {
@@ -4300,6 +4387,54 @@ pub fn lapack_dgesvd(mat: &Matrix) -> Option<DGESVD> {
     }
 }
 
+#[allow(non_snake_case)]
+#[cfg(feature = "O3")]
+pub fn lapack_dpotrf(mat: &Matrix, UPLO: UPLO) -> Option<DPOTRF> {
+    match mat.shape {
+        Row => lapack_dpotrf(&mat.change_shape(), UPLO),
+        Col => {
+            let lda = mat.row as i32;
+            let N = mat.col as i32;
+            let mut A = mat.clone();
+            let mut info = 0i32;
+            let uplo = match UPLO {
+                UPLO::Upper => b'U',
+                UPLO::Lower => b'L'
+            };
+
+            unsafe {
+                dpotrf(
+                    uplo,
+                    N,
+                    &mut A.data,
+                    lda,
+                    &mut info,
+                )
+            }
+
+            if info == 0 {
+                Some(
+                    DPOTRF {
+                        fact_mat: matrix(A.data, mat.row, mat.col, Col),
+                        uplo: UPLO,
+                        status: POSITIVE_STATUS::Success
+                    }
+                )
+            } else if info > 0 {
+                Some(
+                    DPOTRF {
+                        fact_mat: matrix(A.data, mat.row, mat.col, Col),
+                        uplo: UPLO,
+                        status: POSITIVE_STATUS::Failed(info)
+                    }
+                )
+            } else {
+                None
+            }
+        }
+    }
+}
+
 #[allow(non_snake_case)]
 #[cfg(feature = "O3")]
 impl DGETRF {
@@ -4396,6 +4531,42 @@ impl DGEQRF {
     }
 }
 
+#[allow(non_snake_case)]
+impl DPOTRF {
+    pub fn get_U(&self) -> Option<Matrix> {
+        if self.uplo == UPLO::Lower {
+            return None;
+        }
+
+        let mat = &self.fact_mat;
+        let n = mat.col;
+        let mut result = matrix(vec![0f64; n.pow(2)], n, n, mat.shape);
+        for i in 0 .. n {
+            for j in i .. n {
+                result[(i, j)] = mat[(i, j)];
+            }
+        }
+        Some(result)
+    }
+
+    pub fn get_L(&self) -> Option<Matrix> {
+        if self.uplo == UPLO::Upper {
+            return None;
+        }
+
+        let mat = &self.fact_mat;
+        let n = mat.col;
+        let mut result = matrix(vec![0f64; n.pow(2)], n, n, mat.shape);
+
+        for i in 0 .. n {
+            for j in 0 .. i+1 {
+                result[(i, j)] = mat[(i, j)];
+            }
+        }
+        Some(result)
+    }
+}
+
 #[allow(non_snake_case)]
 pub fn gen_householder(a: &Vec<f64>) -> Matrix {
     let mut v = a.fmap(|t| t / (a[0] + a.norm(Norm::L2) * a[0].signum()));

src/structure/sparse.rs

@@ -7,6 +7,7 @@ use crate::traits::math::LinearOp;
 //use crate::traits::math::{InnerProduct, LinearOp, Norm, Normed, Vector};
 use crate::util::non_macro::zeros;
 use std::ops::Mul;
+use crate::fuga::UPLO;
 
 #[derive(Debug, Clone)]
 pub struct SPMatrix {
@@ -188,6 +189,14 @@ impl LinearAlgebra for SPMatrix {
     fn svd(&self) -> SVD {
         unimplemented!()
     }
+
+    fn cholesky(&self, uplo: UPLO) -> Matrix {
+        unimplemented!()
+    }
+
+    fn is_symmetric(&self) -> bool {
+        unimplemented!()
+    }
 }
 
 /// Matrix multiplication with vector