matrix.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include "classreg.h"

#include <stdio.h>
#include <math.h>
#include "ocmatrix.h"
#include "oc2iv.h"
#include "parse.hpp"
#include "ivocvect.h"

#define EPS hoc_epsilon
static Symbol* smat_;

extern int hoc_return_type_code;

extern double hoc_scan(FILE*);
extern "C" FILE* hoc_obj_file_arg(int i);
extern Object** hoc_temp_objptr(Object*);

#if 0
    extern void install_matrix_method(const char* name, double (*)(...));
    extern void* matrix_arg(int);
    extern double* matrix_pelm(void*, int i, int j);
    extern int matrix_nrow(void*);
    extern int matrix_ncol(void*);
    extern int matrix_type(void*);  // FULL 1, SPARSE 2, BAND 3
    extern MAT* matrix_full(void*); // hoc_execerror if void* not right type
    extern SPMAT* matrix_sparse(void*);
#endif

static void check_domain(int i, int j) {
    if (i > j || i < 0) {
        char buf[256];
        sprintf(buf, "index=%d  max_index=%d\n", i, j);
        hoc_execerror("Matrix index out of range:", buf);
    }
}

static void check_capac(int i, int j) {
    if (i != j) {
        hoc_execerror("wrong size for Matrix or Vector operation", 0);
    }
}

Matrix* matrix_arg(int i) {
    Object* ob = *hoc_objgetarg(i);
    if (!ob || ob->ctemplate != smat_->u.ctemplate) {
        check_obj_type(ob, "Matrix");
    }
    return (Matrix*)(ob->u.this_pointer);
}

Object** Matrix::temp_objvar() {
    Matrix* m = (Matrix*)this;
    Object** po;
    if (m->obj_) {
        po = hoc_temp_objptr(m->obj_);
    }else{
        po = hoc_temp_objvar(smat_, (void*)m);
        obj_ = *po;
    }
    return po;
}

static double m_nrow(void* v) {
    hoc_return_type_code = 1; // integer
    Matrix* m = (Matrix*)v;
    return (double)m->nrow();
}

static double m_ncol(void* v) {
    hoc_return_type_code = 1; // integer
    Matrix* m = (Matrix*)v;
    return (double)m->ncol();
}

static double m_setval(void* v) {
    Matrix* m = (Matrix*)v;
    int i, j;
    double val, *pval;
    i = (int)chkarg(1, 0, m->nrow()-1);
    j = (int)chkarg(2, 0, m->ncol()-1);
    val = *getarg(3);
    pval =m->mep(i, j);
    *pval = val;
    return val;
}

static double m_getval(void* v) {
    Matrix* m = (Matrix*)v;
    int i, j;
    i = (int)chkarg(1, 0, m->nrow()-1);
    j = (int)chkarg(2, 0, m->ncol()-1);
    return m->getval(i, j);
}

static double m_sprowlen(void* v) {
    Matrix* m = (Matrix*)v;
    int i;
    hoc_return_type_code = 1; // integer
    i = (int)chkarg(1, 0, m->nrow()-1);
    return double(m->sprowlen(i));
}

static double m_spgetrowval(void* v) {
    Matrix* m = (Matrix*)v;
    int i, jx, j;
    double x;
    i = (int)chkarg(1, 0, m->nrow()-1);
    jx = (int)chkarg(2, 0, m->sprowlen(i)-1);
    x = m->spgetrowval(i, jx, &j);
    if (ifarg(3)) {
        *hoc_pgetarg(3) = double(j);
    }
    return x;
}

static double m_printf(void* v) {
    Matrix* m = (Matrix*)v;
    int i, j, nrow = m->nrow(), ncol = m->ncol();
    const char* f1 = " %-8.3g";
    const char* f2 = "\n";
    if (ifarg(1)) {
        f1 = gargstr(1);
    }
    if (ifarg(2)) {
        f2 = gargstr(2);
    }
    for (i=0; i < nrow; ++i) {
        for (j=0; j < ncol; ++j) {
            Printf(f1, m->getval(i, j));
        }
        Printf("%s", f2);
    }
    return 0.;  
}

