rusolver/
matrix.rs

1// SPDX-License-Identifier: Apache-2.0
2use crate::SolverError;
3use crate::numerics::{
4    finite, zeros
5};
6/// Owned finite row-major FP64 data. Empty axes are rejected in this first release.
7#[derive(Clone, Debug, PartialEq)]
8pub struct Matrix {
9    pub(crate) rows: usize, pub(crate) cols: usize, pub(crate) data: Vec<f64>
10}
11impl Matrix {
12    pub fn new(rows: usize, cols: usize, data: Vec<f64>) -> Result<Self, SolverError> {
13        let size = checked_shape(rows, cols)?;
14        if data.len() != size {
15            return Err(SolverError::Shape("rows * columns != data length"));
16        }
17        finite(&data)?;
18        Ok(Self {
19            rows, cols, data
20        })
21    }
22    /// Explicit FP32 -> FP64 host conversion; no device transfer occurs.
23    pub fn from_f32(rows: usize, cols: usize, data: &[f32]) -> Result<Self, SolverError> {
24        let size = checked_shape(rows, cols)?;
25        if data.len() != size {
26            return Err(SolverError::Shape("FP32 input length"));
27        }
28        let mut out = zeros(size)?;
29        for (i, &x) in data.iter().enumerate() {
30            out[i] = f64::from(x);
31        }
32        Self::new(rows, cols, out)
33    }
34    pub fn zeros(rows: usize, cols: usize) -> Result<Self, SolverError> {
35        let size = checked_shape(rows, cols)?;
36        Ok(Self {
37            rows, cols, data: zeros(size)?
38        })
39    }
40    pub fn identity(n: usize) -> Result<Self, SolverError> {
41        let mut result = Self::zeros(n, n)?;
42        for i in 0..n {
43            result.data[i*n+i] = 1.0;
44        }
45        Ok(result)
46    }
47    pub fn rows(&self) -> usize {
48        self.rows
49    }
50    pub fn columns(&self) -> usize {
51        self.cols
52    }
53    pub fn values(&self) -> &[f64] {
54        &self.data
55    }
56    pub fn into_values(self) -> Vec<f64> {
57        self.data
58    }
59    pub fn view(&self) -> MatrixView<'_> {
60        MatrixView {
61            data: &self.data, rows: self.rows, cols: self.cols,
62            offset: 0, row_stride: self.cols, col_stride: 1
63        }
64    }
65    pub fn get(&self, row: usize, column: usize) -> Option<f64> {
66        self.view().get(row, column)
67    }
68    pub fn transpose(&self) -> Result<Self, SolverError> {
69        self.view().transpose().to_owned()
70    }
71}
72/// Immutable matrix view. Constructor validates reachable indices and values;
73/// unrelated padding elements may contain arbitrary values. Overlapping READS
74/// are allowed; the view cannot be used to create mutable aliases.
75#[derive(Clone, Copy, Debug)]
76pub struct MatrixView<'a> {
77    data: &'a [f64], rows: usize, cols: usize,
78    offset: usize, row_stride: usize, col_stride: usize,
79}
80impl<'a> MatrixView<'a> {
81    pub fn contiguous(rows: usize, cols: usize, data: &'a [f64]) -> Result<Self, SolverError> {
82        if checked_shape(rows, cols)? != data.len() {
83            return Err(SolverError::Shape("contiguous data length"));
84        }
85        Self::strided(rows, cols, data, 0, cols, 1)
86    }
87    pub fn strided(rows: usize, cols: usize, data: &'a [f64], offset: usize,
88    row_stride: usize, col_stride: usize) -> Result<Self, SolverError> {
89        checked_shape(rows, cols)?;
90        if row_stride == 0 || col_stride == 0 {
91            return Err(SolverError::Shape("strides must be positive"));
92        }
93        let last = (rows-1).checked_mul(row_stride)
94        .and_then(|r| (cols-1).checked_mul(col_stride).and_then(|c| r.checked_add(c)))
95        .and_then(|x| offset.checked_add(x)).ok_or(SolverError::SizeOverflow)?;
96        if last >= data.len() {
97            return Err(SolverError::Shape("strided view exceeds backing slice"));
98        }
99        let view = Self {
100            data, rows, cols, offset, row_stride, col_stride
101        };
102        for i in 0..rows {
103            for j in 0..cols {
104                if !view.at(i, j).is_finite() {
105                    return Err(SolverError::NonFinite {
106                        index: i*cols+j
107                    });
108                }
109            }
110        }
111        Ok(view)
112    }
113    pub fn rows(self) -> usize {
114        self.rows
115    }
116    pub fn columns(self) -> usize {
117        self.cols
118    }
119    pub fn get(self, row: usize, column: usize) -> Option<f64> {
120        (row < self.rows && column < self.cols).then(|| self.at(row, column))
121    }
122    pub fn transpose(self) -> Self {
123        Self {
124            rows: self.cols, cols: self.rows, row_stride: self.col_stride,
125            col_stride: self.row_stride, ..self
126        }
127    }
128    pub fn to_owned(self) -> Result<Matrix, SolverError> {
129        let mut out = Matrix::zeros(self.rows, self.cols)?;
130        for i in 0..self.rows {
131            for j in 0..self.cols {
132                out.data[i*self.cols+j] = self.at(i, j);
133            }
134        }
135        Ok(out)
136    }
137    pub(crate) fn at(self, row: usize, col: usize) -> f64 {
138        self.data[self.offset + row*self.row_stride + col*self.col_stride]
139    }
140    pub(crate) fn max_abs(self) -> f64 {
141        let mut max = 0.0f64;
142        for i in 0..self.rows {
143            for j in 0..self.cols {
144                max = max.max(self.at(i, j).abs());
145            }
146        }
147        max
148    }
149    pub(crate) fn square(self) -> Result<usize, SolverError> {
150        if self.rows != self.cols {
151            Err(SolverError::Shape("square matrix required"))
152        } else {
153            Ok(self.rows)
154        }
155    }
156}
157fn checked_shape(rows: usize, cols: usize) -> Result<usize, SolverError> {
158    if rows == 0 || cols == 0 {
159        return Err(SolverError::Shape("zero axes are not supported"));
160    }
161    rows.checked_mul(cols).filter(|x| *x <= isize::MAX as usize / std::mem::size_of::<f64>())
162    .ok_or(SolverError::SizeOverflow)
163}
rusolver/matrix.rs

rusolver/
matrix.rs
