SF_linalg_utils.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: Apache-2.0
 
#ifndef _SF_LINALG_H
#define _SF_LINALG_H
 
#include "SF_vector.h"
#include "SF_sort.h"
 
namespace SF {
 
template<class T, class S>
class sparse_multiply_transpose
{
public:
 
    inline void operator()(const vector<T> &_acnt,
                    const vector<T> &_acol,
                    const vector<S> &_aele,
                    const vector<T> &_bcnt,
                    const vector<T> &_bcol,
                    const vector<S> &_bele,
                    const T _csize,
                    vector<T> &_ccnt,
                    vector<T> &_ccol,
                    vector<S> &_cele)
    {
        int arows = (int)_acnt.size();
        const T *acnt = &_acnt[0];
        const T *acol = &_acol[0];
        const S *aele = &_aele[0];
        const T *bcnt = &_bcnt[0];
        const T *bcol = &_bcol[0];
        const S *bele = &_bele[0];
 
        int brows = (int)_bcnt.size();
        vector<T> _bdsp(brows+1);
        dsp_from_cnt(_bcnt, _bdsp);
 
        _ccnt.resize(_csize);
        _ccnt.zero();
 
        vector<T> _clst(_csize, -1);
        T *ccnt = _ccnt.data();
        T *clst = _clst.data();
        T *bdsp = _bdsp.data();
        int csize = 0;
        for(int m = 0, i = 0; i < arows; i++)
        {
            for(int j = 0; j < acnt[i]; j++)
            {
                T c = acol[m];
                m++;
                for(int k = bdsp[c]; k < bdsp[c] + bcnt[c]; k++)
                {
                    T d = bcol[k];
                    if(clst[d] != i)
                    {
                        clst[d] = i;
                        ccnt[d]++;
                        csize++;
                    }
                }
            }
        }
 
        _clst.assign(size_t(_csize), -1);
        clst = _clst.data();
 
        vector<T> _cdsp(_csize+1);
        dsp_from_cnt(_ccnt, _cdsp);
 
        _ccol.resize(csize);
        _cele.resize(csize);
        T *ccol = _ccol.data();
        S *cele = _cele.data();
        T *cdsp = _cdsp.data();
        for(int m = 0, i = 0; i < arows; i++)
        {
            for(int j = 0; j < acnt[i]; j++)
            {
                T c = acol[m];
                S s = aele[m];
                m++;
                for(int k = bdsp[c]; k < bdsp[c] + bcnt[c]; k++)
                {
                    T d = bcol[k];
                    T e = cdsp[d];
                    S t = s * bele[k];
                    if(clst[d] != i)
                    {
                        clst[d] = i;
                        ccol[e] = i;
                        cele[e] = t;
                        cdsp[d] = e + 1;
                    }
                    else
                    {
                        cele[e - 1] += t;
                    }
                }
            }
        }
    }
};
 
 
template<class T>
inline void multiply_connectivities(const vector<T> & a_cnt,
                             const vector<T> & a_con,
                             const vector<T> & b_cnt,
                             const vector<T> & b_con,
                             vector<T> & c_cnt,
                             vector<T> & c_con)
{
  int numnodes = a_cnt.size();
 
  vector<T> a_ele, b_ele, c_ele;
  a_ele.assign(a_con.size(), 1);
  b_ele.assign(b_con.size(), 1);
 
  sparse_multiply_transpose<T, T> multiply;
  multiply(a_cnt, a_con, a_ele, b_cnt, b_con, b_ele, numnodes, c_cnt, c_con, c_ele);
}
 
template<class T>
inline void transpose_connectivity(const vector<T> & a_cnt,
                            const vector<T> & a_con,
                            vector<T> & b_cnt,
                            vector<T> & b_con)
{
  if(a_con.size() == 0) return;
 
  // get the largest node index to determine number of nodes
  T numnodes = *(std::max_element(a_con.begin(), a_con.end())) + 1;
  vector<T> b_row;
 
  // we compute bcol := arow
  b_con.resize(a_con.size());
  for(size_t i=0, k=0; i < a_cnt.size(); i++)
    for(int j=0; j < a_cnt[i]; j++, k++) b_con[k] = i;
 
  b_row.assign(a_con.begin(), a_con.end());
  binary_sort_sort(b_row, b_con);
 
  b_cnt.resize(numnodes); b_cnt.zero();
  count(b_row, b_cnt);
}
 
}
 
#endif