static double m_fprint(void* v) {
    Matrix* m = (Matrix*)v;
    int i, j, nrow = m->nrow(), ncol = m->ncol();
    int ia = 1;
    bool pr_size = true;
    const char* f1 = " %-8.3g";
    const char* f2 = "\n";
    if (hoc_is_double_arg(ia)) {
        pr_size = ((int)chkarg(ia, 0, 1) == 1) ? true : false;
        ++ia;
    }
    FILE* f = hoc_obj_file_arg(ia);
    if (ifarg(ia+1)) {
        f1 = gargstr(ia+1);
    }
    if (ifarg(ia+2)) {
        f2 = gargstr(ia+2);
    }
    if (pr_size) {
        fprintf(f, "%d %d\n", nrow, ncol);
    }
    for (i=0; i < nrow; ++i) {
        for (j=0; j < ncol; ++j) {
            fprintf(f, f1, m->getval(i, j));
        }
        fprintf(f, "%s", f2);
    }
    return 0.;  
}

static double m_scanf(void* v) {
    // file assumed to be an array of numbers. Numbers in rows
    // are contiguous in the stream.
    // first two numbers are nrow and ncol unless
    // arguments 2 and 3 specify them
    Matrix* m = (Matrix*)v;
    FILE* f = hoc_obj_file_arg(1);
    int i, j, nrow, ncol;
    if (ifarg(2)) {
        nrow = (int)chkarg(2, 1, 1e9);
        ncol = (int)chkarg(3, 1, 1e9);
    }else{
        nrow = (int)hoc_scan(f);
        ncol = (int)hoc_scan(f);
    }
    m->resize(nrow, ncol);
    for (i=0; i < nrow; ++i) for (j=0; j < ncol; ++j) {
        *(m->mep(i, j)) = hoc_scan(f);
    }
    return 0.;
}

static Object** m_resize(void* v) {
    Matrix* m = (Matrix*)v;
    m->resize((int)(chkarg(1, 1., 1e9) + EPS), (int)(chkarg(2, 1., 1e9) + EPS));
    return m->temp_objvar();
}

static Object** m_mulv(void* v) {
    Matrix* m = (Matrix*)v;
    Vect* vin = vector_arg(1);
    Vect* vout;
    bool f = false;
    if (ifarg(2)) {
        vout = vector_arg(2);
    }else{
#ifdef WIN32
        vout = vector_new1(m->nrow());
#else
        vout = new Vect(m->nrow());
#endif
    }
    if (vin == vout) {
        f = true;
#ifdef WIN32
        vin = vector_new2(vin);
#else
        vin = new Vect(*vin);
#endif
    }
#ifdef WIN32
    check_capac(vector_capacity(vin), m->ncol());
    vector_resize(vout, m->nrow());
#else
    check_capac(vin->size(), m->ncol());
    vout->resize(m->nrow());
#endif
    m->mulv(vin, vout);
    if (f) {
#ifdef WIN32
        vector_delete(vin);
#else
        delete vin;
#endif
    }
#ifdef WIN32
    return vector_temp_objvar(vout);
#else
    return vout->temp_objvar();
#endif
}


static Matrix* get_out_mat(Matrix* mat, int n, int m, int i, const char* mes = NULL);

static Matrix* get_out_mat(Matrix* mat, int n, int m, int i, const char* mes) {
    Matrix* out;
    if (ifarg(i)) {
        out = matrix_arg(i);
    }else{
        out = Matrix::instance(n, m, Matrix::MFULL);
        out->obj_ = NULL;
    }
    if (mat == out && mes) {
        hoc_execerror(mes, " matrix operation cannot be done in place");
    }
    return out;
}

static Matrix* get_out_mat(Matrix* m, int i , const char* mes = NULL) {
    return get_out_mat(m, m->nrow(), m->ncol(), i, mes);
}

static Object** m_add(void* v) {
    Matrix* m = (Matrix*)v;
    Matrix *out;
    out = m;
    if (ifarg(2)) {
        out = matrix_arg(2);
    }
    m->add(matrix_arg(1), out);
    return out->temp_objvar();
}

static Object** m_bcopy(void* v) {
    Matrix* m = (Matrix*)v;
    Matrix* out;
    int i0, j0, m0, n0, i1, j1, i;
    i0 = (int)chkarg(1, 0, m->nrow()-1);
    j0 = (int)chkarg(2, 0, m->ncol()-1);
    m0 = (int)chkarg(3, 1, m->nrow()-i0);
    n0 = (int)chkarg(4, 1, m->ncol()-j0);
    if (ifarg(5) && hoc_is_double_arg(5)) {
        i1 = (int)chkarg(5, 0, 1e9);
        j1 = (int)chkarg(6, 0, 1e9);
        i = 7;
    }else{
        i1 = 0;
        j1 = 0;
        i = 5;
    }
    out = get_out_mat(m, m0, n0, i);
    m->bcopy(out, i0, j0, m0, n0, i1, j1);
    return out->temp_objvar();
}

