extern crate ohsl;
pub use ohsl::vector::{Vector, Vec64};
pub use ohsl::matrix::{Matrix, Mat64};
pub use ohsl::complex::{Complex, Cmplx};
fn main() {
println!( "------------------ Linear system ------------------" );
println!( " * Solving the linear system Ax = b, for x, where" );
let mut a = Mat64::new( 3, 3, 0.0 );
a[(0,0)] = 1.0; a[(0,1)] = 1.0; a[(0,2)] = 1.0;
a[(1,0)] = 0.0; a[(1,1)] = 2.0; a[(1,2)] = 5.0;
a[(2,0)] = 2.0; a[(2,1)] = 5.0; a[(2,2)] = - 1.0;
println!( " * A ={}", a );
let b = Vec64::create( vec![ 6.0, -4.0, 27.0 ] );
println!( " * b^T ={}", b );
let x = a.clone().solve_basic( &b );
println!( " * gives the solution vector");
println!( " * x^T ={}", x );
println!( "---------------------------------------------------" );
println!( " * Lets solve the complex linear system Cy = d where" );
let mut c = Matrix::<Cmplx>::new( 2, 2, Cmplx::new( 1.0, 1.0 ) );
c[(0,1)] = Cmplx::new( -1.0, 0.0 );
c[(1,0)] = Cmplx::new( 1.0, -1.0 );
println!( " * C ={}", c );
let mut d = Vector::<Cmplx>::new( 2, Cmplx::new( 0.0, 1.0 ) );
d[1] = Cmplx::new( 1.0, 0.0 );
println!( " * d^T ={}", d );
println!( " * Using LU decomposition we find that ");
let y = c.clone().solve_lu( &d );
println!( " * y^T ={}", y );
println!( "---------------------------------------------------" );
println!( " * We may also find the determinant of a matrix" );
println!( " * |A| = {}", a.determinant() );
println!( " * |C| = {}", c.determinant() );
println!( " * or the inverse of a matrix" );
let inverse = a.inverse();
println!( " * A^-1 =\n{}", inverse );
println!( " * C^-1 =\n{}", c.inverse() );
println!( " * We can check this by multiplication" );
println!( " * I = A * A^-1 =\n{}", a * inverse );
println!( " * which is sufficiently accurate for our purposes.");
println!( "-------------------- FINISHED ---------------------" );
}