home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / NUMPAS.ZIP / RMATINV.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-06-20  |  4.2 KB  |  164 lines
  1. program rmatinv(input, output);
  2. (* Program calculates the inverse of real matrix A by gaussian *)
  3. (* elimination and solves the matrix equation Ax = f if desired *)
  4.  
  5. type mat = array[1..10,1..10] of real;
  6.      mat1 = array[1..10,1..20] of real;
  7. var a,b,c,d,f: mat;
  8.     q: char;
  9.     n: integer;
  10.  
  11. procedure matread(var a1: mat; m1,n1: integer);
  12. (* Reads in the coefficients of matrix A1 *)
  13.  
  14. var i,j: integer;
  15. begin
  16.    for i := 1 to m1 do
  17.       begin
  18.          for j := 1 to n1 do
  19.             begin
  20.                write('Input the ',j,' element of row ',i,':  ');
  21.                readln(a1[i,j]);
  22.             end;
  23.          writeln;
  24.       end;
  25. end;
  26.  
  27. procedure matwrite(a1: mat; m1,n1: integer);
  28. (* Prints out matrix A1 *)
  29.  
  30. var i,j: integer;
  31. begin
  32.    writeln;
  33.    for i := 1 to m1 do
  34.       begin
  35.          for j := 1 to n1 do
  36.             begin
  37.                write(a1[i,j]:10,' ':2);
  38.             end;
  39.          writeln;
  40.       end;
  41. end;
  42.  
  43.  
  44. procedure matinv(a1: mat; var b2: mat; n1: integer);
  45. (* Calculates the inverse of matrix A1 and designates it as matrix C *)
  46.  
  47. label 10,20,30;
  48. var d: real;
  49.     b1: mat1;
  50.     i,j,k,m: integer;
  51.  
  52. begin  (* Set up a matrix [A U] for Gauss elimination procedure *)
  53.    for i := 1 to n1 do
  54.       begin
  55.          for j := 1 to n1 do
  56.             begin
  57.                b1[i,n1+j] := 0;
  58.                b1[i,j] := a1[i,j];
  59.             end;
  60.          b1[i,n1+i] := 1;
  61.       end;
  62.    for k := 1 to n1 do  (* Find the largest element in the k column *)
  63.       begin
  64.          if k = n1 then goto 10;
  65.          m := k;
  66.          for i := k+1 to n1 do
  67.             begin
  68.                if abs(b1[i,k]) > abs(b1[m,k]) then
  69.                m := i;
  70.             end;
  71.          if m = k then goto 10;
  72.          for j := k to 2*n1 do
  73.             begin  (* Exchange k and m rows. *)
  74.                d := b1[k,j];
  75.                b1[k,j] := b1[m,j];
  76.                b1[m,j] := d;
  77.             end;
  78.       10:for j := k+1 to 2*n1 do
  79.             begin  (* Divide row k by the k-k element. *)
  80.                b1[k,j] := b1[k,j]/b1[k,k];
  81.             end;
  82.          if k = 1 then goto 20;
  83.          for i := 1 to k-1 do  (* Clear out above main diagonal *)
  84.             begin
  85.                for j := k+1 to 2*n1 do
  86.                   begin
  87.                      b1[i,j] := b1[i,j] - b1[i,k]*b1[k,j];
  88.                   end;
  89.             end;
  90.          if k = n1 then goto 30;
  91.       20:for i := k+1 to n1 do
  92.             begin
  93.                for j := k+1 to 2*n1 do  (* Clear out below main diagonal *)
  94.                   begin
  95.                      b1[i,j] := b1[i,j] - b1[i,k]*b1[k,j];
  96.                   end;
  97.             end;
  98.          end;
  99. 30:for i := 1 to n1 do    (* Copy A-inverse to the left side of B. *)
  100.       begin
  101.          for j := 1 to n1 do
  102.             begin
  103.                b2[i,j] := b1[i,n1+j];
  104.             end;
  105.       end;
  106. end;
  107.  
  108. procedure matmult(a1,b1: mat; var c1: mat; m1,n1,n2: integer);
  109. (* Calculates the matrix product [A1][B1] = C1 *)
  110.  
  111. var i,j,k: integer;
  112.  
  113. begin
  114.    for i := 1 to m1 do
  115.       begin
  116.          for j := 1 to n2 do
  117.             begin
  118.                c1[i,j] := 0.0;
  119.                for k := 1 to n1 do
  120.                   begin
  121.                      c1[i,j] := c1[i,j] + a1[i,k]*b1[k,j];
  122.                   end;
  123.             end;
  124.       end;
  125. end;
  126.  
  127.  
  128. begin  (* main *)
  129.    writeln('This program calculates the inverse of a matrix.');
  130.    writeln;
  131.    write('Input the order of the square matrix A: ');
  132.    readln(n);
  133.    writeln;
  134.    writeln('Input the elements of A as prompted: ');
  135.    matread(a,n,n);
  136.    writeln('The matrix A is: ');
  137.    writeln;
  138.    matwrite(a,n,n);
  139.    writeln;
  140.    matinv(a,b,n);
  141.    writeln('The inverse of A is: ');
  142.    writeln;
  143.    matwrite(b,n,n);
  144.    matmult(a,b,c,n,n,n);
  145.    writeln;
  146.    writeln('The product of A and A-inverse is:');
  147.    writeln;
  148.    matwrite(c,n,n);
  149.    writeln;
  150.    writeln('Solve Ax = f? (y,n) ');
  151.    readln(q);
  152.    if q = 'y' then
  153.       begin
  154.          writeln('Input the elements of vector f as directed by prompt:');
  155.          writeln;
  156.          matread(f,n,1);
  157.          writeln;
  158.          writeln('The solution vector is:');
  159.          writeln;
  160.          matmult(b,f,d,n,n,1);
  161.          matwrite(d,n,1);
  162.       end;
  163. end.
  164.