static Object** m_mulm(void* v) {
    Matrix* m = (Matrix*)v;
    Matrix* in, *out;
    in = matrix_arg(1);
    if (ifarg(2)) {
        out = matrix_arg(2);
    }else{
        out = Matrix::instance(m->nrow(), in->ncol(), Matrix::MFULL);
    }
    if (in == out || m == out) {
        hoc_execerror("matrix multiplication cannot be done in place", 0);
    }
    out->resize(m->nrow(), in->ncol());
    check_domain(m->ncol(), in->nrow());
    m->mulm(in, out);
    return out->temp_objvar();
}

static Object** m_c(void* v) {
    Matrix* m = (Matrix*)v;
    Matrix* out = get_out_mat(m, 1);
    m->copy(out);
    return out->temp_objvar();
}

static Object** m_transpose(void* v) {
    Matrix* m = (Matrix*)v;
    Matrix* out = get_out_mat(m, 1);
    out->resize(m->ncol(), m->nrow());
    m->transpose(out);
    return out->temp_objvar();
}

static Object** m_symmeig(void* v) {
    Matrix* m = (Matrix*)v;
    Matrix* out = matrix_arg(1);
    Object** p;
    out->resize(m->nrow(), m->ncol());
    Vect* vout;
#ifdef WIN32
    vout = vector_new1(m->nrow());
    p = vector_temp_objvar(vout);
#else
    vout = new Vect(m->nrow());
    p = vout->temp_objvar();
#endif
    m->symmeigen(out, vout);
    return p;
}

static Object** m_svd(void* vv) {
    Matrix* m = (Matrix*)vv;
    Matrix *u=NULL, *v=NULL;
    if (ifarg(2)) {
        u = matrix_arg(1);
        v = matrix_arg(2);
        u->resize(m->nrow(), m->nrow());
        v->resize(m->ncol(), m->ncol());
    }
    Object** p;
    Vect* d;
    int dsize = m->nrow() < m->ncol() ? m->nrow() : m->ncol();
#ifdef WIN32
    d = vector_new1(dsize);
    p = vector_temp_objvar(d);
#else
    d = new Vect(dsize);
    p = d->temp_objvar();
#endif
    m->svd1(u, v, d);
    return p;
}

static Object** m_muls(void* v) {
    Matrix* m = (Matrix*)v;
    Matrix *out;
    out = m;
    if (ifarg(2)) {
        out = matrix_arg(2);
    }
//  I believe meschach does this for us
//  if (out != m) {
//      out->resize(...
//  }
    m->muls(*getarg(1), out);
    return out->temp_objvar();
}

static Object** m_getrow(void* v) {
    Matrix* m = (Matrix*)v;
    int k = (int)chkarg(1, 0, m->nrow()-1); 
    Vect* vout;
    if (ifarg(2)) {
        vout = vector_arg(2);
#ifdef WIN32
        vector_resize(vout, m->ncol());
#else
        vout->resize(m->ncol());
#endif
    }else{
#ifdef WIN32
        vout = vector_new1(m->ncol());
#else
        vout = new Vect(m->ncol());
#endif
    }
    m->getrow(k, vout);
#ifdef WIN32
    return vector_temp_objvar(vout);
#else
    return vout->temp_objvar();
#endif
}

static Object** m_getcol(void* v) {
    Matrix* m = (Matrix*)v;
    int k = (int)chkarg(1, 0, m->ncol()-1); 
    Vect* vout;
    if (ifarg(2)) {
        vout = vector_arg(2);
#ifdef WIN32
        vector_resize(vout, m->nrow());
#else
        vout->resize(m->nrow());
#endif
    }else{
#ifdef WIN32
        vout = vector_new1(m->nrow());
#else
        vout = new Vect(m->nrow());
#endif
    }
    m->getcol(k, vout);
#ifdef WIN32
    return vector_temp_objvar(vout);
#else
    return vout->temp_objvar();
#endif
}

static Object** m_setrow(void* v) {
    Matrix* m = (Matrix*)v;
    int k = (int)chkarg(1, 0, m->nrow()-1); 
    if (hoc_is_double_arg(2)) {
        m->setrow(k, *getarg(2));
    }else{
        Vect* in = vector_arg(2);
#ifdef WIN32
        check_domain(vector_capacity(in), m->ncol());
#else
        check_domain(in->size(), m->ncol());
#endif
        m->setrow(k, in);
    }
    return m->temp_objvar();
}

