matrixmap.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include "matrixmap.h"
#include <vector>
using std::vector;

#include "spmatrix.h"
extern spREAL *spGetElement(char*, int ,int);

MatrixMap::MatrixMap(Matrix& mat)
    : m_(mat), plen_(0), ptree_(NULL),  pm_(NULL) {
}

MatrixMap::MatrixMap(Matrix* mat)
    : m_(*mat), plen_(0), ptree_(NULL), pm_(NULL) {
}

MatrixMap::~MatrixMap() {
    mmfree();
}

void MatrixMap::mmfree() {
    // safe to delete NULL pointers
    delete [] ptree_;
    delete [] pm_;
    ptree_ = NULL;
    pm_ = NULL;
}

void MatrixMap::add(double fac) {
    for (int i=0 ; i < plen_; ++i) {
//printf("i=%d %g += %g * %g\n", i, *ptree_[i], fac, *pm_[i]);
        *ptree_[i] += fac * (*pm_[i]);
    }
}

void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
    NrnThread* _nt = nrn_threads;
    mmfree();
    // how many elements
    int nrow = m_.nrow();
    int ncol = m_.ncol();
//printf("MatrixMap::alloc nrow=%d ncol=%d\n", nrow, ncol);

    plen_ = 0;
    vector<int> nonzero_i, nonzero_j;
    m_.nonzeros(nonzero_i, nonzero_j);
    pm_ = new double*[nonzero_i.size()];
    ptree_ = new double*[nonzero_i.size()];
    for (int k = 0; k < nonzero_i.size(); k++) {
        const int i = nonzero_i[k];
        const int j = nonzero_j[k];
        int it;
        if (i < nnode) {
            it = nodes[i] -> eqn_index_ + layer[i];
//printf("i=%d it=%d area=%g\n", i, it, nodes[i]->area);
            if (layer[i] > 0 && !nodes[i] -> extnode) {
                it = 0;
            }
        } else {
            it = start + i - nnode;
        }
        int jt;
        pm_[plen_] = m_.mep(i, j);
        if (j < nnode) {
            jt = nodes[j] -> eqn_index_ + layer[j];
            if (layer[j] > 0 && !nodes[j] -> extnode) {
                jt = 0;
            }
        } else {
            jt = start + j - nnode;
        }
//printf("MatrixMap::alloc getElement(%d,%d)\n", it, jt);
        ptree_[plen_] = spGetElement(_nt->_sp13mat, it, jt);
        ++plen_;
    }
}