tensor4all_core/matrixaca.rs
1//! Adaptive Cross Approximation (ACA) implementation.
2//!
3//! [`MatrixACA`] provides the ACA algorithm for incrementally building a
4//! low-rank approximation of a matrix. It stores the approximation in
5//! factored form `A ~ U * diag(alpha) * V` and supports adding pivots
6//! one at a time.
7
8use crate::error::{MatrixCIError, Result};
9use crate::scalar::Scalar;
10use crate::traits::AbstractMatrixCI;
11use tensor4all_tensorbackend::{submatrix_argmax, Matrix};
12
13fn matrix_row<T: Scalar>(m: &Matrix<T>, i: usize) -> Vec<T> {
14 (0..m.ncols()).map(|j| m[[i, j]]).collect()
15}
16
17fn matrix_col<T: Scalar>(m: &Matrix<T>, j: usize) -> Vec<T> {
18 (0..m.nrows()).map(|i| m[[i, j]]).collect()
19}
20
21fn append_col<T: Scalar>(m: &Matrix<T>, col: &[T]) -> Matrix<T> {
22 assert_eq!(col.len(), m.nrows());
23
24 let mut result = Matrix::zeros(m.nrows(), m.ncols() + 1);
25 for i in 0..m.nrows() {
26 for j in 0..m.ncols() {
27 result[[i, j]] = m[[i, j]];
28 }
29 result[[i, m.ncols()]] = col[i];
30 }
31 result
32}
33
34fn append_row<T: Scalar>(m: &Matrix<T>, row: &[T]) -> Matrix<T> {
35 assert_eq!(row.len(), m.ncols());
36
37 let mut result = Matrix::zeros(m.nrows() + 1, m.ncols());
38 for i in 0..m.nrows() {
39 for j in 0..m.ncols() {
40 result[[i, j]] = m[[i, j]];
41 }
42 }
43 for j in 0..m.ncols() {
44 result[[m.nrows(), j]] = row[j];
45 }
46 result
47}
48
49/// Adaptive Cross Approximation representation.
50///
51/// Stores a low-rank approximation `A ~ U * diag(alpha) * V`, where
52/// `alpha[k] = 1/delta[k]`. Implements [`AbstractMatrixCI`] for
53/// element-wise evaluation and submatrix extraction.
54///
55/// # Examples
56///
57/// ```
58/// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
59/// use tensor4all_tensorbackend::from_vec2d;
60///
61/// let m = from_vec2d(vec![
62/// vec![1.0_f64, 2.0, 3.0],
63/// vec![4.0, 5.0, 6.0],
64/// vec![7.0, 8.0, 9.0],
65/// ]);
66///
67/// // Build ACA from a starting pivot
68/// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (1, 1)).unwrap();
69/// assert_eq!(aca.rank(), 1);
70///
71/// // Add more pivots
72/// aca.add_pivot(&m, (0, 0)).unwrap();
73/// assert_eq!(aca.rank(), 2);
74///
75/// // Evaluate the approximation
76/// let approx_11 = aca.evaluate(1, 1);
77/// assert!((approx_11 - 5.0).abs() < 1e-10);
78/// ```
79#[derive(Debug, Clone)]
80pub struct MatrixACA<T: Scalar> {
81 /// Row indices (I set)
82 row_indices: Vec<usize>,
83 /// Column indices (J set)
84 col_indices: Vec<usize>,
85 /// U matrix: u_k(x) for all x, shape (nrows, npivots)
86 u: Matrix<T>,
87 /// V matrix: v_k(y) for all y, shape (npivots, ncols)
88 v: Matrix<T>,
89 /// Alpha values: 1/delta for each pivot
90 alpha: Vec<T>,
91}
92
93impl<T: Scalar> MatrixACA<T> {
94 /// Create an empty MatrixACA for a matrix of given size
95 ///
96 /// # Examples
97 ///
98 /// ```
99 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
100 ///
101 /// let aca = MatrixACA::<f64>::new(3, 4);
102 /// assert_eq!(aca.nrows(), 3);
103 /// assert_eq!(aca.ncols(), 4);
104 /// assert_eq!(aca.rank(), 0);
105 /// assert!(aca.is_empty());
106 /// ```
107 pub fn new(nr: usize, nc: usize) -> Self {
108 Self {
109 row_indices: Vec::new(),
110 col_indices: Vec::new(),
111 u: Matrix::zeros(nr, 0),
112 v: Matrix::zeros(0, nc),
113 alpha: Vec::new(),
114 }
115 }
116
117 /// Create a MatrixACA from a matrix with an initial pivot
118 ///
119 /// Returns an error if the pivot value is zero or near-zero.
120 ///
121 /// # Errors
122 ///
123 /// Returns [`MatrixCIError::IndexOutOfBounds`] when `first_pivot` is outside
124 /// `a`, or [`MatrixCIError::SingularMatrix`] when its value is near-zero.
125 ///
126 /// # Examples
127 ///
128 /// ```
129 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
130 /// use tensor4all_tensorbackend::from_vec2d;
131 ///
132 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
133 /// let aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
134 /// assert_eq!(aca.rank(), 1);
135 /// assert_eq!(aca.row_indices(), &[0]);
136 /// assert_eq!(aca.col_indices(), &[0]);
137 /// ```
138 pub fn from_matrix_with_pivot(a: &Matrix<T>, first_pivot: (usize, usize)) -> Result<Self> {
139 let (i, j) = first_pivot;
140 if i >= a.nrows() || j >= a.ncols() {
141 return Err(MatrixCIError::IndexOutOfBounds {
142 row: i,
143 col: j,
144 nrows: a.nrows(),
145 ncols: a.ncols(),
146 });
147 }
148 let pivot_val = a[[i, j]];
149
150 if pivot_val.abs_sq() < f64::EPSILON * f64::EPSILON {
151 return Err(MatrixCIError::SingularMatrix);
152 }
153
154 // u = A[:, j]
155 let mut u = Matrix::zeros(a.nrows(), 1);
156 for r in 0..a.nrows() {
157 u[[r, 0]] = a[[r, j]];
158 }
159
160 // v = A[i, :]
161 let mut v = Matrix::zeros(1, a.ncols());
162 for c in 0..a.ncols() {
163 v[[0, c]] = a[[i, c]];
164 }
165
166 let alpha = vec![T::one() / pivot_val];
167
168 Ok(Self {
169 row_indices: vec![i],
170 col_indices: vec![j],
171 u,
172 v,
173 alpha,
174 })
175 }
176
177 /// Number of pivots
178 ///
179 /// # Examples
180 ///
181 /// ```
182 /// use tensor4all_core::MatrixACA;
183 /// use tensor4all_tensorbackend::from_vec2d;
184 ///
185 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
186 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
187 /// assert_eq!(aca.npivots(), 1);
188 /// aca.add_pivot(&m, (1, 1)).unwrap();
189 /// assert_eq!(aca.npivots(), 2);
190 /// ```
191 pub fn npivots(&self) -> usize {
192 self.alpha.len()
193 }
194
195 /// Get reference to U matrix
196 ///
197 /// # Examples
198 ///
199 /// ```
200 /// use tensor4all_core::MatrixACA;
201 /// use tensor4all_tensorbackend::from_vec2d;
202 ///
203 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
204 /// let aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
205 /// let u = aca.u();
206 /// assert_eq!(u.nrows(), 2); // nrows of original matrix
207 /// assert_eq!(u.ncols(), 1); // one pivot
208 /// ```
209 pub fn u(&self) -> &Matrix<T> {
210 &self.u
211 }
212
213 /// Get reference to V matrix
214 ///
215 /// # Examples
216 ///
217 /// ```
218 /// use tensor4all_core::MatrixACA;
219 /// use tensor4all_tensorbackend::from_vec2d;
220 ///
221 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
222 /// let aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
223 /// let v = aca.v();
224 /// assert_eq!(v.nrows(), 1); // one pivot
225 /// assert_eq!(v.ncols(), 2); // ncols of original matrix
226 /// ```
227 pub fn v(&self) -> &Matrix<T> {
228 &self.v
229 }
230
231 /// Get reference to alpha values
232 ///
233 /// # Examples
234 ///
235 /// ```
236 /// use tensor4all_core::MatrixACA;
237 /// use tensor4all_tensorbackend::from_vec2d;
238 ///
239 /// let m = from_vec2d(vec![vec![2.0_f64, 1.0], vec![1.0, 3.0]]);
240 /// let aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
241 /// let alpha = aca.alpha();
242 /// assert_eq!(alpha.len(), 1);
243 /// // alpha = 1 / pivot_value = 1 / m[0,0] = 0.5
244 /// assert!((alpha[0] - 0.5).abs() < 1e-10);
245 /// ```
246 pub fn alpha(&self) -> &[T] {
247 &self.alpha
248 }
249
250 /// Compute u_k(x) for all x (the k-th column of U after adding a new pivot column)
251 fn compute_uk(&self, a: &Matrix<T>, yk: usize) -> Result<Vec<T>> {
252 let k = self.col_indices.len();
253
254 // result = A[:, yk]
255 let mut result: Vec<T> = matrix_col(a, yk);
256
257 // Subtract contribution from previous pivots
258 for l in 0..k {
259 let xl = self.row_indices[l];
260 let v_l_yk = self.v[[l, yk]];
261 let u_xl_l = self.u[[xl, l]];
262 if u_xl_l.abs_sq() < f64::EPSILON * f64::EPSILON {
263 return Err(MatrixCIError::SingularMatrix);
264 }
265 let factor = v_l_yk / u_xl_l;
266
267 for (i, res) in result.iter_mut().enumerate() {
268 *res = *res - factor * self.u[[i, l]];
269 }
270 }
271
272 Ok(result)
273 }
274
275 /// Compute v_k(y) for all y (the k-th row of V after adding a new pivot row)
276 fn compute_vk(&self, a: &Matrix<T>, xk: usize) -> Result<Vec<T>> {
277 let k = self.row_indices.len();
278
279 // result = A[xk, :]
280 let mut result: Vec<T> = matrix_row(a, xk);
281
282 // Subtract contribution from previous pivots
283 for l in 0..k {
284 let xl = self.row_indices[l];
285 let u_xk_l = self.u[[xk, l]];
286 let u_xl_l = self.u[[xl, l]];
287 if u_xl_l.abs_sq() < f64::EPSILON * f64::EPSILON {
288 return Err(MatrixCIError::SingularMatrix);
289 }
290 let factor = u_xk_l / u_xl_l;
291
292 for (j, res) in result.iter_mut().enumerate() {
293 *res = *res - factor * self.v[[l, j]];
294 }
295 }
296
297 Ok(result)
298 }
299
300 /// Add a pivot column
301 ///
302 /// # Errors
303 ///
304 /// Returns [`MatrixCIError::DimensionMismatch`] when `a` differs from this
305 /// approximation's shape, [`MatrixCIError::IndexOutOfBounds`] when
306 /// `col_index` is outside the column count, or
307 /// [`MatrixCIError::SingularMatrix`] when the update hits a zero pivot.
308 ///
309 /// # Examples
310 ///
311 /// ```
312 /// use tensor4all_core::MatrixACA;
313 /// use tensor4all_tensorbackend::from_vec2d;
314 ///
315 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
316 /// let mut aca = MatrixACA::<f64>::new(2, 2);
317 /// aca.add_pivot_col(&m, 0).unwrap();
318 /// assert_eq!(aca.u().ncols(), 1);
319 /// ```
320 pub fn add_pivot_col(&mut self, a: &Matrix<T>, col_index: usize) -> Result<()> {
321 if a.nrows() != self.nrows() || a.ncols() != self.ncols() {
322 return Err(MatrixCIError::DimensionMismatch {
323 expected_rows: self.nrows(),
324 expected_cols: self.ncols(),
325 actual_rows: a.nrows(),
326 actual_cols: a.ncols(),
327 });
328 }
329 if col_index >= self.ncols() {
330 return Err(MatrixCIError::IndexOutOfBounds {
331 row: 0,
332 col: col_index,
333 nrows: self.nrows(),
334 ncols: self.ncols(),
335 });
336 }
337
338 let uk = self.compute_uk(a, col_index)?;
339 self.col_indices.push(col_index);
340 self.u = append_col(&self.u, &uk);
341
342 Ok(())
343 }
344
345 /// Add a pivot row
346 ///
347 /// # Errors
348 ///
349 /// Returns [`MatrixCIError::DimensionMismatch`] when `a` differs from this
350 /// approximation's shape, [`MatrixCIError::IndexOutOfBounds`] when
351 /// `row_index` is outside the row count, or
352 /// [`MatrixCIError::SingularMatrix`] when the update hits a zero pivot.
353 ///
354 /// # Examples
355 ///
356 /// ```
357 /// use tensor4all_core::MatrixACA;
358 /// use tensor4all_tensorbackend::from_vec2d;
359 ///
360 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
361 /// let mut aca = MatrixACA::<f64>::new(2, 2);
362 /// aca.add_pivot_col(&m, 0).unwrap();
363 /// aca.add_pivot_row(&m, 0).unwrap();
364 /// assert_eq!(aca.npivots(), 1);
365 /// assert_eq!(aca.v().nrows(), 1);
366 /// ```
367 pub fn add_pivot_row(&mut self, a: &Matrix<T>, row_index: usize) -> Result<()> {
368 if a.nrows() != self.nrows() || a.ncols() != self.ncols() {
369 return Err(MatrixCIError::DimensionMismatch {
370 expected_rows: self.nrows(),
371 expected_cols: self.ncols(),
372 actual_rows: a.nrows(),
373 actual_cols: a.ncols(),
374 });
375 }
376 if row_index >= self.nrows() {
377 return Err(MatrixCIError::IndexOutOfBounds {
378 row: row_index,
379 col: 0,
380 nrows: self.nrows(),
381 ncols: self.ncols(),
382 });
383 }
384
385 let vk = self.compute_vk(a, row_index)?;
386
387 // Compute alpha before mutating any state.
388 let xk = row_index;
389 let u_xk_last = self.u[[xk, self.u.ncols() - 1]];
390 if u_xk_last.abs_sq() < f64::EPSILON * f64::EPSILON {
391 return Err(MatrixCIError::SingularMatrix);
392 }
393 let alpha = T::one() / u_xk_last;
394 self.row_indices.push(row_index);
395 self.v = append_row(&self.v, &vk);
396 self.alpha.push(alpha);
397
398 Ok(())
399 }
400
401 /// Add a pivot at the given position
402 ///
403 /// # Errors
404 ///
405 /// Returns the same [`MatrixCIError`] variants as
406 /// [`MatrixACA::add_pivot_col`] and [`MatrixACA::add_pivot_row`], and leaves
407 /// this approximation unchanged on failure.
408 ///
409 /// # Examples
410 ///
411 /// ```
412 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
413 /// use tensor4all_tensorbackend::from_vec2d;
414 ///
415 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
416 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
417 /// aca.add_pivot(&m, (1, 1)).unwrap();
418 /// assert_eq!(aca.rank(), 2);
419 /// // Exact reconstruction at all positions
420 /// for i in 0..2 {
421 /// for j in 0..2 {
422 /// assert!((aca.evaluate(i, j) - m[[i, j]]).abs() < 1e-10);
423 /// }
424 /// }
425 /// ```
426 pub fn add_pivot(&mut self, a: &Matrix<T>, pivot: (usize, usize)) -> Result<()> {
427 let mut staged = self.clone();
428 staged.add_pivot_col(a, pivot.1)?;
429 staged.add_pivot_row(a, pivot.0)?;
430 *self = staged;
431 Ok(())
432 }
433
434 /// Add a pivot that maximizes the error using ACA heuristic
435 ///
436 /// # Errors
437 ///
438 /// Returns [`MatrixCIError::DimensionMismatch`] when `a` differs from this
439 /// approximation's shape, [`MatrixCIError::FullRank`] when no available row
440 /// or column remains, or [`MatrixCIError::SingularMatrix`] when the update
441 /// hits a zero pivot.
442 ///
443 /// # Examples
444 ///
445 /// ```
446 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
447 /// use tensor4all_tensorbackend::from_vec2d;
448 ///
449 /// let m = from_vec2d(vec![
450 /// vec![1.0_f64, 2.0, 3.0],
451 /// vec![4.0, 5.0, 6.0],
452 /// vec![7.0, 8.0, 10.0],
453 /// ]);
454 /// let mut aca = MatrixACA::<f64>::new(3, 3);
455 /// let (r, c) = aca.add_best_pivot(&m).unwrap();
456 /// assert!(r < 3);
457 /// assert!(c < 3);
458 /// assert_eq!(aca.rank(), 1);
459 /// ```
460 pub fn add_best_pivot(&mut self, a: &Matrix<T>) -> Result<(usize, usize)> {
461 let mut staged = self.clone();
462 let pivot = staged.add_best_pivot_inner(a)?;
463 *self = staged;
464 Ok(pivot)
465 }
466
467 fn add_best_pivot_inner(&mut self, a: &Matrix<T>) -> Result<(usize, usize)> {
468 if a.nrows() != self.nrows() || a.ncols() != self.ncols() {
469 return Err(MatrixCIError::DimensionMismatch {
470 expected_rows: self.nrows(),
471 expected_cols: self.ncols(),
472 actual_rows: a.nrows(),
473 actual_cols: a.ncols(),
474 });
475 }
476 if self.is_empty() {
477 // Find global maximum
478 let (i, j, _) = submatrix_argmax(a, 0..self.nrows(), 0..self.ncols());
479 self.add_pivot(a, (i, j))?;
480 return Ok((i, j));
481 }
482
483 // ACA heuristic: find max in last row of V, then max in corresponding column of U
484 let avail_cols = self.available_cols();
485 if avail_cols.is_empty() {
486 return Err(MatrixCIError::FullRank);
487 }
488
489 // Find column with max |v[last, j]| among available columns
490 let last_row = self.v.nrows() - 1;
491 let mut max_val = self.v[[last_row, avail_cols[0]]].abs_sq();
492 let mut best_col = avail_cols[0];
493 for &c in &avail_cols {
494 let val = self.v[[last_row, c]].abs_sq();
495 if val > max_val {
496 max_val = val;
497 best_col = c;
498 }
499 }
500
501 self.add_pivot_col(a, best_col)?;
502
503 // Find row with max |u[i, last]| among available rows
504 let avail_rows = self.available_rows();
505 if avail_rows.is_empty() {
506 return Err(MatrixCIError::FullRank);
507 }
508
509 let last_col = self.u.ncols() - 1;
510 let mut max_val = self.u[[avail_rows[0], last_col]].abs_sq();
511 let mut best_row = avail_rows[0];
512 for &r in &avail_rows {
513 let val = self.u[[r, last_col]].abs_sq();
514 if val > max_val {
515 max_val = val;
516 best_row = r;
517 }
518 }
519
520 self.add_pivot_row(a, best_row)?;
521
522 Ok((best_row, best_col))
523 }
524
525 /// Set columns with new pivot rows and permutation
526 ///
527 /// Updates the column indices and V matrix according to the given
528 /// permutation and new pivot row data.
529 ///
530 /// # Errors
531 /// Returns [`MatrixCIError::InvalidArgument`] when the pivot-row shape or
532 /// permutation is incompatible with the current approximation.
533 ///
534 /// # Examples
535 ///
536 /// ```
537 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
538 /// use tensor4all_tensorbackend::from_vec2d;
539 ///
540 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
541 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
542 /// // Identity permutation keeps columns the same
543 /// let pivot_rows = from_vec2d(vec![vec![1.0_f64, 2.0]]);
544 /// aca.set_cols(&pivot_rows, &[0, 1]).unwrap();
545 /// assert_eq!(aca.col_indices(), &[0]);
546 /// ```
547 pub fn set_cols(&mut self, new_pivot_rows: &Matrix<T>, permutation: &[usize]) -> Result<()> {
548 let mut sorted_permutation = permutation.to_vec();
549 sorted_permutation.sort_unstable();
550 if new_pivot_rows.nrows() != self.v.nrows()
551 || permutation.len() != self.v.ncols()
552 || new_pivot_rows.ncols() < permutation.len()
553 || sorted_permutation
554 .windows(2)
555 .any(|window| window[0] == window[1])
556 || permutation
557 .iter()
558 .any(|&index| index >= new_pivot_rows.ncols())
559 || self
560 .col_indices
561 .iter()
562 .any(|&index| index >= permutation.len())
563 {
564 return Err(MatrixCIError::InvalidArgument {
565 message: "set_cols received incompatible pivot rows or permutation".to_string(),
566 });
567 }
568 let mut staged = self.clone();
569 staged.set_cols_inner(new_pivot_rows, permutation);
570 *self = staged;
571 Ok(())
572 }
573
574 fn set_cols_inner(&mut self, new_pivot_rows: &Matrix<T>, permutation: &[usize]) {
575 // Permute column indices
576 self.col_indices = self.col_indices.iter().map(|&c| permutation[c]).collect();
577
578 // Permute V matrix columns
579 let mut temp_v = Matrix::zeros(self.v.nrows(), new_pivot_rows.ncols());
580 for i in 0..self.v.nrows() {
581 for (new_j, &old_j) in permutation.iter().enumerate() {
582 if old_j < self.v.ncols() {
583 temp_v[[i, new_j]] = self.v[[i, old_j]];
584 }
585 }
586 }
587 self.v = temp_v;
588
589 // Insert new elements
590 let new_indices: Vec<usize> = (0..new_pivot_rows.ncols())
591 .filter(|j| !permutation.contains(j))
592 .collect();
593
594 for k in 0..new_pivot_rows.nrows() {
595 for &j in &new_indices {
596 let mut val = new_pivot_rows[[k, j]];
597 for l in 0..k {
598 let factor = self.u[[self.row_indices[k], l]] * self.alpha[l];
599 val = val - self.v[[l, j]] * factor;
600 }
601 self.v[[k, j]] = val;
602 }
603 }
604 }
605
606 /// Set rows with new pivot columns and permutation
607 ///
608 /// Updates the row indices and U matrix according to the given
609 /// permutation and new pivot column data.
610 ///
611 /// # Errors
612 /// Returns [`MatrixCIError::InvalidArgument`] when the pivot-column shape or
613 /// permutation is incompatible with the current approximation.
614 ///
615 /// # Examples
616 ///
617 /// ```
618 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
619 /// use tensor4all_tensorbackend::from_vec2d;
620 ///
621 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
622 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
623 /// // Identity permutation keeps rows the same
624 /// let pivot_cols = from_vec2d(vec![vec![1.0_f64], vec![3.0]]);
625 /// aca.set_rows(&pivot_cols, &[0, 1]).unwrap();
626 /// assert_eq!(aca.row_indices(), &[0]);
627 /// ```
628 pub fn set_rows(&mut self, new_pivot_cols: &Matrix<T>, permutation: &[usize]) -> Result<()> {
629 let mut sorted_permutation = permutation.to_vec();
630 sorted_permutation.sort_unstable();
631 if new_pivot_cols.ncols() != self.u.ncols()
632 || permutation.len() != self.u.nrows()
633 || new_pivot_cols.nrows() < permutation.len()
634 || sorted_permutation
635 .windows(2)
636 .any(|window| window[0] == window[1])
637 || permutation
638 .iter()
639 .any(|&index| index >= new_pivot_cols.nrows())
640 || self
641 .row_indices
642 .iter()
643 .any(|&index| index >= permutation.len())
644 {
645 return Err(MatrixCIError::InvalidArgument {
646 message: "set_rows received incompatible pivot columns or permutation".to_string(),
647 });
648 }
649 let mut staged = self.clone();
650 staged.set_rows_inner(new_pivot_cols, permutation);
651 *self = staged;
652 Ok(())
653 }
654
655 fn set_rows_inner(&mut self, new_pivot_cols: &Matrix<T>, permutation: &[usize]) {
656 // Permute row indices
657 self.row_indices = self.row_indices.iter().map(|&r| permutation[r]).collect();
658
659 // Permute U matrix rows
660 let mut temp_u = Matrix::zeros(new_pivot_cols.nrows(), self.u.ncols());
661 for (new_i, &old_i) in permutation.iter().enumerate() {
662 if old_i < self.u.nrows() {
663 for j in 0..self.u.ncols() {
664 temp_u[[new_i, j]] = self.u[[old_i, j]];
665 }
666 }
667 }
668 self.u = temp_u;
669
670 // Insert new elements
671 let new_indices: Vec<usize> = (0..new_pivot_cols.nrows())
672 .filter(|i| !permutation.contains(i))
673 .collect();
674
675 for k in 0..new_pivot_cols.ncols() {
676 for &i in &new_indices {
677 let mut val = new_pivot_cols[[i, k]];
678 for l in 0..k {
679 let factor = self.v[[l, self.col_indices[k]]] * self.alpha[l];
680 val = val - self.u[[i, l]] * factor;
681 }
682 self.u[[i, k]] = val;
683 }
684 }
685 }
686}
687
688impl<T: Scalar> AbstractMatrixCI<T> for MatrixACA<T> {
689 fn nrows(&self) -> usize {
690 self.u.nrows()
691 }
692
693 fn ncols(&self) -> usize {
694 self.v.ncols()
695 }
696
697 fn rank(&self) -> usize {
698 self.row_indices.len()
699 }
700
701 fn row_indices(&self) -> &[usize] {
702 &self.row_indices
703 }
704
705 fn col_indices(&self) -> &[usize] {
706 &self.col_indices
707 }
708
709 fn evaluate(&self, i: usize, j: usize) -> T {
710 if self.is_empty() {
711 return T::zero();
712 }
713
714 let mut sum = T::zero();
715 for k in 0..self.rank() {
716 sum = sum + self.u[[i, k]] * self.alpha[k] * self.v[[k, j]];
717 }
718 sum
719 }
720
721 fn submatrix(&self, rows: &[usize], cols: &[usize]) -> Matrix<T> {
722 if self.is_empty() {
723 return Matrix::zeros(rows.len(), cols.len());
724 }
725
726 let r = self.rank();
727 let mut result = Matrix::zeros(rows.len(), cols.len());
728 for (j_out, &col) in cols.iter().enumerate() {
729 for (i_out, &row) in rows.iter().enumerate() {
730 let mut sum = T::zero();
731 for k in 0..r {
732 sum = sum + self.u[[row, k]] * self.alpha[k] * self.v[[k, col]];
733 }
734 result[[i_out, j_out]] = sum;
735 }
736 }
737 result
738 }
739}
740
741#[cfg(test)]
742mod tests;