diffsol_la/matrix/
sparsity.rs

1use crate::{
2    error::LaError,
3    matrix_error,
4    scalar::IndexType,
5    vector::{Vector, VectorIndex},
6    ColMajBlock, VectorCommon,
7};
8
9use super::{Matrix, MatrixCommon};
10
11pub trait MatrixSparsity<M: Matrix>: Sized + Clone {
12    fn nrows(&self) -> IndexType;
13    fn ncols(&self) -> IndexType;
14    fn is_sparse() -> bool;
15    fn try_from_indices(
16        nrows: IndexType,
17        ncols: IndexType,
18        indices: Vec<(IndexType, IndexType)>,
19    ) -> Result<Self, LaError>;
20    fn indices(&self) -> Vec<(IndexType, IndexType)>;
21    fn union(self, other: M::SparsityRef<'_>) -> Result<M::Sparsity, LaError>;
22    fn new_diagonal(n: IndexType) -> Self;
23    fn as_ref(&self) -> M::SparsityRef<'_>;
24    fn get_index(
25        &self,
26        indices: &[(IndexType, IndexType)],
27        ctx: <M::V as VectorCommon>::C,
28    ) -> <M::V as Vector>::Index;
29}
30
31pub trait MatrixSparsityRef<'a, M: Matrix> {
32    fn nrows(&self) -> IndexType;
33    fn ncols(&self) -> IndexType;
34    fn is_sparse() -> bool;
35    fn indices(&self) -> Vec<(IndexType, IndexType)>;
36    fn to_owned(&self) -> M::Sparsity;
37    fn split(
38        &self,
39        algebraic_indices: &<M::V as Vector>::Index,
40    ) -> [(M::Sparsity, <M::V as Vector>::Index); 4];
41}
42
43#[derive(Clone)]
44pub struct Dense<M: Matrix> {
45    nrows: IndexType,
46    ncols: IndexType,
47    _phantom: std::marker::PhantomData<M>,
48}
49
50pub struct DenseRef<'a, M: Matrix> {
51    dense: &'a Dense<M>,
52}
53
54impl<M: Matrix> Dense<M> {
55    pub fn new(nrows: IndexType, ncols: IndexType) -> Self {
56        Dense {
57            nrows,
58            ncols,
59            _phantom: std::marker::PhantomData,
60        }
61    }
62    pub fn nrows(&self) -> IndexType {
63        self.nrows
64    }
65    pub fn ncols(&self) -> IndexType {
66        self.ncols
67    }
68    pub(crate) fn split(
69        &self,
70        algebraic_indices: &<M::V as Vector>::Index,
71    ) -> [(Self, <<M as MatrixCommon>::V as Vector>::Index); 4] {
72        let (ul_blk, ur_blk, ll_blk, lr_blk) =
73            ColMajBlock::split(self.nrows, self.ncols, algebraic_indices);
74        let ul = Dense::new(ul_blk.nrows, ul_blk.ncols);
75        let ur = Dense::new(ur_blk.nrows, ur_blk.ncols);
76        let ll = Dense::new(ll_blk.nrows, ll_blk.ncols);
77        let lr = Dense::new(lr_blk.nrows, lr_blk.ncols);
78        [
79            (ul, ul_blk.src_indices),
80            (ur, ur_blk.src_indices),
81            (ll, ll_blk.src_indices),
82            (lr, lr_blk.src_indices),
83        ]
84    }
85}
86
87impl<'a, M: Matrix> DenseRef<'a, M> {
88    pub fn new(dense: &'a Dense<M>) -> Self {
89        DenseRef { dense }
90    }
91}
92
93impl<M> MatrixSparsity<M> for Dense<M>
94where
95    for<'a> M: Matrix<Sparsity = Dense<M>, SparsityRef<'a> = DenseRef<'a, M>>,
96{
97    fn union(self, other: M::SparsityRef<'_>) -> Result<M::Sparsity, LaError> {
98        if self.nrows() != other.nrows() || self.ncols() != other.ncols() {
99            return Err(matrix_error!(UnionIncompatibleShapes));
100        }
101        Ok(Self::new(self.nrows(), self.ncols()))
102    }
103
104    fn as_ref(&self) -> M::SparsityRef<'_> {
105        DenseRef::new(self)
106    }
107
108    fn nrows(&self) -> IndexType {
109        self.nrows
110    }
111
112    fn ncols(&self) -> IndexType {
113        self.ncols
114    }
115
116    fn is_sparse() -> bool {
117        false
118    }
119
120    fn try_from_indices(
121        nrows: IndexType,
122        ncols: IndexType,
123        _indices: Vec<(IndexType, IndexType)>,
124    ) -> Result<Self, LaError> {
125        if nrows == 0 || ncols == 0 {
126            return Err(matrix_error!(MatrixShapeError));
127        }
128        Ok(Dense::new(nrows, ncols))
129    }
130
131    fn indices(&self) -> Vec<(IndexType, IndexType)> {
132        Vec::new()
133    }
134
135    fn new_diagonal(n: IndexType) -> Self {
136        Dense::new(n, n)
137    }
138    fn get_index(
139        &self,
140        indices: &[(IndexType, IndexType)],
141        ctx: <M::V as VectorCommon>::C,
142    ) -> <M::V as Vector>::Index {
143        let indices: Vec<_> = indices
144            .iter()
145            .map(|(i, j)| {
146                if i >= &self.nrows() || j >= &self.ncols() {
147                    panic!("Index out of bounds")
148                }
149                j * self.nrows() + i
150            })
151            .collect();
152        <M::V as Vector>::Index::from_vec(indices, ctx)
153    }
154}
155
156impl<'a, M> MatrixSparsityRef<'a, M> for DenseRef<'a, M>
157where
158    M: Matrix<Sparsity = Dense<M>, SparsityRef<'a> = Self>,
159{
160    fn to_owned(&self) -> M::Sparsity {
161        Dense::new(self.nrows(), self.ncols())
162    }
163
164    fn split(
165        &self,
166        indices: &<M::V as Vector>::Index,
167    ) -> [(M::Sparsity, <<M as MatrixCommon>::V as Vector>::Index); 4] {
168        self.dense.split(indices)
169    }
170
171    fn nrows(&self) -> IndexType {
172        self.dense.nrows
173    }
174
175    fn ncols(&self) -> IndexType {
176        self.dense.ncols
177    }
178
179    fn is_sparse() -> bool {
180        false
181    }
182
183    fn indices(&self) -> Vec<(IndexType, IndexType)> {
184        Vec::new()
185    }
186}
diffsol_la/matrix/sparsity.rs

diffsol_la/matrix/
sparsity.rs
