module fide; % FInite difference method for partial Differential Eqn
% systems.
% Author: Richard Liska.
% Version: 1.1.2 for REDUCE 3.6, May 29, 1995.
%***********************************************************************
%** (C) 1991-1995, Richard Liska **
%** Faculty of Nuclear Science and Physical Engineering **
%** Technical University of Prague **
%** Brehova 7 **
%** 115 19 Prague 1 **
%** Czech Republic **
%** Email: Richard Liska <liska@siduri.fjfi.cvut.cz> **
%** This package can be distributed through REDUCE Network Library. **
%***********************************************************************
% The FIDE package consists of the following modules:
%
% EXPRES for transforming PDES into any orthogonal coordinate system.
% IIMET for discretization of PDES by integro-interpolation method.
% APPROX for determining the order of approximation of difference
% scheme
% CHARPOL for calculation of amplification matrix and characteristic
% polynomial of difference scheme, which are needed in Fourier
% stability analysis.
% HURWP for polynomial roots locating necessary in verifying the von
% Neumann stability condition.
% LINBAND for generating the block of FORTRAN code, which solves a
% system of linear algebraic equations with band matrix
% appearing quite often in difference schemes.
%
% Changes since version 1.1:
% Patches in SIMPINTERPOL and SIMPINTT 13/06/91
% Patch in TDIFPAIR 08/07/91
% Two FEXPR routines F2VAL, FPLUS changed to MACROs 17/03/92
% Patches in IIM1, AMPMAT, HURW, CHARPOL for 3.5 01/11/93
% Version 1.1.1 of the FIDE package is the result of porting the FIDE
% package version 1.1.0 to REDUCE 3.5.
% Infix uses of NOT removed (not a memq b -> not(a memq b) ) ACH
% MAP* functions replaced by FOR EACH syntax 16/03/94
% Version 1.1.2 of the FIDE package is the result of porting the FIDE
% package version 1.1.1 to REDUCE 3.6.
create!-package('(fide approx charpol hurwp linband), '(contrib fide));
load!-package 'matrix;
load!-package 'fide1; % We need this loaded.
algebraic;
%***********************************************************************
%***** *****
%***** D A T A FOR DISCRETIZATION (CHANGE IF YOU NEED) *****
%***** *****
%***********************************************************************
difmatch all,1,
0,1$
difmatch all,u,
u=one,0,
u(i),
u=half,0,
(u(i-1/2)+u(i+1/2))/2$
difmatch all,diff(u,x),
u=one,2,
(u(i+1)-u(i-1))/(dip1+dim1),
u=half,0,
(u(i+1/2)-u(i-1/2))/di$
difmatch all,diff(u,x,2),
u=one,0,
((u(i+1)-u(i))/dip1-(u(i)-u(i-1))/dim1)/di,
u=half,2,
((u(i+3/2)-u(i+1/2))/dip2-(u(i-1/2)-u(i-3/2))/dim2)/(dip1+dim1)$
difmatch all,u*v,
u=one,v=one,0,
u(i)*v(i),
u=one,v=half,0,
u(i)*(v(i-1/2)+v(i+1/2))/2,
u=half,v=one,0,
(u(i-1/2)+u(i+1/2))/2*v(i),
u=half,v=half,0,
(u(i-1/2)*v(i-1/2)+u(i+1/2)*v(i+1/2))/2$
difmatch all,u**n,
u=one,0,
u(i)**n,
u=half,0,
(u(i-1/2)**n+u(i+1/2)**n)/2$
difmatch all,u*v**n,
u=one,v=one,0,
u(i)*v(i)**n,
u=one,v=half,0,
u(i)*(v(i-1/2)**n+v(i+1/2)**n)/2,
u=half,v=one,0,
(u(i-1/2)+u(i+1/2))/2*v(i)**n,
u=half,v=half,0,
(u(i-1/2)*v(i-1/2)**n+u(i+1/2)*v(i+1/2)**n)/2$
difmatch all,u*v*w,
u=one,v=one,w=one,0,
u(i)*v(i)*w(i),
u=one,v=one,w=half,0,
u(i)*v(i)*(w(i+1/2)+w(i-1/2))/2,
u=one,v=half,w=one,0,
u(i)*(v(i-1/2)+v(i+1/2))/2*w(i),
u=one,v=half,w=half,0,
u(i)*(v(i-1/2)*w(i-1/2)+v(i+1/2)*w(i+1/2))/2,
u=half,v=one,w=one,0,
(u(i-1/2)+u(i+1/2))/2*v(i)*w(i),
u=half,v=one,w=half,0,
(u(i-1/2)*w(i-1/2)+u(i+1/2)*w(i+1/2))/2*v(i),
u=half,v=half,w=one,0,
(u(i-1/2)*v(i-1/2)+u(i+1/2)*v(i+1/2))/2*w(i),
u=half,v=half,w=half,0,
(u(i-1/2)*v(i-1/2)*w(i-1/2)+u(i+1/2)*v(i+1/2)*w(i+1/2))/2$
difmatch all,v*diff(u,x),
u=one,v=one,2,
v(i)*(u(i+1)-u(i-1))/(dip1+dim1),
u=one,v=half,2,
(v(i+1/2)+v(i-1/2))/2*(u(i+1)-u(i-1))/(dip1+dim1),
u=half,v=one,0,
v(i)*(u(i+1/2)-u(i-1/2))/di,
u=half,v=half,0,
(v(i+1/2)+v(i-1/2))/2*(u(i+1/2)-u(i-1/2))/di$
difmatch all,v*w*diff(u,x),
u=one,v=one,w=one,2,
v(i)*w(i)*(u(i+1)-u(i-1))/(dip1+dim1),
u=one,v=one,w=half,2,
v(i)*(w(i-1/2)+w(i+1/2))/2*(u(i+1)-u(i-1))/(dip1+dim1),
u=one,v=half,w=one,2,
(v(i+1/2)+v(i-1/2))/2*w(i)*(u(i+1)-u(i-1))/(dip1+dim1),
u=one,v=half,w=half,2,
(v(i+1/2)*w(i+1/2)+v(i-1/2)*w(i-1/2))/2*(u(i+1)-u(i-1))/(dip1+dim1),
u=half,v=one,w=one,0,
v(i)*w(i)*(u(i+1/2)-u(i-1/2))/di,
u=half,v=one,w=half,0,
v(i)*(w(i-1/2)+w(i+1/2))/2*(u(i+1/2)-u(i-1/2))/di,
u=half,v=half,w=one,0,
(v(i+1/2)+v(i-1/2))/2*w(i)*(u(i+1/2)-u(i-1/2))/di,
u=half,v=half,w=half,0,
(v(i+1/2)*w(i+1/2)+v(i-1/2)*w(i-1/2))/2*(u(i+1/2)-u(i-1/2))/di$
difmatch all,x*u,
u=one,0,
x(i)*u(i),
u=half,1,
(x(i-1/2)*u(i-1/2)+x(i+1/2)*u(i+1/2))/2$
difmatch all,u/x**n,
u=one,0,
u(i)/x(i)**n,
u=half,0,
(u(i-1/2)/x(i-1/2)**n+u(i+1/2)/x(i+1/2)**n)/2$
difmatch all,u*v/x**n,
u=one,v=one,0,
u(i)*v(i)/x(i)**n,
u=one,v=half,0,
u(i)*(v(i-1/2)+v(i+1/2))/2/x(i)**n,
u=half,v=one,0,
(u(i-1/2)+u(i+1/2))/2*v(i)/x(i)**n,
u=half,v=half,0,
(u(i-1/2)*v(i-1/2)/x(i-1/2)**n+u(i+1/2)*v(i+1/2)/x(i+1/2)**n)/2$
difmatch all,diff(x**n*u,x)/x**n,
u=one,2,
(x(i+1)**n*u(i+1)-x(i-1)**n*u(i-1))/x(i)**n/(dim1+dip1),
u=half,0,
(x(i+1/2)**n*u(i+1/2)-x(i-1/2)**n*u(i-1/2))/di/x(i)**n$
difmatch all,diff(u*v,x),
u=one,v=one,4,
(u(i+1)*v(i+1)-u(i-1)*v(i-1))/(dim1+dip1),
u=one,v=half,2,
((u(i+1)+u(i))/2*v(i+1/2)-(u(i-1)+u(i))/2*v(i-1/2))/di,
u=half,v=one,2,
((v(i+1)+v(i))/2*u(i+1/2)-(v(i-1)+v(i))/2*u(i-1/2))/di,
u=half,v=half,0,
(u(i+1/2)*v(i+1/2)-u(i-1/2)*v(i-1/2))/di$
difmatch all,diff(u*v,x)/x**n,
u=one,v=one,4,
(u(i+1)*v(i+1)-u(i-1)*v(i-1))/x(i)**n/(dim1+dip1),
u=one,v=half,2,
((u(i+1)+u(i))/2*v(i+1/2)-(u(i-1)+u(i))/2*v(i-1/2))/x(i)**n/di,
u=half,v=one,2,
((v(i+1)+v(i))/2*u(i+1/2)-(v(i-1)+v(i))/2*u(i-1/2))/x(i)**n/di,
u=half,v=half,0,
(u(i+1/2)*v(i+1/2)-u(i-1/2)*v(i-1/2))/x(i)**n/di$
difmatch all,diff(u*diff(v,x),x)/x**n,
u=half,v=one,0,
(u(i+1/2)*(v(i+1)-v(i))/dip1-u(i-1/2)*(v(i)-v(i-1))/dim1)/di/x(i)**n,
u=half,v=half,2,
(u(i+1/2)*(v(i+3/2)-v(i-1/2))/(di+dip2)-u(i-1/2)*(v(i+1/2)-
v(i-3/2))/(di+dim2))/di/x(i)**n,
u=one,v=one,2,
((u(i+1)+u(i))/2*(v(i+1)-v(i))/dip1-(u(i)+u(i-1))/2*(v(i)-v(i-1))
/dim1)/di/x(i)**n,
u=one,v=half,4,
((u(i+1)+u(i))/2*(v(i+3/2)-v(i-1/2))/(di+dip2)-
(u(i)+u(i-1))/2*(v(i+1/2)-v(i-3/2))/(di+dim2))/di/x(i)**n$
%***********************************************************************
%***** E N D OF D A T A FOR DISCRETIZATION *****
%***********************************************************************
endmodule;
module approx;
% Author: R. Liska$
% Version REDUCE 3.6 05/1991$
fluid '(!*prapprox)$
switch prapprox$
!*prapprox:=nil$
global '(cursym!* coords!* icoords!* functions!* hipow!* lowpow!*)$
% Implicitely given indices
icoords!*:='(i j k l m n i1 j1 k1 l1 m1 n1)$
algebraic$
procedure fact(n)$
if n=0 then 1
else n*fact(n-1)$
procedure taylor(fce,var,step,ord)$
if step=0 or ord=0 then fce
else fce+for j:=1:ord sum step**j/fact(j)*df(fce,var,j)$
symbolic$
procedure maxorder u$
begin
scalar movar,var$
a:movar:=car u$
if not eqexpr movar then return errpri2(movar,'hold)$
movar:=cdr movar$
var:=car movar$
movar:=reval cadr movar$
if not atom var or not(var memq coords!*) then return msgpri(
" Parameter ",var," must be coordinate",nil,'hold)
else if not fixp movar then return msgpri(
" Parameter ", movar," must be integer",nil,'hold)
else put(var,'maxorder,movar)$
u:=cdr u$
if u then go to a$
return nil
end$
put('maxorder,'stat,'rlis)$
procedure center u$
begin
scalar movar,var$
a:movar:=car u$
if not eqexpr movar then return errpri2(movar,'hold)$
movar:=cdr movar$
var:=car movar$
movar:=reval cadr movar$
if not atom var or not(var memq coords!*) then return msgpri(
" Parameter ",var," must be coordinate",nil,'hold)
else if not(fixp movar or (eqcar(movar,'quotient) and
(fixp cadr movar or
(eqcar(cadr movar,'minus) and fixp cadadr movar))
and fixp caddr movar)) then return msgpri(
" Parameter ", movar," must be integer or rational number",nil,
'hold)
else put(var,'center,movar)$
u:=cdr u$
if u then go to a$
return nil
end$
put('center,'stat,'rlis)$
procedure functions u$
<<functions!* := u$
for each a in u do put(a,'simpfn,'simpiden) >>$
put('functions,'stat,'rlis)$
procedure simptaylor u$
begin
scalar ind,var,movar,step,fce,ifce$
fce:=car u$
if null cdr u then return simp fce$
ifce:=cadr u$
if cddr u then fce:= fce . cddr u$
ind:=mvar numr simp ifce$
var:=tcar get(ind,'coord)$
step:=reval list('difference,
ifce,
list('plus,
if (movar:=get(var,'center)) then movar
else 0,
ind))$
step:=list('times,
step,
get(var,'gridstep))$
movar:=if (movar:=get(var,'maxorder)) then movar
else 3$
return simp list('taylor,
fce,
var,
step,
movar)
end$
algebraic$
procedure approx difsch$
begin
scalar ldifsch,rdifsch,nrcoor,coors,rest,ldifeq,rdifeq,alglist!*$
symbolic
<<for each a in functions!* do
<<put(a,'simpfn,'simptaylor)$
eval list('depend,mkquote (a . coords!*)) >>$
flag(functions!*,'full)$
for each a in coords!* do put(a,'gridstep, intern compress append
(explode 'h,explode a))$
nrcoor:=length coords!* - 1$
eval list('array,
mkquote list('steps . add1lis list(nrcoor)))$
coors:=coords!*$
for j:=0:nrcoor do
<<setel(list('steps,j),aeval get(car coors,'gridstep))$
coors:=cdr coors >> >>$
ldifsch:=lhs difsch$
rdifsch:=rhs difsch$
ldifeq:=ldifsch$
rdifeq:=rdifsch$
ldifeq:=substeps(ldifeq)$
rdifeq:=substeps(rdifeq)$
rest:=ldifsch-ldifeq-rdifsch+rdifeq$
for j:=0:nrcoor do
steps(j):=steps(j)**minorder(rest,steps(j))$
write " Difference scheme approximates differential equation ",
ldifeq=rdifeq$
write " with orders of approximation:"$
on div$
for j:=0:nrcoor do write steps(j)$
off div$
symbolic if !*prapprox
then algebraic write " Rest of approximation : ",rest$
symbolic
<<for each a in functions!* do
<<put(a,'simpfn,'simpiden)$
eval list('nodepend,mkquote (a . coords!*)) >>$
remflag(functions!*,'full)>>$
clear steps
end$
procedure substeps u$
begin
scalar step,nu,du$
nu:=num u$
du:=den u$
symbolic for each a in coords!* do
<<step:=get(a,'gridstep)$
flag(list step,'used!*)$
put(step,'avalue,'(scalar 0)) >>$
symbolic rmsubs()$
nu:=nu$
du:=du$
symbolic for each a in coords!* do
<<step:=get(a,'gridstep)$
remflag(list step,'used!*)$
remprop(step,'avalue) >>$
symbolic rmsubs()$
if du=0 then <<write
" Reformulate difference scheme, grid steps remain in denominators"$
u:=0 >>
else u:=nu/du$
return u
end$
procedure minorder(pol,var)$
begin
scalar lcofs,mord$
coeff(den pol,var)$
mord:=-hipow!*$
lcofs := rest coeff(num pol,var)$
if not(mord=0) then return (mord+lowpow!*)$
mord:=1$
a:if lcofs={} then return 0
else if first lcofs=0 then lcofs:=rest lcofs
else return mord$
mord:=mord+1$
go to a
end$
endmodule;
module charpol;
% Author: R. Liska.
% Version REDUCE 3.6 05/1991.
fluid '(!*exp !*gcd !*prfourmat)$
switch prfourmat$
!*prfourmat:=t$
procedure coefc1 uu$
begin
scalar lco,l,u,v,a$
u:=car uu$
v:=cadr uu$
a:=caddr uu$
lco:=aeval list('coeff,u,v)$
lco:=cdr lco$
l:=length lco - 1$
for i:=0:l do
<<setel(list(a,i),car lco)$
lco:=cdr lco >>$
return (l . 1)
end$
deflist('((coefc1 coefc1)),'simpfn)$
global '(cursym!* coords!* icoords!* unvars!*)$
icoords!*:='(i j k l m n i1 j1 k1 l1 m1 n1)$
flag('(tcon unit charmat ampmat denotemat),'matflg)$
put('unit,'rtypefn,'getrtypecar)$
put('charmat,'rtypefn,'getrtypecar)$
put('ampmat,'rtypefn,'getrtypecar)$
put('denotemat,'rtypefn,'getrtypecar)$
procedure unit u$
generateident length matsm u$
procedure charmat u$
matsm list('difference,list('times,'lam,list('unit,u)),u)$
procedure charpol u$
begin
scalar x,complexx;
complexx:=!*complex;
algebraic on complex;
x:=simp list('det,list('charmat,carx(u,'charpol)))$
if null complexx then algebraic off complex;
return x
end;
put('charpol,'simpfn,'charpol)$
algebraic$
korder lam$
procedure re(x)$
sub(i=0,x)$
procedure im(x)$
(x-re(x))/i$
procedure con(x)$
sub(i=-i,x)$
procedure complexpol x$
begin
scalar y$
y:=troot1 x$
return if im y=0 then y
else y*con y
end$
procedure troot1 x$
begin
scalar y$
y:=x$
while not(sub(lam=0,y)=y) and sub(lam=1,y)=0 do y:=y/(lam-1)$
x:=sub(lam=1,y)$
if not(numberp y or (numberp num y and numberp den y)) then
write " If ",re x," = 0 and ",im x,
" = 0 , a root of the polynomial is equal to 1."$
return y
end$
procedure hurw(x)$
% X is a polynomial in LAM, all its roots are |LAMI|<1 <=> for all roots
% of the polynomial HURW(X) holds RE(LAMI)<0.
(lam-1)**deg(num x,lam)*sub(lam=(lam+1)/(lam-1),x)$
symbolic$
procedure unfunc u$
<<unvars!*:=u$
for each a in u do put(a,'simpfn,'simpiden) >>$
put('unfunc,'stat,'rlis)$
global '(denotation!* denotid!*)$
denotation!*:=nil$
denotid!*:='a$
procedure denotid u$
<<denotid!*:=car u$nil>>$
put('denotid,'stat,'rlis)$
procedure cleardenot$
denotation!*:=nil$
put('cleardenot,'stat,'endstat)$
flag('(cleardenot),'eval)$
algebraic$
array cofpol!*(20)$
procedure denotepol u$
begin
scalar nco,dco$
dco:=den u$
u:=num u$
nco:=coefc1 (u,lam,cofpol!*)$
for j:=0:nco do cofpol!*(j):=cofpol!*(j)/dco$
denotear nco$
u:=for j:=0:nco sum lam**j*cofpol!*(j)$
return u
end$
symbolic$
put('denotear,'simpfn,'denotear)$
procedure denotear u$
begin
scalar nco,x$
nco:=car u$
for i:=0:nco do
<<x:=list('cofpol!*,i)$
setel(x,mk!*sq denote(getel x,0,i)) >>$
return (nil .1)
end$
procedure denotemat u$
begin
scalar i,j,x$
i:=0$
x:=for each a in matsm u collect
<<i:=i+1$
j:=0$
for each b in a collect
<<j:=j+1$
denote(mk!*sq b,i,j) >> >>$
return x
end$
procedure denote(u,i,j)$
% U is prefix form, I,J are integers
begin
scalar reu,imu,ireu,iimu,eij,fgcd$
if atom u then return simp u$
fgcd:=!*gcd$
!*gcd:=t$
reu:=simp!* list('re,u)$
imu:=simp!* list('im,u)$
!*gcd:=fgcd$
eij:=append(explode i,explode j)$
ireu:=intern compress append(append(explode denotid!* ,'(r)),eij)$
iimu:=intern compress append(append(explode denotid!* ,'(i)),eij)$
if car reu then insdenot(ireu,reu)$
if car imu then insdenot(iimu,imu)$
return simp list('plus,
if car reu then ireu
else 0,
list('times,
'i,
if car imu then iimu
else 0))
end$
procedure insdenot(iden,u)$
denotation!*:=(u . iden) . denotation!*$
procedure prdenot$
for each a in reverse denotation!* do
assgnpri(list('!*sq,car a,t),list cdr a,'only)$
put('prdenot,'stat,'endstat)$
flag('(prdenot),'eval)$
procedure ampmat u$
begin
scalar x,i,h1,h0,un,rh1,rh0,ru,ph1,ph0,!*exp,!*gcd,complexx$
complexx:=!*complex;
!*exp:=t$
fouriersubs()$
u:=car matsm u$
x:=for each a in coords!* collect if a='t then 0
else list('times,
tcar get(a,'index),
get(a,'wave),
get(a,'step))$
x:=list('expp,'plus . x)$
x:=simp x$
u:=for each a in u collect resimp quotsq(a,x)$
gonsubs()$
algebraic on complex;
u:=for each a in u collect resimp a$
remfourier()$
a:if null u then go to d$
ru:=caar u$
un:=unvars!*$
i:=1$
b:if un then go to c$
rh1:=reverse rh1$
rh0:=reverse rh0$
h1:=rh1 . h1$
h0:=rh0 . h0$
rh0:=rh1:=nil$
u:=cdr u$
go to a$
c:rh1:=coefck(ru,list('u1!*,i)) . rh1$
rh0:=negsq coefck(ru,list('u0!*,i)) . rh0$
un:=cdr un$
i:=i+1$
go to b$
d:h1:=reverse h1$
h0:=reverse h0$
if !*prfourmat then
<<ph1:=for each a in h1 collect
for each b in a collect prepsq!* b$
setmatpri('h1,nil . ph1)$
ph1:=nil$
ph0:=for each a in h0 collect
for each b in a collect prepsq!* b$
setmatpri('h0,nil . ph0)$
ph0:=nil >>$
!*gcd:=t;
x:=if length h1=1 then list list quotsq(caar h0,caar h1)
else lnrsolve(h1,h0)$
if null complexx then algebraic off complex;
return x
end$
procedure coefck(x,y)$
% X is standard form, Y is prefix form, returns coefficient of Y
% appearing in X, i.e. X contains COEFCK(X,Y)*Y
begin
scalar ky,xs$
ky:=!*a2k y$
xs:=car subf(x,list(ky . 0))$
xs:=addf(x,negf xs)$
if null xs then return(nil . 1)$
xs:=quotf1(xs,!*k2f ky)$
return if null xs then <<msgpri
("COEFCK failed on ",list('sq!*,x . 1,nil)," with ",y,'hold)$
xread nil$
!*f2q xs>>
else !*f2q xs
end$
procedure simpfour u$
begin
scalar nrunv,x,ex,arg,mv,cor,incr,lcor$
nrunv:=get(car u,'nrunvar)$
a:u:=cdr u$
if null u then go to r$
arg:=simp car u$
mv:=mvar car arg$
if not atom mv or not numberp cdr arg then return msgpri
("Bad index ",car u,nil,nil,'hold)$
cor:=tcar get(mv,'coord)$
if not(cor member coords!*) then return msgpri
("Term ",car u," contains non-coordinate ",mv,'hold)$
if cor member lcor then return msgpri
("Term ",car u," means second appearance of coordinate ",cor,
'hold)$
if not(cor='t) and cdr arg>get(cor,'maxden)
then put(cor,'maxden,cdr arg)$
lcor:=cor . lcor$
incr:=addsq(arg,negsq !*k2q mv)$
if not flagp(cor,'uniform) then return lprie
("Non-uniform grids not yet supported")$
if cor='t then go to ti$
ex:=list('times,car u,get(cor,'step),get(cor,'wave)) . ex$
go to a$
ti:if null car incr then x:=list('u0!*,nrunv)
else if incr= 1 . 1 then x:=list('u1!*,nrunv)
else return lprie "Scheme is not twostep in time"$
go to a$
r:for each a in setdiff(coords!*,lcor) do
if a='t then x:=list('u0!*,nrunv)
else ex:=list('times,tcar get(a,'index),get(a,'step),get(a,'wave))
. ex$
return simp list('times,x,list('expp,'plus . ex))
end$
procedure fouriersubs$
begin
scalar x,i$
for each a in '(expp u1!* u0!*) do put(a,'simpfn,'simpiden)$
x:=unvars!*$
i:=1$
a:if null x then go to b$
put(car x,'nrunvar,i)$
i:=i+1$
x:=cdr x$
go to a$
b:flag(unvars!*,'full)$
for each a in unvars!* do put(a,'simpfn,'simpfour)$
for each a in coords!* do
if not(a='t) then
<<x:=intern compress append(explode 'h,explode a)$
put(a,'step,x)$
if not flagp(a,'uniform) then put(x,'simpfn,'simpiden)$
x:=intern compress append(explode 'k,explode a)$
put(a,'wave,x)$
x:=intern compress append(explode 'a,explode a)$
put(a,'angle,x)$
put(a,'maxden,1) >>$
algebraic for all z,y,v let
expp((z+y)/v)=expp(z/v)*expp(y/v),
expp(z+y)=expp z*expp y
end$
procedure gonsubs$
begin
scalar xx$
algebraic for all z,y,v clear expp((z+y)/v),expp(z+y)$
for each a in coords!* do
if not(a='t) then
<<xx:=list('quotient,
list('times,
get(a,'maxden),
get(a,'angle)),
get(a,'step))$
setk(get(a,'wave),aeval xx)$
xx:=list('times,
get(a,'wave),
get(a,'step))$
mathprint list('setq,
get(a,'angle),
if get(a,'maxden)=1 then xx
else list('quotient,
xx,
get(a,'maxden))) >>$
algebraic for all x let expp x=cos x+i*sin x$
algebraic for all x,n such that numberp n and n>1 let
sin(n*x)=sin x*cos((n-1)*x)+cos x*sin((n-1)*x),
cos(n*x)=cos x*cos((n-1)*x)-sin x*sin((n-1)*x)$
for each a in unvars!* do
<<put(a,'simpfn,'simpiden)$
remprop(a,'nrunvar) >>
end$
procedure remfourier$
<<algebraic for all x clear expp x$
algebraic for all x,n such that numberp n and n>1 clear
sin(n*x),cos(n*x)$
for each a in coords!* do
if not(a='t) then
<<remprop(a,'step)$
remprop(remprop(a,'wave),'avalue)$
remprop(a,'maxden) >> >>$
operator numberp$
endmodule;
module hurwp;
% Author: R. Liska.
% Version REDUCE 3.6 05/1991
global '(ofl!* mlist!*)$
fluid '(!*exp !*gcd)$
flag('(tcon),'matflg)$
put('tcon,'msimpfn,'tcon)$
put('tcon,'rtypefn,'getrtypecar)$
procedure tcon u$
% Calculates complex conjugate and transpose matrix
begin
scalar v,b$
v:=matsm list('tp,u)$
for each a in v do
<<b:=a$
while b do
<<rplaca(b,quotsq(subf(numr car b, '((i minus i))),
!*f2q denr car b))$
b:=cdr b >> >>$
return v
end$
algebraic$
korder lam$
symbolic$
global '(positive!* userpos!* userneg!* !*pfactor)$
!*pfactor:=nil$
procedure positivep u$
% U is prefix form. Procedure tests if U>0, eventually writes this
% condition and puts U into POSITIVE!*. If U<=0 then returns NIL,
% if U>0 then T, in other cases 'COND.
% If it does not know if U>0 and program is running in interactive
% mode it asks user if U>0 and return value is based on user reply.
if numberp u then
if u>0 then t
else nil
else if eqcar(u,'!*sq) and fixp caadr u and fixp cdadr u then
if caadr u*cdadr u>0 then t
else nil
else
begin
scalar x,exp$
exp:=!*exp$
if !*pfactor and
member('factor,mlist!*) then
<<!*exp:=nil$
u:=aeval list('ppfactor,u) >>$
u:=prepsq!* simp u$
!*exp:=exp$
x:=if terminalp() and null ofl!* then
begin
scalar y,z$
prin2!* "Is it true, that "$
maprin u$
prin2!* " > 0 ?"$
a:prin2!* " Reply (Y/N/?)"$
terpri!* t$
y:=read()$
if y eq 'y then <<z:=t$ userpos!*:=u . userpos!* >>
else if y eq 'n
then <<z:=nil$ userneg!*:=u . userneg!*>>
else if y eq '? then z:='cond
else go to a$
return z
end
else 'cond$
if x eq 'cond then
if null positive!* then positive!*:= list (1 . u)
else positive!* := ((caar positive!* + 1) . u) . positive!*$
return x
end$
global'(hconds!*)$
algebraic$
array cof(20),fcof(20)$
share hconds!*$
procedure ppfactor x$
begin
scalar d,n,de$
d:=factorize(num x)$
n:=for each a in d product a$
if den x=1 then return n$
d:=factorize(den x)$
de:=for each a in d product a$
return (n/de)
end$
procedure hurwitzp u$
% U is a polynomial in LAM. Procedure tests if it is Hurwitz polynomial
% i.e. for all its rools LAMI holds RE(LAMI)<0.
% Returned values: YES - definitely yes
% NO - definitely no
% COND - if conditions holds (all members of POSITIVE!*
% are >0)
if im u=0 then rehurwp u
else cohurwp u$
symbolic$
procedure coef1(u,v,a)$
begin
scalar lco,l$
lco:=aeval list('coeff,u,v)$
lco:=cdr lco$
l:=length lco - 1$
for i:=0:l do
<<setel(list(a,i),car lco)$
lco:=cdr lco >>$
return l
end$
procedure rehurwp u$
begin
scalar deg,hurp,gcd$
gcd:=!*gcd$
!*gcd:=t$
deg:=coef1(car u,'lam,'cof)$
if deg=0 then return typerr(u,"Polynomial in LAM")$
positive!* := userpos!* := userneg!* := nil$
if deg <= 2 then
<<for i:=0:deg do setel(list('cof,i),
aeval list('quotient,
getel list('cof,i),
getel list('cof,deg)))$
if deg=1 then hurp:=positivep getel list('cof,0)
else if deg=2 then hurp:=
if positivep getel list('cof,0) and
positivep getel list('cof,1) then
if positive!* then 'cond
else t
else if positive!* then 'cond % added 08/08/91
else nil$
printcond(nil) >>
else hurp:=rehurwp1 deg$
!*gcd:=gcd$
return rethurp hurp
end$
procedure rethurp hurp$
<<hconds!*:= 'list . if positive!* then
for each a in positive!* collect cdr a
else nil$
!*k2q(if null hurp then 'no
else if null positive!* then 'yes
else 'cond) >>$
put('rehurwp,'simpfn,'rehurwp)$
procedure cohurwp u$
begin
scalar deg$
u:=reval list('sub,'(equal lam (times i lam)),car u)$
deg:=coef1(u,'lam,'cof)$
if deg=0 then return typerr(u,"Polynomial in LAM")$
positive!* := userpos!* := userneg!* :=nil$
if aeval list('im,getel list('cof,deg))=0 then
for j:= 0:deg do setel(list('cof,j),
aeval list('times,'i,getel list('cof,j)))$
return rethurp cohurwp1 (deg)
end$
put('cohurwp,'simpfn,'cohurwp)$
procedure rehurwp1 deg$
begin
scalar i,bai,bdi,x,lich,sud,bsud,matr,hmat,csud,clich,dsud,dlich$
a:i:=deg$
csud:=clich:=nil$
bsud:=t$
b:x:=positivep getel list('cof,i)$
if null x then go to c
else if x eq t then bai:=t
else if x eq 'cond then
if i=deg or i=0 then <<csud:=caar positive!* . csud$
clich:=caar positive!* . clich >>
else if bsud then csud:=caar positive!* . csud
else clich:=caar positive!* . clich$
i:=i-1$
bsud:=not bsud$
if i>=0 then go to b$
go to d$
% Change of sign AI = - AI
c:if bai or bdi then go to n
else bai:=t$
for i:=0:deg do setel(list('cof,i),
aeval list('minus,getel list('cof,i)))$
go to a$
% Checking DI > 0 - Hurwitz determinants
% Splitting to odd and even coeffs. AI, A0 is coeff. by LAM**DEG
d:bsud:=t$
for i:=deg step -1 until 0 do
<<x:=simp getel list('cof,i)$
if bsud then sud:=x . sud
else lich:=x . lich$
bsud:=not bsud >>$
sud:=reverse sud$
lich:=reverse lich$
% Filling of SUD and LICH on the length DEG by zeroes from right
sud:=filzero(sud,deg)$
lich:=filzero(lich,deg)$
dsud:=dlich:=nil$
matr:=nil$
i:=1$
bsud:=nil$
d1:matr:=nconc(matr,list lich)$
lich:=(nil . 1) . lich$
d2:hmat:=cutmat(matr,i)$
x:=mk!*sq detq hmat$
x:=positivep x$ % Necessary to add storing of odd and even DIs
if null x then
if bsud then go to n
else go to c
else if x eq t and not bsud then bdi:=t
else if x eq 'cond then
if bsud then dsud:=caar positive!* . dsud
else dlich:=caar positive!* . dlich$
i:=i+1$
bsud:=not bsud$
if i>deg then go to k$
if not bsud then go to d1$
matr:=nconc(matr,list sud)$
sud:=(nil . 1) . sud$
go to d2$
n:return nil$
k:if null positive!* or ((null csud or null clich) and
(null dsud or null dlich))
then return <<printuser()$ t>>$
prin2t "If we denote:"$
printcond(t)$
printdef('c1,clich:=reverse clich)$
printdef('c2,csud:=reverse csud)$
printdef('d1,dlich:=reverse dlich)$
printdef('d2,dsud:=reverse dsud)$
prin2t "Necessary and sufficient conditions are:"$
prin2t if null csud or null clich then " (D1) OR (D2)"
else if null dsud or null dlich then " (C1) OR (C2)"
else " ( (C1) OR (C2) ) AND ( (D1) OR (D2) )"$
printuser()$
return 'cond
end$
procedure printcond(x)$
<<if not x then
prin2t "Necessary and sufficient conditions are:"$
positive!*:=reverse positive!*$
for each a in positive!* do
<<if x then <<prin2!* " ("$
prin2!* car a$
prin2!* ") " >>$
maprin cdr a$
prin2!* " > 0"$
terpri!* t >>$
if not x then printuser() >>$
procedure printuser()$
if userpos!* or userneg!* then
<<prin2t"You have explicitly stated:"$
for each a in userpos!* do <<maprin a$
prin2!* " > 0"$
terpri!* t >>$
for each a in userneg!* do <<maprin a$
prin2!* " <= 0"$
terpri!* t >> >>$
procedure printdef(x,y)$
if y then
<<prin2!* " ("$
prin2!* x$
prin2!* ") ("$
prin2!* car y$
prin2!* ")"$
if cdr y then for each a in cdr y do
<<prin2!* " AND ("$
prin2!* a$
prin2!* ")" >>$
terpri!* t >>$
procedure filzero(x,n)$
% Adds zeros (in S.Q. form) to the list X from right on the length N
begin
scalar y,i$
y:=x$
i:=1$
if null x then return typerr(x,"Empty list")$
while cdr y do
<<y:=cdr y$
i:=i+1>>$
while i<n do
<<rplacd(y,list(nil . 1))$
y:=cdr y$
i:=i+1 >>$
return x
end$
procedure cutmat(x,n)$
% From each member of list X, i.e. row of a matrix, remains
% the first N elements
for each a in x collect cutrow(a,n)$
procedure cutrow(y,n)$
begin
scalar i,z,zz$
i:=1$
z:=list car y$
zz:=z$
y:=cdr y$
while i<n do
<<rplacd(zz,list car y)$
y:=cdr y$
zz:=cdr zz$
i:=i+1 >>$
return z
end$
procedure cohurwp1 (deg)$
begin
scalar k,x,y,ak,bk,akk,bkk,matr,hmat$
% Splitting on RE and IM part
for j:=0:deg do
<<x:=getel list('cof,j)$
y:=simp list('re,x)$
x:=resimp simp list('quotient,list('difference,x,mk!*sq y),'i)$
ak:=x . ak$
bk:=y . bk >>$
% Construction of coeffs. AI, BI
positive!*:=userpos!*:=userneg!*:=nil$
akk:=filzero(ak,2*deg)$
bkk:=filzero(bk,2*deg)$
k:=2$
d1:matr:=nconc(matr,list akk)$
matr:=nconc(matr,list bkk)$
akk:=(nil . 1) . akk$
bkk:=(nil . 1) . bkk$
hmat:=cutmat(matr,k)$
x:=mk!*sq detq hmat$
x:=positivep x$
if null x then go to n$
if k=2*deg then go to ko$
k:=k+2$
go to d1$
n:return nil$
ko:printcond(nil)$
return t
end$
endmodule;
module linband;
% Author: R. Liska
% Version REDUCE 3.6 05/1991
% GENTRAN package has to be loaded prior to this module
global'(fortcurrind!* genstmtnum!* genstmtincr!*)$
fluid'(!*period)$ % declaration for 3.4
fluid'(!*imsl !*nag !*essl)$
switch imsl,nag,essl$
!*imsl:=nil$
!*nag:=nil$
!*essl:=nil$
procedure ison x$
if eval x then 1
else 0$
operator ison$
if null getd 'gentempst then
procedure gentempst$
list('gentemp,xread t)$
global'(temp!*)$
temp!*:=nil$
procedure gentemp u$
<<temp!* := ((!*period . fortcurrind!*) . u) . temp!*$ nil>>$
put('gentemp,'stat,'gentempst)$
put('gentemp,'formfn,'formgentran)$
load!-package 'gentran;
procedure outtemp$
begin
scalar period,fortind$
period:=!*period$
fortind:=fortcurrind!*$
for each a in reverse temp!* do
<<!*period:=caar a$
fortcurrind!*:=cdar a$
eval list('gentran,mkquote cdr a,nil)>>$
temp!* := nil$
!*period:=period$
fortcurrind!*:=fortind$
return nil
end$
put('outtemp,'stat,'endstat)$
flag('(outtemp),'eval)$
algebraic$
procedure genlinbandsol(nlc,nuc,system)$
% Generates FORTRAN program for solving of linear algebraic system
% of equations with band matrix with NLC lower codiagonals and NUC
% upper codiagonals.
begin
scalar pvars,svars,vareq,fveq$
% PVARS - list of variables before actual variable
% SVARS - list of variables after actual variable
% VAREQ - actual v-equation (list {variable equation})
symbolic
<<put('list,'evfno,get('list,'evfn))$
put('list,'evfn,'listnoeval)$
put('equal,'psopfno,get('equal,'psopfn))$
put('equal,'psopfn,'equalaeval)>>$
system:=expanddo(nlc,nuc,system)$
vareq:=first system$
pvars:={}$
svars:=findsvars(nuc,vareq,system)$
off period$
gentran n:=1$
gentemp n:=1$
on period$
ncol!*:=nlc+nuc+1$
for i:=1:nlc do
<<genvareq(pvars,svars,vareq,nlc,nlc-i+1,pfix!*)$
pvars:=append(pvars,first vareq . {})$
system:=nextvareqsys(vareq,system)$
vareq:=first system$
system:=rest system$
gennp1()$
svars:=findsvars(nuc,vareq,system) >>$
while length svars=nuc do
<<genvareq(pvars,svars,vareq,nlc,0,0)$
pvars:=append(rest pvars,first vareq . {})$
fveq:=first system$
system:=nextvareqsys(vareq,system)$
vareq:=first system$
system:=rest system$
% Get in and get out of loop
if (ffst system=do) and (first vareq=first frrfst system and
rest vareq=rest frrfst system) then
pvars:=findpvars(nlc,first system)
else if first fveq=do and not(ffst system=do) then
pvars:=lastvars(nlc,fveq)$
gennp1()$
svars:=findsvars(nuc,vareq,system) >>$
for i:=1:nuc do
<<genvareq(pvars,svars,vareq,nlc,i,sfix!*)$
pvars:=append(rest pvars,first vareq . {})$
system:=nextvareqsys(vareq,system)$
vareq:=first system$
system:=rest system$
if not(svars={}) then
<<svars:=rest svars$
gennp1() >> >>$
off period$
if ison !*imsl = 1 then pvars:=gencall!-imsl(nlc,nuc)
else if ison !*nag = 1 then pvars:=gencall!-nag(nlc,nuc)
else if ison !*essl= 1 then pvars:=gencall!-essl(nlc,nuc)
else pvars:=gencall!-linpack(nlc,nuc)$
on period$
outtemp$
symbolic <<put('list,'evfn,remprop('list,'evfno))$
put('equal,'psopfn,remprop('equal,'psopfno))>>
end$
procedure gencall!-imsl (nlc,nuc)$
gentran
<<literal tab!*,"call leqt1b(acof,n,",eval nlc,",",eval nuc,
",iacof,arhs,1,iarhs,0,xl,ier)",cr!*$
literal "c iacof is actual 1-st dimension of the acof array",cr!*$
literal "c iarhs is actual 1-st dimension of the arhs array",cr!*$
literal "c xl is working array with size n*(nlc+1)",cr!*$
literal
"c where n is number of equations nlc number of lower",cr!*$
literal "c codiagonals",cr!*$
literal
"c if ier=129( .ne.0) matrix acof is algorithmically singular",
cr!*$
literal tab!*,"if(ier.ne.0) call errout",cr!*>>$
procedure gencall!-linpack(nlc,nuc)$
if ncol!*=3 and nlc=1 then gencall!-linpack!-trid(nlc,nuc)
else gentran
<<literal tab!*,"call dgbfa(acof,iacof,n,",eval nlc,",",eval nuc,
",ipvt,ier)",cr!*$
literal "c n is number of equations",cr!*$
literal "c acof is array of dimension (iacof,p), p >= n",cr!*$
literal "c iacof >= ",eval(nlc+ncol!*),cr!*$
literal "c ipvt is array of dimension at least (n)",cr!*$
literal "c if (ier.ne.0) matrix acof is algorithmically singular",
cr!*$
literal tab!*,"if(ier.ne.0) call errout",cr!*$
literal tab!*,"call dgbsl(acof,iacof,n,",eval nlc,",",eval nuc,
",ipvt,arhs,0)",cr!*>>$
procedure gencall!-linpack!-trid(nlc,nuc)$
gentran
<<literal tab!*,"call dgtsl(n,alcd,ad,aucd,arhs,ier)",cr!*$
literal "c n is number of equations",cr!*$
literal
"c alcd,ad,aucd,arhs are arrays of dimension at least (n)",cr!*$
literal "c if (ier.ne.0) matrix acof is algorithmically singular",
cr!*$
literal tab!*,"if(ier.ne.0) call errout",cr!* >>$
procedure gencall!-essl(nlc,nuc)$
if ncol!*=3 and nlc=1 then gencall!-essl!-trid(nlc,nuc)
else gentran
<<literal tab!*,"call dgbf(acof,iacof,n,",eval nlc,",",eval nuc,
",ipvt)",cr!*$
literal "c n is number of equations",cr!*$
literal "c acof and arhs are double precision type",cr!*$
literal "c for single precision change dgbf to sgbf and ",
"dgbs to sgbs",cr!*$
literal "c acof is array of dimension (iacof,p), p >= n",cr!*$
literal "c iacof >= ",eval(nlc+ncol!*+15),cr!*$
literal "c arhs is array of dimension at least (n)",cr!*$
literal "c ipvt is integer array of dimension at least (n)",cr!*$
literal tab!*,"call dgbs(acof,iacof,n,",eval nlc,",",eval nuc,
",ipvt,arhs)",cr!*>>$
procedure gencall!-essl!-trid(nlc,nuc)$
gentran
<<literal tab!*,"call dgtf(n,alcd,ad,aucd,af,ipvt)",cr!*$
literal "c n is number of equations",cr!*$
literal
"c alcd,ad,aucd,af,arhs are arrays of dimension at least (n)",cr!*$
literal "c these arrays are double precision type",cr!*$
literal "c for single precision change dgtf to sgtf and ",
"dgts to sgts",cr!*$
literal
"c ipvt is integer array of dimension at least (n+3)/8",cr!*$
literal tab!*,"call dgts(n,alcd,ad,aucd,af,ipvt,arhs)",cr!* >>$
procedure gencall!-nag(nlc,nuc)$
if ncol!*=3 and nlc=1 then gencall!-nag!-trid(nlc,nuc)
else gentran
<<ier:=0$
literal tab!*,"call f01lbf(n,",eval nlc,",",eval nuc,
",acof,iacof,al,ial,in,iv,ier)",cr!*$
literal "c n is number of equations",cr!*$
literal "c acof is array of dimension (iacof,p), p >= n",cr!*$
literal "c iacof >= min(n,",eval ncol!*,")",cr!*$
literal "c al is array of dimension (ial,p), p >= n",cr!*$
literal "c ial >= max(1,",eval nlc,")",cr!*$
literal "c in is integer array of dimension at least (n)",cr!*$
literal tab!*,"if(ier.ne.0) call errout",cr!*$
literal tab!*,"call f04ldf(n,",eval nlc,",",eval nuc,
",1,acof,iacof,al,ial,in,arhs,iarhs,ier)",cr!*$
literal "c arhs is array of dimension (iarhs), iarhs >= n",cr!*$
literal tab!*,"if(ier.ne.0) call errout",cr!* >>$
procedure gencall!-nag!-trid(nlc,nuc)$
gentran
<<ier:=0$
literal tab!*,
"call f01lef(n,ad,0.,aucd,alcd,1.e-10,au2cd,in,ier)",cr!*$
literal "c n is number of equations",cr!*$
literal
"c alcd,ad,aucd,au2cd,arhs are arrays of dimension at least (n)",cr!*$
literal "c in is integer array of dimension at least (n)",cr!*$
literal tab!*,"if(ier.ne.0 .or. in(n).ne.0) call errout",cr!*$
literal tab!*,
"call f04lef(1,n,ad,aucd,alcd,au2cd,in,arhs,0.,ier)",cr!*$
literal tab!*,"if(ier.ne.0) call errout",cr!* >>$
procedure gennp1$
<<off period$
gentran n:=n+1$
gentemp n:=n+1$
on period >>$
% Definition of operator SUBE
symbolic$
symbolic procedure simpsube u$
begin
scalar x$
a:if null cdr u then go to d
else if null eqexpr car u then errpri2(car u,t)$
x:=list('equal,reval cadar u,caddar u) . x$
u:=cdr u$
go to a$
d:x:=reverse(car u . x)$
x:=subeval x$
return x
end$
symbolic put('sube,'psopfn,'simpsube)$
algebraic$
% Procedures FFRRST etc.
procedure ffst u$
first first u$
procedure frst u$
first rest u$
procedure rfst u$
rest first u$
procedure rrst u$
rest rest u$
procedure fffst u$
first ffst u$
procedure ffrst u$
first frst u$
procedure frfst u$
first rfst u$
procedure frrst u$
first rrst u$
procedure rffst u$
rest ffst u$
procedure rfrst u$
rest frst u$
procedure rrfst u$
rest rfst u$
procedure rrrst u$
rest rrst u$
procedure ffffst u$
ffst ffst u$
procedure fffrst u$
ffst frst u$
procedure ffrfst u$
ffst rfst u$
procedure ffrrst u$
ffst rrst u$
procedure frffst u$
frst ffst u$
procedure frfrst u$
frst frst u$
procedure frrfst u$
frst rfst u$
procedure frrrst u$
frst rrst u$
procedure rfffst u$
rfst ffst u$
procedure rffrst u$
rfst frst u$
procedure rfrfst u$
rfst rfst u$
procedure rfrrst u$
rfst rrst u$
procedure rrffst u$
rrst ffst u$
procedure rrfrst u$
rrst frst u$
procedure rrrfst u$
rrst rfst u$
procedure rrrrst u$
rrst rrst u$
procedure findsvars(nuc,vareq,system)$
% Looks for NUC next variables in SYSTEM
% VAREQ is actual v-equation
if ffst system=do then findsvarsdo(nuc,vareq,first system)
else findsvars1(nuc,rest system)$
procedure findsvars1(nuc,system)$
% Substitutes values for loop variable
if nuc=0 or system={} then {}
else if ffst system=do then fsvars1do(nuc,first system)
else ffst system . findsvars1(nuc-1,rest system)$
procedure fsvars1do(nuc,cykl)$
% Substitutes into the loop CYKL
begin
scalar id,from,step,syst,x,y$
cykl:=rest cykl$
syst:=first cykl$
id:=first syst$
from:=frst syst$
step:=frrrst syst$
syst:=rest cykl$
x:={}$
a:y:=sube(id=from,ffst syst)$
x:=y . x$
nuc:=nuc-1$
if nuc=0 then go to r$
syst:=rest syst$
if not(syst={}) then go to a$
syst:=rest cykl$
from:=from+step$
go to a$
r:x:=reverse x$
return x
end$
procedure findsvarsdo(nuc,vareq,cykl)$
% Does not substitute for loop variable, only increases it
% by STEP if it is necessary
begin
scalar id,add1,step,syst,x,y$
cykl:=rest cykl$
syst:=first cykl$
id:=first syst$
step:=frrrst syst$
syst:=rest cykl$
while not(first vareq=ffst syst and rest vareq=rfst syst)
do syst:=rest syst$
syst:=rest syst$
add1:=0$
x:={}$
a:if syst={} then go to b$
y:=sube(id=id+add1,ffst syst)$
x:=y . x$
nuc:=nuc-1$
if nuc=0 then go to r$
syst:=rest syst$
go to a$
b:syst:=rest cykl$
add1:=add1+step$
go to a$
r:x:=reverse x$
return x
end$
procedure expanddo(nlc,nuc,system)$
% Every loop in SYSTEM is expanded so that more than or equal to
% NLC first elements and more than or equal NUC last elements are
% excluded from the loop, and changes the parameters of loop so
% that its meaning remains the same
begin
scalar x$
x:={}$
a:if system={} then go to r$
if ffst system=do then x:=append(expddo(nlc,nuc,first system),x)
else x:=first system . x$
system:=rest system$
go to a$
r:x:=reverse x$
return x
end$
procedure expddo(nlc,nuc,cykl)$
% Performs the expansion of the loop CYKL - returns reverse list
begin
scalar id,from,to1,step,syst,lsyst,ns,x,y,bn$
cykl:=rest cykl$
syst:=first cykl$
id:=first syst$
from:=frst syst$
to1:=frrst syst$
step:=frrrst syst$
syst:=rest cykl$
lsyst:=length syst$
ns:=quotient1(nlc,lsyst)$
if nlc>ns*lsyst then ns:=ns+1$
bn:=0$
x:={}$
a:y:=sube(id=from,ffst syst) . sube(id=from,frfst syst) . {}$
x:=y . x$
syst:=rest syst$
if not(syst={}) then go to a$
ns:=ns-1$
from:=from+step$
if ns=0 then go to b$
syst:=rest cykl$
go to a$
b:if bn=1 then go to r$
syst:=rest cykl$
ns:=quotient1(nuc,lsyst)$
if nuc>ns*lsyst then ns:=ns+1$
to1:=to1-ns*step$
y:=do . (id . from . to1 . step . {}) . syst$
x:=y . x$
from:=to1+step$
bn:=1$
go to a$
r:return x
end$
symbolic procedure quotient1(u,v)$
quotient(u,v)$
symbolic operator quotient1$
operator acof,arhs$
procedure genvareq(pvars,svars,vareq,nlc,nzero,mode)$
if ison !*imsl = 1 then
genvareq!-imsl(pvars,svars,vareq,nlc,nzero,mode)
else if ison !*nag = 1 then
genvareq!-nag(pvars,svars,vareq,nlc,nzero,mode)
else genvareq!-linpack(pvars,svars,vareq,nlc,nzero,mode)$
procedure genvareq!-imsl(pvars,svars,vareq,nlc,nzero,mode)$
% Generates N-th row of coeff. matrix ACOF and right hand side ARHS
% according to the v-equation VAREQ.
% NZERO is number of zeroes before or after (according to MODE).
% Matrix ACOF is transformed to IMSL band matrix storage.
begin
integer j$
scalar var,rhside,lhside,x,y$
if not(length pvars + length svars+1+nzero=ncol!*) then return
write" Unconsistent PVARS:",pvars," SVARS:",svars," NZERO:",nzero$
var:=first vareq$
vareq:=frst vareq$
rhside:=rhs vareq$
lhside:=lhs vareq$
j:=1$
x:=0$
if mode=pfix!* then while j<=nzero do
<<gentran acof(n,eval j):=0$
j:=j+1 >>$
for each a in pvars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran acof(n,eval j):=:y$
j:=j+1>>$
y:=lincof(lhside,var)$
x:=x+var*y$
gentran acof(n,eval j):=:y$
j:=j+1$
for each a in svars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran acof(n,eval j):=:y$
j:=j+1>>$
if mode=sfix!* then while j<=ncol!* do
<<gentran acof(n,eval j):=0$
j:=j+1 >>$
gentran arhs(n):=:rhside$
gentemp eval(var):=arhs(n)$
if not(x-lhside=0) then write " For equation ",vareq," given only ",
"variables ",pvars,svars,var$
return
end$
procedure genvareq!-linpack(pvars,svars,vareq,nlc,nzero,mode)$
% Generates N-th row of coeff. matrix ACOF and right hand side ARHS
% according to the v-equation VAREQ.
% NZERO is number of zeroes before or after (according to MODE).
% Matrix ACOF is transformed to LINPACK band matrix storage.
% NCOL!* is the band width.
begin
integer j,jj,nn$
scalar var,rhside,lhside,x,y$
if not(length pvars + length svars+1+nzero=ncol!*) then return
write" Unconsistent PVARS:",pvars," SVARS:",svars," NZERO:",nzero$
if nlc=1 and ncol!*=3 then return
genvareq!-linpack!-trid(pvars,svars,vareq,nlc,nzero,mode)$
var:=first vareq$
vareq:=frst vareq$
rhside:=rhs vareq$
lhside:=lhs vareq$
j:=n-nlc$
jj:=1$
nn:=ncol!*+nlc$
x:=0$
if mode=pfix!* then while jj<=nzero do
<<nn:=nn-1$
jj:=jj+1$
j:=j+1 >>$
for each a in pvars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran acof(nn,j)::=:y$
nn:=nn-1$
j:=j+1>>$
y:=lincof(lhside,var)$
x:=x+var*y$
gentran acof(nn,j)::=:y$
nn:=nn-1$
j:=j+1$
for each a in svars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran acof(nn,j)::=:y$
nn:=nn-1$
j:=j+1>>$
gentran arhs(n):=:rhside$
gentemp eval(var):=arhs(n)$
if not(x-lhside=0) then write " For equation ",vareq," given only ",
"variables ",pvars,svars,var$
return
end$
procedure genvareq!-linpack!-trid(pvars,svars,vareq,nlc,nzero,mode)$
begin
scalar var,rhside,lhside,x,y$
var:=first vareq$
vareq:=frst vareq$
rhside:=rhs vareq$
lhside:=lhs vareq$
x:=0$
for each a in pvars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran alcd(n):=:y >>$
y:=lincof(lhside,var)$
x:=x+var*y$
gentran ad(n):=:y$
for each a in svars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran aucd(n):=:y >>$
gentran arhs(n):=:rhside$
gentemp eval(var):=arhs(n)$
if not(x-lhside=0) then write " For equation ",vareq," given only ",
"variables ",pvars,svars,var$
return
end$
procedure genvareq!-nag(pvars,svars,vareq,nlc,nzero,mode)$
% Generates N-th row of coeff. matrix ACOF and right hand side ARHS
% according to the v-equation VAREQ.
% NZERO is number of zeroes before or after (according to MODE).
% Matrix ACOF is transformed to NAG band matrix storage.
% NCOL!* is the band width.
begin
integer j$
scalar var,rhside,lhside,x,y$
if not(length pvars + length svars+1+nzero=ncol!*) then return
write" Unconsistent PVARS:",pvars," SVARS:",svars," NZERO:",nzero$
if nlc=1 and ncol!*=3 then return
genvareq!-nag!-trid(pvars,svars,vareq,nlc,nzero,mode)$
var:=first vareq$
vareq:=frst vareq$
rhside:=rhs vareq$
lhside:=lhs vareq$
j:=1$
x:=0$
for each a in pvars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran acof(eval j,n):=:y$
j:=j+1>>$
y:=lincof(lhside,var)$
x:=x+var*y$
gentran acof(eval j,n):=:y$
j:=j+1$
for each a in svars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran acof(eval j,n):=:y$
j:=j+1>>$
gentran arhs(n):=:rhside$
gentemp eval(var):=arhs(n)$
if not(x-lhside=0) then write " For equation ",vareq," given only ",
"variables ",pvars,svars,var$
return
end$
procedure genvareq!-nag!-trid(pvars,svars,vareq,nlc,nzero,mode)$
begin
scalar var,rhside,lhside,x,y$
var:=first vareq$
vareq:=frst vareq$
rhside:=rhs vareq$
lhside:=lhs vareq$
x:=0$
for each a in pvars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran alcd(n):=:y >>$
y:=lincof(lhside,var)$
x:=x+var*y$
gentran ad(n):=:y$
for each a in svars do
<<y:=lincof(lhside,a)$
x:=x+a*y$
gentran aucd(n+1):=:y >>$
gentran arhs(n):=:rhside$
gentemp eval(var):=arhs(n)$
if not(x-lhside=0) then write " For equation ",vareq," given only ",
"variables ",pvars,svars,var$
return
end$
procedure lincof(expre,ker)$
% Expression EXPRE is linear in kernel KER.
% Returns coeff. of KER in EXPRE.
(expre-sube(ker=0,expre))/ker$
stackdolabel!*:={}$
procedure nextvareqsys(vareq,system)$
% Looks for the next v-equation. Returns the new v-equation . SYSTEM.
% During get into the loop generates the beginning of the loop,
% during get out of the loop generates end of the loop.
if rest system={} then {} . {}
else if ffst system=do then nextvesdo(vareq,system)
else if ffrst system=do then nextvesdofst(rest system)
else frst system . rest system$
procedure nextvesdofst(system)$
% Get into the loop
begin
scalar id,from,to1,step$
id:=frfst system$
from:=frst id$
to1:=frrst id$
step:=frrrst id$
id:=first id$
genstmtnum!*:=genstmtnum!*+genstmtincr!*$
gentran literal tab!*,"do ",eval(genstmtnum!*)," ",eval(id),"=",
eval(from),",",eval(to1),",",eval(step),cr!*$
stackdolabel!*:=genstmtnum!* . stackdolabel!*$
genstmtnum!*:=genstmtnum!*+genstmtincr!*$
gentemp <<literal tab!*,"do ",eval(genstmtnum!*)," ",
eval(id),"=",eval(from),
",",eval(to1),",",eval(step),cr!*>>$
fortcurrind!*:=fortcurrind!* + 4$
stackdolabel!*:=genstmtnum!* . stackdolabel!*$
id:=frrfst system . system$
return id
end$
procedure nextvesdo(vareq,system)$
% SYSTEM begins with a loop - test on the end of loop.
% Suppose that after the loop cannot be another loop, which
% follows from EXPANDDO.
begin
scalar vareqs$
vareqs:=rrfst system$
while not(first vareq=ffst vareqs and rest vareq=rfst vareqs) and
not(rest vareqs={}) do vareqs:=rest vareqs$
vareqs:=rest vareqs$
if vareqs={} then
% end of loop
<<fortcurrind!* := fortcurrind!* - 4$
gentemp<<literal eval first stackdolabel!*,tab!*,"continue",
cr!*>>$
stackdolabel!*:=rest stackdolabel!*$
gentran literal eval first stackdolabel!*,tab!*,"continue",cr!*$
stackdolabel!*:=rest stackdolabel!*$
vareqs:=frst system . rest system >>
else vareqs:=first vareqs . system$
return vareqs
end$
procedure findpvars(nlc,cykl)$
% Looks for NLC previous variables during geting into the loop
begin
scalar id,step$
id:=frst cykl$
step:=frrrst id$
id:=first id$
cykl:=reverse rrst cykl$
id:=reverse fsvars1do(nlc,
do . (id . (id-step) . 0 . (-step) . {}) . cykl)$
return id
end$
procedure lastvars(nlc,cykl)$
% Looks for the NLC last variables of the loop CYKL
begin
scalar id,step,to1$
id:=frst cykl$
to1:=frrst id$
step:=frrrst id$
id:=first id$
cykl:=reverse rrst cykl$
id:=reverse fsvars1do(nlc,do . (id . to1 . 0 . (-step) . {}) . cykl)$
return id
end$
symbolic$
flag('(ffst frst rfst rrst fffst ffrst frfst frrst rffst rfrst rrfst
rrrst ffffst fffrst ffrfst ffrrst frffst frfrst frrfst frrrst
rfffst rffrst rfrfst rfrrst rrffst rrfrst rrrfst rrrrst
findsvars findsvars1 fsvars1do findsvarsdo expanddo expddo
genvareq nextvareqsys nextvesdofst nextvesdo findpvars lastvars),
'noval)$
procedure equalaeval u$
'equal . aevlis u$
procedure aevlis u$
for each a in u collect aeval a$
procedure listnoeval(u,v)$
if atom u then listnoeval(cadr get(u,'avalue),v)
else u$
endmodule;
end;