parall.p

早期freebsd实现

(*
 * Copyright (c) 1980, 1993
 *	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	8.1 (Berkeley) 6/6/93
 *)

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;