Return to parall.p CVS log

Up to [CSRG BSD Unix] / 43BSD / ucb / pascal / pdx / test

Annotation of 43BSD/ucb/pascal/pdx/test/parall.p, revision 1.1

1.1     ! root        1: program Parall(input,output);
        !             2: 
        !             3: {Declare constants for unit conversions, convergence tests, etc.}
        !             4: 
        !             5: const  SQRT2 = 1.4142136;      {Square root of 2}
        !             6:        TWOPI = 6.2831853;      {Two times pi}
        !             7:        MINALPHA = 0.001;       {Minimum y-step size}
        !             8:        ROUGHLYZERO = 0.001;    {Approximation to zero for convergence}
        !             9:        YTHRESHOLD = 40.0;      {Heuristic constant for thresholding}
        !            10:        N = 8;                  {Vector and matrix dimension}
        !            11: 
        !            12: 
        !            13: {Declare types for circuit elements, vectors, and matrices.}
        !            14: 
        !            15: type   VSOURCE = record
        !            16:                        ampl: real;      {Volts (peak)}
        !            17:                        freq: real;      {Radians/second}
        !            18:                        xindex: integer;        {Index for x value}
        !            19:                        yindex: integer;        {Index for y value}
        !            20:                  end;
        !            21: 
        !            22:        RLPAIR = record
        !            23:                        r: real;         {Ohms}
        !            24:                        l: real;         {Henries}
        !            25:                        islope: real;    {Amps/second}
        !            26:                        invariant: real; {Trapezoidal invariant}
        !            27:                        lasttime: real;  {Most recent time}
        !            28:                        xindex: array [1..2] of integer;  {x value indices}
        !            29:                        yindex: array [1..2] of integer;  {y value indices}
        !            30:                 end;
        !            31: 
        !            32:        CAPACITOR = record
        !            33:                        c: real;         {Farads}
        !            34:                        vslope: real;    {Volts/second}
        !            35:                        invariant: real; {Trapezoidal invariant}
        !            36:                        lasttime: real;  {Most recent time}
        !            37:                        xindex: array [1..2] of integer;  {x value indices}
        !            38:                        yindex: array [1..2] of integer;  {y value indices}
        !            39:                    end;
        !            40: 
        !            41:        VECTOR = array [1..N] of real;
        !            42: 
        !            43:        MATRIX = array [1..N,1..N] of real;
        !            44: 
        !            45: 
        !            46: {Declare global variables for central simulation information.}
        !            47: 
        !            48: var    ground: VSOURCE;        {Ground -- a fake source}
        !            49:        itcount: integer;       {Main routine iteration count}
        !            50:        update: integer;        {Update loop counter for main}
        !            51:        optimcount: integer;    {Number of optimization steps}
        !            52:        newtoncount: integer;   {Number of Newton steps}
        !            53:        v1: VSOURCE;            {Voltage source V1}
        !            54:        rl1: RLPAIR;            {R1/L1 resistor-inductor pair}
        !            55:        rl2: RLPAIR;            {R2/L2 resistor-inductor pair}
        !            56:        c1: CAPACITOR;          {Capacitor C1}
        !            57:        a: MATRIX;              {Global matrix A}
        !            58:        b: MATRIX;              {Global matrix B}
        !            59:        jac: MATRIX;            {Global Jacobian matrix}
        !            60:        x: VECTOR;              {Global vector of dependents}
        !            61:        xnext: VECTOR;          {Next x-vector for simulation}
        !            62:        y: VECTOR;              {Global vector of independents}
        !            63:        ynext: VECTOR;          {Next y-vector for simulation}
        !            64:        gradient: VECTOR;       {Global gradient vector}
        !            65:        h: real;                {Time step value}
        !            66:        time: real;             {Current time}
        !            67:        lastychange: real;      {YStep's most recent y-change}
        !            68:        timestep: integer;      {Current timestep number}
        !            69:        maxsteps: integer;      {Number of time steps to run}
        !            70:         oldxnorm: real;                {Old one-norm of x-vector}
        !            71:        newxnorm: real;         {New one-norm of x-vector}
        !            72:        closenough: boolean;    {Flag to indicate convergence}
        !            73: 
        !            74: 
        !            75: 
        !            76: 
        !            77: {The following procedure initializes everything for the program based
        !            78:  on the little test circuit suggested by Sarosh.  The user is asked
        !            79:  to specify the simulation and circuit parameters, and then the matrix
        !            80:  and vector values are set up.}
        !            81: 
        !            82: procedure InitializeEverything;
        !            83: var i,j: integer;
        !            84:     rtemp: real;
        !            85: begin
        !            86: 
        !            87: {Ready the input and output files (almost nil for Berkeley).}
        !            88: writeln(output);
        !            89: writeln(output,'*** Simulation Output Record ***');
        !            90: writeln(output);
        !            91: writeln(output);
        !            92: 
        !            93: {Initialize convergence test/indication variables.}
        !            94: oldxnorm := 0.0;
        !            95: newxnorm := 0.0;
        !            96: lastychange := 0.0;
        !            97: 
        !            98: {Get desired time step size for simulation.}
        !            99: readln(input,h);
        !           100: writeln(output,'h (Seconds): ',h:12:7);
        !           101: 
        !           102: {Get desired number of time steps for simulation.}
        !           103: readln(input,maxsteps);
        !           104: writeln(output,'maxsteps: ',maxsteps:4);
        !           105: 
        !           106: {Get parameters for source V1 and initialize the source.}
        !           107: with v1 do
        !           108:    begin
        !           109:    readln(input,rtemp);
        !           110:    writeln(output,'V1 (Volts RMS): ',rtemp:9:3);
        !           111:    ampl := rtemp * SQRT2;
        !           112:    readln(input,rtemp);
        !           113:    writeln(output,'f (Hertz): ',rtemp:9:3);
        !           114:    freq := rtemp * TWOPI;
        !           115:    xindex := 1;
        !           116:    yindex := 1;
        !           117:    end;
        !           118: 
        !           119: {Get parameters for R1/L1 pair and initialize the pair.}
        !           120: with rl1 do
        !           121:    begin
        !           122:    readln(input,r);
        !           123:    writeln(output,'R1 (Ohms): ',r:9:3);
        !           124:    readln(input,l);
        !           125:    writeln(output,'L1 (Henries): ',l:12:7);
        !           126:    islope := 0.0;
        !           127:    invariant := 0.0;
        !           128:    lasttime := -1.0;        {Negative to force first update}
        !           129:    xindex[1] := 2;
        !           130:    yindex[1] := 2;
        !           131:    xindex[2] := 3;
        !           132:    yindex[2] := 3;
        !           133:    end;
        !           134: 
        !           135: {Get parameters for R2/L2 pair and initialize the pair.}
        !           136: with rl2 do
        !           137:    begin
        !           138:    readln(input,r);
        !           139:    writeln(output,'R2 (Ohms): ',r:9:3);
        !           140:    readln(input,l);
        !           141:    writeln(output,'L2 (Henries): ',l:12:7);
        !           142:    islope := 0.0;
        !           143:    invariant := 0.0;
        !           144:    lasttime := -1.0;      {Negative to force first update}
        !           145:    xindex[1] := 4;
        !           146:    yindex[1] := 4;
        !           147:    xindex[2] := 5;
        !           148:    yindex[2] := 5;
        !           149:    end;
        !           150: 
        !           151: {Get parameters for capacitor C1 and initialize the device.}
        !           152: with c1 do
        !           153:    begin
        !           154:    readln(input,c);
        !           155:    writeln(output,'C1 (Farads): ',c:12:7);
        !           156:    vslope := 0.0;
        !           157:    invariant := 0.0;
        !           158:    lasttime := -1.0;      {Negative to force first update}
        !           159:    xindex[1] := 6;
        !           160:    yindex[1] := 6;
        !           161:    xindex[2] := 7;
        !           162:    yindex[2] := 7;
        !           163:    end;
        !           164: 
        !           165: {Initialize the ground node.}
        !           166: with ground do
        !           167:    begin
        !           168:    ampl := 0.0;
        !           169:    freq := 0.0;
        !           170:    xindex := 8;
        !           171:    yindex := 8;
        !           172:    end;
        !           173: 
        !           174: {Zero out all the vectors and matrices.}
        !           175: for i := 1 to N do
        !           176:    begin
        !           177:    x[i] := 0.0;
        !           178:    y[i] := 0.0;
        !           179:    for j := 1 to N do
        !           180:       begin
        !           181:       a[i,j] := 0.0;
        !           182:       b[i,j] := 0.0;
        !           183:       jac[i,j] := 0.0;
        !           184:       end;
        !           185:    end;
        !           186: 
        !           187: {Initialize the A matrix.}
        !           188: a[1,2] := -1.0;
        !           189: a[2,3] := 1.0;
        !           190: a[2,4] := -1.0;
        !           191: a[3,5] := 1.0;
        !           192: a[4,7] := 1.0;
        !           193: a[5,1] := 1.0;
        !           194: a[7,6] := 1.0;
        !           195: a[8,8] := 1.0;
        !           196: 
        !           197: {Initialize the B matrix.}
        !           198: b[1,1] := 1.0;
        !           199: b[3,7] := -1.0;
        !           200: b[4,8] := -1.0;
        !           201: b[5,2] := 1.0;
        !           202: b[6,3] := 1.0;
        !           203: b[6,4] := 1.0;
        !           204: b[7,5] := 1.0;
        !           205: b[8,6] := 1.0;
        !           206: 
        !           207: {Initialize the Jacobian matrix.}
        !           208: rtemp := h / (2.0 * rl1.l  +  rl1.r * h);
        !           209: jac[2,2] := rtemp;
        !           210: jac[3,3] := rtemp;
        !           211: jac[2,3] := -rtemp;
        !           212: jac[3,2] := -rtemp;
        !           213: rtemp := h / (2.0 * rl2.l  +  rl2.r * h);
        !           214: jac[4,4] := rtemp;
        !           215: jac[5,5] := rtemp;
        !           216: jac[4,5] := -rtemp;
        !           217: jac[5,4] := -rtemp;
        !           218: jac[6,6] := -1.0;
        !           219: jac[7,7] := 1.0;
        !           220: jac[7,6] := h / (2.0 * c1.c);
        !           221: end;
        !           222: 
        !           223: 
        !           224: 
        !           225: 
        !           226: {The following procedure solves the equation Ax=b for an N x N system
        !           227:  of linear equations, where A is the coefficient matrix, b is the
        !           228:  right-hand-side vector, and x is the vector of unknowns.  Gaussian
        !           229:  elimination with maximal column pivots is used.                     }
        !           230: 
        !           231: procedure Solve(a: MATRIX; b: VECTOR; var x: VECTOR);
        !           232: var y,z: real;
        !           233:     i,j,k,k1: integer;
        !           234: begin
        !           235: for k := 1 to N-1 do
        !           236:    begin
        !           237:    y := abs(a[k,k]);
        !           238:    j := k;
        !           239:    k1 := k + 1;
        !           240:    for i := k1 to N do
        !           241:       if abs(a[i,k]) > y then
        !           242:          begin
        !           243:          j := i;
        !           244:          y := abs(a[i,k]);
        !           245:          end;
        !           246:    for i := k to N do
        !           247:       begin
        !           248:       y := a[k,i];
        !           249:       a[k,i] := a[j,i];
        !           250:       a[j,i] := y;
        !           251:       end;
        !           252:    y := b[k];
        !           253:    b[k] := b[j];
        !           254:    b[j] := y;
        !           255:    z := a[k,k];
        !           256:    for i := k1 to N do
        !           257:       begin
        !           258:       y := a[i,k] / z;
        !           259:       a[i,k] := y;
        !           260:       for j := k1 to N do a[i,j] := a[i,j]  -  y * a[k,j];
        !           261:       b[i] := b[i] - y * b[k];
        !           262:       end;
        !           263:    end;
        !           264: x[N] := b[N] / a[N,N];
        !           265: for i := N-1 downto 1 do
        !           266:    begin
        !           267:    y := b[i];
        !           268:    for j := i+1 to N do y := y  -  a[i,j] * x[j];
        !           269:    x[i] := y / a[i,i];
        !           270:    end;
        !           271: end;
        !           272: 
        !           273: 
        !           274: {The following procedure computes and returns a vector called "deltay",
        !           275:  which is the change in the y-vector prescribed by the Newton-Rhapson
        !           276:  algorithm.}
        !           277: 
        !           278: procedure NewtonStep(var deltay: VECTOR);
        !           279: var phi: VECTOR;
        !           280:     m: MATRIX;
        !           281:     i,j,k: integer;
        !           282: begin
        !           283: for i := 1 to N do
        !           284:    begin
        !           285:    phi[i] := 0.0;
        !           286:    for j := 1 to N do
        !           287:       begin
        !           288:       phi[i] := phi[i]  +  a[i,j] * y[j]  +  b[i,j] * x[j];
        !           289:       m[i,j] := -a[i,j];
        !           290:       for k := 1 to N do m[i,j] := m[i,j] - b[i,k] * jac[k,j];
        !           291:       end;
        !           292: 
        !           293:    end;
        !           294: Solve(m,phi,deltay);
        !           295: end;
        !           296: 
        !           297: 
        !           298: 
        !           299: 
        !           300: {The following function computes the value of theta, the objective
        !           301:  function, given the x and y vectors.}
        !           302: 
        !           303: function ThetaValue(x,y: VECTOR): real;
        !           304: var i,j: integer;
        !           305:     phielem: real;
        !           306:     theta: real;
        !           307: begin
        !           308: theta := 0.0;
        !           309: for i:= 1 to N do
        !           310:    begin
        !           311:    phielem := 0.0;
        !           312:    for j := 1 to N do
        !           313:       phielem := phielem  +  a[i,j] * y[j]  +  b[i,j] * x[j];
        !           314:    theta := theta  +  phielem * phielem;
        !           315:    end;
        !           316: ThetaValue := theta;
        !           317: end;
        !           318: 
        !           319: 
        !           320: {The following function computes the theta value associated with a
        !           321:  proposed step of size alpha in the direction of the gradient.}
        !           322: 
        !           323: function Theta(alpha: real): real;
        !           324: var ythere: VECTOR;
        !           325:     i: integer;
        !           326: begin
        !           327: for i := 1 to N do
        !           328:    ythere[i] := y[i]  -  alpha * gradient[i];
        !           329: Theta := ThetaValue(x,ythere);
        !           330: end;
        !           331: 
        !           332: 
        !           333: {The following procedure computes the gradient of the objective 
        !           334:  function (theta) with respect to the vector y.}
        !           335: 
        !           336: procedure ComputeGradient;
        !           337: var i,j,k: integer;
        !           338:     m: MATRIX;
        !           339:     v: VECTOR;
        !           340: begin
        !           341: {Compute v = Ay + Bx and M = A' + J'B'.}
        !           342: for i := 1 to N do
        !           343:    begin
        !           344:    v[i] := 0.0;
        !           345:    for j := 1 to N do
        !           346:       begin
        !           347:       v[i] := v[i]  +  a[i,j] * y[j]  +  b[i,j] * x[j];
        !           348:       m[i,j] := a[j,i];
        !           349:       for k := 1 to N do
        !           350:          m[i,j] := m[i,j]  +  jac[k,i] * b[j,k];
        !           351:       end;
        !           352:    end;
        !           353: {Compute gradient = 2Mv.}
        !           354: for i := 1 to N do
        !           355:    begin
        !           356:    gradient[i] := 0.0;
        !           357:    for j := 1 to N do
        !           358:       gradient[i] := gradient[i]  +  m[i,j] * v[j];
        !           359:    gradient[i] := 2.0 * gradient[i];
        !           360:    end;
        !           361: end;
        !           362: 
        !           363: 
        !           364: {The following procedure computes the bounds on alpha, the step size
        !           365:  to take in the gradient direction.  The bounding is done according
        !           366:  to an algorithm suggested in S.W.Director's text on circuits.}
        !           367: 
        !           368: procedure GetAlphaBounds(var lower,upper: real);
        !           369: var alpha: real;
        !           370:     oldtheta,newtheta: real;
        !           371: begin
        !           372: if ThetaValue(x,y) = 0.0
        !           373:    then
        !           374:       begin
        !           375:       lower := 0;
        !           376: 
        !           377:       upper := 0;
        !           378:       end
        !           379:    else
        !           380:       begin
        !           381:       lower := MINALPHA;
        !           382:       oldtheta := Theta(lower);
        !           383:       upper := MINALPHA * 2.0;
        !           384:       newtheta := Theta(upper);
        !           385:       if newtheta <= oldtheta then
        !           386:          begin
        !           387:          alpha := upper;
        !           388:          repeat
        !           389:             begin
        !           390:             oldtheta := newtheta;
        !           391:             alpha := alpha * 2.0;
        !           392:             newtheta := Theta(alpha);
        !           393:             end
        !           394:          until (newtheta > oldtheta);
        !           395:          lower := alpha / 4.0;
        !           396:          upper := alpha;
        !           397:          end;
        !           398:       end;
        !           399: end;
        !           400: 
        !           401: 
        !           402: {The following function computes the best value of alpha for the step
        !           403:  in the gradient direction.  This best value is the value that minimizes
        !           404:  theta along the gradient-directed path.}
        !           405: 
        !           406: function BestAlpha(lower,upper: real): real;
        !           407: var alphaL,alphaU,alpha1,alpha2: real;
        !           408:     thetaL,thetaU,theta1,theta2: real;
        !           409: 
        !           410: begin
        !           411: alphaL := lower;
        !           412: thetaL := Theta(alphaL);
        !           413: alphaU := upper;
        !           414: thetaU := Theta(alphaU);
        !           415: alpha1 := 0.381966 * alphaU  +  0.618034 * alphaL;
        !           416: theta1 := Theta(alpha1);
        !           417: alpha2 := 0.618034 * alphaU  +  0.381966 * alphaL;
        !           418: theta2 := Theta(alpha2);
        !           419: repeat
        !           420:    if theta1 < theta2
        !           421:       then
        !           422:          begin
        !           423:          alphaU := alpha2;
        !           424:          thetaU := theta2;
        !           425:          alpha2 := alpha1;
        !           426:          theta2 := theta1;
        !           427:          alpha1 := 0.381966 * alphaU  +  0.618034 * alphaL;
        !           428:          theta1 := Theta(alpha1);
        !           429:          end
        !           430:       else
        !           431:          begin
        !           432:          alphaL := alpha1;
        !           433:          thetaL := theta1;
        !           434:          alpha1 := alpha2;
        !           435:          theta1 := theta2;
        !           436:          alpha2 := 0.618034 * alphaU  +  0.381966 * alphaL;
        !           437:          theta2 := Theta(alpha2);
        !           438:          end
        !           439: until abs(thetaU - thetaL) <= ROUGHLYZERO;
        !           440: BestAlpha := (alphaL + alphaU) / 2.0;
        !           441: end;
        !           442: 
        !           443: 
        !           444: {The following procedure computes and returns a vector called "deltay",
        !           445:  which is the change in the y-vector prescribed by the steepest-descent
        !           446:  approach.}
        !           447: 
        !           448: procedure OptimizationStep(var deltay: VECTOR);
        !           449: var lower,upper: real;
        !           450:     alpha: real;
        !           451:     i: integer;
        !           452: begin
        !           453: ComputeGradient;
        !           454: GetAlphaBounds(lower,upper);
        !           455: if lower <> upper then
        !           456:    begin
        !           457:    alpha := BestAlpha(lower,upper);
        !           458:    for i:= 1 to N do deltay[i] := - alpha * gradient[i];
        !           459:    end;
        !           460: end;
        !           461: 
        !           462: 
        !           463: 
        !           464: 
        !           465: {The following function computes the one-norm of a vector argument.
        !           466:  The length of the argument vector is assumed to be N.}
        !           467: 
        !           468: function OneNorm(vec: VECTOR): real;
        !           469: var sum: real;
        !           470:     i: integer;
        !           471: begin
        !           472: sum := 0;
        !           473: for i := 1 to N do  sum := sum + abs(vec[i]);
        !           474: OneNorm := sum;
        !           475: end;
        !           476: 
        !           477: 
        !           478: {The following procedure takes a y-step, using the optimization
        !           479: approach when far from the solution and the Newton-Rhapson approach
        !           480: when fairly close to the solution.}
        !           481: 
        !           482: procedure YStep;
        !           483: var deltay: VECTOR;
        !           484:     ychange: real;
        !           485:     scale: real;
        !           486:     i: integer;
        !           487: begin
        !           488: NewtonStep(deltay);
        !           489: ychange := OneNorm(deltay);
        !           490: if ychange > YTHRESHOLD
        !           491:    then
        !           492: {
        !           493:       begin
        !           494:       OptimizationStep(deltay);
        !           495:       ychange := OneNorm(deltay);
        !           496:       if ychange > YTHRESHOLD then
        !           497: }
        !           498:          begin
        !           499:          scale := YTHRESHOLD/ychange;
        !           500:          for i := 1 to N do deltay[i] := scale * deltay[i];
        !           501:         optimcount := optimcount + 1;
        !           502:         end            {;}
        !           503: {
        !           504:       optimcount := optimcount + 1;
        !           505:       end
        !           506: }
        !           507:    else
        !           508:       begin
        !           509:       newtoncount := newtoncount + 1;
        !           510:       end;
        !           511: for i := 1 to N do ynext[i] := y[i] + deltay[i];
        !           512: end;
        !           513: 
        !           514: 
        !           515: 
        !           516: 
        !           517: {The following procedure updates the output of a voltage source
        !           518:  given the current time.}
        !           519: 
        !           520: procedure VsourceStep(vn: VSOURCE);
        !           521: begin
        !           522: with vn do
        !           523:    xnext[xindex] := ampl * sin(freq * time);
        !           524: end;
        !           525: 
        !           526: 
        !           527: {The following procedure updates the outputs of a resistor-inductor
        !           528:  pair given the time step to take...that is, this routine takes a
        !           529:  time step for resistor-inductor pairs.  The new outputs are found
        !           530:  by applying the trapezoidal rule.}
        !           531: 
        !           532: procedure RLPairStep(var rln: RLPAIR);
        !           533: begin
        !           534: with rln do
        !           535:    begin
        !           536:    if (time > lasttime) then
        !           537:       begin
        !           538:       lasttime := time;
        !           539:       invariant := xnext[xindex[1]]  +  (h / 2.0) * islope;
        !           540:       end;
        !           541:    islope := (y[yindex[1]]  -  y[yindex[2]]  -  r * xnext[xindex[1]]) / l;
        !           542:    xnext[xindex[1]] := invariant  +  (h / 2.0) * islope;
        !           543:    xnext[xindex[2]] := - xnext[xindex[1]];
        !           544:    end;
        !           545: end;
        !           546: 
        !           547: 
        !           548: {The following procedure updates the outputs of a capacitor given the 
        !           549:  time step...it takes the time step using a trapezoidal rule iteration.}
        !           550: 
        !           551: procedure CapacitorStep(var cn: CAPACITOR);
        !           552: var v: real;
        !           553: begin
        !           554: with cn do
        !           555:    begin
        !           556:    if (time > lasttime) then
        !           557:       begin
        !           558:       lasttime := time;
        !           559:       v := xnext[xindex[2]]  -  y[yindex[2]];
        !           560:       invariant := v  +  (h / 2.0) * vslope;
        !           561:       end;
        !           562:    vslope := y[yindex[1]] / c;
        !           563:    v := invariant  +  (h / 2.0) * vslope;
        !           564:    xnext[xindex[1]] := - y[yindex[1]];
        !           565:    xnext[xindex[2]] := y[yindex[2]] + v;
        !           566:    end;
        !           567: end;
        !           568: 
        !           569: 
        !           570: {The following procedure controls the overall x-step for the 
        !           571:  specific circuit to be simulated.}
        !           572: 
        !           573: procedure XStep;
        !           574: begin
        !           575: VsourceStep(v1);
        !           576: RLPairStep(rl1);
        !           577: RLPairStep(rl2);
        !           578: CapacitorStep(c1);
        !           579: VsourceStep(ground);
        !           580: end;
        !           581: 
        !           582: 
        !           583: 
        !           584: 
        !           585: {The following procedure prints out the values of all the voltages and
        !           586:  currents for the circuit at each time step.}
        !           587: 
        !           588: procedure PrintReport;
        !           589: begin
        !           590: writeln(output);
        !           591: writeln(output);
        !           592: writeln(output,'REPORT at Time = ',time:8:5,' seconds');
        !           593: writeln(output,'Number of iterations used: ',itcount:2);
        !           594: writeln(output,'Number of truncations: ',optimcount:2);
        !           595: writeln(output,'Number of Newton y-steps: ',newtoncount:2);
        !           596: writeln(output,'The voltages and currents are:');
        !           597: writeln(output,'  V1 = ',x[1]:12:7,'   I1 = ',y[1]:12:7);
        !           598: writeln(output,'  V2 = ',y[2]:12:7,'   I2 = ',x[2]:12:7);
        !           599: writeln(output,'  V3 = ',y[3]:12:7,'   I3 = ',x[3]:12:7);
        !           600: writeln(output,'  V4 = ',y[4]:12:7,'   I4 = ',x[4]:12:7);
        !           601: writeln(output,'  V5 = ',y[5]:12:7,'   I5 = ',x[5]:12:7);
        !           602: writeln(output,'  V6 = ',x[7]:12:7,'   I6 = ',y[6]:12:7);
        !           603: writeln(output,'  V7 = ',y[7]:12:7,'   I7 = ',x[6]:12:7);
        !           604: writeln(output,'  V8 = ',x[8]:12:7,'   I8 = ',y[8]:12:7);
        !           605: end;
        !           606: 
        !           607: 
        !           608: 
        !           609: 
        !           610: {Finally, the main routine controls the whole simulation process.}
        !           611: 
        !           612: begin
        !           613: InitializeEverything;
        !           614: PrintReport;
        !           615: for timestep := 1 to maxsteps do
        !           616:    begin
        !           617:    itcount := 0;
        !           618:    optimcount := 0;
        !           619:    newtoncount := 0;
        !           620:    closenough := FALSE;
        !           621:    time := h * timestep;
        !           622:    repeat
        !           623:       begin
        !           624:       itcount := itcount + 1;
        !           625:       YStep;
        !           626:       XStep;
        !           627:       for update := 1 to N do
        !           628:         begin
        !           629:         x[update] := xnext[update];
        !           630:         y[update] := ynext[update];
        !           631:         end;
        !           632:       oldxnorm := newxnorm;
        !           633:       newxnorm := OneNorm(x);
        !           634:       if abs(newxnorm - oldxnorm) <= ROUGHLYZERO 
        !           635:          then closenough := TRUE;
        !           636:       end;
        !           637:       if itcount < 4 then closenough := FALSE;
        !           638:    until (itcount = 25) or closenough;
        !           639:    PrintReport;
        !           640:    end;
        !           641: end.
        !           642: 

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.