Task

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.



See also

   the Wikipedia entry:   Gaussian elimination


 
module checkit {
      Dim Base 1, a(6, 6), b(6)
      a(1,1)= 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
      \\ remove \\ to feed next array
      \\ a(1,1)=1.1,0.12,0.13,0.12,0.14,-0.12,1.21,0.63,0.39,0.25,0.16,0.1,1.03,1.26,1.58,1.98,2.49,3.13, 1.06,1.88,3.55,6.7,12.62,23.8, 1.12,2.51,6.32,15.88,39.9,100.28,1.16,3.14,9.87,31.01,97.41,306.02
      for i=1 to 6 : for j=1 to 6 : a(i,j)=val(a(i,j)->Decimal) :Next j:Next i
      b(1)=-0.01, 0.61, 0.91, 0.99, 0.60, 0.02
      for i=1 to 6 : b(i)=val(b(i)->Decimal) :Next i
      function GaussJordan(a(), b()) {
            cols=dimension(a(),1)
            rows=dimension(a(),2)
            \\ make augmented matrix
            Dim Base 1, a(cols, rows)
            \\ feed array with rationals
            Dim Base 1, b(Len(b()))
            for diag=1 to rows {
                        max_row=diag
                        max_val=abs(a(diag, diag))
                        if diag<rows Then {
                              for ro=diag+1 to rows {
                                    d=abs(a(ro, diag))
                                    if d>max_val then max_row=ro : max_val=d
                              }
                        }
            \\         SwapRows diag, max_row
                        if diag<>max_row then {
                              for i=1 to cols {
                                    swap a(diag, i), a(max_row, i)
                              }
                              swap b(diag), b(max_row)
                        }
                        invd= a(diag, diag)
                        if diag<=cols then {
                              for col=diag to cols {
                                    a(diag, col)/=invd
                              }
                        }
                        b(diag)/=invd
                        for ro=1 to rows {
                              d1=a(ro,diag)
                              d2=d1*b(diag)
                              if ro<>diag Then {
                                         for col=diag to cols {a(ro, col)-=d1*a(diag, col)}
                                          b(ro)-=d2  
                              }
                        }        
                  }
            =b()
      }
      Function ArrayLines$(a(), leftmargin=6, maxwidth=8,decimals$="") {
            \\ defualt  no set  decimals, can show any number
            ex$={
            }
           const way$=", {0:"+decimals$+":-"+str$(maxwidth,"")+"}"
            if dimension(a())=1 then {
                  m=each(a())
                  while m {ex$+=format$(way$,array(m))}
                  Insert 3, 2  ex$=string$(" ", leftmargin)
                  =ex$ :    Break
            } 
            for i=1 to dimension(a(),1)  {
                  ex1$=""
                  for j=1 to dimension(a(),2 ) {
                              ex1$+=format$(way$,a(i,j))
                  }
                  Insert 1,2  ex1$=string$(" ", leftmargin)
                  ex$+=ex1$+{
                  }
            }
            =ex$
      }
      mm=GaussJordan(a(), b())
            c=each(mm)
            while c {
                  print array(c)
            }
      \\ check accuracy
      link mm to r()
      \\ prepare output document
      Document out$={Algorithm using decimals
                  }+"Matrix A:"+ArrayLines$(a(),,,"2")+{
                  }+"Vector B:"+ArrayLines$(b(),,,"2")+{
                  }+"Solution: "+{
                  }
      acc=25
      for i=1 to  dimension(a(),1)
            sum=a(1,1)-a(1,1)
            For j=1 to dimension(a(),2)
                  sum+=r(j)*a(i,j)
            next j
            p$=format$("Coef. {0::-2},  rounding to {1} decimal, compare {2:-5}, solution: {3}", i, acc, round(sum-b(i),acc)=0@, r(i) )
            Print p$
            Out$=p$+{
            }
      next i
      Report out$
      clipboard out$
}
checkit
 