kryst 3.2.1

Krylov subspace and preconditioned iterative solvers for dense and sparse linear systems, with shared and distributed memory parallelism.
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

use crate::algebra::prelude::*;
use crate::error::KError;

/// Column-major dense block vector storage used by block Krylov variants.
#[derive(Debug, Clone, Default)]
pub struct BlockVec {
    data: Vec<S>,
    n: usize,
    p: usize,
}

impl BlockVec {
    /// Create a new block vector with `n` rows and `p` columns.
    pub fn new(n: usize, p: usize) -> Self {
        Self {
            data: vec![S::zero(); n.saturating_mul(p)],
            n,
            p,
        }
    }

    /// Resize the block vector to `n` rows and `p` columns, zero-filling new entries.
    pub fn resize(&mut self, n: usize, p: usize) {
        let need = n.saturating_mul(p);
        if self.data.len() != need {
            self.data.resize(need, S::zero());
        }
        self.n = n;
        self.p = p;
    }

    /// Number of rows in the block vector.
    #[inline]
    pub fn nrows(&self) -> usize {
        self.n
    }

    /// Number of columns in the block vector.
    #[inline]
    pub fn ncols(&self) -> usize {
        self.p
    }

    /// Immutable view into the `j`-th column.
    #[inline]
    pub fn col(&self, j: usize) -> &[S] {
        debug_assert!(
            j < self.p,
            "BlockVec::col: column {} out of range {}",
            j,
            self.p
        );
        let offset = j * self.n;
        &self.data[offset..offset + self.n]
    }

    /// Mutable view into the `j`-th column.
    #[inline]
    pub fn col_mut(&mut self, j: usize) -> &mut [S] {
        debug_assert!(
            j < self.p,
            "BlockVec::col_mut: column {} out of range {}",
            j,
            self.p
        );
        let offset = j * self.n;
        &mut self.data[offset..offset + self.n]
    }

    /// Immutable view into the raw column-major storage.
    #[inline]
    pub fn as_slice(&self) -> &[S] {
        &self.data
    }

    /// Mutable view into the raw column-major storage.
    #[inline]
    pub fn as_mut_slice(&mut self) -> &mut [S] {
        &mut self.data
    }
}

impl BlockVec {
    /// Fill the block vector with zeros.
    pub fn fill_zero(&mut self) {
        for v in &mut self.data {
            *v = S::zero();
        }
    }
}

/// Convenience helper for verifying block dimensions.
#[allow(dead_code)]
pub(crate) fn assert_block_dims(expected_rows: usize, vec: &BlockVec) -> Result<(), KError> {
    if vec.nrows() != expected_rows {
        return Err(KError::InvalidInput(format!(
            "BlockVec has {} rows but expected {}",
            vec.nrows(),
            expected_rows
        )));
    }
    Ok(())
}