parall.p

早期freebsd实现

    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.