source: anuga_work/development/sudi/sw_1d/periodic_waves/johns/gaussPivot.py @ 7837

Last change on this file since 7837 was 7837, checked in by mungkasi, 14 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

File size: 1.1 KB
Line 
1## module gaussPivot
2''' x = gaussPivot(a,b,tol=1.0e-9).
3    Solves [a]{x} = {b} by Gauss elimination with
4    scaled row pivoting
5'''   
6from Numeric import *
7import swap 
8import error
9
10def gaussPivot(a,b,tol=1.0e-12):
11    n = len(b)
12   
13  # Set up scale factors
14    s = zeros((n),Float)
15    for i in range(n):
16        s[i] = max(abs(a[i,:]))
17           
18    for k in range(0,n-1):
19       
20      # Row interchange, if needed
21        p = int(argmax(abs(a[k:n,k])/s[k:n])) + k
22        if abs(a[p,k]) < tol: error.err('Matrix is singular')
23        if p != k:
24            swap.swapRows(b,k,p)
25            swap.swapRows(s,k,p)
26            swap.swapRows(a,k,p)
27           
28      # Elimination
29        for i in range(k+1,n):
30            if a[i,k] != 0.0:
31                lam = a[i,k]/a[k,k]
32                a[i,k+1:n] = a [i,k+1:n] - lam*a[k,k+1:n]
33                b[i] = b[i] - lam*b[k]
34    if abs(a[n-1,n-1]) < tol: error.err('Matrix is singular')
35                   
36  # Back substitution
37    for k in range(n-1,-1,-1):
38        b[k] = (b[k] - dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
39    return b
40
41
42
43       
Note: See TracBrowser for help on using the repository browser.