SF_linalg_utils.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#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