      program tgauss
      integer n1,m1,L1
      integer NP,MP,SR
      INTEGER i,j,k
      INTEGER fp
      PARAMETER(NP=50,MP=50)
      real ma1(NP,MP),TA1(NP,MP),XMM(NP,MP), b(NP),x(NP)

      SR=1
      open (unit=8,file="linsys1.dat")
      fp=8

      read (fp,*) 
      read (fp,*) n1,m1   
      call readmatrix(MA1,n1,m1,fp)

      read (fp,*) 
      read (fp,*) L1 
      call readvektor(b,L1,fp)
      do 13 i=1,n1
         do 11 j=1,m1
            TA1(i,j)=MA1(i,j)
 11      continue

 13   continue

      do  k=1,L1
         x(k)=b(k)
      End Do

      write(*,*) ' Matrix A : '
      call displaymatrix(TA1,N1,M1)
         write(*,*) ' L, and vector ',L1
      call displayvektor(x,L1)
      call gaussj(TA1,n1,NP,x,1,MP)
      write(*,*) 'Inverse of Matrix A : '
      call displaymatrix(TA1,n1,m1)
      call multMatrixMatrix(TA1,MA1,XMM,n1,m1,n1,m1)
       write(*,*) 'Inverse of Matrix A * A: '
      call displaymatrix(XMM,n1,m1)
       print *, "Solution is: "
       call displayvektor(x,L1)
      end program
