Gaussian elimination: Difference between revisions

{{trans|C}}

 

I r1 != r2

swap(&a[r1], &a[r2])

V b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02]

 

 

{{out}}

=={{header|360 Assembly}}==

{{trans|PL/I}}

GAUSSEL CSECT

USING GAUSSEL,R13 base register

PG DC CL91' ' buffer

REGEQU

{{out}}

<pre>

 

=={{header|Ada}}==

with Ada.Numerics.Generic_Real_Arrays;

 

begin

Put (X);

 

{{out}}

{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-2.4.1 algol68g-2.4.1].}}

{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted]}}

COMMENT PROVIDES

MODE FIXED; INT fixed exception, unfixed exception;

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

 

COMMENT PRELUDE REQUIRES

MODE SCAL = REAL;

);

 

COMMENT PRELUDE REQUIRES

MODE SCAL = REAL,

);

 

COMMENT POSTLUDE PROIVIDES

PROC VOID exception too many iterations, exception value error;

exception too many iterations:

exception value error:

# -*- coding: utf-8 -*- #

 

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

 

<pre>

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

</pre>

 

=={{header|C}}==

This modifies A and b in place, which might not be quite desirable.

#include <stdlib.h>

#include <math.h>

 

return 0;

-0.01

1.60279

=={{header|C sharp|C#}}==

This modifies A and b in place, which might not be quite desirable.

using System;

 

}

}

using System;

 

}

}

<pre>{{out}}

-0.0597391027501976

=={{header|C++}}==

{{trans|Go}}

#include <cassert>

#include <cmath>

}

return 0;

 

{{out}}

 

=={{header|Common Lisp}}==

