Net2/usr/src/usr.bin/pascal/pdx/test/parall.p
(*
* Copyright (c) 1980 The Regents of the University of California.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. All advertising materials mentioning features or use of this software
* must display the following acknowledgement:
* This product includes software developed by the University of
* California, Berkeley and its contributors.
* 4. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* @(#)parall.p 5.1 (Berkeley) 4/16/91
*)
program Parall(input,output);
{Declare constants for unit conversions, convergence tests, etc.}
const SQRT2 = 1.4142136; {Square root of 2}
TWOPI = 6.2831853; {Two times pi}
MINALPHA = 0.001; {Minimum y-step size}
ROUGHLYZERO = 0.001; {Approximation to zero for convergence}
YTHRESHOLD = 40.0; {Heuristic constant for thresholding}
N = 8; {Vector and matrix dimension}
{Declare types for circuit elements, vectors, and matrices.}
type VSOURCE = record
ampl: real; {Volts (peak)}
freq: real; {Radians/second}
xindex: integer; {Index for x value}
yindex: integer; {Index for y value}
end;
RLPAIR = record
r: real; {Ohms}
l: real; {Henries}
islope: real; {Amps/second}
invariant: real; {Trapezoidal invariant}
lasttime: real; {Most recent time}
xindex: array [1..2] of integer; {x value indices}
yindex: array [1..2] of integer; {y value indices}
end;
CAPACITOR = record
c: real; {Farads}
vslope: real; {Volts/second}
invariant: real; {Trapezoidal invariant}
lasttime: real; {Most recent time}
xindex: array [1..2] of integer; {x value indices}
yindex: array [1..2] of integer; {y value indices}
end;
VECTOR = array [1..N] of real;
MATRIX = array [1..N,1..N] of real;
{Declare global variables for central simulation information.}
var ground: VSOURCE; {Ground -- a fake source}
itcount: integer; {Main routine iteration count}
update: integer; {Update loop counter for main}
optimcount: integer; {Number of optimization steps}
newtoncount: integer; {Number of Newton steps}
v1: VSOURCE; {Voltage source V1}
rl1: RLPAIR; {R1/L1 resistor-inductor pair}
rl2: RLPAIR; {R2/L2 resistor-inductor pair}
c1: CAPACITOR; {Capacitor C1}
a: MATRIX; {Global matrix A}
b: MATRIX; {Global matrix B}
jac: MATRIX; {Global Jacobian matrix}
x: VECTOR; {Global vector of dependents}
xnext: VECTOR; {Next x-vector for simulation}
y: VECTOR; {Global vector of independents}
ynext: VECTOR; {Next y-vector for simulation}
gradient: VECTOR; {Global gradient vector}
h: real; {Time step value}
time: real; {Current time}
lastychange: real; {YStep's most recent y-change}
timestep: integer; {Current timestep number}
maxsteps: integer; {Number of time steps to run}
oldxnorm: real; {Old one-norm of x-vector}
newxnorm: real; {New one-norm of x-vector}
closenough: boolean; {Flag to indicate convergence}
{The following procedure initializes everything for the program based
on the little test circuit suggested by Sarosh. The user is asked
to specify the simulation and circuit parameters, and then the matrix
and vector values are set up.}
procedure InitializeEverything;
var i,j: integer;
rtemp: real;
begin
{Ready the input and output files (almost nil for Berkeley).}
writeln(output);
writeln(output,'*** Simulation Output Record ***');
writeln(output);
writeln(output);
{Initialize convergence test/indication variables.}
oldxnorm := 0.0;
newxnorm := 0.0;
lastychange := 0.0;
{Get desired time step size for simulation.}
readln(input,h);
writeln(output,'h (Seconds): ',h:12:7);
{Get desired number of time steps for simulation.}
readln(input,maxsteps);
writeln(output,'maxsteps: ',maxsteps:4);
{Get parameters for source V1 and initialize the source.}
with v1 do
begin
readln(input,rtemp);
writeln(output,'V1 (Volts RMS): ',rtemp:9:3);
ampl := rtemp * SQRT2;
readln(input,rtemp);
writeln(output,'f (Hertz): ',rtemp:9:3);
freq := rtemp * TWOPI;
xindex := 1;
yindex := 1;
end;
{Get parameters for R1/L1 pair and initialize the pair.}
with rl1 do
begin
readln(input,r);
writeln(output,'R1 (Ohms): ',r:9:3);
readln(input,l);
writeln(output,'L1 (Henries): ',l:12:7);
islope := 0.0;
invariant := 0.0;
lasttime := -1.0; {Negative to force first update}
xindex[1] := 2;
yindex[1] := 2;
xindex[2] := 3;
yindex[2] := 3;
end;
{Get parameters for R2/L2 pair and initialize the pair.}
with rl2 do
begin
readln(input,r);
writeln(output,'R2 (Ohms): ',r:9:3);
readln(input,l);
writeln(output,'L2 (Henries): ',l:12:7);
islope := 0.0;
invariant := 0.0;
lasttime := -1.0; {Negative to force first update}
xindex[1] := 4;
yindex[1] := 4;
xindex[2] := 5;
yindex[2] := 5;
end;
{Get parameters for capacitor C1 and initialize the device.}
with c1 do
begin
readln(input,c);
writeln(output,'C1 (Farads): ',c:12:7);
vslope := 0.0;
invariant := 0.0;
lasttime := -1.0; {Negative to force first update}
xindex[1] := 6;
yindex[1] := 6;
xindex[2] := 7;
yindex[2] := 7;
end;
{Initialize the ground node.}
with ground do
begin
ampl := 0.0;
freq := 0.0;
xindex := 8;
yindex := 8;
end;
{Zero out all the vectors and matrices.}
for i := 1 to N do
begin
x[i] := 0.0;
y[i] := 0.0;
for j := 1 to N do
begin
a[i,j] := 0.0;
b[i,j] := 0.0;
jac[i,j] := 0.0;
end;
end;
{Initialize the A matrix.}
a[1,2] := -1.0;
a[2,3] := 1.0;
a[2,4] := -1.0;
a[3,5] := 1.0;
a[4,7] := 1.0;
a[5,1] := 1.0;
a[7,6] := 1.0;
a[8,8] := 1.0;
{Initialize the B matrix.}
b[1,1] := 1.0;
b[3,7] := -1.0;
b[4,8] := -1.0;
b[5,2] := 1.0;
b[6,3] := 1.0;
b[6,4] := 1.0;
b[7,5] := 1.0;
b[8,6] := 1.0;
{Initialize the Jacobian matrix.}
rtemp := h / (2.0 * rl1.l + rl1.r * h);
jac[2,2] := rtemp;
jac[3,3] := rtemp;
jac[2,3] := -rtemp;
jac[3,2] := -rtemp;
rtemp := h / (2.0 * rl2.l + rl2.r * h);
jac[4,4] := rtemp;
jac[5,5] := rtemp;
jac[4,5] := -rtemp;
jac[5,4] := -rtemp;
jac[6,6] := -1.0;
jac[7,7] := 1.0;
jac[7,6] := h / (2.0 * c1.c);
end;
{The following procedure solves the equation Ax=b for an N x N system
of linear equations, where A is the coefficient matrix, b is the
right-hand-side vector, and x is the vector of unknowns. Gaussian
elimination with maximal column pivots is used. }
procedure Solve(a: MATRIX; b: VECTOR; var x: VECTOR);
var y,z: real;
i,j,k,k1: integer;
begin
for k := 1 to N-1 do
begin
y := abs(a[k,k]);
j := k;
k1 := k + 1;
for i := k1 to N do
if abs(a[i,k]) > y then
begin
j := i;
y := abs(a[i,k]);
end;
for i := k to N do
begin
y := a[k,i];
a[k,i] := a[j,i];
a[j,i] := y;
end;
y := b[k];
b[k] := b[j];
b[j] := y;
z := a[k,k];
for i := k1 to N do
begin
y := a[i,k] / z;
a[i,k] := y;
for j := k1 to N do a[i,j] := a[i,j] - y * a[k,j];
b[i] := b[i] - y * b[k];
end;
end;
x[N] := b[N] / a[N,N];
for i := N-1 downto 1 do
begin
y := b[i];
for j := i+1 to N do y := y - a[i,j] * x[j];
x[i] := y / a[i,i];
end;
end;
{The following procedure computes and returns a vector called "deltay",
which is the change in the y-vector prescribed by the Newton-Rhapson
algorithm.}
procedure NewtonStep(var deltay: VECTOR);
var phi: VECTOR;
m: MATRIX;
i,j,k: integer;
begin
for i := 1 to N do
begin
phi[i] := 0.0;
for j := 1 to N do
begin
phi[i] := phi[i] + a[i,j] * y[j] + b[i,j] * x[j];
m[i,j] := -a[i,j];
for k := 1 to N do m[i,j] := m[i,j] - b[i,k] * jac[k,j];
end;
end;
Solve(m,phi,deltay);
end;
{The following function computes the value of theta, the objective
function, given the x and y vectors.}
function ThetaValue(x,y: VECTOR): real;
var i,j: integer;
phielem: real;
theta: real;
begin
theta := 0.0;
for i:= 1 to N do
begin
phielem := 0.0;
for j := 1 to N do
phielem := phielem + a[i,j] * y[j] + b[i,j] * x[j];
theta := theta + phielem * phielem;
end;
ThetaValue := theta;
end;
{The following function computes the theta value associated with a
proposed step of size alpha in the direction of the gradient.}
function Theta(alpha: real): real;
var ythere: VECTOR;
i: integer;
begin
for i := 1 to N do
ythere[i] := y[i] - alpha * gradient[i];
Theta := ThetaValue(x,ythere);
end;
{The following procedure computes the gradient of the objective
function (theta) with respect to the vector y.}
procedure ComputeGradient;
var i,j,k: integer;
m: MATRIX;
v: VECTOR;
begin
{Compute v = Ay + Bx and M = A' + J'B'.}
for i := 1 to N do
begin
v[i] := 0.0;
for j := 1 to N do
begin
v[i] := v[i] + a[i,j] * y[j] + b[i,j] * x[j];
m[i,j] := a[j,i];
for k := 1 to N do
m[i,j] := m[i,j] + jac[k,i] * b[j,k];
end;
end;
{Compute gradient = 2Mv.}
for i := 1 to N do
begin
gradient[i] := 0.0;
for j := 1 to N do
gradient[i] := gradient[i] + m[i,j] * v[j];
gradient[i] := 2.0 * gradient[i];
end;
end;
{The following procedure computes the bounds on alpha, the step size
to take in the gradient direction. The bounding is done according
to an algorithm suggested in S.W.Director's text on circuits.}
procedure GetAlphaBounds(var lower,upper: real);
var alpha: real;
oldtheta,newtheta: real;
begin
if ThetaValue(x,y) = 0.0
then
begin
lower := 0;
upper := 0;
end
else
begin
lower := MINALPHA;
oldtheta := Theta(lower);
upper := MINALPHA * 2.0;
newtheta := Theta(upper);
if newtheta <= oldtheta then
begin
alpha := upper;
repeat
begin
oldtheta := newtheta;
alpha := alpha * 2.0;
newtheta := Theta(alpha);
end
until (newtheta > oldtheta);
lower := alpha / 4.0;
upper := alpha;
end;
end;
end;
{The following function computes the best value of alpha for the step
in the gradient direction. This best value is the value that minimizes
theta along the gradient-directed path.}
function BestAlpha(lower,upper: real): real;
var alphaL,alphaU,alpha1,alpha2: real;
thetaL,thetaU,theta1,theta2: real;
begin
alphaL := lower;
thetaL := Theta(alphaL);
alphaU := upper;
thetaU := Theta(alphaU);
alpha1 := 0.381966 * alphaU + 0.618034 * alphaL;
theta1 := Theta(alpha1);
alpha2 := 0.618034 * alphaU + 0.381966 * alphaL;
theta2 := Theta(alpha2);
repeat
if theta1 < theta2
then
begin
alphaU := alpha2;
thetaU := theta2;
alpha2 := alpha1;
theta2 := theta1;
alpha1 := 0.381966 * alphaU + 0.618034 * alphaL;
theta1 := Theta(alpha1);
end
else
begin
alphaL := alpha1;
thetaL := theta1;
alpha1 := alpha2;
theta1 := theta2;
alpha2 := 0.618034 * alphaU + 0.381966 * alphaL;
theta2 := Theta(alpha2);
end
until abs(thetaU - thetaL) <= ROUGHLYZERO;
BestAlpha := (alphaL + alphaU) / 2.0;
end;
{The following procedure computes and returns a vector called "deltay",
which is the change in the y-vector prescribed by the steepest-descent
approach.}
procedure OptimizationStep(var deltay: VECTOR);
var lower,upper: real;
alpha: real;
i: integer;
begin
ComputeGradient;
GetAlphaBounds(lower,upper);
if lower <> upper then
begin
alpha := BestAlpha(lower,upper);
for i:= 1 to N do deltay[i] := - alpha * gradient[i];
end;
end;
{The following function computes the one-norm of a vector argument.
The length of the argument vector is assumed to be N.}
function OneNorm(vec: VECTOR): real;
var sum: real;
i: integer;
begin
sum := 0;
for i := 1 to N do sum := sum + abs(vec[i]);
OneNorm := sum;
end;
{The following procedure takes a y-step, using the optimization
approach when far from the solution and the Newton-Rhapson approach
when fairly close to the solution.}
procedure YStep;
var deltay: VECTOR;
ychange: real;
scale: real;
i: integer;
begin
NewtonStep(deltay);
ychange := OneNorm(deltay);
if ychange > YTHRESHOLD
then
{
begin
OptimizationStep(deltay);
ychange := OneNorm(deltay);
if ychange > YTHRESHOLD then
}
begin
scale := YTHRESHOLD/ychange;
for i := 1 to N do deltay[i] := scale * deltay[i];
optimcount := optimcount + 1;
end {;}
{
optimcount := optimcount + 1;
end
}
else
begin
newtoncount := newtoncount + 1;
end;
for i := 1 to N do ynext[i] := y[i] + deltay[i];
end;
{The following procedure updates the output of a voltage source
given the current time.}
procedure VsourceStep(vn: VSOURCE);
begin
with vn do
xnext[xindex] := ampl * sin(freq * time);
end;
{The following procedure updates the outputs of a resistor-inductor
pair given the time step to take...that is, this routine takes a
time step for resistor-inductor pairs. The new outputs are found
by applying the trapezoidal rule.}
procedure RLPairStep(var rln: RLPAIR);
begin
with rln do
begin
if (time > lasttime) then
begin
lasttime := time;
invariant := xnext[xindex[1]] + (h / 2.0) * islope;
end;
islope := (y[yindex[1]] - y[yindex[2]] - r * xnext[xindex[1]]) / l;
xnext[xindex[1]] := invariant + (h / 2.0) * islope;
xnext[xindex[2]] := - xnext[xindex[1]];
end;
end;
{The following procedure updates the outputs of a capacitor given the
time step...it takes the time step using a trapezoidal rule iteration.}
procedure CapacitorStep(var cn: CAPACITOR);
var v: real;
begin
with cn do
begin
if (time > lasttime) then
begin
lasttime := time;
v := xnext[xindex[2]] - y[yindex[2]];
invariant := v + (h / 2.0) * vslope;
end;
vslope := y[yindex[1]] / c;
v := invariant + (h / 2.0) * vslope;
xnext[xindex[1]] := - y[yindex[1]];
xnext[xindex[2]] := y[yindex[2]] + v;
end;
end;
{The following procedure controls the overall x-step for the
specific circuit to be simulated.}
procedure XStep;
begin
VsourceStep(v1);
RLPairStep(rl1);
RLPairStep(rl2);
CapacitorStep(c1);
VsourceStep(ground);
end;
{The following procedure prints out the values of all the voltages and
currents for the circuit at each time step.}
procedure PrintReport;
begin
writeln(output);
writeln(output);
writeln(output,'REPORT at Time = ',time:8:5,' seconds');
writeln(output,'Number of iterations used: ',itcount:2);
writeln(output,'Number of truncations: ',optimcount:2);
writeln(output,'Number of Newton y-steps: ',newtoncount:2);
writeln(output,'The voltages and currents are:');
writeln(output,' V1 = ',x[1]:12:7,' I1 = ',y[1]:12:7);
writeln(output,' V2 = ',y[2]:12:7,' I2 = ',x[2]:12:7);
writeln(output,' V3 = ',y[3]:12:7,' I3 = ',x[3]:12:7);
writeln(output,' V4 = ',y[4]:12:7,' I4 = ',x[4]:12:7);
writeln(output,' V5 = ',y[5]:12:7,' I5 = ',x[5]:12:7);
writeln(output,' V6 = ',x[7]:12:7,' I6 = ',y[6]:12:7);
writeln(output,' V7 = ',y[7]:12:7,' I7 = ',x[6]:12:7);
writeln(output,' V8 = ',x[8]:12:7,' I8 = ',y[8]:12:7);
end;
{Finally, the main routine controls the whole simulation process.}
begin
InitializeEverything;
PrintReport;
for timestep := 1 to maxsteps do
begin
itcount := 0;
optimcount := 0;
newtoncount := 0;
closenough := FALSE;
time := h * timestep;
repeat
begin
itcount := itcount + 1;
YStep;
XStep;
for update := 1 to N do
begin
x[update] := xnext[update];
y[update] := ynext[update];
end;
oldxnorm := newxnorm;
newxnorm := OneNorm(x);
if abs(newxnorm - oldxnorm) <= ROUGHLYZERO
then closenough := TRUE;
end;
if itcount < 4 then closenough := FALSE;
until (itcount = 25) or closenough;
PrintReport;
end;
end.