nonlin.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include <stdlib.h>
#include <math.h>
#include "hoc.h"
#include "parse.hpp"
#include "hocparse.h"
#include "equation.h"
#include "lineq.h"
#include "code.h"




int do_equation;    /* switch for determining access to dep vars */
int *hoc_access;    /* links to next accessed variables */
int  var_access;    /* variable number as pointer into access array */


static double   **varble;   /* pointer to symbol values */
typedef struct elm *Elm;

#define diag(s) hoc_execerror(s, (char*)0);

void dep_make(void)/* tag the variable as dependent with a variable number */
{
#if !OCSMALL
    Symbol *sym;
    unsigned    *numpt=0;
#if defined(__TURBOC__)
    Inst *pcsav=pc;
#endif

    sym = spop();
#if defined(__TURBOC__)
    pc = pcsav;
#endif
    switch (sym->type)
    {

    case UNDEF:
        hoc_execerror(sym->name, "undefined in dep_make");
        sym->type = VAR;
        OPVAL(sym) = (double *)emalloc(sizeof(double));
        *(OPVAL(sym)) = 0.;
    case VAR:
        if (sym->subtype != NOTUSER) {
            execerror(sym->name, "can't be a dependent variable");
        }
if (!ISARRAY(sym)) {
            numpt = &(sym->s_varn);
}else{
        Arrayinfo *aray = OPARINFO(sym);
        if (sym->s_varn == 0)   /* allocate varnum array */
        {
            int total = 1;
            int i;
            for (i=0; i < aray->nsub; i++)
                total *= (aray->sub)[i];
aray->a_varn = (unsigned *)ecalloc((unsigned)total, sizeof(unsigned));
            sym->s_varn = (unsigned)total;  /* set_varble() uses this */
        }
        numpt = &((aray->a_varn)[araypt(sym, OBJECTVAR)]);
}
        break;

    default:
        execerror(sym->name, "can't be a dependent variable");
    }
    if (*numpt > 0)
        execerror(sym->name, "made dependent twice");
    *numpt = ++neqn;
#endif
}


void init_access(void) /* zero the access array */
{
#if !OCSMALL
    if (hoc_access != (int *)0)
        free((char *)hoc_access);
    hoc_access = (int *)ecalloc((neqn+1), sizeof(int));
    var_access = -1;
#endif
}

static void eqn_space(void);    /* reallocate space for matrix */
static void set_varble(void);   /* set up varble array by searching for tags */
static void eqn_side(int lhs);
static unsigned row;
static unsigned maxeqn;

void eqn_name(void) /* save row number for lhs and/or rhs */
{
#if !OCSMALL

    if (maxeqn != neqn) /* discard equations and reallocate space */
    {
        eqn_space();
        set_varble();
    }
    init_access();
    do_equation = 1;
    eval();
    do_equation = 0;
    if (var_access < 1)
        execerror("illegal equation name",(pc - 2)->sym->name);
    row = var_access;
    nopop();
#endif
}

static void set_varble(void) {  /* set up varble array by searching for tags */
#if !OCSMALL
    Symbol  *sp;

    for (sp = symlist->first; sp != (Symbol *)0; sp = sp->next)
    {
        if (sp->s_varn != 0)
        {
            if (sp->type == VAR) {          
                if (!ISARRAY(sp)) {
                    varble[sp->s_varn] = OPVAL(sp);
                }else{
                    int i;
                    Arrayinfo   *aray = OPARINFO(sp);
                    for (i = 0; i < sp->s_varn; i++)
                        if ((aray->a_varn)[i] > 0)
varble[(aray->a_varn)[i]] = OPVAL(sp) + i;
                }
            }
        }
    }
#endif
}

static double Delta = .001; /* variable variation */

void eqinit(void)   /* built in function to initialize equation solver */
{
#if !OCSMALL
    Symbol  *sp;

    if (ifarg(1))
        Delta = *getarg(1);
    for (sp = symlist->first; sp != (Symbol *)0; sp = sp->next)
    {
        if (sp->s_varn != 0)
        {
            if (ISARRAY(sp)) {
                if (OPARINFO(sp)->a_varn != (unsigned *)0)
                    free((char *)(OPARINFO(sp)->a_varn));
            }
            sp->s_varn = 0;
        }
    }
    neqn = 0;
    eqn_space();
#endif
    ret();
    pushx(0.);
}

