pub fn ldl_numeric<N, I, PStorage>(mat: CsMatViewI<N, I>,
                                   l_colptr: &[I],
                                   parents: linalg::etree::ParentsView,
                                   perm: &Permutation<I, PStorage>,
                                   l_nz: &mut [I],
                                   l_indices: &mut [I],
                                   l_data: &mut [N],
                                   diag: &mut [N],
                                   y_workspace: &mut [N],
                                   pattern_workspace: &mut DStack<I>,
                                   flag_workspace: &mut [I])
where N: Clone + Copy + PartialEq + Num + PartialOrd,
      I: SpIndex,
      PStorage: Deref<Target = [I]>
{
    let outer_it = mat.outer_iterator_perm(perm.view());
    for (k, (_, vec)) in outer_it.enumerate() {

        // compute the nonzero pattern of the kth row of L
        // in topological order

        flag_workspace[k] = I::from_usize(k); // this node is visited
        y_workspace[k] = N::zero();
        l_nz[k] = I::zero();
        pattern_workspace.clear_right();

        for (inner_ind, &val) in vec.iter_perm(perm.inv())
                                    .filter(|&(i, _)| i <= k) {
            y_workspace[inner_ind] = y_workspace[inner_ind] + val;
            let mut i = inner_ind;
            pattern_workspace.clear_left();
            while flag_workspace[i].index() != k {
                pattern_workspace.push_left(I::from_usize(i));
                flag_workspace[i] = I::from_usize(k);
                i = parents.get_parent(i).expect("enforced by ldl_symbolic");
            }
            pattern_workspace.push_left_on_right();
        }

        // use a sparse triangular solve to compute the values
        // of the kth row of L
        diag[k] = y_workspace[k];
        y_workspace[k] = N::zero();
        'pattern: for &i in pattern_workspace.iter_right() {
            let i = i.index();
            let yi = y_workspace[i];
            y_workspace[i] = N::zero();
            let p2 = l_colptr[i] + l_nz[i];
            for p in l_colptr[i].index()..p2.index() {
                // we cannot go inside this loop before something has actually
                // be written into l_indices[l_colptr[i]..p2] so this
                // read is actually not into garbage
                // actually each iteration of the 'pattern loop adds writes the
                // value in l_indices that will be read on the next iteration
                // TODO: can some design change make this fact more obvious?
                let y_index = l_indices[p].index();
                y_workspace[y_index] = y_workspace[y_index] - l_data[p] * yi;
            }
            let l_ki = yi / diag[i];
            diag[k] = diag[k] - l_ki * yi;
            l_indices[p2.index()] = I::from_usize(k);
            l_data[p2.index()] = l_ki;
            l_nz[i] += I::one();
        }
        if diag[k] == N::zero() {
            panic!("Matrix is singular");
        }
    }
}
/// Perform numeric LDLT decomposition
///
/// pattern_workspace is a DStack of capacity n
Perform numeric LDLT decomposition
pattern_workspace is a DStack of capacity n
Perform numeric LDLT decomposition pattern_workspace is a DStack of capacity n
pub fn ldl_numeric<N, I, PStorage>(mat: CsMatViewI<N, I>,
                                   l_colptr: &[I],
                                   parents: linalg::etree::ParentsView,
                                   perm: &Permutation<I, PStorage>,
                                   l_nz: &mut [I],
                                   l_indices: &mut [I],
                                   l_data: &mut [N],
                                   diag: &mut [N],
                                   y_workspace: &mut [N],
                                   pattern_workspace: &mut DStack<I>,
                                   flag_workspace: &mut [I])
