[comp.os.vms] GAMMA.PAS

ksimpson%sfsu1.hepnet@LBL.GOV (Kenneth Simpson) (05/28/88)

PROGRAM FixedPointGamma(input,output);

	{ Evaluates the zero of a polynominal by the fixed point method.}

        CONST
          zero = 1E-10;
          maxdim = 1000;

        TYPE
          states  = (start,iterating,done,exceeded_maxit,divisor_zero);
          dim = 0..maxdim;
          vector  = ARRAY[dim] OF double;
          ivector = ARRAY[dim] OF integer;
 
	VAR
          in_file,out_file           : text; 
          y,xinit,xroot              : vector;
          count                      : ivector;
          tolerence                  : double;
          max_iter                : integer;
          output_dat                 : boolean;
          status                     : states;
     
{==========================================================================}

PROCEDURE Open_IO;

  { Open files.} 

  BEGIN

    open ( in_file ,  'FP.DAT', OLD );{ Contains initialization data } 
    open ( out_file , 'FPD.DAT' );    { Output file.               } 
    reset( in_file );
    rewrite( out_file );
  END { Open_IO };


{============================================================================}

PROCEDURE Close_IO;

 { Close all data files.}

  BEGIN
    close ( in_file );
    close ( out_file );
  END { Close_IO };

{============================================================================}


PROCEDURE Initialize ( VAR xinit       : Vector;{ x is initial data point. }   
                       VAR tolerence   : double;  { tolerence is on root.} 
                       VAR max_iter    : integer; { maxint maximum iterations.}
                       VAR output_dat  : boolean );

VAR
  i   : INTEGER;

  { Reads initial x, tolerence of root and maximum number of iterations.}

  BEGIN
    FOR i := 1 TO 2 DO
      readln( in_file, xinit[i] );
    readln( in_file, tolerence );
    readln( in_file, max_iter );
    readln( in_file, output_dat );
  END { Initialize };

{=========================================================================}

PROCEDURE Error_Messages;

VAR 
   i  : INTEGER;

  BEGIN
     writeln( out_file );
    CASE status OF
      divisor_zero:
        BEGIN
          writeln( out_file,' ERROR detected:');
          writeln( out_file );
          writeln( out_file,' ===> DIVISION BY ZERO <===');
          writeln( out_file );
          writeln( out_file, ' ZERO defined to be =',zero:20:14);
          writeln( out_file);
        END;
      exceeded_maxit:
        BEGIN
          writeln( out_file,' ERROR detected:');
          writeln( out_file );
          write( out_file,' ===> MAX ITERATIONS OF ',max_iter:4);
          write( out_file,' EXCEEDED <===');
          writeln( out_file );
        END;
      done: 
        BEGIN
          writeln( out_file,' DONE. No errors detected.');
          writeln( out_file);
          FOR i:=1 TO 2 DO
            BEGIN
              writeln( out_file,' Root x[',i:2,'] =',xroot[i]:20:14);
              writeln( out_file,' Iterations  =',count[i]:6);
            END;
        END;
    END { CASE }
      
  END { Error_Messages };

{=========================================================================}

PROCEDURE Print_Header;

 { Prints header information and in out_file data file. }

VAR 
  i   : INTEGER;

 BEGIN             
   writeln(out_file );
   FOR i := 1 TO 2 DO
     writeln(out_file ,' Initial value for x[',i:2,'] = ',xinit[i]:20:14);
   writeln(out_file ,' Tolerance = ',tolerence:16:14 );
   writeln(out_file ,' Maximum possible iterations = ',max_iter:6);
   IF output_dat THEN
     BEGIN
       writeln(out_file );
       write(out_file,'======================================================');
       writeln(out_file, '================');
       write(out_file ,'   N     ROOT        Xn               G(Xn)           ');
       writeln(out_file,'  |Xn-G(Xn)|');
       write(out_file,'======================================================');
       writeln(out_file, '================');
       writeln(out_file );
     END
   ELSE
     BEGIN
       writeln(out_file,' ( ===> Data output shut off. <=== )');
       writeln(out_file);
     END;
 END { PrintHeader };


{============================================================================}

  PROCEDURE Gamma( VAR y : Vector;
                       x : Vector ); 
                      
 
  VAR
     sum  : double;
 
  BEGIN
      sum  := x[1] + x[2];
      y[1] := sin(sum);
      y[2] := cos(sum);
  END;


{==========================================================================}

PROCEDURE MVFP;                  

 { Finds roots of a vector valued function by fixed point method.}

	VAR
          x,xnext,differ     : Vector;
          i                  : INTEGER;          
                   
  BEGIN
    status:= start;
    FOR i := 1 TO 2 DO
      BEGIN
        x[i] :=xinit[i]; 
        count[i] := 0;
      END;
    Print_Header;
    REPEAT
      BEGIN
        Gamma(y,x);
        FOR i := 1 TO 2 DO
          BEGIN
            xnext[i]  := y[i];
            differ[i] := xnext[i]-x[i];
            IF ( status <> divisor_zero ) THEN
              IF ( count[i] >= max_iter ) THEN
                status:= exceeded_maxit
              ELSE IF ( abs(differ[i]) < tolerence ) THEN
                BEGIN
                  status:= done;
                  xroot[i] := xnext[i];
                  count[i] := count[i] + 1;
                  IF ( output_dat ) THEN
                    writeln(out_file,count[i]:4,i:6,x[i]:20:14,
                            xnext[i]:20:14,abs(differ[i]):20:14);
                END
              ELSE
                BEGIN
                  status := iterating;
                  count[i] := count[i] + 1;
                  IF ( output_dat ) THEN
                    writeln(out_file,count[i]:4,i:6,x[i]:20:14,
                                     xnext[i]:20:14,abs(differ[i]):20:14);
                  x[i] := xnext[i];
                END;
            END;
      END;
    UNTIL (status <> iterating);  
    Error_Messages;
  END { Fix_Point };

{==========================================================================}
{**************************************************************************}
{==========================================================================}

BEGIN { Program body. }
   Open_IO;
   Initialize( xinit, tolerence, max_iter, output_dat );
   MVFP;
   Close_IO
END.