void eqn_init(void) /* initialize equation row */
{
#if !OCSMALL
    struct elm *el;

    for (el = rowst[row]; el != (struct elm *)0; el = el->c_right)
        el->value = 0.;
    rhs[row] = 0.;
#endif
}

void eqn_lhs(void)  /* add terms to left hand side */
{
    eqn_side(1);
}

void eqn_rhs(void)  /* add terms to right hand side */
{
    eqn_side(0);
}



static void eqn_side(int lhs) {
#if !OCSMALL
    int i;
    struct elm *el;
    double f0, f1;
    Inst *savepc = pc;
#if defined(__TURBOC__)
    Inst *pc1;
#endif

    init_access();
    do_equation = 1;
    execute(savepc);
#if defined(__TURBOC__)
    pc1 = pc;
#endif
    do_equation = 0;

    if (lhs)
    {
        f0 = xpop();
    }
    else
    {
        f0 = -xpop();
    }

    rhs[row] -= f0;
    for (i = var_access; i > 0; i = hoc_access[i])
    {
        *varble[i] += Delta;
        execute(savepc);
        *varble[i] -= Delta;
        if (lhs)
        {
            f1 = xpop();
        }
        else
        {
            f1 = -xpop();
        }
        el = getelm((struct elm *)0, row, (unsigned)i);
        el->value += (f1 - f0)/Delta;
    }
#if defined(__TURBOC__)
    pc=pc1;
#endif
    pc++;
#endif
}

static void eqn_space(void) {   /* reallocate space for matrix */
#if !OCSMALL
    int i;
    struct elm *el;

    if (maxeqn > 0 && rowst == (Elm *)0)
        diag("matrix coefficients cannot be released");
    for(i=1 ; i <= maxeqn ; i++)
        for(el = rowst[i] ; el; el = el->c_right)
            free((char *)el);
    maxeqn = neqn;
    if (varble)
        free((char *)varble);
    if (rowst)
        free((char *)rowst);
    if (colst)
        free((char *)colst);
    if (eqord)
        free((char *)eqord);
    if (varord)
        free((char *)varord);
    if (rhs)
        free((char *)rhs);
    varble = (double **)0;
    rowst = colst = (Elm *)0;
    eqord = varord = (unsigned *)0;
    rhs = (double *)0;
    rowst = (Elm *)ecalloc((maxeqn+1),sizeof(struct elm *));
    varble = (double **)emalloc((maxeqn+1)*sizeof(double *));
    colst = (Elm *)ecalloc((maxeqn+1),sizeof(struct elm *));
    eqord = (unsigned *)emalloc((maxeqn+1)*sizeof(unsigned));
    varord = (unsigned *)emalloc((maxeqn+1)*sizeof(unsigned));
    rhs = (double *)emalloc((maxeqn+1)*sizeof(double));
    for (i=1 ; i<= maxeqn ; i++)
    {
        eqord[i] = i;
        varord[i] = i;
    }
#endif
}

void hoc_Prmat(void)
{
#if !OCSMALL
    prmat();
#endif
    ret();
    pushx(1.);
}

void solve(void)
{
#if !OCSMALL
    /* Sum is a measure of the dependent variable accuracy
       and how well the equations are solved */

    int i;
    double sum;
    struct elm  *el;

    sum = 0.;
    for (i=1 ; i <= neqn ; i++)
        sum += fabs(rhs[i]);
    if (!matsol())
        diag("indeterminate system");
    for (i=1 ; i<= neqn ; i++)
    {
        *varble[varord[i]] += rhs[eqord[i]];
        sum += fabs(rhs[i]);
    }
    /* free all elements */
    for (i=1; i <= neqn; i++)
    {
        struct elm *el2;
        for (el = rowst[i]; el != (struct elm *)0; el = el2) {
            el2 = el->c_right;
            free(el);
        }
        rowst[i] = colst[i] = (struct elm *)0;
    }
#else
    double sum = 0;
#endif
    ret();
    pushx(sum);
}