where N: Clone + Copy + PartialEq + Num + PartialOrd,
      I: SpIndex,
      PStorage: Deref<Target = [I]>
{
    let outer_it = mat.outer_iterator_perm(perm.view());
    for (k, (_, vec)) in outer_it.enumerate() {
        // compute the nonzero pattern of the kth row of L
        // in topological order
        flag_workspace[k] = I::from_usize(k); // this node is visited
        y_workspace[k] = N::zero();
        l_nz[k] = I::zero();
        pattern_workspace.clear_right();
        for (inner_ind, &val) in vec.iter_perm(perm.inv())
                                    .filter(|&(i, _)| i <= k) {
            y_workspace[inner_ind] = y_workspace[inner_ind] + val;
            let mut i = inner_ind;
            pattern_workspace.clear_left();
            while flag_workspace[i].index() != k {
                pattern_workspace.push_left(I::from_usize(i));
                flag_workspace[i] = I::from_usize(k);
                i = parents.get_parent(i).expect("enforced by ldl_symbolic");
            }
            pattern_workspace.push_left_on_right();
        }
        // use a sparse triangular solve to compute the values
        // of the kth row of L
        diag[k] = y_workspace[k];
        y_workspace[k] = N::zero();
        'pattern: for &i in pattern_workspace.iter_right() {
            let i = i.index();
            let yi = y_workspace[i];
            y_workspace[i] = N::zero();
            let p2 = l_colptr[i] + l_nz[i];
            for p in l_colptr[i].index()..p2.index() {
                // we cannot go inside this loop before something has actually
                // be written into l_indices[l_colptr[i]..p2] so this
                // read is actually not into garbage
                // actually each iteration of the 'pattern loop adds writes the
                // value in l_indices that will be read on the next iteration
                // TODO: can some design change make this fact more obvious?
                let y_index = l_indices[p].index();
                y_workspace[y_index] = y_workspace[y_index] - l_data[p] * yi;
            }
            let l_ki = yi / diag[i];
            diag[k] = diag[k] - l_ki * yi;
            l_indices[p2.index()] = I::from_usize(k);
            l_data[p2.index()] = l_ki;
            l_nz[i] += I::one();
        }
        if diag[k] == N::zero() {
            panic!("Matrix is singular");
        }
    }
}
pub fn ldl_numeric < N , I , PStorage > ( mat : CsMatViewI < N , I > , l_colptr : & [ I ] , parents : linalg :: etree :: ParentsView , perm : & Permutation < I , PStorage > , l_nz : & mut [ I ] , l_indices : & mut [ I ] , l_data : & mut [ N ] , diag : & mut [ N ] , y_workspace : & mut [ N ] , pattern_workspace : & mut DStack < I > , flag_workspace : & mut [ I ] ) where N : Clone + Copy + PartialEq + Num + PartialOrd , I : SpIndex , PStorage : Deref < Target = [ I ] > { let outer_it = mat . outer_iterator_perm ( perm . view ( ) ) ; for ( k , ( _ , vec ) ) in outer_it . enumerate ( ) { // compute the nonzero pattern of the kth row of L // in topological order flag_workspace [ k ] = I :: from_usize ( k ) ; // this node is visited y_workspace [ k ] = N :: zero ( ) ; l_nz [ k ] = I :: zero ( ) ; pattern_workspace . clear_right ( ) ; for ( inner_ind , & val ) in vec . iter_perm ( perm . inv ( ) ) . filter ( | & ( i , _ ) | i <= k ) { y_workspace [ inner_ind ] = y_workspace [ inner_ind ] + val ; let mut i = inner_ind ; pattern_workspace . clear_left ( ) ; while flag_workspace [ i ] . index ( ) != k { pattern_workspace . push_left ( I :: from_usize ( i ) ) ; flag_workspace [ i ] = I :: from_usize ( k ) ; i = parents . get_parent ( i ) . expect ( " " ) ; } pattern_workspace . push_left_on_right ( ) ; } // use a sparse triangular solve to compute the values // of the kth row of L diag [ k ] = y_workspace [ k ] ; y_workspace [ k ] = N :: zero ( ) ; ' pattern : for & i in pattern_workspace . iter_right ( ) { let i = i . index ( ) ; let yi = y_workspace [ i ] ; y_workspace [ i ] = N :: zero ( ) ; let p2 = l_colptr [ i ] + l_nz [ i ] ; for p in l_colptr [ i ] . index ( ) .. p2 . index ( ) { // we cannot go inside this loop before something has actually // be written into l_indices[l_colptr[i]..p2] so this // read is actually not into garbage // actually each iteration of the 'pattern loop adds writes the // value in l_indices that will be read on the next iteration // TODO: can some design change make this fact more obvious? let y_index = l_indices [ p ] . index ( ) ; y_workspace [ y_index ] = y_workspace [ y_index ] - l_data [ p ] * yi ; } let l_ki = yi / diag [ i ] ; diag [ k ] = diag [ k ] - l_ki * yi ; l_indices [ p2 . index ( ) ] = I :: from_usize ( k ) ; l_data [ p2 . index ( ) ] = l_ki ; l_nz [ i ] += I :: one ( ) ; } if diag [ k ] == N :: zero ( ) { panic ! ( " " ) ; } } }
Perform numeric LDLT decomposition
pattern_workspace is a DStack of capacity n
Perform numeric LDLT decomposition pattern_workspace is a DStack of capacity n
mat
CsMatViewI<N, I>
l_colptr
&[I]
parents
linalg::etree::ParentsView
perm
&Permutation<I, PStorage>
l_nz
&mut [I]
l_indices
&mut [I]
l_data
&mut [N]
diag
&mut [N]
y_workspace
&mut [N]
pattern_workspace
&mut DStack<I>
flag_workspace
&mut [I]
mat
CsMatViewI<N, I>
l_colptr
&[I]
parents
linalg::etree::ParentsView
perm
&Permutation<I, PStorage>
l_nz
&mut [I]
l_indices
&mut [I]
l_data
&mut [N]
diag
&mut [N]
y_workspace
&mut [N]
pattern_workspace
&mut DStack<I>
flag_workspace
&mut [I]