/* allocate.c  93.12.19
 * Copyright 1983-1992   Albert Davis
 * Allocates space for the admittance matrix.  Allocation below the
 * diagonal is by row, above the diagonal is by column, and stored backwards.
 * This broken vector is stored.  The length of it increases with increasing
 * position.  The maximum length of the nth vector is 2n-1.  For a band matrix
 * only those elements that are non-zero or are nearer to the diagonal than a
 * non-zero element are stored.
 */
#include "ecah.h"
#include "branch.h"
#include "error.h"
#include "mode.h"
#include "nodestat.h"
#include "options.h"
#include "status.h"
#include "declare.h"
/*--------------------------------------------------------------------------*/
	void	 allocate(int);
static	void	 count_nodes(void);
static	void	 expand_list(branch_t*);
static	void	 map_nodes(void);
static	void	 order_reverse(void);
static	void	 order_forward(void);
static	void	 order_auto(void);
static	void	 map_list(branch_t*);
static	void	 build_basenode_table(void);
static	void	 addto_basenode_table(const branch_t*);
static	void	 lownode(int,int);
static	void	 alloc_hold_vectors(void);
static	unsigned size_arrays(void);
static	void	 alloc_matrix(unsigned);
static	void	 alloc_constants(void);
static	void	 alloc_pointer_table(void);
static	void	 fill_pointer_table(void);
/*--------------------------------------------------------------------------*/
extern const struct options opt;
extern       struct status stats;
extern       struct nodestuff ns;
extern const char e_om[];
extern unsigned aspace;		    /* count of non-zero array elements	    */
extern double **rerow, **imrow;	    /* array of row pointers		    */
extern double **recol, **imcol;	    /* array of column pointers		    */
extern double **redia, **imdia;	    /* array of diagonal pointers	    */
extern double *reals, *imags;	    /* big array allocation bases	    */
extern struct nodestat *nstat;	    /* node status flags		    */
/*--------------------------------------------------------------------------*/
/* allocate: memory allocation for simulation
 */
void allocate(int mode)
{
 if (mode != sSTATUS){
    time_reset(&(stats.load));
    time_reset(&(stats.lu));
    time_reset(&(stats.back));
    time_reset(&(stats.review));
    time_reset(&(stats.output));

    time_zstart(&(stats.setup));
 }
 if (!nstat){		    	    	/* if "nstat" not allocated,	    */
    dealloc(YES);		    	/* must allocate everything	    */
    count_nodes();		 	/* and start solution over.	    */
    expand_list(firstbranch_dev());
    map_nodes();
    alloc_hold_vectors();
    build_basenode_table();
    aspace = size_arrays();
 }
 if (mode==sSTATUS)
    return;						/* done if status   */
 if (!exists(firstbranch_all()))
    error(bERROR, "no circuit\n");

				/* lost between commands.  re-- every time  */
 alloc_matrix(aspace);
 alloc_constants();
 alloc_pointer_table();
 fill_pointer_table();
 if (mode != sSTATUS)
    time_stop(&(stats.setup));
}
/*--------------------------------------------------------------------------*/
/* count_nodes: count nodes in main ckt (not subckts)
 * update the variables "stats.total_nodes" and "stats.user_nodes"
 * zeros "stats.subckt_nodes" and "stats.model_nodes"
 */
static void count_nodes(void)
{
 const branch_t *stop;

 stats.user_nodes = 0;
 stop = firstbranch_dev();
 if (isdevice(stop)){
    const branch_t *brh;
    brh = stop;
    do {
       int ii;
       for (ii = 0;  brh->n[ii].e != INVALIDNODE;  ii++){
          if (brh->n[ii].e > stats.user_nodes)
	     stats.user_nodes = brh->n[ii].e;
       }
    } while (brh = nextbranch_dev(brh),  brh != stop);
 }
 stats.total_nodes = stats.user_nodes;
 stats.subckt_nodes = stats.model_nodes = 0;
}
/*--------------------------------------------------------------------------*/
/* expand_list: expand (flatten) a list of components (subckts)
 * Scan component list.  Expand each subckt: create actual elements
 * for flat representation to use for simulation.
 * Recursive to allow for nested subckts.
 */
static void expand_list(branch_t *stop)
{
 if (isdevice(stop)){
    branch_t *brh;
    brh = stop;
    do {
       expand_branch(brh);
       if (brh->subckt){
          expand_list(brh->subckt);
       }
    } while (brh=nextbranch_dev(brh),  brh != stop);
 }
}
/*--------------------------------------------------------------------------*/
/* map_nodes: map intermediate node number to internal node number.
 * Ideally, this function would find some near-optimal order
 * and squash out gaps.
 */
