/*
extern cpoly csolvlin(int cwnin,complexmp cwwmin[cmaxsizepad][cmaxsizepad],complexscalar cwwrin[cmaxsizepad]);
*/

extern cpoly csolvlin(int cwnin,complexmp cwwmin[cmaxsizepad][cmaxsizepad],complexscalar cwwrin[cmaxsizepad]);
extern int csolvlin(int cwnin,int * * * *cwmin,int * * *cwrin);

int cwnin;
int cwninmax;
int cwn;

complexvector cwx_ram;unsigned int * * *cwx=cwx_ram.cv;
complexvector cwr_ram;unsigned int * * *cwr=cwr_ram.cv;
complexvector cwrin_ram;unsigned int * * *cwrin=cwrin_ram.cv;

complexmatrix cwm_ram;unsigned int * * * *cwm=cwm_ram.cm; 
//complexmatrix cwmin_ram;unsigned int * * * *cwmin=cwmin_ram.cm;

unsigned int * * * *cwmin=cwm_ram.cm; //this can replace line above but
                                     //destroys input matrix--no longer
                                     //relevant with overloads.
cpoly czpolyout;


complexmp cwpolyout[cmaxsizepad];
cwpoly cwzpolyout; 

cwpoly csolvlin(int cwnin,complexmp cwwmin[cmaxsizepad][cmaxsizepad],complexmp cwwrin[cmaxsizepad])
{int cwi,cwj,cwk,cwkk,cwl;
 cwn=cwnin;
 for (cwk=0;cwk<prec+2;cwk++) 
   {for (cwkk=0;cwkk<2;cwkk++)
     {for (cwi=1;cwi<cwn+1;cwi++)
        {cwr[cwi][cwkk][cwk]=cwwrin[cwi].gv(cwkk,cwk);
	   for (cwj=1;cwj<cwn+1;cwj++)
             {cwm[cwi][cwj][cwkk][cwk]=cwwmin[cwi][cwj].gv(cwkk,cwk);}
	}
      }
   }
                 //insert old code
    for (cwk=2;cwk<cwn+1;cwk++)
     {//cwl=cwk-1;
      for (cwi=cwk;cwi<cwn+1;cwi++)
        {cwl=cwk-1;                      //moved up
	 cmult(cwm[cwi][cwl],cwr[cwl]);
	 cdiv(cz,cwm[cwl][cwl]);
         csub(cz,cwr[cwi]);
	 cchs(cz);
	 cstore(cwr[cwi]);
	 for (cwj=cwk;cwj<cwn+1;cwj++)
	   {cmult(cwm[cwi][cwl],cwm[cwl][cwj]);
            cdiv(cz,cwm[cwl][cwl]);
            csub(cz,cwm[cwi][cwj]);
	    cchs(cz);
	    cstore(cwm[cwi][cwj]);
	   }
	}
     }
  cdiv(cwr[cwn],cwm[cwn][cwn]);
  cstore(cwx[cwn]);
  for (cwi=2;cwi<cwn+1;cwi++)
     {cwk=cwn-cwi+1;
      crecall(cwr[cwk]);
      cstore(cwx[cwk]);
      for (cwj=1;cwj<cwi;cwj++)
	{cwl=cwn-cwj+1;
	 cmult(cwm[cwk][cwl],cwx[cwl]);
	 csub(cz,cwx[cwk]);
	 cchs(cz);
	 cstore(cwx[cwk]);
	}
      cdiv(cwx[cwk],cwm[cwk][cwk]);
      cstore(cwx[cwk]);
     }
 for (cwk=0;cwk<prec+2;cwk++) 
   {for (cwkk=0;cwkk<2;cwkk++)
     {for (cwi=1;cwi<cwn+1;cwi++)
        {cwpolyout[cwi].sv(cwkk,cwk,cwx[cwi][cwkk][cwk]);}
     }
   }
  return cwzpolyout;
}


               //old code kept for backward compatibility
int csolvlin(int cwnin,unsigned int * * * *cwmin,unsigned int * * *cwrin)
{int cwi,cwj,cwk,cwkk,cwl;
 cwn=cwnin;
 for (cwk=0;cwk<prec+2;cwk++)                      //oops
   {for (cwkk=0;cwkk<2;cwkk++)
     {for (cwi=1;cwi<cwn+1;cwi++)
	{cwr[cwi][cwkk][cwk]=cwrin[cwi][cwkk][cwk];
	   for (cwj=1;cwj<cwn+1;cwj++)
	     {cwm[cwi][cwj][cwkk][cwk]=cwmin[cwi][cwj][cwkk][cwk];}
	}
      }
   }
    for (cwk=2;cwk<cwn+1;cwk++)
     {//cwl=cwk-1;
      for (cwi=cwk;cwi<cwn+1;cwi++)
        {cwl=cwk-1;                      //moved up
	 cmult(cwm[cwi][cwl],cwr[cwl]);
	 cdiv(cz,cwm[cwl][cwl]);
         csub(cz,cwr[cwi]);
	 cchs(cz);
	 cstore(cwr[cwi]);
	 for (cwj=cwk;cwj<cwn+1;cwj++)
	   {cmult(cwm[cwi][cwl],cwm[cwl][cwj]);
            cdiv(cz,cwm[cwl][cwl]);          //36,36 is 0
            csub(cz,cwm[cwi][cwj]);
	    cchs(cz);
	    cstore(cwm[cwi][cwj]);
	   }
	}
     }
  cdiv(cwr[cwn],cwm[cwn][cwn]);
  cstore(cwx[cwn]);
  for (cwi=2;cwi<cwn+1;cwi++)
     {cwk=cwn-cwi+1;
      crecall(cwr[cwk]);
      cstore(cwx[cwk]);
      for (cwj=1;cwj<cwi;cwj++)
	{cwl=cwn-cwj+1;
	 cmult(cwm[cwk][cwl],cwx[cwl]);
	 csub(cz,cwx[cwk]);
	 cchs(cz);
	 cstore(cwx[cwk]);
	}
      cdiv(cwx[cwk],cwm[cwk][cwk]);
      cstore(cwx[cwk]);
     }
  return 0;
}