ldlt.f90

这是一个用fortran编写的“第一类边界条件的有限元算法”

!***************************************************************************!
!                  SUBROUTINE LDLT(a,n,iw,p,ie)                             !
!      功能    对称带型线性方程组的系数矩阵的下三角部分被定带宽地储存在矩型 !
!              数组a中，利用a数组来解方程组。                               ! 
!      参数说明                                                             ! 
!        a-----n*iw的二维实数组,输入参数,存放对称带型线性方程组的系数矩阵的 !
!              下三角部分。                                                 !
!        n-----整型变量,输入参数,方程组的阶数。                             !
!        iw----整型变量,输入参数,对称带型线性方程组的半带宽。               !
!        p-----n个元素的一维实数组，输入、输出参数。开始存放方程组的右侧列  !
!              向量，工作结束时，存放解向量。                               !
!        ie----整型变量，输出参数，标志。ie=0时表示程序正常结束；           !
!              ie=1时，表示系数矩阵奇异。                                   !
!***************************************************************************! 
 SUBROUTINE LDLT(a,n,iw,p,ie)
   DIMENSION a(n,iw),p(n)

    do 15 i=1,n
	  if(i.le.iw) goto 20
	    it=i-iw+1
        goto 30
20	    it=1
30      k=i-1
        if(i.eq.1) goto 40
      do 25 l=it,k
	     il=l+iw-1
		 b=a(i,il)
		 a(i,il)=b/a(i,iw)
		 p(i)=p(i)-a(i,il)*p(l)
		 mi=l+1
      do  25 j=mi,i
	     ij=j+iw-i
		 jl=l+iw-j
25    a(i,ij)=a(i,ij)-a(j,jl)*b
40      if(a(i,iw).eq.0) goto 100
15   continue 
       do 45 j=1,n
	     if(j.le.iw) goto 60
		   it=n-j+iw
		   goto 70
60       it=n
70       i=n-j+1
         p(i)=p(i)/a(i,iw)
		 if(j.eq.1) go to 45
		   k=i+1
		   do 65 mj=k,it
		     ij=i-mj+iw
65         p(i)=p(i)-p(mj)*a(mj,ij)
45         continue 
           ie=0
		   go to 110
100        ie=1
110     return
end subroutine LDLT