DataMuseum.dk

top - metrics - download

    Length: 3072 (0xc00)
    Types: TextFile
    Names: »gausstxt«

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦this⟧ »gausstxt« 
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
    └─⟦4334b4c0b⟧ 
        └─⟦this⟧ »gausstxt« 

TextFile

detgauss=set 5
permanent detgauss.5
detgauss=algol
procedure for solving linear equations, gauss.
07 03 73 19 00 00
external


real procedure Det_Gauss(n,m,A,exit);
value n,m;
integer n,m;
array A;
label exit;
comment Det Gauss solves m systems of linear equations with n
        unknowns and computes the determinant of the coeffi-
        cient matrix. The parameters are:
        
        n    : Number of unknowns.

        m    : Number of right hand sides. m=0 implies that
               only the determinant will be computed (cf. section 2).

        A    : Has the dimension (1:n,1:n+m) and must on entry
               contain the coefficient matrix in A(1:n,1:n) and
               the right hand sides in A(1:n,n+1:n+m). The solu-
               tion vectors are on exit stored as columns in
               A(1:n,n+1:n+m).

        exit : Label to which Det Gauss goes when the system is
               singular;

begin integer i,j,k,i0,j0;
  real factor,max,detfac,twofac;
  integer array permute(1:n);
  m:= n+m;
  detfac:= 1;
  for i:= 1 step 1 until n do
  begin
    max:= 0;
    for j:= 1 step 1 until n do
    max:= max+ A(i,j)**2;
    if max>1 or max<.25 then
    begin
      twofac:= 2**(-entier(ln(max)/1.3863+1));
      for j:= 1 step 1 until m do
      A(i,j) := A(i,j)*twofac;
      detfac:= detfac/twofac
    end
  end equilibration;
  for k:= 1 step 1 until n do
  begin
    max:= 0;
    for i:= k step 1 until n do
    for j:= k step 1 until n do
    begin
      factor:= abs(A(i,j));
      if max <= factor then
      begin
        max:= factor;
        i0:= i;
        j0:= j
      end
    end searching_for_pivotal_element;
    if max <= '-8 then goto exit;
    max:= A(i0,j0);detfac:= detfac*max;
    if i0<>k then
    begin
      detfac:= -detfac;
      for j:= k step 1 until m do
      begin
        factor:= A(k,j);
        A(k,j) :=A(i0,j);
        A(i0,j) := factor
      end
    end interchange_of_rows;
    permute(k) := k;
    if j0<>k then
    begin
      detfac:= -detfac;
      permute(k) := j0;
      for i:= 1 step 1 until n do
      begin
        factor:= A(i,k);
        A(i,k) := A(i,j0);
        A(i,j0) := factor
      end
    end interchange_of_columns;
    for i:= k+1 step 1 until n do
    begin
      factor:= A(i,k)/max;
      for j:= k+1 step 1 until m do
      A(i,j) := A(i,j) - A(k,j)*factor
    end reduction
  end for_k;
  for k:= n+1 step 1 until m do
  for i:= n step -1 until 1 do
  begin
    factor:= A(i,k);
    for j:= i+1 step 1 until n do
    factor:= factor- A(i,j)*A(j,k);
    A(i,k) := factor/A(i,i)
  end solving;
  if m<>n then
  begin
    for i:= n-1 step -1 until 1 do
    begin
      i0:= permute(i);
      if i0<>i then
      for k:= n+1 step 1 until m do
      begin
        factor:= A(i,k);
        A(i,k) := A(i0,k);
        A(i0,k) := factor
      end
    end
  end interchange_of_solutions;
  Det_Gauss:= detfac
end;
end Det_Gauss;          \f

\f



\f

▶EOF◀