/* xsolve.c  93.12.21
 * Copyright 1983-1992   Albert Davis
 * LU decomposition and fwd-back substitution for complex numbers
 * 2 modules:  xlu() and xsolve()
 * xlu() does LU decomposition using the Crout algorithm
 * xsolve() does forward and back substitution after calling xlu() if necessary
 * For a readable version of this, see the real number version: solve
 */
#include "ecah.h"
#include "array.h"
#include "error.h"
#include "nodestat.h"
#include "options.h"
#include "status.h"
#include "declare.h"
/*--------------------------------------------------------------------------*/
	void	  xsolve(void);
	void	  xlu(void);
static	complex_t dotprod(int,int,int);
/*--------------------------------------------------------------------------*/
extern const struct options opt;
extern       struct status stats;
extern       struct nodestuff ns;
/*--------------------------------------------------------------------------*/
void xsolve(void)
{
 int jj, ii;

 time_start(&(stats.back));
 for (ii = 1; ii <= stats.total_nodes; ++ii){	    /* forward substitution */
    complex_t acc;
    acc.x = ns.acx[ii];
    acc.y = ns.acy[ii];
    for ( jj=ns.basnode[ii]; jj<ii; ++jj ){
       /* con[ii] -= Lower[ii][jj] * con[jj] */
       acc = csub(acc, cxmul(*Rel(ii,jj),*Iml(ii,jj), ns.acx[jj],ns.acy[jj]));
    }

    /*  con[ii] /= *xy(ii,ii)   */
    acc = cxdiv(acc.x,acc.y, *Red(ii,ii),*Imd(ii,ii));
    ns.acx[ii] = acc.x;
    ns.acy[ii] = acc.y;
 }
 for (jj = stats.total_nodes; jj > 1; --jj)	    /* back substitution    */
    for (ii = ns.basnode[jj];   ii < jj;   ++ii){
       /* con[ii] -= *Upper[ii][jj] * con[jj] */
       complex_t acc;
       acc = cxmul(*Reu(ii,jj),*Imu(ii,jj), ns.acx[jj],ns.acy[jj]);
       acc = cxsub(ns.acx[ii],ns.acy[ii], acc.x,acc.y);
       ns.acx[ii] = acc.x;
       ns.acy[ii] = acc.y;
    }
 time_stop(&(stats.back));
}
/*--------------------------------------------------------------------------*/
void xlu(void)
{
 complex_t acc;
 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){
       /* *Upper(bn,mm) /= *Diag(bn,bn); */
       acc = cxdiv(*Reu(bn,mm),*Imu(bn,mm), *Red(bn,bn),*Imd(bn,bn));
       *Reu(bn,mm) = acc.x;
       *Imu(bn,mm) = acc.y;

       for (ii =bn+1;  ii<mm;  ii++){
	  acc = dotprod(ii,mm,ii);
	  acc = cxsub(*Reu(ii,mm),*Imu(ii,mm), acc.x,acc.y);
	  acc = cxdiv(acc.x,acc.y, *Red(ii,ii),*Imd(ii,ii));
	  *Reu(ii,mm) = acc.x;
	  *Imu(ii,mm) = acc.y;
       }
       for (jj = bn+1;  jj<mm;  jj++){
	  acc = dotprod(mm,jj,jj);
	  acc = cxsub(*Rel(mm,jj),*Iml(mm,jj), acc.x,acc.y);
	  *Rel(mm,jj) = acc.x;
	  *Iml(mm,jj) = acc.y;
       }
       /* jj == mm */{
	  acc = dotprod(mm,mm,mm);
	  acc = cxsub(*Red(mm,mm),*Imd(mm,mm), acc.x,acc.y);
	  if (acc.x == 0.  &&  acc.y == 0.){
	     error(bWARNING, "open circuit: node %u\n", mm);
	     acc.x = opt.gmin;
	  }
	  *Red(mm,mm) = acc.x;
	  *Imd(mm,mm) = acc.y;
       }
    }else{    /* bn == mm */
       if (*Red(mm,mm)==0. && *Imd(mm,mm)==0.){
	  if (ns.acx[mm]!=0. || ns.acy[mm]!=0.)
	     error(bWARNING, "open circuit: node %u\n", mm);
	  *Red(mm,mm) = opt.gmin;
       }
    }
 }
 time_stop(&(stats.lu));
}
/*--------------------------------------------------------------------------*/
static complex_t dotprod(int ii, int jj, int mm)
{
 int ind, len;
 int kk;
 double *xrow, *xcol;
 double *yrow, *ycol;
 complex_t dot;

 dot.x = dot.y = 0.;
 kk = (ns.basnode[ii]>ns.basnode[jj]) ? ns.basnode[ii] : ns.basnode[jj];
 len = mm-kk;
 if (len > 0){
    xrow = Rel(ii,kk);
    xcol = Reu(kk,jj);
    yrow = Iml(ii,kk);
    ycol = Imu(kk,jj);
    /* for (i = kk;   i < mm;   i++) */
    for (ind = 0;   ind < len;   ind++){
       dot.x += (xrow[-ind] * xcol[ind]) - (yrow[-ind] * ycol[ind]);
       dot.y += (xrow[-ind] * ycol[ind]) + (yrow[-ind] * xcol[ind]);
       /*dot = cadd(dot, cxmul(xrow[-ind],yrow[-ind],xcol[ind],ycol[ind]));*/
    }
 }
 return dot;
}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
