/* solve.c  94.02.11
 * Copyright 1983-1992   Albert Davis
 * LU decomposition and fwd-back substitution for real numbers
 * 2 modules:  lu() and solve()
 * lu() does LU decomposition using the Crout algorithm (modified)
 * solve() does forward and back substitution after calling lu() if necessary
 */
#include "ecah.h"
#include "array.h"
#include "error.h"
#include "mode.h"
#include "nodestat.h"
#include "options.h"
#include "status.h"
#include "declare.h"
/*--------------------------------------------------------------------------*/
	void	solve(void);
static	void	lu(void);
static	void	lu_partial(void);
static	void	lu_full(void);
static	double	dotprod(int,int,int);
/*--------------------------------------------------------------------------*/
extern const struct options opt;
extern       struct status stats;
extern       struct nodestuff ns;
extern struct nodestat *nstat;
extern const double trtime0;
extern const int inc_mode;	    /* make incremental changes		    */
/*--------------------------------------------------------------------------*/
void solve(void)
{
 int jj, ii;

 lu();
 time_start(&(stats.back));
 for (ii = 1;  ii <= stats.total_nodes;  ++ii){	    /* forward substitution */
    ns.fw[ii] = ns.i[ii];
    for (jj = ns.basnode[ii];   jj < ii;   ++jj){
       ns.fw[ii] -= *LUl(ii,jj) * ns.fw[jj];
    }
    ns.fw[ii] /= *LUd(ii,ii);
 }
 if (opt.foooo!=4){
    for (ii = stats.total_nodes;  ii >= 1;  --ii){
	  ns.v0[ii] = ns.fw[ii];
	  nstat[ii].aiter = stats.iter[iTOTAL];
    }
    for (jj = stats.total_nodes;  jj > 1;  --jj){   /* back substitution    */
       for (ii = ns.basnode[jj];   ii < jj;   ++ii){
	   ns.v0[ii] -= *LUu(ii,jj) * ns.v0[jj];
       }
    }
 }else{
    for (ii = stats.total_nodes;  ii >= 1;  --ii){
       /*if (nstat[ii].mode != mGATE){*/
	  ns.v0[ii] = ns.fw[ii];
	  for (jj = ii+1;  jj <= stats.total_nodes;  ++jj){
	     ns.v0[ii] -= *im(ii,jj) * ns.v0[jj];
	  }
	  nstat[ii].aiter = stats.iter[iTOTAL];
	  /*nstat[ii].mode = mLU;*/
	  /*nstat[ii].lastchange = trtime0;*/
       /*}*/
    }
 }
 time_stop(&(stats.back));
 if (opt.foooo == 1){
    for (ii=1; ii<=stats.total_nodes; ++ii){
       for (jj=1; jj<=stats.total_nodes; ++jj){
	  printf("%9.2e ", *re(ii,jj));
       }
       printf("%9.2e\n", ns.i[ii]);
    }
 }
}
/*--------------------------------------------------------------------------*/
static void lu(void)
{
 if (opt.lubypass && inc_mode){
    lu_partial();
 }else{
    lu_full();
 }
}
/*--------------------------------------------------------------------------*/
static void lu_partial(void)
{
 int ii, jj, mm;
 int bn;
 int prop = 0;
 
 time_start(&(stats.lu));
 for (mm = 1;   mm <= stats.total_nodes;   ++mm){
    bn = ns.basnode[mm];
    if (nstat[mm].needslu  ||  bn <= prop){
       nstat[mm].needslu = NO;
       prop = mm;
       if (bn < mm){
	  prop = mm;
	  *LUu(bn,mm) = *AAu(bn,mm) / *LUd(bn,bn);
	  for (ii = bn+1;  ii<mm;  ii++){
	     *LUu(ii,mm) = (*AAu(ii,mm) - dotprod(ii,mm,ii)) / *LUd(ii,ii);
	  }
	  *LUl(mm,bn) = *AAl(mm,bn);
	  for (jj = bn+1;  jj<mm;  jj++){
	     *LUl(mm,jj) = *AAl(mm,jj) - dotprod(mm,jj,jj);
	  }
	  /* jj == mm */{
	     *LUd(mm,mm) = *AAd(mm,mm) - dotprod(mm,mm,mm);
	     if (opt.gmin > *LUd(mm,mm)  &&  *LUd(mm,mm) > -opt.gmin){
		error(bWARNING, "open circuit: internal node %u\n", mm);
		*LUd(mm,mm) = opt.gmin;
	     }
	  }
       }else{    /* bn == mm */
	  *LUd(mm,mm) = *AAd(mm,mm);
	  if (opt.gmin > *LUd(mm,mm)  &&  *LUd(mm,mm) > -opt.gmin){
	     error(bWARNING, "open circuit: internal node %u\n", mm);
	     *LUd(mm,mm) = opt.gmin;
	  }
       }
    }
 }
 time_stop(&(stats.lu));
}
/*--------------------------------------------------------------------------*/
static void lu_full(void)
{
 int ii, jj, mm;
 int bn;

 time_start(&(stats.lu));
 for (mm = 1;   mm <= stats.total_nodes;   ++mm){
    bn = ns.basnode[mm];
    if (bn < mm){
       *LUu(bn,mm) = *AAu(bn,mm) / *LUd(bn,bn);
       for (ii = bn+1;  ii<mm;  ii++){
	  *LUu(ii,mm) = (*AAu(ii,mm) - dotprod(ii,mm,ii)) / *LUd(ii,ii);
       }
       *LUl(mm,bn) = *AAl(mm,bn);
       for (jj = bn+1;  jj<mm;  jj++){
	  *LUl(mm,jj) = *AAl(mm,jj) - dotprod(mm,jj,jj);
       }
       /* jj == mm */{
	  *LUd(mm,mm) = *AAd(mm,mm) - dotprod(mm,mm,mm);
	  if (opt.gmin > *LUd(mm,mm)  &&  *LUd(mm,mm) > -opt.gmin){
	     error(bWARNING, "open circuit: internal node %u\n", mm);
	     *LUd(mm,mm) = opt.gmin;
	  }
       }
    }else{    /* bn == mm */
       *LUd(mm,mm) = *AAd(mm,mm);
       if (opt.gmin > *LUd(mm,mm)  &&  *LUd(mm,mm) > -opt.gmin){
	  error(bWARNING, "open circuit: internal node %u\n", mm);
	  *LUd(mm,mm) = opt.gmin;
       }
    }
 }
 time_stop(&(stats.lu));
}
/*--------------------------------------------------------------------------*/
static double dotprod(int ii, int jj, int mm)
{
 int ind, len;
 int kk;
 double *row, *col;
 double dot = 0.;

 kk = MAX(ns.basnode[ii], ns.basnode[jj]);
 len = mm - kk;
 if (len > 0){
    row = LUl(ii,kk);
    col = LUu(kk,jj);
    /* for (i = kk;   i < mm;   i++) */
    for (ind = 0;   ind < len;   ind++)
        dot += row[-ind] * col[ind];
 }
 return dot;
}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
