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 an error when the matrix dimensions are invalid (a shape mismatch)
124 /// /// or the initial pivot is out of range (an out of bounds failure).
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 an error when the pivot column is out of range (an out of bounds
305 /// /// failure) or the ACA update fails.
306 ///
307 /// # Examples
308 ///
309 /// ```
310 /// use tensor4all_core::MatrixACA;
311 /// use tensor4all_tensorbackend::from_vec2d;
312 ///
313 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
314 /// let mut aca = MatrixACA::<f64>::new(2, 2);
315 /// aca.add_pivot_col(&m, 0).unwrap();
316 /// assert_eq!(aca.u().ncols(), 1);
317 /// ```
318 pub fn add_pivot_col(&mut self, a: &Matrix<T>, col_index: usize) -> Result<()> {
319 if a.nrows() != self.nrows() || a.ncols() != self.ncols() {
320 return Err(MatrixCIError::DimensionMismatch {
321 expected_rows: self.nrows(),
322 expected_cols: self.ncols(),
323 actual_rows: a.nrows(),
324 actual_cols: a.ncols(),
325 });
326 }
327 if col_index >= self.ncols() {
328 return Err(MatrixCIError::IndexOutOfBounds {
329 row: 0,
330 col: col_index,
331 nrows: self.nrows(),
332 ncols: self.ncols(),
333 });
334 }
335
336 let uk = self.compute_uk(a, col_index)?;
337 self.col_indices.push(col_index);
338 self.u = append_col(&self.u, &uk);
339
340 Ok(())
341 }
342
343 /// Add a pivot row
344 ///
345 /// # Errors
346 ///
347 /// Returns an error when the pivot row is out of range (an out of bounds
348 /// /// failure) or the ACA update fails.
349 ///
350 /// # Examples
351 ///
352 /// ```
353 /// use tensor4all_core::MatrixACA;
354 /// use tensor4all_tensorbackend::from_vec2d;
355 ///
356 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
357 /// let mut aca = MatrixACA::<f64>::new(2, 2);
358 /// aca.add_pivot_col(&m, 0).unwrap();
359 /// aca.add_pivot_row(&m, 0).unwrap();
360 /// assert_eq!(aca.npivots(), 1);
361 /// assert_eq!(aca.v().nrows(), 1);
362 /// ```
363 pub fn add_pivot_row(&mut self, a: &Matrix<T>, row_index: usize) -> Result<()> {
364 if a.nrows() != self.nrows() || a.ncols() != self.ncols() {
365 return Err(MatrixCIError::DimensionMismatch {
366 expected_rows: self.nrows(),
367 expected_cols: self.ncols(),
368 actual_rows: a.nrows(),
369 actual_cols: a.ncols(),
370 });
371 }
372 if row_index >= self.nrows() {
373 return Err(MatrixCIError::IndexOutOfBounds {
374 row: row_index,
375 col: 0,
376 nrows: self.nrows(),
377 ncols: self.ncols(),
378 });
379 }
380
381 let vk = self.compute_vk(a, row_index)?;
382
383 // Compute alpha before mutating any state.
384 let xk = row_index;
385 let u_xk_last = self.u[[xk, self.u.ncols() - 1]];
386 if u_xk_last.abs_sq() < f64::EPSILON * f64::EPSILON {
387 return Err(MatrixCIError::SingularMatrix);
388 }
389 let alpha = T::one() / u_xk_last;
390 self.row_indices.push(row_index);
391 self.v = append_row(&self.v, &vk);
392 self.alpha.push(alpha);
393
394 Ok(())
395 }
396
397 /// Add a pivot at the given position
398 ///
399 /// # Errors
400 ///
401 /// Returns an error when the pivot is out of range (an out of bounds failure)
402 /// /// or the ACA update fails.
403 ///
404 /// # Examples
405 ///
406 /// ```
407 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
408 /// use tensor4all_tensorbackend::from_vec2d;
409 ///
410 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
411 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
412 /// aca.add_pivot(&m, (1, 1)).unwrap();
413 /// assert_eq!(aca.rank(), 2);
414 /// // Exact reconstruction at all positions
415 /// for i in 0..2 {
416 /// for j in 0..2 {
417 /// assert!((aca.evaluate(i, j) - m[[i, j]]).abs() < 1e-10);
418 /// }
419 /// }
420 /// ```
421 pub fn add_pivot(&mut self, a: &Matrix<T>, pivot: (usize, usize)) -> Result<()> {
422 let mut staged = self.clone();
423 staged.add_pivot_col(a, pivot.1)?;
424 staged.add_pivot_row(a, pivot.0)?;
425 *self = staged;
426 Ok(())
427 }
428
429 /// Add a pivot that maximizes the error using ACA heuristic
430 ///
431 /// # Errors
432 ///
433 /// Returns an error when no valid pivot exists (a singular failure) or the
434 /// /// ACA update fails.
435 ///
436 /// # Examples
437 ///
438 /// ```
439 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
440 /// use tensor4all_tensorbackend::from_vec2d;
441 ///
442 /// let m = from_vec2d(vec![
443 /// vec![1.0_f64, 2.0, 3.0],
444 /// vec![4.0, 5.0, 6.0],
445 /// vec![7.0, 8.0, 10.0],
446 /// ]);
447 /// let mut aca = MatrixACA::<f64>::new(3, 3);
448 /// let (r, c) = aca.add_best_pivot(&m).unwrap();
449 /// assert!(r < 3);
450 /// assert!(c < 3);
451 /// assert_eq!(aca.rank(), 1);
452 /// ```
453 pub fn add_best_pivot(&mut self, a: &Matrix<T>) -> Result<(usize, usize)> {
454 let mut staged = self.clone();
455 let pivot = staged.add_best_pivot_inner(a)?;
456 *self = staged;
457 Ok(pivot)
458 }
459
460 fn add_best_pivot_inner(&mut self, a: &Matrix<T>) -> Result<(usize, usize)> {
461 if a.nrows() != self.nrows() || a.ncols() != self.ncols() {
462 return Err(MatrixCIError::DimensionMismatch {
463 expected_rows: self.nrows(),
464 expected_cols: self.ncols(),
465 actual_rows: a.nrows(),
466 actual_cols: a.ncols(),
467 });
468 }
469 if self.is_empty() {
470 // Find global maximum
471 let (i, j, _) = submatrix_argmax(a, 0..self.nrows(), 0..self.ncols());
472 self.add_pivot(a, (i, j))?;
473 return Ok((i, j));
474 }
475
476 // ACA heuristic: find max in last row of V, then max in corresponding column of U
477 let avail_cols = self.available_cols();
478 if avail_cols.is_empty() {
479 return Err(MatrixCIError::FullRank);
480 }
481
482 // Find column with max |v[last, j]| among available columns
483 let last_row = self.v.nrows() - 1;
484 let mut max_val = self.v[[last_row, avail_cols[0]]].abs_sq();
485 let mut best_col = avail_cols[0];
486 for &c in &avail_cols {
487 let val = self.v[[last_row, c]].abs_sq();
488 if val > max_val {
489 max_val = val;
490 best_col = c;
491 }
492 }
493
494 self.add_pivot_col(a, best_col)?;
495
496 // Find row with max |u[i, last]| among available rows
497 let avail_rows = self.available_rows();
498 if avail_rows.is_empty() {
499 return Err(MatrixCIError::FullRank);
500 }
501
502 let last_col = self.u.ncols() - 1;
503 let mut max_val = self.u[[avail_rows[0], last_col]].abs_sq();
504 let mut best_row = avail_rows[0];
505 for &r in &avail_rows {
506 let val = self.u[[r, last_col]].abs_sq();
507 if val > max_val {
508 max_val = val;
509 best_row = r;
510 }
511 }
512
513 self.add_pivot_row(a, best_row)?;
514
515 Ok((best_row, best_col))
516 }
517
518 /// Set columns with new pivot rows and permutation
519 ///
520 /// Updates the column indices and V matrix according to the given
521 /// permutation and new pivot row data.
522 ///
523 /// # Errors
524 /// Returns [`MatrixCIError::InvalidArgument`] when the pivot-row shape or
525 /// permutation is incompatible with the current approximation.
526 ///
527 /// # Examples
528 ///
529 /// ```
530 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
531 /// use tensor4all_tensorbackend::from_vec2d;
532 ///
533 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
534 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
535 /// // Identity permutation keeps columns the same
536 /// let pivot_rows = from_vec2d(vec![vec![1.0_f64, 2.0]]);
537 /// aca.set_cols(&pivot_rows, &[0, 1]).unwrap();
538 /// assert_eq!(aca.col_indices(), &[0]);
539 /// ```
540 pub fn set_cols(&mut self, new_pivot_rows: &Matrix<T>, permutation: &[usize]) -> Result<()> {
541 let mut sorted_permutation = permutation.to_vec();
542 sorted_permutation.sort_unstable();
543 if new_pivot_rows.nrows() != self.v.nrows()
544 || permutation.len() != self.v.ncols()
545 || new_pivot_rows.ncols() < permutation.len()
546 || sorted_permutation
547 .windows(2)
548 .any(|window| window[0] == window[1])
549 || permutation
550 .iter()
551 .any(|&index| index >= new_pivot_rows.ncols())
552 || self
553 .col_indices
554 .iter()
555 .any(|&index| index >= permutation.len())
556 {
557 return Err(MatrixCIError::InvalidArgument {
558 message: "set_cols received incompatible pivot rows or permutation".to_string(),
559 });
560 }
561 let mut staged = self.clone();
562 staged.set_cols_inner(new_pivot_rows, permutation);
563 *self = staged;
564 Ok(())
565 }
566
567 fn set_cols_inner(&mut self, new_pivot_rows: &Matrix<T>, permutation: &[usize]) {
568 // Permute column indices
569 self.col_indices = self.col_indices.iter().map(|&c| permutation[c]).collect();
570
571 // Permute V matrix columns
572 let mut temp_v = Matrix::zeros(self.v.nrows(), new_pivot_rows.ncols());
573 for i in 0..self.v.nrows() {
574 for (new_j, &old_j) in permutation.iter().enumerate() {
575 if old_j < self.v.ncols() {
576 temp_v[[i, new_j]] = self.v[[i, old_j]];
577 }
578 }
579 }
580 self.v = temp_v;
581
582 // Insert new elements
583 let new_indices: Vec<usize> = (0..new_pivot_rows.ncols())
584 .filter(|j| !permutation.contains(j))
585 .collect();
586
587 for k in 0..new_pivot_rows.nrows() {
588 for &j in &new_indices {
589 let mut val = new_pivot_rows[[k, j]];
590 for l in 0..k {
591 let factor = self.u[[self.row_indices[k], l]] * self.alpha[l];
592 val = val - self.v[[l, j]] * factor;
593 }
594 self.v[[k, j]] = val;
595 }
596 }
597 }
598
599 /// Set rows with new pivot columns and permutation
600 ///
601 /// Updates the row indices and U matrix according to the given
602 /// permutation and new pivot column data.
603 ///
604 /// # Errors
605 /// Returns [`MatrixCIError::InvalidArgument`] when the pivot-column shape or
606 /// permutation is incompatible with the current approximation.
607 ///
608 /// # Examples
609 ///
610 /// ```
611 /// use tensor4all_core::{AbstractMatrixCI, MatrixACA};
612 /// use tensor4all_tensorbackend::from_vec2d;
613 ///
614 /// let m = from_vec2d(vec![vec![1.0_f64, 2.0], vec![3.0, 4.0]]);
615 /// let mut aca = MatrixACA::from_matrix_with_pivot(&m, (0, 0)).unwrap();
616 /// // Identity permutation keeps rows the same
617 /// let pivot_cols = from_vec2d(vec![vec![1.0_f64], vec![3.0]]);
618 /// aca.set_rows(&pivot_cols, &[0, 1]).unwrap();
619 /// assert_eq!(aca.row_indices(), &[0]);
620 /// ```
621 pub fn set_rows(&mut self, new_pivot_cols: &Matrix<T>, permutation: &[usize]) -> Result<()> {
622 let mut sorted_permutation = permutation.to_vec();
623 sorted_permutation.sort_unstable();
624 if new_pivot_cols.ncols() != self.u.ncols()
625 || permutation.len() != self.u.nrows()
626 || new_pivot_cols.nrows() < permutation.len()
627 || sorted_permutation
628 .windows(2)
629 .any(|window| window[0] == window[1])
630 || permutation
631 .iter()
632 .any(|&index| index >= new_pivot_cols.nrows())
633 || self
634 .row_indices
635 .iter()
636 .any(|&index| index >= permutation.len())
637 {
638 return Err(MatrixCIError::InvalidArgument {
639 message: "set_rows received incompatible pivot columns or permutation".to_string(),
640 });
641 }
642 let mut staged = self.clone();
643 staged.set_rows_inner(new_pivot_cols, permutation);
644 *self = staged;
645 Ok(())
646 }
647
648 fn set_rows_inner(&mut self, new_pivot_cols: &Matrix<T>, permutation: &[usize]) {
649 // Permute row indices
650 self.row_indices = self.row_indices.iter().map(|&r| permutation[r]).collect();
651
652 // Permute U matrix rows
653 let mut temp_u = Matrix::zeros(new_pivot_cols.nrows(), self.u.ncols());
654 for (new_i, &old_i) in permutation.iter().enumerate() {
655 if old_i < self.u.nrows() {
656 for j in 0..self.u.ncols() {
657 temp_u[[new_i, j]] = self.u[[old_i, j]];
658 }
659 }
660 }
661 self.u = temp_u;
662
663 // Insert new elements
664 let new_indices: Vec<usize> = (0..new_pivot_cols.nrows())
665 .filter(|i| !permutation.contains(i))
666 .collect();
667
668 for k in 0..new_pivot_cols.ncols() {
669 for &i in &new_indices {
670 let mut val = new_pivot_cols[[i, k]];
671 for l in 0..k {
672 let factor = self.v[[l, self.col_indices[k]]] * self.alpha[l];
673 val = val - self.u[[i, l]] * factor;
674 }
675 self.u[[i, k]] = val;
676 }
677 }
678 }
679}
680
681impl<T: Scalar> AbstractMatrixCI<T> for MatrixACA<T> {
682 fn nrows(&self) -> usize {
683 self.u.nrows()
684 }
685
686 fn ncols(&self) -> usize {
687 self.v.ncols()
688 }
689
690 fn rank(&self) -> usize {
691 self.row_indices.len()
692 }
693
694 fn row_indices(&self) -> &[usize] {
695 &self.row_indices
696 }
697
698 fn col_indices(&self) -> &[usize] {
699 &self.col_indices
700 }
701
702 fn evaluate(&self, i: usize, j: usize) -> T {
703 if self.is_empty() {
704 return T::zero();
705 }
706
707 let mut sum = T::zero();
708 for k in 0..self.rank() {
709 sum = sum + self.u[[i, k]] * self.alpha[k] * self.v[[k, j]];
710 }
711 sum
712 }
713
714 fn submatrix(&self, rows: &[usize], cols: &[usize]) -> Matrix<T> {
715 if self.is_empty() {
716 return Matrix::zeros(rows.len(), cols.len());
717 }
718
719 let r = self.rank();
720 let mut result = Matrix::zeros(rows.len(), cols.len());
721 for (j_out, &col) in cols.iter().enumerate() {
722 for (i_out, &row) in rows.iter().enumerate() {
723 let mut sum = T::zero();
724 for k in 0..r {
725 sum = sum + self.u[[row, k]] * self.alpha[k] * self.v[[k, col]];
726 }
727 result[[i_out, j_out]] = sum;
728 }
729 }
730 result
731 }
732}
733
734#[cfg(test)]
735mod tests;