home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / NUMPAS.ZIP / CMATINV.PAS next >
Encoding:
Pascal/Delphi Source File  |  1985-06-24  |  5.7 KB  |  216 lines
  1. program cmatinv(input, output);
  2. (* Program calculates the inverse of complex matrix A by gauss elimination, *)
  3. (* 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,f: mat;
  8.     n: integer;
  9.     q: char;
  10.  
  11. procedure cmatread(var a1: mat; m1,n1: integer);
  12. (* Reads in coefficients of n-order matrix A1 from the keyboard and *)
  13. (* forms the 2n-order equivalent real matrix. *)
  14.  
  15. var aa,ab: real;
  16.     i,j: integer;
  17.  
  18. begin
  19.    writeln('Enter matrix as directed by prompt.');
  20.    writeln;
  21.    for i := 1 to m1 do
  22.       begin
  23.          for j := 1 to n1 do
  24.             begin
  25.                writeln;
  26.                writeln('Enter complex ',i,'-',j, ' element: ');
  27.                readln(aa,ab);
  28.                a1[2*i-1,2*j-1] := aa;
  29.                a1[2*i-1,2*j] := ab;
  30.                a1[2*i,2*j-1] := -ab;
  31.                a1[2*i,2*j] := aa;
  32.             end;
  33.          writeln;
  34.       end;
  35. end;
  36.  
  37. procedure cmatwrite(a1: mat; m1,n1: integer);
  38. (* Prints out the complex matrix A1 from the equivalent real matrix. *)
  39.  
  40. var i,j: integer;
  41.  
  42. begin
  43.    writeln;
  44.    for i := 1 to m1 do
  45.       begin
  46.          for j := 1 to n1 do
  47.             begin
  48.                write('(',a1[2*i-1,2*j-1]:9,', ',a1[2*i-1,2*j]:9,')','  ');
  49.             end;
  50.          writeln;
  51.       end;
  52.    writeln;
  53. end;
  54.  
  55.  
  56. procedure matinv(a1: mat; var b2: mat; n1: integer);
  57. (* Calculates the inverse of matrix A1 and designates it as B2. *)
  58.  
  59. label 10,20,30;
  60. var d: real;
  61.     b1: mat1;
  62.     i,j,k,m: integer;
  63.  
  64. begin  (* Set up a matrix [A U] for Gauss elimination procedure *)
  65.    for i := 1 to n1 do
  66.       begin
  67.          for j := 1 to n1 do
  68.             begin
  69.                b1[i,n1+j] := 0;
  70.                b1[i,j] := a1[i,j];
  71.             end;
  72.          b1[i,n1+i] := 1;
  73.       end;
  74.    for k := 1 to n1 do  (* Find the largest element in the k column *)
  75.       begin
  76.          if k = n1 then goto 10;
  77.          m := k;
  78.          for i := k+1 to n1 do
  79.             begin
  80.                if abs(b1[i,k]) > abs(b1[m,k]) then
  81.                m := i;
  82.             end;
  83.          if m = k then goto 10;
  84.          for j := k to 2*n1 do
  85.             begin  (* Exchange k and m rows. *)
  86.                d := b1[k,j];
  87.                b1[k,j] := b1[m,j];
  88.                b1[m,j] := d;
  89.             end;
  90.       10:for j := k+1 to 2*n1 do
  91.             begin  (* Divide row k by the k-k element. *)
  92.                b1[k,j] := b1[k,j]/b1[k,k];
  93.             end;
  94.          if k = 1 then goto 20;
  95.          for i := 1 to k-1 do  (* Clear out above main diagonal *)
  96.             begin
  97.                for j := k+1 to 2*n1 do
  98.                   begin
  99.                      b1[i,j] := b1[i,j] - b1[i,k]*b1[k,j];
  100.                   end;
  101.             end;
  102.          if k = n1 then goto 30;
  103.       20:for i := k+1 to n1 do
  104.             begin
  105.                for j := k+1 to 2*n1 do  (* Clear out below main diagonal *)
  106.                   begin
  107.                      b1[i,j] := b1[i,j] - b1[i,k]*b1[k,j];
  108.                   end;
  109.             end;
  110.          end;
  111. 30:for i := 1 to n1 do    (* Copy A-inverse to the left side of B. *)
  112.       begin
  113.          for j := 1 to n1 do
  114.             begin
  115.                b2[i,j] := b1[i,n1+j];
  116.             end;
  117.       end;
  118. end;
  119.  
  120. procedure matmult(a1,b1: mat; var c1: mat; m1,n1,n2: integer);
  121. (* Calculates the matrix product [C1] = [A1][B1]. *)
  122.  
  123. var i,j,k: integer;
  124.  
  125. begin
  126.    for i := 1 to m1 do
  127.       begin
  128.          for j := 1 to n2 do
  129.             begin
  130.                c1[i,j] := 0.0;
  131.                for k := 1 to n1 do
  132.                   begin
  133.                      c1[i,j] := c1[i,j] + a1[i,k]*b1[k,j];
  134.                   end;
  135.             end;
  136.       end;
  137. end;
  138.  
  139.  
  140.  
  141. procedure cmatpol(a1: mat; m1: integer);
  142. (* Calculates the polar form of a complex vector in rectangular form, *)
  143. (* and prints out the results. *)
  144.  
  145. var c11,c2,cmag1,angle1,eta: real;
  146.     j,i: integer;
  147.  
  148. begin
  149.    for i := 1 to m1 do
  150.    begin
  151.       j := 2*i-1;
  152.       c11 := a1[j,1];
  153.       c2 := a1[j,2];
  154.       cmag1 := sqrt(c11*c11 + c2*c2);
  155.       if (c11 = 0.0) and (c2 = 0.0) then
  156.          angle1 := 0.0;
  157.       if (c11 = 0.0) and (c2 > 0.0) then
  158.          angle1 := 90.0;
  159.       if (c11 = 0.0) and (c2 < 0.0) then
  160.          angle1 := -90.0;
  161.       if (abs(c11) > 0.0) then
  162.          eta := arctan(c2/c11)*57.296;
  163.       if (c11 > 0) then
  164.          angle1 := eta;
  165.       if (c11 < 0.0) and (c2 >= 0.0) then
  166.          angle1 := eta + 180.0;
  167.       if (c11 < 0.0) and (c2 < 0.0) then
  168.          angle1 := eta - 180.0;
  169.       writeln(cmag1:9,'  ',angle1:5:1);
  170.    end;
  171. end;
  172.  
  173.  
  174. begin  (* main *)
  175.    writeln('This program finds the inverse of a complex matrix.');
  176.    writeln;
  177.    write('Input the order of the square complex matrix A: ');
  178.    readln(n);
  179.    writeln;
  180.    writeln('The elements of A will next be input: ');
  181.    cmatread(a,n,n);
  182.    writeln('The matrix A is: ');
  183.    writeln;
  184.    cmatwrite(a,n,n);
  185.    writeln;
  186.    matinv(a,b,2*n);
  187.    writeln('The inverse of A is: ');
  188.    writeln;
  189.    cmatwrite(b,n,n);
  190.    matmult(a,b,c,2*n,2*n,2*n);
  191.    writeln;
  192.    writeln('The product of A and A-inverse is:');
  193.    writeln;
  194.    cmatwrite(c,n,n);
  195.    writeln;
  196.    writeln('Solve Ax = f ?  (y,n)');
  197.    readln(q);
  198.    if q = 'y' then
  199.       begin
  200.          writeln('Input the vector f as directed by prompt:');
  201.          writeln;
  202.          cmatread(f,n,1);
  203.          matmult(b,f,c,2*n,2*n,2);
  204.          writeln;
  205.          writeln('The complex vector solution is:');
  206.          writeln;
  207.          writeln('In rectangular form:');
  208.          writeln;
  209.          cmatwrite(c,n,1);
  210.          writeln;
  211.          writeln('In polar form (magnitude, and angle in degrees):');
  212.          writeln;
  213.          cmatpol(c,n);
  214.       end;
  215. end.
  216.