1module bareiss; % Inversion routines using the Bareiss 2-step method.
2
3% Author: Anthony C. Hearn.
4% Modifications by: David Hartley.
5
6% Redistribution and use in source and binary forms, with or without
7% modification, are permitted provided that the following conditions are met:
8%
9%    * Redistributions of source code must retain the relevant copyright
10%      notice, this list of conditions and the following disclaimer.
11%    * Redistributions in binary form must reproduce the above copyright
12%      notice, this list of conditions and the following disclaimer in the
13%      documentation and/or other materials provided with the distribution.
14%
15% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
16% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
17% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
18% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
19% CONTRIBUTORS
20% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
26% POSSIBILITY OF SUCH DAMAGE.
27%
28
29
30% This module is rife with essential references to RPLAC-based
31% functions.
32
33fluid '(!*exp asymplis!* subfg!* wtl!* !*trsparse powlis!* powlis1!*
34        bareiss!-step!-size!*);   % !*solveinconsistent
35
36global '(assumptions requirements);
37
38bareiss!-step!-size!* := 2;     % Seems fastest on average.
39
40symbolic procedure matinverse u;
41   lnrsolve(u,generateident length u);
42
43symbolic procedure lnrsolve(u,v);
44   %U is a matrix standard form, V a compatible matrix form.
45   %Value is U**(-1)*V.
46   begin scalar temp,vlhs,vrhs,ok,
47                !*exp,!*solvesingular;
48   if !*ncmp then return clnrsolve(u,v);
49   !*exp := t;
50   if asymplis!* or wtl!* then
51    <<temp := asymplis!* . wtl!*;
52      asymplis!* := wtl!* := nil>>;
53   vlhs := for i:=1:length car u collect intern gensym();
54   vrhs := for i:=1:length car v collect intern gensym();
55   u := car normmat augment(u,v);
56   v := append(vlhs,vrhs);
57   ok := setkorder v;
58   u := foreach r in u collect prsum(v,r);
59   v := errorset!*({function solvebareiss, mkquote u,mkquote vlhs},t);
60   if caar v memq {'singular,'inconsistent} then
61      <<setkorder ok; rerror(matrix,13,"Singular matrix")>>;
62   v := pair(cadr s,car s) where s = cadar v;
63   u := foreach j in vlhs collect
64           coeffrow(negf numr q,vrhs,denr q) where q = cdr atsoc(j,v);
65   setkorder ok;
66   if temp then <<asymplis!* := car temp; wtl!* := cdr temp>>;
67   return for each j in u collect
68             for each k in j collect
69                if temp then resimp k else cancel k;
70   end;
71
72symbolic procedure prsum(kl,cl);
73   % kl: list of kernel, cl: list of sf -> prsum: sf
74   % kl and cl assumed to have same length
75   if null kl then nil
76   else if null car cl then prsum(cdr kl,cdr cl)
77   else car kl .** 1 .* car cl .+ prsum(cdr kl,cdr cl);
78
79symbolic procedure solvebareiss(exlis,varlis);
80   % exlis: list of sf, varlis: list of kernel
81   % -> solvebareiss: tagged solution list
82   % Solve linear system exlis for variables in varlis using multi-step
83   % Bareiss elimination and fraction-free back-substitution.  The
84   % equations in exlis are not converted to a matrix, but kept as
85   % (sparse) standard forms.
86   begin
87   % if asymplis!* or wtl!* then
88   %   <<temp := asymplis!* . wtl!*; asymplis!* := wtl!* := nil>>;
89   exlis := sparse_bareiss(exlis,varlis,bareiss!-step!-size!*);
90   if car exlis = 'inconsistent then return 'inconsistent . nil;
91   exlis := cdr exlis;
92   if not !*solvesingular and length exlis < length varlis then
93      return 'singular . nil;
94   if !*trsparse then
95      solvesparseprint("Reduced system",reverse exlis,varlis);
96   exlis := sparse_backsub(exlis,varlis);
97   varlis := foreach p in exlis collect car p;
98   exlis := foreach p in exlis collect cdr p;
99   % if temp then <<asymplis!* := car temp; wtl!* := cdr temp>>;
100   exlis := for each ex in exlis collect resimp subs2!* ex;
101   return t . {{exlis,varlis,1}};
102   end;
103
104symbolic procedure coeffrow(u,v,d);
105   % u:sf, v:list of kernel, d:sf -> coeffrow: list of sq
106   % u is linear homogeneous in the kernels in v
107   if null v then nil
108   else if null u or mvar u neq car v
109    then (nil ./ 1) . coeffrow(u,cdr v,d)
110   else (lc u ./ d) . coeffrow(red u,cdr v,d);
111
112
113symbolic procedure augment(u,v);
114   if null u then nil else append(car u,car v) . augment(cdr u,cdr v);
115
116symbolic procedure generateident n;
117  %returns matrix canonical form of identity matrix of order N.
118   begin scalar u,v;
119        for i := 1:n do
120         <<u := nil;
121           for j := 1:n do u := ((if i=j then 1 else nil) . 1) . u;
122           v := u . v>>;
123        return v
124   end;
125
126symbolic procedure normmat u;
127   %U is a matrix standard form.
128   %Value is dotted pair of matrix polynomial form and factor.
129   begin scalar x,y,z;
130      x := 1;
131      for each v in u do
132         <<y := 1;
133           for each w in v do y := lcm(y,denr w);
134           z := (for each w in v
135                    collect multf(numr w,quotf(y,denr w)))
136              . z;
137           x := multf(y,x)>>;
138      return reverse z . x
139   end;
140
141symbolic procedure sparse_bareiss(u,v,k);
142   % u: list of sf, v: list of kernel, k: posint
143   % -> sparse_bareiss: (t|'inconsistent) . list of sf
144   % Multi-step Bareiss elimination using exterior multiplication to
145   % calculate and organise determinants efficiently. Individual blocks
146   % are solved using Cramer's rule. Exterior forms are decomposed into
147   % {constant,linear} parts in non-pivot variables (non-linear part is
148   % not needed). The leading coefficient of the first expression
149   % returned is the determinant of the system.
150   begin scalar p,d,w,pivs,s,asymplis!*,powlis!*,powlis1!*,wtl!*;
151   d := 1;
152   u := foreach f in u join if f then {!*sf2ex(f,v)};
153   while p := choose_pivot_rows(u,v,k,d) do
154      begin
155      u := car p; v := cadr p;  % throws out free vars as well
156      p := cddr p;
157      pivs := lpow car p;       % pivot variables
158      u := foreach r in u join  % multi-step elim. on remaining rows
159              begin
160              r := splitup(r,v);
161              r := extadd(extmult(cadr r,car p),extmult(car r,cadr p));
162              if null (r := subs2chkex r) then return nil;
163              r := innprodpex(pivs,quotexf!*(r,d));
164              % since we did r := r^pivs and then r := pivs _| r,
165              % sign has changed if degree(pivs) is odd
166              if not evenp length pivs then r := negex r;
167              return {r};
168              end;
169      d := lc car p;            % update divisor
170      assumptions := 'list . mk!*sq !*f2q d .
171                            (pairp assumptions and cdr assumptions);
172      p := extadd(car p,cadr p);% recombine pivot rows
173      s := evenp length pivs;
174      foreach x in pivs do      % Cramer's rule on pivot rows
175         w := if (s := not s) then innprodpex(delete(x,pivs),p) . w
176              else negex innprodpex(delete(x,pivs),p) . w;
177      end;
178   foreach f in u do % inconsistent system
179      requirements := 'list . mk!*sq !*f2q !*ex2sf f .
180                          (pairp requirements and cdr requirements);
181   return if u then 'inconsistent . nil
182          else t . foreach f in w collect !*ex2sf f;
183   end;
184
185symbolic procedure choose_pivot_rows(u,v,k,d);
186   % u: list of ex, v: list of kernel, k: posint, d: sf
187   % -> choose_pivot_rows: nil or (list of ex).(list of kernel).ex
188   % Choose pivots in the first k variables from v (or the first k-1
189   % variables from the first pivot variable in v). If k pivots can't be
190   % found, don't waste time looking in further columns (so number of
191   % pivot rows is <= k).  If pivots found, return remaining rows,
192   % remaining variables and decomposed exterior product of pivot rows.
193   if null u or null v then nil else
194   begin scalar w,s,ss,p,x,y,rows,pivots;
195   w := u;
196   for i:=1:k do if v then v := cdr v;
197   while k neq 0 do
198      if null u then % ran out of rows before finding k pivots
199         if null v or null w or pivots then k := 0
200         else % skip k more variables and reset everything
201          << for i:=1:k do if v then v := cdr v; s := nil; u := w>>
202      else if car(x := splitup(car u,v)) and
203              (y := if null pivots then car x
204                    else subs2chkex extmult(car x,car pivots)) then
205         begin % found one
206         rows := x . rows;
207         pivots := (if null pivots then y else quotexf!*(y,d)) . pivots;
208         % if # rows skipped is odd, then reverse sign
209         if s then ss := not ss;
210         w := delete(car u,w); u := cdr u; k := k - 1;
211         end
212      else <<u := cdr u; s := not s>>; % skip row
213   if null pivots then return; % couldn't find any pivots
214   % next line adjusts sign to return row1^...^rowk
215   % instead of rowk^...^row1
216   if remainder(length lpow car pivots,4) member {2,3} then
217      ss := not ss;
218   rows := reverse rows;       % calculate dets along pivot rows
219   pivots := reverse pivots;
220   p := car rows;
221   foreach r in cdr rows do
222      p := {car(pivots := cdr pivots),
223            quotexf!*(extadd(extmult(cadr r,car p),
224                             extmult(car r,cadr p)),d)};
225   return w . v . if ss then {negex car p,negex cadr p} else p;
226   end;
227
228symbolic procedure sparse_backsub(exlis,varlis);
229   % exlis: list of sf, varlis: list of kernel
230   % -> sparse_backsub: list of kernel.sq
231   % Fraction-free back-substitution for exlis, where reverse exlis is
232   % a list of rows of an upper-triangular matrix wrt varlis.  Since
233   % exlis has been produced in a fraction-free way, the leading
234   % coefficient of the first row is the determinant of the system.
235   begin scalar d,z,c;
236   if null exlis then return nil;       % trivial case
237   d := lc car exlis;                  % determinant
238   foreach x in exlis do % Almost redundant for first x.
239      begin scalar s,p,v,r;
240      p := lc x;   % pivot.
241      v := mvar x;
242      x := red x;
243      while not domainp x and mvar x member varlis do
244          <<if (c := atsoc(mvar x,z)) then % back-substitute
245               s := addf(multf(lc x,cdr c),s)
246            else % move free variable terms to rhs
247               r := addf(!*t2f lt x,r);
248            x := red x>>;
249      % Used to be quotf!*, but that could give a terminal error.
250      s := negf quotff(addf(multf(addf(r,x),d),s),p);
251      z := (v . s) . z;
252      end;
253   for each p in z do cdr p := cancel(cdr p ./ d);
254   return z
255   end;
256
257symbolic procedure quotff(u,v);
258   % We do the rationalizesq step to allow for surd divisors.
259   if null u then nil
260    else (if x then x
261           else (if denr y = 1 then numr y
262                  else rederr "Invalid division in backsub")
263                 where y=rationalizesq(u ./ v))
264          where x=quotf(u,v);
265
266endmodule;
267
268end;
269