lu_fac.f

二维热能方程的求解

c=======================================================================
      subroutine lufac(nband,neqn,level,A)
c=======================================================================
      implicit none
      integer nband,neqn,level
      real*8 A(-nband:nband,neqn)

      integer i,j,k,kb,kpi
      real*8  zero,pick,mult

      if(level.ge.1) write(*,'("LUFAC: level=",i1)') level
      if(level.ge.2) print*,"bandwidth=",nband," neqn=",neqn

      zero=0.0d0

       do 10 k=1,neqn
          kb = min(neqn-k,nband)
          do 20 i= 1, kb
             kpi = k+i
             pick= A(-i,kpi)
             if(pick.eq.zero) goto 20
             mult= pick/A(0,k)
             A(-i,kpi) = mult
             do j= 1, kb
                A(j-i,kpi)= A(j-i,kpi)-mult*A(j,k)
             end do
 20       continue
 10    continue

      return
      end