clarabel/solver/chordal/decomp/
reverse_standard.rs

1// -----------------------------------
2//  reverse the standard decomposition
3// -----------------------------------
4
5use crate::solver::chordal::ChordalInfo;
6use crate::solver::DefaultVariables;
7use crate::{
8    algebra::*,
9    solver::SupportedConeT::{self},
10};
11use std::iter::zip;
12
13impl<T> ChordalInfo<T>
14where
15    T: FloatT,
16{
17    pub(crate) fn decomp_reverse_standard(
18        &self,
19        new_vars: &mut DefaultVariables<T>,
20        old_vars: &DefaultVariables<T>,
21        _old_cones: &[SupportedConeT<T>],
22    ) {
23        let H = &self.H.as_ref().unwrap();
24        let (_, m) = new_vars.dims();
25
26        H.gemv(&mut new_vars.s, &old_vars.s[m..], T::one(), T::zero());
27        H.gemv(&mut new_vars.z, &old_vars.z[m..], T::one(), T::zero());
28
29        //  to remove the overlaps we take the average of the values for
30        //  each overlap by dividing by the number of blocks that overlap
31        // in a particular entry, i.e. number of 1s in each row of H
32
33        let (rows, nnzs) = number_of_overlaps_in_rows(H);
34
35        for (ri, nnz) in zip(rows, nnzs) {
36            new_vars.z[ri] /= nnz;
37        }
38    }
39}
40
41fn number_of_overlaps_in_rows<T>(A: &CscMatrix<T>) -> (Vec<usize>, Vec<T>)
42where
43    T: FloatT,
44{
45    // sum the entries row-wise
46    let mut n_overlaps: Vec<T> = vec![T::zero(); A.m];
47    A.row_sums(&mut n_overlaps);
48    let ri = n_overlaps.iter().position_all(|&x| *x > T::one());
49
50    //n_overlaps <- n_overlaps[ri]
51    let n_overlaps: Vec<T> = ri.iter().map(|&i| n_overlaps[i]).collect();
52
53    (ri, n_overlaps)
54}