static void map_nodes(void)
{
 ns.nm = (int*)calloc((unsigned)stats.total_nodes+1, sizeof(int));
 if (!ns.nm)				/* node map table		    */
    error(bERROR, e_om, "nm");
 time_zstart(&(stats.order));
 switch (opt.order){
    default:
       error(bWARNING, "invalid order spec: %d\n", opt.order);
       /* fall through to order_reverse. */
    case oAUTO:
       order_auto();
       break;
    case oREVERSE:
       order_reverse();
       break;
    case oFORWARD:
       order_forward();
       break;
 }
 time_stop(&(stats.order));
 map_list(firstbranch_dev());
}
/*--------------------------------------------------------------------------*/
/* order_reverse: force ordering to reverse of user ordering
 *  subcircuits at beginning, results on border at the bottom
 */
static void order_reverse(void)
{
 int node;
 ns.nm[0] = 0;
 for (node=1;  node<=stats.total_nodes;  ++node)
    ns.nm[node] = stats.total_nodes - node + 1;
}
/*--------------------------------------------------------------------------*/
/* order_forward: use user ordering, with subcircuits added to end
 * results in border at the top (worst possible if lots of subcircuits)
 */
static void order_forward(void)
{
 int node;
 ns.nm[0] = 0;
 for (node=1;  node<=stats.total_nodes;  ++node)
    ns.nm[node] = node;
}
/*--------------------------------------------------------------------------*/
/* order_auto: full automatic ordering
 * reverse, for now
 */
static void order_auto(void)
{
 int node;
 ns.nm[0] = 0;
 for (node=1;  node<=stats.total_nodes;  ++node)
    ns.nm[node] = stats.total_nodes - node + 1;
}
/*--------------------------------------------------------------------------*/
/* map_list: map intermediate node to internal (working) node number.
 * recursive: actual mapping done here
 */
static void map_list(branch_t *stop)
{
 if (isdevice(stop)){
    branch_t *brh;
    brh = stop;
    do {
       int ii;
       for (ii = 0;  brh->n[ii].t != INVALIDNODE;  ii++){
	  brh->n[ii].m = ns.nm[brh->n[ii].t];
       }
       brh->n[ii].m = INVALIDNODE;
       if (brh->subckt)
          map_list(brh->subckt);
    } while (brh = nextbranch_dev(brh),  brh != stop);
 }
}
/*--------------------------------------------------------------------------*/
/* build_basenode_table: build a table containing the lowest node number
 * connected to each node.  This determines the bandwidth of the matrix,
 * or the length of the vectors (in the matrix) corresponding to each node.
 * Actual work is done by addto_basenode_table.
 */
static void build_basenode_table(void)
{
 int node;

 ns.basnode = (int*)calloc((unsigned)stats.total_nodes+1, sizeof(int));
 if (!ns.basnode)		    /* base-node, (vector size) table	    */
    error(bERROR, e_om, "basnode");

 for (node = 0;  node <= stats.total_nodes;  ++node){
    ns.basnode[node] = node;	    /* initialize each to itself	    */
 }
 addto_basenode_table(firstbranch_dev());
}
/*--------------------------------------------------------------------------*/
/* addto_basenode_table: actual work of build_basenode_table.
 * Find and store "base-node" for each node.  This is the lowest numbered
 * node that it connects to.  We don't do anything with higher numbered nodes.
 * (Unnecessary because of symmetry.)   Recursive, for subckts.
 */
static void addto_basenode_table(const branch_t *stop)
{
 if (isdevice(stop)){
    const branch_t *brh;
    brh = stop;
    do {
       if (brh->subckt){
          addto_basenode_table(brh->subckt);
       }else{
	  int ii,jj;
	  for (ii = 0;  brh->n[ii].m != INVALIDNODE;  ii++){
	     if (brh->n[ii].m != 0){
		for (jj = 0;  jj < ii ;  jj++)
		   lownode(brh->n[ii].m,brh->n[jj].m);
		++nstat[brh->n[ii].m].needsanalog;
	     }
	  }
       }
    } while (brh = nextbranch_dev(brh), brh != stop);
 }
}
/*--------------------------------------------------------------------------*/
/* lownode: find the lowest node that each node connects to.
 * (connects to == has distance 1 from)
 */
