c    11.LU$BJ,2rK!(B
c      $B?9@n(B $BLwBg(B , 22980138
c
      program ex11
c
      implicit none
      integer :: i,j,n
      real :: a(2,3),b(3,4),c(3,4)
c
      data (a(1,j),j=1,3)/ 4.0e0, 3.0e0, 0.0e0/
      data (a(2,j),j=1,3)/ 3.0e0, 2.0e0, 1.0e0/

      data (b(1,j),j=1,4)/ 2.0e0, 5.0e0,-3.0e0, 3.0e0/
      data (b(2,j),j=1,4)/ 4.0e0, 5.0e0,-3.0e0, 5.0e0/
      data (b(3,j),j=1,4)/-1.0e0,-5.0e0, 4.0e0, 1.0e0/

      data (c(1,j),j=1,4)/ 1.0e0, 1.0e0, 1.0e0, 3.0e0/
      data (c(2,j),j=1,4)/ 2.0e0, 3.0e0, 4.0e0, 5.0e0/
      data (c(3,j),j=1,4)/ 1.0e0, 2.0e0, 3.0e0, 2.0e0/

      write(6,*)'11.LU$BJ,2rK!(B'
      write(6,*)''
      write(6,*) '(1)'
      write(6,*)'$B=i4|CM(B'
      do i=1,2
         write(6,100) i,'   ',(a(i,j),j=1,3)
      enddo
      n=2
      call LU(a,n)

      write(6,*)''
      write(6,*) '(2)'
      write(6,*)'$B=i4|CM(B'
      do i=1,3
         write(6,100) i,'   ',(b(i,j),j=1,4)
      enddo
      n=3
      call LU(b,n)

      write(6,*)''
      write(6,*) '(3)'
      write(6,*)'$B=i4|CM(B'
      do i=1,3
         write(6,100) i,'   ',(c(i,j),j=1,4)
      enddo
      n=3
      call LU(c,n)

      stop
  100 format(1x,i2,1x,a3,1x,f7.4,1x,f7.4,1x,f7.4,1x,f7.4) 
      end

c-----LU$BJ,2rK!$N%5%V%k!<%A%s(B
      subroutine LU(a,n)
      implicit none
      integer :: i,j,k,n,m,l,maxi
      real :: maxa
      real :: a(n,n+1),aa(n,n),bb(n),ll(n,n),
     -        uu(n,n),da(n),db,ma(3),x(n),y(n)
      
      do i=1,n
       do j=1,n
        aa(i,j)=a(i,j)
        bb(i)=a(i,n+1)
       enddo
      enddo
      write(6,*)''
      do m=1,n
        if (n==2) then
        write(6,200) m,'   ',(aa(m,l),l=1,n),'   ',bb(m)
        elseif (n==3) then
        write(6,300) m,'   ',(aa(m,l),l=1,n),'   ',bb(m)
        endif
      enddo
      do k=1,1
          maxi=k
          do i=k,n
             if ( abs(a(i,k)) > abs(a(maxi,k)) ) then
                maxi=i
             endif
          enddo
          do j=k,n
             da(j)     =aa(maxi,j)
             aa(maxi,j)=aa(k,j)
             aa(k,j)   =da(j)
             db        =bb(maxi)
             bb(maxi)  =bb(k)
             bb(k)     =db
          enddo
          if (i/=k) then
            write(6,*)''
            do m=1,n
             if (n==2) then
               write(6,200) m,'   ',(aa(m,l),l=1,n),'   ',bb(m)
             elseif (n==3) then
               write(6,300) m,'   ',(aa(m,l),l=1,n),'   ',bb(m)
             endif
            enddo
          end if
c         do i=k,n
c           ma(i)=abs(aa(i,k))
c         enddo
c         maxa = max(ma(1),ma(2),ma(3))
c         do i=k,n
c          if (ma(i) == maxa) then
c             do j=k,n
c                da(j)  =aa(i,j)
c                aa(i,j)=aa(k,j)
c                aa(k,j)=da(j)
c                db     =bb(i)
c                bb(i)  =bb(k)
c                bb(k)  =db
c             enddo
c             if (i/=k) then
c               write(6,*)''
c               do m=1,n
c                if (n==2) then
c                  write(6,200) m,'   ',(aa(m,l),l=1,n),'   ',bb(m)
c                elseif (n==3) then
c                  write(6,300) m,'   ',(aa(m,l),l=1,n),'   ',bb(m)
c                endif
c               enddo
c             end if
c          end if
c         enddo
      enddo
      do k=1,n
         do j=1,n
          if (j <= k) then
            ll(k,j)=aa(k,j)
             do m=1,j-1
               ll(k,j)=ll(k,j)-ll(k,m)*uu(m,j)
             enddo
          elseif (j > k) then
            uu(k,j)=aa(k,j)/ll(k,k)
             do m=1,k-1
                uu(k,j)=uu(k,j)-(ll(k,m)*uu(m,j))/ll(k,k)
             enddo
          endif
         enddo
      enddo
      do i=1,n
       do j=i,n
          aa(i,j)=ll(i,j)+uu(i,j)
       enddo
      enddo
      write(6,*)''
      write(6,*)'LU$BJ,2r(B:'
      do m=1,n
        if (n==2) then
        write(6,200) m,'   ',(aa(m,l),l=1,n),'   ',bb(m)
        elseif (n==3) then
        write(6,300) m,'   ',(aa(m,l),l=1,n),'   ',bb(m)
        endif
      enddo
      if (aa(n,n)==0.0e0) go to 10
      write(6,*)''
      write(6,*)'Ly=b:'
      do i=1,n
         y(i)=bb(i)/ll(i,i)
         do j=1,i-1
            y(i)=y(i)-(ll(i,j)*y(j))/ll(i,i)
         enddo
      enddo
      write(6,'(a4,f7.4)')'y1=',y(1)
      write(6,'(a4,f7.4)')'y2=',y(2)
      if (n==3) then
      write(6,'(a4,f7.4)')'y3=',y(3)
      endif

      write(6,*)''
      write(6,*)'Ux=y:'
      do i=n,1,-1
         x(i)=y(i)
         do j=n,i+1,-1
            x(i)=x(i)-uu(i,j)*x(j)
         enddo
      enddo
      write(6,'(a4,f7.4)')'x1=',x(1)
      write(6,'(a4,f7.4)')'x2=',x(2)
      if (n==3) then
      write(6,'(a4,f7.4)')'x3=',x(3)
      endif

      write(6,*)''
      write(6,*)'$B2r(B'
      write(6,'(a3,f7.4)')'x=',x(1)
      write(6,'(a3,f7.4)')'y=',x(2)
      if (n==3) then
      write(6,'(a3,f7.4)')'z=',x(3)
      endif
      go to 30
   10 write(6,*)''
      write(6,*)'$B2r$OITDj(B'  
      go to 30
   20 write(6,*)'$B2r$OITG=(B'
   30 continue  
      return
  100 format(1x,i2,1x,a3,1x,f7.4,1x,f7.4,1x,f7.4,1x,f7.4) 
  200 format(1x,i2,1x,a3,1x,f7.4,1x,f7.4,
     -       1x,a3,1x,f7.4) 
  300 format(1x,i2,1x,a3,1x,f7.4,1x,f7.4,1x,f7.4,
     -       1x,a3,1x,f7.4) 
      end