static Object** m_setcol(void* v) {
    Matrix* m = (Matrix*)v;
    int k = (int)chkarg(1, 0, m->ncol()-1); 
    if (hoc_is_double_arg(2)) {
        m->setcol(k, *getarg(2));
    }else{
        Vect* in = vector_arg(2);
#ifdef WIN32
        check_domain(vector_capacity(in), m->nrow());
#else
        check_domain(in->size(), m->nrow());
#endif
        m->setcol(k, in);
    }
    return m->temp_objvar();
}

static Object** m_setdiag(void* v) {
    Matrix* m = (Matrix*)v;
    int k = (int)chkarg(1, -(m->nrow() - 1), m->ncol()-1);  
    if (hoc_is_double_arg(2)) {
        m->setdiag(k, *getarg(2));
    }else{
        Vect* in = vector_arg(2);
#ifdef WIN32
        check_domain(vector_capacity(in), m->nrow());
#else
        check_domain(in->size(), m->nrow());
#endif
        m->setdiag(k, in);
    }
    return m->temp_objvar();
}

static Object** m_getdiag(void* v) {
    Matrix* m = (Matrix*)v;
    int k = (int)chkarg(1, -(m->nrow()-1), m->ncol()-1);    
    Vect* vout;
    if (ifarg(2)) {
        vout = vector_arg(2);
#ifdef WIN32
        vector_resize(vout, m->nrow());
#else
        vout->resize(m->nrow());
#endif
    }else{
#ifdef WIN32
        vout = vector_new1(m->nrow());
#else
        vout = new Vect(m->nrow());
#endif
    }
    m->getdiag(k, vout);
#ifdef WIN32
    return vector_temp_objvar(vout);
#else
    return vout->temp_objvar();
#endif
}

static Object** m_zero(void* v) {
    Matrix* m = (Matrix*)v;
    m->zero();
    return m->temp_objvar();
}

static Object** m_ident(void* v) {
    Matrix* m = (Matrix*)v;
    m->ident();
    return m->temp_objvar();
}

static Object** m_exp(void* v) {
    Matrix* m = (Matrix*)v;
    Matrix* out = get_out_mat(m, 1, "exponentiation");
    m->exp(out);
    return out->temp_objvar();
}

static Object** m_pow(void* v) {
    Matrix* m = (Matrix*)v;
    int k = (int)chkarg(1, 0., 100.);
    Matrix* out = get_out_mat(m, 2, "raising to a power");
    m->pow(k, out);
    return out->temp_objvar();
}

static Object** m_inverse(void* v) {
    Matrix* m = (Matrix*)v;
    Matrix* out = get_out_mat(m, 1);
    m->inverse(out);
    return out->temp_objvar();
}

static double m_det(void* v) {
    Matrix* m = (Matrix*)v;
    int e;
    double a =  m->det(&e);
    double* pe = hoc_pgetarg(1);
    *pe = double(e);
    return a;
}

static Object** m_solv(void* v) {
    Matrix* m = (Matrix*)v;
    check_capac(m->nrow(), m->ncol());
    Vect* vin = vector_arg(1);
#ifdef WIN32
    check_capac(vector_capacity(vin), m->ncol());
#else
    check_capac(vin->size(), m->ncol());
#endif
    Vect* vout = NULL;
    bool f = false;
    bool use_lu = false;
    // args 2 and 3 are optional [vout, use previous LU factorization]
    // and in either order
    for (int i=2; i <=3; ++i) {
        if (ifarg(i)) {
            if (hoc_is_object_arg(i)) {
                if (vout) {
                }
                vout = vector_arg(i);
            }else{
                use_lu = ((int)(*getarg(i))) ? true : false;
            }
        }
    }
    if (!vout) {
#ifdef WIN32
        vout = vector_new1(m->nrow());
#else
        vout = new Vect(m->nrow());
#endif
    }
#ifdef WIN32
    vector_resize(vout, m->ncol());
#else
    vout->resize(m->ncol());
#endif
    if (vin == vout) {
        f = true;
#ifdef WIN32
        vin = vector_new2(vin);
#else
        vin = new Vect(*vin);
#endif
    }
    m->solv(vin, vout, use_lu);
    if (f) {
#ifdef WIN32
        vector_delete(vin);
#else
        delete vin;
#endif
    }
#ifdef WIN32
    return vector_temp_objvar(vout);
#else
    return vout->temp_objvar();
#endif
}

static Object** m_set(void* v) {
    Matrix* m = (Matrix*)v;
    int i, j, nrow = m->nrow(), ncol = m->ncol();
    int k;
    for (k=0, i=0; i < nrow; ++i) {
        for (j=0; j < ncol; ++j) {
            *(m->mep(i, j)) = *getarg(++k);
        }
    }
    return m->temp_objvar();
}