static void lownode(int node1, int node2)
{
 if (node1 == 0  ||  node2 == 0)	/* node 0 is ground, and doesn't    */
    ;					/* count as a connection	    */
 else if ( node1 < ns.basnode[node2] ) 
    ns.basnode[node2]=node1;
 else if ( node2 < ns.basnode[node1] ) 
    ns.basnode[node1]=node2;
}
/*--------------------------------------------------------------------------*/
/* alloc_hold_vectors: allocate space to hold data between commands.
 * store voltage between commands, for restart and convergence assistance
 */
static void alloc_hold_vectors(void)
{
 nstat = (struct nodestat*)calloc((unsigned)stats.total_nodes+2,
    			sizeof(struct nodestat));
 if (!nstat)
    error(bERROR, e_om, "nstat");

 ns.vdc = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
 if (!ns.vdc)
    error(bERROR, e_om, "vdc");
}
/*--------------------------------------------------------------------------*/
/* size_arrays: determine memory needed for main matrix
 * (count non-zero array elememts)
 */
static unsigned size_arrays(void)
{
 int node;
 unsigned space;

 space = 0;
 for (node = 0;   node <= stats.total_nodes;   ++node)
    space += 2 * ( node-ns.basnode[node] ) + 1;
 return space;
}
/*--------------------------------------------------------------------------*/
/* alloc_matrix: allocate main matrix
 * (the thing we do LU decomposition to, to solve)
 */
static void alloc_matrix(unsigned space)
{
 if (space > SEGSIZ/sizeof(double))	/* attempt to trap allocation of    */
    error(bERROR, e_om, "seg size");	/* single matrix large enough to    */
    					/* cause segmentation problems	    */
					/* also traps space*sizeof(double) */
					/* overflow inside calloc.	    */
 reals = (double*)calloc(space, sizeof(double));
 if (!reals)
    error(bERROR, e_om, "reals");
 imags = (double*)calloc(space, sizeof(double));
 if (!imags)
    error(bERROR, e_om, "imags");
}
/*--------------------------------------------------------------------------*/
/* alloc_constants: allocate constant vectors
 * (constant is not the right word...)  allocate space for the right-side
 * vector (initially current sources, on solution becomes voltages) and copies
 * used for one-time-ago and convergence checking
 */
static void alloc_constants(void)
{
 ns.i = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
 if (!ns.i)
    error(bERROR, e_om, "i,acx");
 ns.acx = ns.i;

 ns.v0 = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
 if (!ns.v0)
    error(bERROR, e_om, "v0,acy");
 ns.acy = ns.v0;

 ns.vt1 = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
 if (!ns.vt1)
    error(bERROR, e_om, "vt1");

 ns.fw = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
 if (!ns.fw)
    error(bERROR, e_om, "fw");
}
/*--------------------------------------------------------------------------*/
/* alloc_pointer_table: allocate tables for sparse matrix indexing
 */
static void alloc_pointer_table(void)
{
 redia = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
 if (!redia)		/* pointers to diagonal */
    error(bERROR, e_om, "redia");
 rerow = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
 if (!rerow)		/* row pointers: points to col 0 of each row */
    error(bERROR, e_om, "rerow");
 recol = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
 if (!recol)		/* column pointers: points to row 0 */
    error(bERROR, e_om, "recol");

 imdia = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
 if (!imdia)
    error(bERROR, e_om, "imdia");
 imrow = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
 if (!imrow)
    error(bERROR, e_om, "imrow");
 imcol = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
 if (!imcol)
    error(bERROR, e_om, "imcol");
}
/*--------------------------------------------------------------------------*/
/* fill_pointer_table: build tables for sparse matrix access
 * there are 3 tables: row, column and diagonal
 * row pointer points to each row, where column 0 would be if it existed,
 * even though col 0 doesn't actually exist.
 * Matrix is stored by row below the diagonal, and by column, backwards, above,
 * in such a way that the diagonal is in both the row and column.
 */
static void fill_pointer_table(void)
{
 double *repnt, *impnt;		    /* running location of place in array   */
 int node;

 repnt = reals;					/* fill in pointer tables   */
 impnt = imags;
 for (node = 0;   node <= stats.total_nodes;   ++node){
    recol[node] = repnt - ns.basnode[node];	/* pointers to the bogus    */
    rerow[node] = recol[node] + 2*node;		/* row or column 0	    */
    redia[node]	= recol[node] + node;
    repnt = &repnt[2*(node-ns.basnode[node])+1];/* Next free slot.	    */

    imcol[node] = impnt - ns.basnode[node];
    imrow[node] = imcol[node] + 2*node;
    imdia[node] = imcol[node] + node;
    impnt = &impnt[2*(node-ns.basnode[node])+1];
 }
}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
