Gaussian elimination: Difference between revisions

From Rosetta Code
Content added Content deleted
(+ D entry)
(No more a draft task)
Line 1: Line 1:
{{draft task}} [[Category:Matrices]]
{{task}}[[Category:Matrices]]
Problem: Solve Ax=b using Gaussian elimination then backwards substitution. A being an n by n matrix. Also, x and b are n by 1 vectors. To improve accuracy, please use partial pivoting and scaling.
Problem: Solve Ax=b using Gaussian elimination then backwards substitution. A being an n by n matrix. Also, x and b are n by 1 vectors. To improve accuracy, please use partial pivoting and scaling.


Revision as of 10:23, 1 January 2014

Task
Gaussian elimination
You are encouraged to solve this task according to the task description, using any language you may know.

Problem: Solve Ax=b using Gaussian elimination then backwards substitution. A being an n by n matrix. Also, x and b are n by 1 vectors. To improve accuracy, please use partial pivoting and scaling.

ALGOL 68

Works with: ALGOL 68 version Revision 1 - extension to language used - "PRAGMA READ" (similar to C's #include directive.)
Works with: ALGOL 68G version Any - tested with release algol68g-2.4.1.

File: prelude_exception.a68<lang algol68># -*- coding: utf-8 -*- # COMMENT PROVIDES 

 MODE FIXED; INT fixed exception, unfixed exception;
 PROC (STRING message) FIXED raise, raise value error

END COMMENT

  1. Note: ℵ indicates attribute is "private", and
       should not be used outside of this prelude #

MODE FIXED = BOOL; # if an exception is detected, can it be fixed "on-site"? # FIXED fixed exception = TRUE, unfixed exception = FALSE;

MODE #ℵ#SIMPLEOUTV = [0]UNION(CHAR, STRING, INT, REAL, BOOL, BITS); MODE #ℵ#SIMPLEOUTM = [0]#ℵ#SIMPLEOUTV; MODE #ℵ#SIMPLEOUTT = [0]#ℵ#SIMPLEOUTM; MODE SIMPLEOUT = [0]#ℵ#SIMPLEOUTT;

PROC raise = (#ℵ#SIMPLEOUT message)FIXED: ( 

 putf(stand error, ($"Exception:"$, $xg$, message, $l$));
 stop

);

PROC raise value error = (#ℵ#SIMPLEOUT message)FIXED: 

 IF raise(message) NE fixed exception THEN exception value error; FALSE FI;

SKIP</lang>File: prelude_mat_lib.a68<lang algol68># -*- coding: utf-8 -*- # COMMENT PRELUDE REQUIRES 

 MODE SCAL = REAL;
 FORMAT scal repr = real repr
 # and various SCAL OPerators #

END COMMENT

COMMENT PRELUDE PROIVIDES 

 MODE VEC, MAT;
 OP :=:, -:=, +:=, *:=, /:=;
 FORMAT sub, sep, bus;
 FORMAT vec repr, mat repr

END COMMENT

  1. Note: ℵ indicates attribute is "private", and
       should not be used outside of this prelude #

INT #ℵ#lwb vec := 1, #ℵ#upb vec := 0; INT #ℵ#lwb mat := 1, #ℵ#upb mat := 0; MODE VEC = [lwb vec:upb vec]SCAL, 

    MAT = [lwb mat:upb mat,lwb vec:upb vec]SCAL;

FORMAT sub := $"( "$, sep := $", "$, bus := $")"$, nl:=$lxx$; FORMAT vec repr := $f(sub)n(upb vec - lwb vec)(f(scal repr)f(sep))f(scal repr)f(bus)$; FORMAT mat repr := $f(sub)n(upb mat - lwb mat)(f( vec repr)f(nl))f( vec repr)f(bus)$;

  1. OPerators to swap the contents of two VECtors #

PRIO =:= = 1; OP =:= = (REF VEC u, v)VOID: 

 FOR i TO UPB u DO SCAL scal=u[i]; u[i]:=v[i]; v[i]:=scal OD;

OP +:= = (REF VEC lhs, VEC rhs)REF VEC: ( 

 FOR i TO UPB lhs DO lhs[i] +:= rhs[i] OD;
 lhs

);

OP -:= = (REF VEC lhs, VEC rhs)REF VEC: ( 

 FOR i TO UPB lhs DO lhs[i] -:= rhs[i] OD;
 lhs

);

OP *:= = (REF VEC lhs, SCAL rhs)REF VEC: ( 

 FOR i TO UPB lhs DO lhs[i] *:= rhs OD;
 lhs

);

OP /:= = (REF VEC lhs, SCAL rhs)REF VEC: ( 

 SCAL inv = 1 / rhs; # multiplication is faster #
 FOR i TO UPB lhs DO lhs[i] *:= inv OD;
 lhs

);

SKIP</lang>File: prelude_gaussian_elimination.a68<lang algol68># -*- coding: utf-8 -*- # COMMENT PRELUDE REQUIRES 

 MODE SCAL = REAL,
 REAL near min scal = min real ** 0.99,
 MODE VEC = []REAL,
 MODE MAT = [,]REAL,
 FORMAT scal repr = real repr,
 and various OPerators of MAT and VEC

END COMMENT

COMMENT PRELUDE PROVIDES 

 PROC(MAT a, b)MAT gaussian elimination;
 PROC(REF MAT a, b)REF MAT in situ gaussian elimination

END COMMENT

  2. using Gaussian elimination, find x where A*x = b #

PROC in situ gaussian elimination = (REF MAT a, b)REF MAT: (

  1. Note: a and b are modified "in situ", and b is returned as x #
 FOR diag TO UPB a-1 DO
   INT pivot row := diag; SCAL pivot factor := ABS a[diag,diag];
   FOR row FROM diag + 1 TO UPB a DO # Full pivoting #
     SCAL abs a diag = ABS a[row,diag];
     IF abs a diag>=pivot factor THEN 
       pivot row := row; pivot factor := abs a diag FI
   OD;
 # now we have the "best" diag to full pivot, do the actual pivot #
   IF diag NE pivot row THEN
   # a[pivot row,] =:= a[diag,]; XXX: unoptimised # #DB#
     a[pivot row,diag:] =:= a[diag,diag:]; # XXX: optimised #
     b[pivot row,] =:= b[diag,] # swap/pivot the diags of a & b #
   FI;

   IF ABS a[diag,diag] <= near min scal THEN 
     raise value error("singular matrix") FI;
   SCAL a diag reciprocal := 1 / a[diag, diag];

   FOR row FROM diag+1 TO UPB a DO
     SCAL factor = a[row,diag] * a diag reciprocal;
   # a[row,] -:= factor * a[diag,] XXX: "unoptimised" # #DB#
     a[row,diag+1:] -:= factor * a[diag,diag+1:];# XXX: "optimised" #
     b[row,] -:= factor * b[diag,]
   OD
 OD;
  1. We have a triangular matrix, at this point we can traverse backwards
 up the diagonal calculating b\A Converting it initial to a diagonal
 matrix, then to the identity.  #

 FOR diag FROM UPB a BY -1 TO 1+LWB a DO

   IF ABS a[diag,diag] <= near min scal THEN 
     raise value error("Zero pivot encountered?") FI;
   SCAL a diag reciprocal = 1 / a[diag,diag];

   FOR row TO diag-1 DO
     SCAL factor = a[row,diag] * a diag reciprocal;
   # a[row,diag] -:= factor * a[diag,diag]; XXX: "unoptimised" so remove # #DB#
     b[row,] -:= factor * b[diag,]
   OD;
  1. Now we have only diagonal elements we can simply divide b
 by the values along the diagonal of A. #
   b[diag,] *:= a diag reciprocal
 OD;

 b # EXIT #

);

PROC gaussian elimination = (MAT in a, in b)MAT: (

  1. Note: a and b are cloned and not modified "in situ" #
 [UPB in a, 2 UPB in a]SCAL a := in a;
 [UPB in b, 2 UPB in b]SCAL b := in b;
 in situ gaussian elimination(a,b)

);

SKIP</lang>File: postlude_exception.a68<lang algol68># -*- coding: utf-8 -*- # COMMENT POSTLUDE PROIVIDES 

 PROC VOID exception too many iterations, exception value error;

END COMMENT

SKIP EXIT exception too many iterations: exception value error: 

 stop</lang>File: test_Gaussian_elimination.a68<lang algol68>#!/usr/bin/algol68g-full --script #
  1. -*- coding: utf-8 -*- #

PR READ "prelude_exception.a68" PR;

  1. define the attributes of the scalar field being used #

MODE SCAL = REAL; FORMAT scal repr = $g(-0,real width)$;

  1. create "near min scal" as is scales better then small real #

SCAL near min scal = min real ** 0.99;

PR READ "prelude_mat_lib.a68" PR; PR READ "prelude_gaussian_elimination.a68" PR;

MAT a =(( 1.00, 0.00, 0.00, 0.00, 0.00, 0.00), 

       ( 1.00, 0.63, 0.39,  0.25,  0.16,   0.10),
       ( 1.00, 1.26, 1.58,  1.98,  2.49,   3.13),
       ( 1.00, 1.88, 3.55,  6.70, 12.62,  23.80),
       ( 1.00, 2.51, 6.32, 15.88, 39.90, 100.28),
       ( 1.00, 3.14, 9.87, 31.01, 97.41, 306.02));

VEC b = (-0.01, 0.61, 0.91, 0.99, 0.60, 0.02);

[UPB b,1]SCAL col b; col b[,1]:= b;

upb vec := 2 UPB a;

printf((vec repr, gaussian elimination(a,col b)));

PR READ "postlude_exception.a68" PR</lang>Output: 

( -.010000000000002, 1.602790394502130, -1.613203059905640, 1.245494121371510, -.490989719584686, .065760696175236)

C

This modifies A and b in place, which might not be quite desirable. <lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <math.h>
  1. define mat_elem(a, y, x, n) (a + ((y) * (n) + (x)))

void swap_row(double *a, double *b, int r1, int r2, int n) { double tmp, *p1, *p2; int i;

if (r1 == r2) return; for (i = 0; i < n; i++) { p1 = mat_elem(a, r1, i, n); p2 = mat_elem(a, r2, i, n); tmp = *p1, *p1 = *p2, *p2 = tmp; } tmp = b[r1], b[r1] = b[r2], b[r2] = tmp; }

void gauss_eliminate(double *a, double *b, double *x, int n) {

  1. define A(y, x) (*mat_elem(a, y, x, n))

int i, j, col; double max, tmp;

for (col = 0; col < n; col++) { j = col, max = A(j, j);

for (i = col + 1; i < n; i++) if ((tmp = fabs(A(i, col))) > max) j = i, max = tmp;

swap_row(a, b, col, j, n);

for (i = col + 1; i < n; i++) { tmp = A(i, col) / A(col, col); for (j = col + 1; j < n; j++) A(i, j) -= tmp * A(col, j); A(i, col) = 0; b[i] -= tmp * b[col]; } } for (col = n - 1; col >= 0; col--) { tmp = b[col]; for (j = n - 1; j > col; j--) tmp -= x[j] * A(col, j); x[col] = tmp / A(col, col); }

  1. undef A

}

int main(void) { double a[] = { 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.63, 0.39, 0.25, 0.16, 0.10, 1.00, 1.26, 1.58, 1.98, 2.49, 3.13, 1.00, 1.88, 3.55, 6.70, 12.62, 23.80, 1.00, 2.51, 6.32, 15.88, 39.90, 100.28, 1.00, 3.14, 9.87, 31.01, 97.41, 306.02 }; double b[] = { -0.01, 0.61, 0.91, 0.99, 0.60, 0.02 }; double x[6]; int i;

gauss_eliminate(a, b, x, 6);

for (i = 0; i < 6; i++) printf("%g\n", x[i]);

return 0;

}</lang>

Output:


-0.01
1.60279
-1.6132
1.24549
-0.49099
0.0657607


D

Translation of: Go

<lang d>import std.stdio, std.math, std.algorithm, std.range, std.numeric;

Tuple!(double[],"x", string,"err") gaussPartial(in double[][] a0, in double[] b0) /*pure nothrow*/ in { 

   assert(a0.length == a0[0].length);
   assert(a0.length == b0.length);
   assert(a0.all!(row => row.length == a0[0].length));

} body { 

   enum eps = 1e-6;
   immutable m = b0.length;

   // Make augmented matrix.
   auto a = a0.zip(b0).map!(c => cast(double[])(c[0] ~ c[1])).array;

   // Wikipedia algorithm from Gaussian elimination page,
   // produces row-eschelon form.
   foreach (immutable k; 0 .. a.length) {
       // Find pivot for column k and swap.
       a[k .. m].minPos!((a, b) => a[k] > b[k]).front.swap(a[k]);
       if (a[k][k].abs < eps)
           return typeof(return)(null, "singular");

       // Do for all rows below pivot.
       foreach (immutable i; k + 1 .. m) {
           // Do for all remaining elements in current row.
           a[i][k+1 .. m+1] -= a[k][k+1 .. m+1] * (a[i][k] / a[k][k]);

           a[i][k] = 0; // Fill lower triangular matrix with zeros.
       }
   }

   // End of WP algorithm. Now back substitute to get result.
   auto x = new double[m];
   foreach_reverse (immutable i; 0 .. m)
       x[i] = (a[i][m] - a[i][i+1 .. m].dotProduct(x[i+1 .. m]))
              / a[i][i];

   return typeof(return)(x, null);

}

void main() { 

   // The test case result is correct to this tolerance.
   enum eps = 1e-13;

   // Common RC example.  Result computed with rational arithmetic
   // then converted to double, and so should be about as close to
   // correct as double represention allows.
   immutable a = [[1.00, 0.00, 0.00, 0.00, 0.00, 0.00],
                  [1.00, 0.63, 0.39, 0.25, 0.16, 0.10],
                  [1.00, 1.26, 1.58, 1.98, 2.49, 3.13],
                  [1.00, 1.88, 3.55, 6.70, 12.62, 23.80],
                  [1.00, 2.51, 6.32, 15.88, 39.90, 100.28],
                  [1.00, 3.14, 9.87, 31.01, 97.41, 306.02]];
   immutable b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02];

   /*immutable*/ const r = gaussPartial(a, b);
   if (!r.err.empty)
       return writeln("Error: ", r.err);
   r.x.writeln;

   immutable result = [-0.01, 1.602790394502114, -1.6132030599055613,
       1.2454941213714368, -0.4909897195846576, 0.065760696175232];
   foreach (immutable i, immutable xi; result)
       if (abs(r.x[i] - xi) > eps)
           return writeln("Out of tolerance: ", r.x[i], " ", xi);

}</lang>

Output:
[-0.01, 1.60279, -1.6132, 1.24549, -0.49099, 0.0657607]

Go

Gaussian elimination with partial pivoting by pseudocode on WP page "Gaussian elimination." Scaling not implemented. <lang go>package main

import ( 

   "errors"
   "fmt"
   "log"
   "math"

)

type testCase struct { 

   a [][]float64
   b []float64
   x []float64

}

var tc = testCase{ 

   // common RC example.  Result x computed with rational arithmetic then
   // converted to float64, and so should be about as close to correct as
   // float64 represention allows.
   a: [][]float64{
       {1.00, 0.00, 0.00, 0.00, 0.00, 0.00},
       {1.00, 0.63, 0.39, 0.25, 0.16, 0.10},
       {1.00, 1.26, 1.58, 1.98, 2.49, 3.13},
       {1.00, 1.88, 3.55, 6.70, 12.62, 23.80},
       {1.00, 2.51, 6.32, 15.88, 39.90, 100.28},
       {1.00, 3.14, 9.87, 31.01, 97.41, 306.02}},
   b: []float64{-0.01, 0.61, 0.91, 0.99, 0.60, 0.02},
   x: []float64{-0.01, 1.602790394502114, -1.6132030599055613,
       1.2454941213714368, -0.4909897195846576, 0.065760696175232},

}

// result from above test case turns out to be correct to this tolerance. const ε = 1e-13

func main() { 

   x, err := GaussPartial(tc.a, tc.b)
   if err != nil {
       log.Fatal(err)
   }
   fmt.Println(x)
   for i, xi := range x {
       if math.Abs(tc.x[i]-xi) > ε {
           log.Println("out of tolerance")
           log.Fatal("expected", tc.x)
       }
   }

}

func GaussPartial(a0 [][]float64, b0 []float64) ([]float64, error) { 

   // make augmented matrix
   m := len(b0)
   a := make([][]float64, m)
   for i, ai := range a0 {
       row := make([]float64, m+1)
       copy(row, ai)
       row[m] = b0[i]
       a[i] = row
   }
   // WP algorithm from Gaussian elimination page
   // produces row-eschelon form
   for k := range a {
       // Find pivot for column k:
       iMax := k
       max := math.Abs(a[k][k])
       for i := k + 1; i < m; i++ {
           if abs := math.Abs(a[i][k]); abs > max {
               iMax = i
               max = abs
           }
       }
       if a[iMax][k] == 0 {
           return nil, errors.New("singular")
       }
       // swap rows(k, i_max)
       a[k], a[iMax] = a[iMax], a[k]
       // Do for all rows below pivot:
       for i := k + 1; i < m; i++ {
           // Do for all remaining elements in current row:
           for j := k + 1; j <= m; j++ {
               a[i][j] -= a[k][j] * (a[i][k] / a[k][k])
           }
           // Fill lower triangular matrix with zeros:
           a[i][k] = 0
       }
   }
   // end of WP algorithm.
   // now back substitute to get result.
   x := make([]float64, m)
   for i := m - 1; i >= 0; i-- {
       x[i] = a[i][m]
       for j := i + 1; j < m; j++ {
           x[i] -= a[i][j] * x[j]
       }
       x[i] /= a[i][i]
   }
   return x, nil

}</lang>

Output:
[-0.01 1.6027903945020987 -1.613203059905494 1.245494121371364 -0.49098971958462834 0.06576069617522803]

J

%. , J's matrix divide verb, directly solves systems of determined and of over-determined linear equations directly. This example J session builds a noisy sine curve on the half circle, fits quintic and quadratic equations, and displays the results of evaluating these polynomials.

<lang J> 

  f=: 6j2&":   NB. formatting verb

  sin=: 1&o.   NB. verb to evaluate circle function 1, the sine

  add_noise=: ] + (* (_0.5 + 0 ?@:#~ #))   NB. AMPLITUDE add_noise SIGNAL

  f RADIANS=: o.@:(%~ i.@:>:)5  NB. monadic circle function is  pi times
 0.00  0.63  1.26  1.88  2.51  3.14

  f SINES=: sin RADIANS
 0.00  0.59  0.95  0.95  0.59  0.00

  f NOISY_SINES=: 0.1 add_noise SINES
_0.01  0.61  0.91  0.99  0.60  0.02

  A=: (^/ i.@:#) RADIANS  NB. A is the quintic coefficient matrix

  NB. display the equation to solve
  (f A) ; 'x' ; '=' ; f@:,. NOISY_SINES

┌────────────────────────────────────┬─┬─┬──────┐ │ 1.00 0.00 0.00 0.00 0.00 0.00│x│=│ _0.01│ │ 1.00 0.63 0.39 0.25 0.16 0.10│ │ │ 0.61│ │ 1.00 1.26 1.58 1.98 2.49 3.13│ │ │ 0.91│ │ 1.00 1.88 3.55 6.70 12.62 23.80│ │ │ 0.99│ │ 1.00 2.51 6.32 15.88 39.90100.28│ │ │ 0.60│ │ 1.00 3.14 9.87 31.01 97.41306.02│ │ │ 0.02│ └────────────────────────────────────┴─┴─┴──────┘ 

  f QUINTIC_COEFFICIENTS=: NOISY_SINES %. A   NB. %. solves the linear system
_0.01  1.71 _1.88  1.48 _0.58  0.08

  quintic=: QUINTIC_COEFFICIENTS&p.  NB. verb to evaluate the polynomial

  NB. %. also solves the least squares fit for overdetermined system
  quadratic=: (NOISY_SINES %. (^/ i.@:3:) RADIANS)&p.  NB. verb to evaluate quadratic.
  quadratic

_0.0200630695393961729 1.26066877804926536 _0.398275112136019516&p. 

  NB. The quintic is agrees with the noisy data, as it should
  f@:(NOISY_SINES ,. sin ,. quadratic ,. quintic) RADIANS
_0.01  0.00 _0.02 _0.01
 0.61  0.59  0.61  0.61
 0.91  0.95  0.94  0.91
 0.99  0.95  0.94  0.99
 0.60  0.59  0.63  0.60
 0.02  0.00  0.01  0.02

  f MID_POINTS=: (+ -:@:(-/@:(2&{.)))RADIANS
_0.31  0.31  0.94  1.57  2.20  2.83

  f@:(sin ,. quadratic ,. quintic) MID_POINTS
_0.31 _0.46 _0.79
 0.31  0.34  0.38
 0.81  0.81  0.77
 1.00  0.98  1.00
 0.81  0.83  0.86
 0.31  0.36  0.27

</lang>

Julia

Using built-in LAPACK-based linear solver (which employs partial-pivoted Gaussian elimination): <lang julia>x = A \ b</lang>

Mathematica

<lang Mathematica>GaussianElimination[A_?MatrixQ, b_?VectorQ] := Last /@ RowReduce[Flatten /@ Transpose[{A, b}]]</lang>

MATLAB

<lang MATLAB> function [ x ] = GaussElim( A, b)

% Ensures A is n by n sz = size(A); if sz(1)~=sz(2) 

   fprintf('A is not n by n\n');
   clear x;
   return;

end

n = sz(1);

% Ensures b is n x 1. if n~=sz(1) 

   fprintf('b is not 1 by n.\n');
   return

end

x = zeros(n,1); aug = [A b]; tempmatrix = aug;

for i=2:sz(1) 


   % Find maximum of row and divide by the maximum
   tempmatrix(1,:) = tempmatrix(1,:)/max(tempmatrix(1,:));
   
   % Finds the maximum in column
   temp = find(abs(tempmatrix) - max(abs(tempmatrix(:,1))));
   if length(temp)>2
       for j=1:length(temp)-1
           if j~=temp(j)
               maxi = j; %maxi = column number of maximum
               break;
           end
       end
   else % length(temp)==2
       maxi=1;
   end
   
   % Row swap if maxi is not 1
   if maxi~=1
       temp = tempmatrix(maxi,:);
       tempmatrix(maxi,:) = tempmatrix(1,:);
       tempmatrix(1,:) = temp;
   end    
   
   % Row reducing
   for j=2:length(tempmatrix)-1
       tempmatrix(j,:) = tempmatrix(j,:)-tempmatrix(j,1)/tempmatrix(1,1)*tempmatrix(1,:);
       if tempmatrix(j,j)==0 || isnan(tempmatrix(j,j)) || abs(tempmatrix(j,j))==Inf
           fprintf('Error: Matrix is singular.\n');
           clear x;
           return
       end
   end
   aug(i-1:end,i-1:end) = tempmatrix;
   
   % Decrease matrix size
   tempmatrix = tempmatrix(2:end,2:end);

end

% Backwards Substitution x(end) = aug(end,end)/aug(end,end-1); for i=n-1:-1:1 

   x(i) = (aug(i,end)-dot(aug(i,1:end-1),x))/aug(i,i);

end

end </lang>

PARI/GP

If A and B have floating-point numbers (t_REALs) then the following uses Gaussian elimination: <lang parigp>matsolve(A,B)</lang>

If the entries are integers, then p-adic lifting (Dixon 1982) is used instead.

PHP

<lang php>function swap_rows(&$a, &$b, $r1, $r2) { 

   if ($r1 == $r2) return;

   $tmp = $a[$r1];
   $a[$r1] = $a[$r2];
   $a[$r2] = $tmp;

   $tmp = $b[$r1];
   $b[$r1] = $b[$r2];
   $b[$r2] = $tmp;

}

function gauss_eliminate($A, $b, $N) { 

   for ($col = 0; $col < $N; $col++)
   {
       $j = $col;
       $max = $A[$j][$j];

       for ($i = $col + 1; $i < $N; $i++)
       {
           $tmp = abs($A[$i][$col]);
           if ($tmp > $max)
           {
               $j = $i;
               $max = $tmp;
           }
       }

       swap_rows($A, $b, $col, $j);

       for ($i = $col + 1; $i < $N; $i++)
       {
           $tmp = $A[$i][$col] / $A[$col][$col];
           for ($j = $col + 1; $j < $N; $j++)
           {
               $A[$i][$j] -= $tmp * $A[$col][$j];
           }
           $A[$i][$col] = 0;
           $b[$i] -= $tmp * $b[$col];
       }
   }
   $x = array();
   for ($col = $N - 1; $col >= 0; $col--)
   {
       $tmp = $b[$col];
       for ($j = $N - 1; $j > $col; $j--)
       {
           $tmp -= $x[$j] * $A[$col][$j];
       }
       $x[$col] = $tmp / $A[$col][$col];
   }
   return $x;

}

function test_gauss() { 

   $a = array(
       array(1.00, 0.00, 0.00,  0.00,  0.00, 0.00),
       array(1.00, 0.63, 0.39,  0.25,  0.16, 0.10),
       array(1.00, 1.26, 1.58,  1.98,  2.49, 3.13),
       array(1.00, 1.88, 3.55,  6.70, 12.62, 23.80),
       array(1.00, 2.51, 6.32, 15.88, 39.90, 100.28),
       array(1.00, 3.14, 9.87, 31.01, 97.41, 306.02)
   );
   $b = array( -0.01, 0.61, 0.91, 0.99, 0.60, 0.02 );

   $x = gauss_eliminate($a, $b, 6);

   ksort($x);
   print_r($x);

}

test_gauss();</lang>

Output:
Array
(
    [0] => -0.01
    [1] => 1.6027903945021
    [2] => -1.6132030599055
    [3] => 1.2454941213714
    [4] => -0.49098971958463
    [5] => 0.065760696175228
)

Racket

<lang racket>

  1. lang racket

(require math/matrix) (define A 

 (matrix [[1.00  0.00  0.00  0.00  0.00   0.00]
          [1.00  0.63  0.39  0.25  0.16   0.10]
          [1.00  1.26  1.58  1.98  2.49   3.13]
          [1.00  1.88  3.55  6.70 12.62  23.80]
          [1.00  2.51  6.32 15.88 39.90 100.28]
          [1.00  3.14  9.87 31.01 97.41 306.02]]))

(define b (col-matrix [-0.01 0.61 0.91 0.99 0.60 0.02]))

(matrix-solve A b) </lang> Output: <lang racket>

  1. <array
 '#(6 1)
 #[-0.01
  1.602790394502109
  -1.613203059905556
  1.2454941213714346
  -0.4909897195846582
  0.06576069617523222]>

</lang>

Tcl

Library: Tcllib (Package: math::linearalgebra)

<lang tcl>package require math::linearalgebra

set A { 

   {1.00  0.00  0.00  0.00  0.00   0.00}
   {1.00  0.63  0.39  0.25  0.16   0.10}
   {1.00  1.26  1.58  1.98  2.49   3.13}
   {1.00  1.88  3.55  6.70 12.62  23.80}
   {1.00  2.51  6.32 15.88 39.90 100.28}
   {1.00  3.14  9.87 31.01 97.41 306.02}

} set b {-0.01 0.61 0.91 0.99 0.60 0.02} puts -nonewline [math::linearalgebra::show [math::linearalgebra::solveGauss $A $b] "%.2f"]</lang>

Output:
-0.01
1.60
-1.61
1.25
-0.49
0.07
Retrieved from "https://rosettacode.org/wiki/Gaussian_elimination?oldid=212379"