(defmacro mapcar-1 (fn n list)

"Maps a function of two parameters where the first one is fixed, over a list"

(let ((m1 (redc (rev (redc m)))))

(reverse (mapcar #'(lambda (r) (let ((pivot (find-if-not #'zerop r))) (if pivot (/ (car (last r)) pivot) (throw 'result 'singular)))) m1)) ))))

 

{{out}}

=={{header|D}}==

{{trans|Go}}

std.typecons;

 

if (abs(r.x[i] - xi) > eps)

return writeln("Out of tolerance: ", r.x[i], " ", xi);

{{out}}

<pre>[-0.01, 1.60279, -1.6132, 1.24549, -0.49099, 0.0657607]</pre>

 

=={{header|Delphi}}==

 

// Modified from:

end.

 

{{out}}

<pre>1.00, 0.00, 0.00, 4.00</pre>

 

=={{header|F_Sharp|F#}}==

===The Function===

// Gaussian Elimination. Nigel Galloway: February 2nd., 2019

let gelim augM=

|_->fN (i+1) (fG e g i@[g])

fN 0 augM

===The Task===

This task uses functionality from [[Continued_fraction/Arithmetic/Construct_from_rational_number#F.23]] and [[Continued_fraction#F.23]]

let test=[[ -6I; -18I; 13I; 6I; -6I; -15I; -2I; -9I; -231I];

[ 2I; 20I; 9I; 2I; 16I; -12I; -18I; -5I; 647I];

let fN (n,g)=cN2S(π(rI2cf n g))

gelim test |> List.map fN |> List.iteri(fun i n->(printfn "x[%d]=%1.14f " (i+1) (snd (Seq.pairwise n|> Seq.find(fun (n,g)->n-g < 0.0000000000001M)))))

{{out}}

<pre>

=={{header|Fortran}}==

Gaussian Elimination with partial pivoting using augmented matrix

program ge

 

 

end program

 

=={{header|FreeBASIC}}==

Gaussian elimination with pivoting.

FreeBASIC version 1.05

 

Sub GaussJordan(matrix() As Double,rhs() As Double,ans() As Double)

next n

sleep

{{out}}

<pre>

</pre>

=={{header|Generic}}==

generic coordinaat

{

}

 

 

=={{header|Go}}==

===Partial pivoting, no scaling===

Gaussian elimination with partial pivoting by [https://en.wikipedia.org/wiki/Gaussian_elimination#Pseudocode pseudocode] on WP page [https://en.wikipedia.org/wiki/Gaussian_elimination Gaussian elimination]."

 

import (

}

return x, nil

{{out}}

<pre>

===Scaled partial pivoting===

Changes from above version noted with comments. For the example data scaling does help a bit.

 

import (

}

return x, nil

{{out}}

<pre>

 

===From scratch===

isMatrix xs = null xs || all ((== (length.head $ xs)).length) xs

 

let b = [[-0.01], [0.61], [0.91], [0.99], [0.60], [0.02]]

mapM_ print $ gauss a b

{{out}}

<pre>

===Another example===

We use Rational numbers for having more precision. a % b is the rational a / b.

mult uss vss = map ((\xs -> if null xs then [] else foldl1 (zipWith (+)) xs). zipWith (\vs u -> map (u*) vs) vss) uss

 

let b = [[-0.01], [0.61], [0.91], [0.99], [0.60], [0.02]]

task a b

{{out}}

<pre>

 

===Determinant and permutation matrix are given===

mult uss vss = map ((\xs -> if null xs then [] else foldl1 (zipWith (+)) xs). zipWith (\vs u -> map (u*) vs) vss) uss

 

let b = [[-0.01], [0.61], [0.91], [0.99], [0.60], [0.02]]

task a b

{{out}}

<pre>

%. , 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.

 

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

 

0.81 0.83 0.86

0.31 0.36 0.27

 

=={{header|Java}}==

Naked implementation, using Java arrays instead of a matrix class.

 

 

public class GaussianElimination {

}

}

 

=={{header|JavaScript}}==

From Numerical Recipes in C:

function lusolve(A, b, update) {

var lu = ludcmp(A, update)

[-0.01, 0.61, 0.91, 0.99, 0.60, 0.02]

)

 

{{output}}

<pre>-0.01000000000000004, 1.6027903945021095, -1.6132030599055475, 1.2454941213714232, -0.4909897195846526, 0.06576069617523138</pre>

 

=={{header|Julia}}==

Using built-in LAPACK-based linear solver (which employs partial-pivoted Gaussian elimination):

 

=={{header|Klong}}==

elim::{[h m];h::*m::x@>*'x;

:[2>#x;x;(,h),0,:\.f({1_x}'{x-h**x%*h}'1_m)]}

{v::v,((*x)-/:[[]~v;[];v*x@1+!#v])%x@1+#v}'||'x;|v}

gauss::{subst(elim(x))}

 

Example, matrix taken from C version:

 

gauss([[1.00 0.00 0.00 0.00 0.00 0.00 -0.01]

[1.00 0.63 0.39 0.25 0.16 0.10 0.61]

-0.490989719584658025

0.0657606961752320591]

 

=={{header|Kotlin}}==

{{trans|Go}}

 

val ta = arrayOf(

}

}

 

{{out}}

 

=={{header|Lambdatalk}}==

{require lib_matrix}

 

->

[-0.01,1.6027903945021094,-1.613203059905548,1.245494121371424,-0.49098971958465304,0.06576069617523143]

 

=={{header|Lobster}}==

{{trans|Go}}

 

// test case from Go version at http://rosettacode.org/wiki/Gaussian_elimination

print("out of tolerance, expected: " + tx[i] + " got: " + xi)

 

{{out}}

<pre>

=={{header|M2000 Interpreter}}==

Faster, with accuracy of 25 decimals

module checkit {

Dim Base 1, a(6, 6), b(6)

}

checkit

 

slower with accuracy of 26 decimals

Module Checkit2 {

Dim Base 1, a(6, 6), b(6)

}

Checkit2

{{out}}

<pre style="height:30ex;overflow:scroll">

 

=={{header|Mathematica}} / {{header|Wolfram Language}}==

 

=={{header|MATLAB}}==

function [ x ] = GaussElim( A, b)

 

 

end

 

=={{header|Modula-3}}==

 

;Matrix interface

 

(*

(* prints the matrix row-by-row; sorry, no special padding to line up columns *)

 

 

;Matrix implementation

 

 

IMPORT IO;

 

BEGIN

 

; interface for the ring of integers modulo an integer

 

 

(*

*)

 

 

;test implementation

It's fairly easy to initialize an array of types in Modula-3, but it can get cumbersome with structured types, so we wrote a procedure to convert an integer matrix to a matrix of integers modulo a number.

 

 

IMPORT IO, ModularRing AS MR, IntMatrix AS IM, ModMatrix AS MM;

IO.PutChar('\n');

 

 

{{out}}

=={{header|Nim}}==

{{trans|Kotlin}}

 

type

echo "Out of tolerance."

echo "Expected values are ", refx

 

{{out}}

=={{header|OCaml}}==

The OCaml stdlib is fairly lean, so these stand-alone solutions often need to include support functions which would be part of a codebase, like these...

module Array = struct

include Array

if n > b then acc else aux (f acc n) (succ n)

in aux init a

The solver:

(* Some less-general support functions for 'solve'. *)

let swap_elem m i j = let x = m.(i) in m.(i) <- m.(j); m.(j) <- x

done;

x

Example data...

let a =

[| [| 1.00; 0.00; 0.00; 0.00; 0.00; 0.00 |];

[| 1.00; 3.14; 9.87; 31.01; 97.41; 306.02 |] |]

let b = [| -0.01; 0.61; 0.91; 0.99; 0.60; 0.02 |]

In the REPL, the solution is:

# let x = solve a b;;

val x : float array =

[|-0.0100000000000000991; 1.60279039450210536; -1.61320305990553226;

1.24549412137140547; -0.490989719584644546; 0.0657606961752301433|]

Further, let's define multiplication and subtraction to check our results...

let mul m v =

Array.mapi (fun i u ->

 

let sub u v = Array.mapi (fun i e -> e -. v.(i)) u

Now 'x' can be plugged into the equation to calculate the residual:

# let residual = sub b (mul a x);;

val residual : float array =

[|9.8879238130678e-17; 1.11022302462515654e-16; 2.22044604925031308e-16;

8.88178419700125232e-16; -5.5511151231257827e-16; 4.26741975090294545e-16|]

 

=={{header|PARI/GP}}==

If A and B have floating-point numbers (<code>t_REAL</code>s) then the following uses Gaussian elimination:

 

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

=={{header|Perl}}==

{{libheader|Math&#58;&#58;Matrix}}

my $a = Math::Matrix->new([0,1,0],

[0,0,1],

[4]);

my $x = $a->concat($b)->solve;

 

<code>Math::Matrix</code> <code>solve()</code> expects the column vector to be an extra column in the matrix, hence <code>concat()</code>. Putting not just a column there but a whole identity matrix (making Nx2N) is how its <code>invert()</code> is implemented. Note that <code>solve()</code> doesn't notice singular matrices and still gives a return when there is in fact no solution to Ax=B.

=={{header|Phix}}==

{{trans|PHP}}

{{out}}

<pre>

 

=={{header|PHP}}==

{

if ($r1 == $r2) return;

}

 

 

{{out}}

 

=={{header|PL/I}}==

 

declare n fixed binary;

end Backward_substitution;

 

{{out}}

<pre>

=={{header|PowerShell}}==

===Gauss===

function gauss($a,$b) {

$n = $a.count

gauss $a $b

 

<b>Output:</b>

<pre>

</pre>

===Gauss-Jordan===

function gauss-jordan($a,$b) {

$n = $a.count

"x ="

gauss-jordan $a $b

<b>Output:</b>

<pre>

 

=={{header|Python}}==

# If 'b' is the identity, then 'x' is the inverse of 'a'.

 

[[Fraction(1, 1), Fraction(0, 1), Fraction(0, 1)],

[Fraction(0, 1), Fraction(1, 1), Fraction(0, 1)],

===Using numpy===

$ python3

Python 3.6.0 |Anaconda custom (64-bit)| (default, Dec 23 2016, 12:22:00)

[ 0.06388889, -0.14444444, 0.14722222]])

>>>

 

=={{header|R}}==

Here 'b' is a matrix. Partial pivoting is used, and the determinant of 'a' is returned as well.

 

n <- nrow(a)

det <- 1

 

a <- matrix(c(2, 9, 4, 7, 5, 3, 6, 1, 8), 3, 3, byrow=T)

 

{{out}}

 

=={{header|Racket}}==

#lang racket

(require math/matrix)

 

(matrix-solve A b)

{{out}}

#<array

'#(6 1)

-0.4909897195846582

0.06576069617523222]>

 

=={{header|Raku}}==

(formerly Perl 6)

Gaussian elimination results in a matrix in row echelon form. Gaussian elimination with back-substitution (also known as Gauss-Jordan elimination) results in a matrix in reduced row echelon form. That being the case, we can reuse much of the code from the [[Reduced row echelon form]] task. Raku stores and does calculations on decimal numbers within its limit of precision using Rational numbers by default, meaning the calculations are exact.

 

@b.kv.map: { @a[$^k].append: $^v };

@a.&rref[*]»[*-1];

}

 

sub rref (@m) {

 

for ^$rows -> $r {

}

@m[$i, $r] = @m[$r, $i] if $r != $i;

for ^$rows -> $n {

next if $n == $r;

}

++$lead;

sub rat-or-int ($num) {

return $num unless $num ~~ Rat;

$num.nude.join: '/';

}

say-it 'x matrix:', (my @gj = gauss-jordan-solve @a, @b), "%16.12f";

say-it 'or, x in exact rationals:', @gj, "%28s";

{{out}}

<pre>A matrix:

=={{header|REXX}}==

===version 1===

* 07.08.2014 Walter Pachl translated from PL/I)

* improved to get integer results for, e.g. this input:

x.j = (b.j - t) / a.j.j

end

{{out}}

<pre>The equations are:

 

Only &nbsp; '''8''' &nbsp; (fractional) decimal digits were used for the output display.

numeric digits 1000 /*heavy─duty decimal digits precision. */

parse arg iFID . /*obtain optional argument from the CL.*/

do o=1 for n; say right('x['o"] = ", 38) left('', x.o>=0) x.o/1

end /*o*/

{{out|input file|text=: &nbsp; &nbsp; <tt> GAUSS_E.DAT </tt>}}

<pre>

Also added was the rounding the residual numbers to zero &nbsp;if&nbsp; the number of significant decimal digits was &nbsp; <big> &le; </big> &nbsp;5%&nbsp; of

<br>the number of significant fractional decimal digits &nbsp; (in this case, &nbsp;5%&nbsp; of &nbsp;1,000&nbsp; digits for the decimal fraction).

numeric digits 1000 /*heavy─duty decimal digits precision. */

parse arg iFID . /*obtain optional argument from the CL.*/

do o=1 for n; say right('x['o"] = ", 38) left('', x.o>=0) x.o/1

end /*o*/

{{out|output|text=&nbsp; when using the same default input file as for version 2:}}

<pre>

 

=={{header|Ruby}}==

require 'bigdecimal/ludcmp'

include LUSolve

one = BigDecimal("1.0")

 

{{Output}}

<pre>

 

=={{header|Rust}}==

// using a Vec<f32> might be a better idea

// for now, let us create a fixed size array

}

}

 

=={{header|Sidef}}==

Uses the '''rref(A)''' function from [https://rosettacode.org/wiki/Reduced_row_echelon_form#Sidef Reduced row echelon form].

{{trans|Raku}}

 

var A = gather {

 

var G = gauss_jordan_solve(a, b)

{{out}}

<pre>

This implementation computes also the determinant of the matrix A, as it requires only a few operations. The matrix B is overwritten with the solution of the system, and A is overwritten with garbage.

 

real scalar i,j,n,s

real vector js

det = det*a[i,i]

}

 

=== LU decomposition and backsubstitution ===

 

real scalar i,j,n,s

real vector js

y[i,.] = y[i,.]/u[i,i]

}

 

=== Example ===

Here we are computing the inverse of a 3x3 matrix (which happens to be a magic square), using both methods.

 

 

: b

2 | .1055555556 .0222222222 -.0611111111 |

3 | .0638888889 -.1444444444 .1472222222 |

 

=={{header|Swift}}==

{{trans|Rust}}

 

var system = sys

 

for (i, f) in sols.enumerated() {

print("X\(i + 1) = \(f)")

 

{{out}}

=={{header|Tcl}}==

{{tcllib|math::linearalgebra}}

 

set A {

}

set b {-0.01 0.61 0.91 0.99 0.60 0.02}

{{out}}

<pre>

On TI-83 or TI-84, another way to solve this task is to use the '''PlySmlt2''' internal apps and choose

"simult equ solver" with 6 equations and 6 unknowns.

[ 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 3.14 9.87 31.01 97.41 306.02 0.02]]→[A]

Matr>List(rref([A]),7,L1)

{{out}}

<pre>

 

=={{header|VBA}}==

Private Function gauss_eliminate(a As Variant, b As Variant) As Variant

Dim n As Integer: n = UBound(b)

t = Join(s, ", ")

Debug.Print t; ")"

<pre>(-0.01, 1.60279039450209, -1.61320305990548, 1.24549412137136, -0.490989719584628, 0.065760696175228)</pre>

 

=={{header|VBScript}}==

const n=6

dim a(6,6),b(6),x(6),ab

buf=buf&right(space(8)&formatnumber(x(i),2),8)&vbcrlf

next

{{out}}

<pre>

=={{header|Wren}}==

{{trans|Kotlin}}

 

var ta = [

}

i = i + 1

 

{{out}}

=={{header|zkl}}==

Using the GNU Scientific Library:

a:=GSL.Matrix(6,6).set(

1.00, 0.00, 0.00, 0.00, 0.00, 0.00,

b:=GSL.VectorFromData(-0.01, 0.61, 0.91, 0.99, 0.60, 0.02);

x:=a.AxEQb(b);

{{out}}

<pre>

Or, using lists:

{{trans|C}}

n:=b.len();

foreach dia in ([0..n-1]){

}

x

List(1.00, 0.63, 0.39, 0.25, 0.16, 0.10,),

List(1.00, 1.26, 1.58, 1.98, 2.49, 3.13,),

List(1.00, 3.14, 9.87, 31.01, 97.41, 306.02) );

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

{{out}}

<pre>L(-0.01,1.60279,-1.6132,1.24549,-0.49099,0.0657607)</pre>