static Object** m_to_vector(void* v) {
    Matrix* m = (Matrix*)v;
    Vect* vout;
    int i, j, k;
    int nrow = m->nrow();
    int ncol = m->ncol();
    if (ifarg(1)) {
        vout = vector_arg(1);
        vector_resize(vout, nrow*ncol);
    }else{
        vout = vector_new1(nrow*ncol);
    }
    k = 0;
    double* ve = vector_vec(vout);
    for (j=0; j < ncol; ++j) for (i=0; i < nrow; ++i) {
        ve[k++] = m->getval(i, j);
    }
    return vector_temp_objvar(vout);
}

static Object** m_from_vector(void* v) {
    Matrix* m = (Matrix*)v;
    Vect* vout;
    int i, j, k;
    int nrow = m->nrow();
    int ncol = m->ncol();
    vout = vector_arg(1);
    check_capac(nrow*ncol, vector_capacity(vout));
    k = 0;
    double* ve = vector_vec(vout);
    for (j=0; j < ncol; ++j) for (i=0; i < nrow; ++i) {
        *(m->mep(i, j)) = ve[k++];
    }
    return m->temp_objvar();
}

static Member_func m_members[] = {
    // returns double scalar
    "x", m_nrow, // will be changed below
    "nrow", m_nrow,
    "ncol", m_ncol,
    "getval", m_getval,
    "setval", m_setval,
    "sprowlen", m_sprowlen,
    "spgetrowval", m_spgetrowval,
    "det", m_det,

    "printf", m_printf,
    "fprint", m_fprint,
    "scanf", m_scanf,
    0, 0
};
    
static Member_ret_obj_func m_retobj_members[] = {
    // returns Vector
    "mulv", m_mulv,
    "getrow", m_getrow,
    "getcol", m_getcol,
    "getdiag", m_getdiag,
    "solv", m_solv,
    "symmeig", m_symmeig,
    "svd", m_svd,
    // returns Matrix
    "c", m_c,
    "add", m_add,
    "bcopy", m_bcopy,
    "resize", m_resize,
    "mulm", m_mulm,
    "muls", m_muls,
    "setrow", m_setrow,
    "setcol", m_setcol,
    "setdiag", m_setdiag,
    "zero", m_zero,
    "ident", m_ident,
    "exp", m_exp,
    "pow", m_pow,
    "inverse", m_inverse,
    "transpose", m_transpose,
    "set", m_set,
    "to_vector", m_to_vector,
    "from_vector", m_from_vector,
    0, 0
};

static void* m_cons(Object* o) {
    int i=1, j=1, storage_type = Matrix::MFULL;
    if (ifarg(1)) i = int(chkarg(1, 1, 1e10) + EPS);
    if (ifarg(2)) j = int(chkarg(2, 1, 1e10) + EPS);
    if (ifarg(3)) storage_type = int(chkarg(3, 1, 3));
    Matrix* m = Matrix::instance(i, j, storage_type);
    m->obj_ = o;
    return m;
}

static void m_destruct(void* v) {
    // supposed to notify freed val array here.
//printf("Matrix deleted\n");
    delete (Matrix*)v;
}

static void steer_x(void* v) {
    Matrix* m = (Matrix*)v;
    int i1, i2;
    Symbol* s = hoc_spop();
    i2 = (int)(hoc_xpop() + EPS);
    i1 = (int)(hoc_xpop() + EPS);
    check_domain(i1, m->nrow()-1);
    check_domain(i2, m->ncol()-1);
    hoc_pushpx(m->mep(i1, i2));
}


#if WIN32 && !USEMATRIX
void Matrix_reg();
#endif

void Matrix_reg() {
    class2oc("Matrix", m_cons, m_destruct, m_members, NULL, m_retobj_members, NULL);
    smat_ = hoc_lookup("Matrix");
    // now make the x variable an actual double
    Symbol* sx = hoc_table_lookup("x", smat_->u.ctemplate->symtable);
    sx->type = VAR;
    sx->arayinfo = (Arrayinfo *)hoc_Emalloc(sizeof(Arrayinfo) + 2*sizeof(int));
    sx->arayinfo->refcount = 1;
    sx->arayinfo->a_varn = NULL;
    sx->arayinfo->nsub = 2;
    sx->arayinfo->sub[0] = 1;
    sx->arayinfo->sub[1] = 1;
    smat_->u.ctemplate->steer = steer_x;
}