mdarray-linalg-faer 0.1.2

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
84
85
86

use dyn_stack::{MemBuffer, MemStack};
use faer_traits::ComplexField;
use mdarray::{DSlice, Layout};
use num_complex::ComplexFloat;

use crate::{into_faer, into_faer_mut};

pub fn qr_faer<
    T: ComplexFloat + ComplexField + Default + 'static,
    La: Layout,
    Lq: Layout,
    Lr: Layout,
>(
    a: &DSlice<T, 2, La>,
    q_mda: Option<&mut DSlice<T, 2, Lq>>,
    r_mda: &mut DSlice<T, 2, Lr>,
) {
    let (m, n) = *a.shape();
    let rank = Ord::min(m, 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 mut r_faer = into_faer_mut(r_mda);
    for i in 0..rank {
        for j in i..n {
            r_faer[(i, j)] = qr_mat[(i, j)];
        }
        for j in 0..i {
            r_faer[(i, j)] = T::zero();
        }
    }
    for i in rank..m {
        for j in 0..n {
            r_faer[(i, j)] = T::zero();
        }
    }

    if let Some(q) = q_mda {
        let mut q_faer = into_faer_mut(q);
        // TODO: check why this is necessary
        for i in 0..m {
            for j in 0..m {
                if i == j {
                    q_faer[(i, j)] = T::one();
                } else {
                    q_faer[(i, j)] = 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,
                            m,
                        )
                    )),
                );
    }
}