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


*        Gaussian elimination      09/02/2019
GAUSSEL  CSECT
         USING  GAUSSEL,R13        base register
         B      72(R15)            skip savearea
         DC     17F'0'             savearea
         SAVE   (14,12)            save previous context
         ST     R13,4(R15)         link backward
         ST     R15,8(R13)         link forward
         LR     R13,R15            set addressability
         LA     R7,1               j=1
       DO WHILE=(C,R7,LE,N)        do j=1 to n
         LA     R9,1(R7)             j+1
         LR     R6,R9                i=j+1
       DO WHILE=(C,R6,LE,N)          do i=j+1 to n
         LR     R1,R7                  j
         MH     R1,=AL2(NN)            *n
         AR     R1,R7                  +j
         BCTR   R1,0                   j*n+j-1
         SLA    R1,2                   ~
         LE     F0,A-(NN*4)(R1)        a(j,j)
         LR     R1,R6                  i
         MH     R1,=AL2(NN)            *n
         AR     R1,R7                  j
         BCTR   R1,0                   i*n+j-1
         SLA    R1,2                   ~
         LE     F2,A-(NN*4)(R1)        a(i,j)
         DER    F0,F2                  a(j,j)/a(i,j)
         STE    F0,W                   w=a(j,j)/a(i,j)
         LR     R8,R9                  k=j+1
       DO WHILE=(C,R8,LE,N)            do k=j+1 to n
         LR     R1,R7                    j
         MH     R1,=AL2(NN)              *n
         AR     R1,R8                    +k
         BCTR   R1,0                     j*n+k-1
         SLA    R1,2                     ~
         LE     F0,A-(NN*4)(R1)          a(j,k)
         LR     R1,R6                    i
         MH     R1,=AL2(NN)              *n
         AR     R1,R8                    +k
         BCTR   R1,0                     i*n+k-1
         SLA    R1,2                     ~
         LE     F2,A-(NN*4)(R1)          a(i,k)
         LE     F6,W                     w
         MER    F6,F2                    *a(i,k)
         SER    F0,F6                    a(j,k)-w*a(i,k)
         STE    F0,A-(NN*4)(R1)          a(i,k)=a(j,k)-w*a(i,k)
         LA     R8,1(R8)                 k=k+1
       ENDDO    ,                      end do k
         LR     R1,R7                  j
         SLA    R1,2                   ~
         LE     F0,B-4(R1)             b(j)
         LR     R1,R6                  i
         SLA    R1,2                   ~
         LE     F2,B-4(R1)             b(i)
         LE     F6,W                   w
         MER    F6,F2                  *b(i)
         SER    F0,F6                  b(j)-w*b(i)
         STE    F0,B-4(R1)             b(i)=b(j)-w*b(i)
         LA     R6,1(R6)               i=i+1
       ENDDO    ,                    end do i
         LA     R7,1(R7)             j=j+1
       ENDDO    ,                  end do j
         L      R2,N               n
         SLA    R2,2               ~
         LE     F0,B-4(R1)         b(n)
         L      R1,N               n
         MH     R1,=AL2(NN)        *n
         A      R1,N               n
         BCTR   R1,0               n*n+n-1
         SLA    R1,2               ~
         LE     F2,A-(NN*4)(R1)    a(n,n)
         DER    F0,F2              b(n)/a(n,n)
         STE    F0,X-4(R2)         x(n)=b(n)/a(n,n)
         L      R7,N               n
         BCTR   R7,0               j=n-1
       DO WHILE=(C,R7,GE,=F'1')    do j=n-1 to 1 by -1
         LE     F0,=E'0'             0
         STE    F0,W                 w=0
         LA     R9,1(R7)             j+1
         LR     R6,R9                i=j+1
       DO WHILE=(C,R6,LE,N)          do i=j+1 to n
         LR     R1,R7                  j
         MH     R1,=AL2(NN)            *n
         AR     R1,R6                  i
         BCTR   R1,0                   j*n+i-1  
         SLA    R1,2                   ~
         LE     F0,A-(NN*4)(R1)        a(j,i)
         LR     R1,R6                  i
         SLA    R1,2                   ~
         LE     F2,X-4(R1)             x(i)
         MER    F0,F2                  a(j,i)*x(i)
         LE     F6,W                   w
         AER    F6,F0                  +a(j,i)*x(i)
         STE    F6,W                   w=w+a(j,i)*x(i)
         LA     R6,1(R6)               i=i+1
       ENDDO    ,                    end do i
         LR     R2,R7                j
         SLA    R2,2                 ~
         LE     F0,B-4(R2)           b(j)
         SE     F0,W                 -w
         LR     R1,R7                j
         MH     R1,=AL2(NN)          *n
         AR     R1,R7                j
         BCTR   R1,0                 j*n+j-1
         SLA    R1,2                 ~
         LE     F2,A-(NN*4)(R1)      a(j,j)
         DER    F0,F2                (b(j)-w)/a(j,j)
         STE    F0,X-4(R2)           x(j)=(b(j)-w)/a(j,j)
         BCTR   R7,0                 j=j-1
       ENDDO    ,                  end do j
         XPRNT  =CL8'SOLUTION',8   print
         MVC    PG,=CL91' '        clear buffer
         LA     R6,1               i=1
       DO WHILE=(C,R6,LE,N)        do i=1 to n
         LR     R1,R6                i
         SLA    R1,2                 ~
         LE     F0,X-4(R1)           x(i)
         LA     R0,5                 number of decimals
         BAL    R14,FORMATF          edit
         MVC    PG(13),0(R1)         output
         XPRNT  PG,L'PG              print
         LA     R6,1(R6)             i=i+1
       ENDDO    ,                  end do i
         L      R13,4(0,R13)       restore previous savearea pointer
         RETURN (14,12),RC=0       restore registers from calling sav
         COPY   plig\$_FORMATF.MLC format F13.n
NN       EQU    (X-B)/4            n
N        DC     A(NN)              n
A        DC  E'1',E'0',E'0',E'0',E'0',E'0'
         DC  E'1',E'0.63',E'0.39',E'0.25',E'0.16',E'0.10'
         DC  E'1',E'1.26',E'1.58',E'1.98',E'2.49',E'3.13'
         DC  E'1',E'1.88',E'3.55',E'6.70',E'12.62',E'23.80'
         DC  E'1',E'2.51',E'6.32',E'15.88',E'39.90',E'100.28'
         DC  E'1',E'3.14',E'9.87',E'31.01',E'97.41',E'306.02'
B        DC  E'-0.01',E'0.61',E'0.91',E'0.99',E'0.60',E'0.02'
X        DS     (NN)E              x(n)
W        DS     E                  w
PG       DC     CL91' '            buffer
         REGEQU
         END    GAUSSEL