mdarray-linalg-faer 0.2.0

Faer backend for mdarray-linalg
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

use dyn_stack::{MemBuffer, MemStack};
use faer_traits::ComplexField;
use mdarray::{Dim, Layout, Shape, Slice};
use num_complex::ComplexFloat;

use crate::{into_faer, into_faer_mut};

pub(super) fn qr_faer<
    T: ComplexFloat + ComplexField + Default,
    D0: Dim,
    D1: Dim,
    D2: Dim,
    La: Layout,
    Lq: Layout,
    Lr: Layout,
>(
    a: &Slice<T, (D0, D1), La>,
    q_mda: Option<&mut Slice<T, (D0, D2), Lq>>,
    r_mda: &mut Slice<T, (D2, D1), Lr>,
) {
    let ash = *a.shape();
    let (m, n) = (ash.dim(0), ash.dim(1));
    let rank = m.min(n);
    let par = faer::get_global_parallelism();
    let block_size = faer::linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);

    let mut qr_mat = into_faer(a).to_owned();
    let mut h_factor = faer::Mat::<T>::zeros(block_size, rank);

    let _ = faer::linalg::qr::no_pivoting::factor::qr_in_place(
        qr_mat.as_mut(),
        h_factor.as_mut(),
        par,
        MemStack::new(&mut MemBuffer::new(
            faer::linalg::qr::no_pivoting::factor::qr_in_place_scratch::<T>(
                m,
                n,
                block_size,
                par,
                faer::prelude::default(),
            ),
        )),
        faer::prelude::default(),
    );

    let r_rows = r_mda.shape().dim(0); 
    let mut r_faer = into_faer_mut(r_mda);
    for i in 0..r_rows {
        for j in 0..n {
            r_faer[(i, j)] = if i < rank && j >= i {
                qr_mat[(i, j)]
            } else {
                T::zero()
            };
        }
    }

    if let Some(q) = q_mda {
        let q_cols = q.shape().dim(1); 
        let mut q_faer = into_faer_mut(q);

        for i in 0..m {
            for j in 0..q_cols {
                q_faer[(i, j)] = if i == j { T::one() } else { T::zero() };
            }
        }

        faer::linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
            qr_mat.as_ref(),
            h_factor.as_ref(),
            faer::Conj::No,
            q_faer,
            par,
            MemStack::new(&mut MemBuffer::new(
                faer::linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
                    m,
                    block_size,
                    q_cols,  
                ),
            )),
        );
    }
}