faer 0.24.0

linear algebra library
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

#![allow(non_snake_case)]
use diol::prelude::*;
use dyn_stack::{MemBuffer, MemStack};
use faer::diag::Diag;
use faer::linalg::cholesky::lblt;
use faer::prelude::*;
use faer::reborrow::ReborrowMut;
use faer::stats::prelude::*;
#[cfg(any(openblas, mkl, blis))]
use lapack_sys as la;
#[cfg(blis)]
extern crate blis_src;
#[cfg(mkl)]
extern crate intel_mkl_src;
#[cfg(openblas)]
extern crate openblas_src;
#[cfg(any(openblas, mkl, blis))]
fn lapack(bencher: Bencher, PlotArg(n): PlotArg) {
	use aligned_vec::avec;
	let rng = &mut StdRng::seed_from_u64(0);
	let A = CwiseMatDistribution {
		nrows: n,
		ncols: n,
		dist: StandardNormal,
	}
	.rand::<Mat<f64>>(rng);
	let A = &A + A.adjoint();
	let mut lblt = A.clone();
	let perm = &mut *vec![0usize; n];
	let perm_inv = &mut *vec![0usize; n];
	let lwork = unsafe {
		let mut lwork = core::mem::zeroed();
		la::dsytrf_(
			&(b'L' as _),
			(&n) as *const _ as *const _,
			lblt.as_ptr_mut() as _,
			(&lblt.col_stride()) as *const _ as *const _,
			perm.as_mut_ptr() as _,
			&mut lwork,
			(&-1isize) as *const _ as *const _,
			(&mut 0usize) as *mut _ as *mut _,
		);
		lwork as usize
	};
	let work = &mut *avec![0.0f64; lwork];
	bencher.bench(|| unsafe {
		lblt.copy_from_triangular_lower(&A);
		perm.fill(0);
		la::dsytrf_(
			&(b'L' as _),
			(&n) as *const _ as *const _,
			lblt.as_ptr_mut() as _,
			(&lblt.col_stride()) as *const _ as *const _,
			perm.as_mut_ptr() as _,
			work.as_mut_ptr() as _,
			(&lwork) as *const _ as *const _,
			(&mut 0usize) as *mut _ as *mut _,
		);
	});
}
fn faer(bencher: Bencher, PlotArg(n): PlotArg) {
	let rng = &mut StdRng::seed_from_u64(0);
	let A = CwiseMatDistribution {
		nrows: n,
		ncols: n,
		dist: StandardNormal,
	}
	.rand::<Mat<f64>>(rng);
	let A = &A + A.adjoint();
	let mut lblt = A.clone();
	let mut subdiag = Diag::zeros(n);
	let perm = &mut *vec![0usize; n];
	let perm_inv = &mut *vec![0usize; n];
	let par = Par::Seq;
	let params = default();
	let scratch =
		lblt::factor::cholesky_in_place_scratch::<usize, f64>(n, par, params);
	let mut mem = MemBuffer::new(scratch);
	let stack = MemStack::new(&mut mem);
	bencher.bench(|| {
		lblt.copy_from_triangular_lower(&A);
		lblt::factor::cholesky_in_place(
			lblt.rb_mut(),
			subdiag.rb_mut(),
			perm,
			perm_inv,
			par,
			stack,
			params,
		);
	});
}
fn main() -> eyre::Result<()> {
	let bench = Bench::from_args()?;
	bench.register_many(
		"lblt",
		{
			let f = list![faer];
			#[cfg(any(openblas, mkl, blis))]
			let f = diol::variadics::Cons {
				head: lapack,
				tail: f,
			};
			f
		},
		[64, 96, 128, 192, 256, 1024, 2048, 4096, 8192].map(PlotArg),
	);
	bench.run()?;
	Ok(())
}