module solvelnr; % Code for solving a general system of linear eqns.

% Authors: Anthony C. Hearn and Eberhard Schruefer.
% Modifications by: David Hartley.

% Based on code by David R. Stoutemyer modified by Donald R. Morrison.

% Copyright (c) 1993 RAND.  All rights reserved.

% The number of equations and the number of unknowns are arbitrary.
% I.e. the system can be under- or overdetermined.

fluid '(!*cramer !*exp !*solvesingular asymplis!* wtl!* 
	!*arbvars !*trsparse !*varopt bareiss!-step!-size!*);

% switch solveinconsistent;

% !*solveinconsistent := t; % Default value.

symbolic procedure solvelnrsys(exlis,varlis);
   % exlis: list of sf, varlis: list of kernel
   % -> solvelnrsys: tagged solution list 
   % Check the system for sparsity, then decide whether to use the
   % Cramer or Bareiss method.  Using the Bareiss method on sparse
   % systems, 4-step elimination seems to be faster than 2-step.
   % The Bareiss code is not good at handling surds at the moment,
   % hence exptexpflistp test.
   begin scalar w,method;
   if w := solvesparsecheck(exlis,varlis) then exlis := w
     else exlis := exlis . varlis;
   if null !*cramer and null exptexpflistp car exlis
      then method := 'solvebareiss
     else method := 'solvecramer;
   exlis := apply2(method,car exlis,cdr exlis)
               where bareiss!-step!-size!* = if w then 4 else 2;
   return solvesyspost(exlis,varlis);
   end;


symbolic procedure exptexpflistp u;
   %  True if any of u contains an expt kernel.
   u and (exptexpfp car u or exptexpflistp cdr u);


symbolic procedure solvesyspost(exlis,varlis);
   % exlis: tagged solution list, varlis: list of kernel
   %  -> solvesyspost: tagged solution list
   % Insert arbitrary constants and present
   % solutions in same order as in varlis.
   % Also reorders expressions to prevailing kernel order.
   car exlis . foreach s in cdr exlis collect
      if car s and null cadr s then s else
      begin scalar arbvars,z;
      if !*arbvars or (null cadr s and length varlis = 1) then
         arbvars := foreach v in setdiff(varlis,cadr s) collect
                       v . mvar makearbcomplex()
      else
         varlis := intersection(varlis,cadr s);
      z := pair(cadr s,sublis(arbvars,car s));
      z := append(z,foreach p in arbvars collect car p . !*k2q cdr p);
      return {foreach v in varlis collect reordsq cdr atsoc(v,z),
              varlis,caddr s};
      end;

symbolic procedure solvecramer(exlis,varlis);
   % exlis: list of sf, varlis: list of kernel
   % -> solvecramer: tagged solution list 
   % Just a different name at the moment.
   glnrsolve(exlis,varlis);

symbolic procedure solvesparsecheck(sys,vl);
   % sys: list of sf, vl: list of kernel
   % -> solvesparsecheck: nil or {list of sf,list of kernel}
   % This program checks for a sparse linear system. If the
   % system is sparse enough, it returns (exlis.varlis) reordered
   % such that a maximum triangular upper diagonal form is
   % established. Otherwise the result is NIL.
   begin scalar vl1,xl,sys1,q,x,y;
         integer sp;    

   % First collect a list vl1 where each variable is followed
   % by the number of equations where it occurs, and then
   % by the number of other variables which occur in these
   % equations (connectivity). At the same time, collect a measure
   % of the sparsity.
   sp:=0;
   vl1:= for each x in vl collect x . 0 . nil;
   foreach q in sys do
      foreach x in (xl := intersection(topkerns q,vl)) do
       <<y := assoc(x,vl1);
         cdr y := (cadr y + 1) . union(xl,cddr y);
         sp := sp + 1>>;
   foreach p in vl1 do
      cddr p := length cddr p - 1; % could drop the -1

   % Drop out if density > 80%
   if sp > length sys * length vl * 0.8 then 
    <<if !*trsparse then prin2t "System is not very sparse";
      return nil>>;

   % If varopt is on, sort variables first by least occurrences and then
   % least connectivity, but allow dependency to override.
   % Reset kernel order and reorder equations.
   if !*trsparse then
     solvesparseprint("Original sparse system",sys,vl);
   
   if !*varopt then
   << vl1 := sort(vl1,function solvevarordp);
      vl1 := foreach x in vl1 collect car x;
      vl1 := solvevaradjust vl1;
      foreach k in reverse vl1 do updkorder k;
      sys := for each q in sys collect reorder q >>
   else vl1 := foreach x in vl1 collect car x;

   % Next sort equations in ascending order of their first variable
   % and then descending order of the next variable.
   sys1:= (nil . nil) . foreach x in vl1 collect x . nil;
   foreach q in sys do
      <<if domainp q or not member(mvar q,vl1) then y := assoc(nil,sys1)
        else y := assoc(mvar q,sys1);
        cdr y := q . cdr y>>;
   foreach p in cdr sys1 do
      if cdr p then cdr p := sort(cdr p, function solvesparsesort);

   % Finally split off a leading diagonal system and push the remaining
   % equations down.
   sys := nconc(foreach p in sys1 join if cdr p then {cadr p},
                reversip foreach p in sys1 join if cdr p then cddr p);
   if !*trsparse then
      solvesparseprint("Variables and/or equations rearranged",sys,vl1);
   return sys . vl1
   end;

symbolic procedure solvevarordp(x,y);
   cadr x < cadr y or
   cadr x = cadr y and cddr x < cddr y;

symbolic procedure solvevaradjust u;
   % u:list of kernel -> solvevaradjust:list of kernel
   % Adjust ordering of u to respect dependency ordering by recursively
   % putting variables with no dependencies last
   begin scalar v,y;
   if null u then return nil;
   v := foreach x in u join
      	<< y := assoc(x,depl!*);
 	   if null y or null xnp(cdr y,u) then {x} >>;
   return nconc(solvevaradjust setdiff(u,v),v)
   end;

symbolic procedure solvesparseprint(text,sys,vl);
   <<terpri(); prin2t text;
     for each e in sys do
       << e := topkerns e;
	  for each x in vl do
	  if memq(x,e) then prin2 "*"  else prin2 "-";
	  terpri()>>>>;

symbolic procedure topkerns u;
   % u:sf -> topkerns:list of kernel
   % kernels in top level of u
   if domainp u then nil else mvar u . topkerns red u;

symbolic procedure solvesparsesort(x,y);
   % x,y: sf -> solvesparsesort: bool
   if domainp x then nil
   else if domainp y then t
   else if mvar x = mvar y then solvesparsesort(red y,red x)
   else if ordop(mvar x,mvar y) then t
   else nil;

endmodule;

end;
