skyslv.f

关于网格自动生成delaunay算法的

************************************************************************
*
*     Modified Cholesky Decomposition
*     of skyline storage al matrix
************************************************************************
      subroutine skydec(msizl,n,al,idiag,ww,eps,ind)
      implicit real*8(a-h,o-z)
      dimension al(msizl),ww(n)
      dimension idiag(0:n)
      if(eps.lt.1.0d-20) then
         eps= 1.0d-20
      endif
 700  format(//5x,75('*'),
     $        /8x,'Message from subroutine skydec',
     $        /8x,'eps (absolute error tolerance) has not been set !',
     $        /8x,'default value eps=1.0d-20 is used',
     $        /5x,75('*'))
      if(n.lt.1) goto 420
*---------------------------------------
*     lu decomposition
*---------------------------------------
      do 20 k= 1,n
         ist= k-idiag(k)+idiag(k-1)+1
         iaki0= idiag(k)-k
         do 40 i= ist,k-1
            ww(i)= al(iaki0+i)
            jst= i-idiag(i)+idiag(i-1)+1
            ijmax= max(ist,jst)
            ilij0= idiag(i)-i
            do 30 j= ijmax,i-1
               ww(i)= ww(i)-ww(j)*al(ilij0+j)
 30         continue
            al(iaki0+i)= ww(i)*al(idiag(i))
 40      continue
         dd= al(idiag(k))
         do 50 i= ist,k-1
            dd= dd-ww(i)*al(iaki0+i)
 50      continue
         if(abs(dd).lt.eps) goto 410
         add= 1.0d0/dd
         al(idiag(k))= add
 20   continue
      ind= 0
      return
 410  ind= k
      write(7777,*) 'Error Message from subroutine skydec 1 !      '
      write(7777,*) 'Decomposition was finished in failure !'
      write(7777,*) 'Singular point was discovered !'
      write(7777,*) 'Degree of freedom =',n
      write(7777,*) 'Result cood = ',ind
      stop
 420  ind= 30000
      write(7777,*) 'Error Message from subroutine skydec 2 !      '
      write(7777,*) 'Irregular Input data is expected !'
      write(7777,*) 'Singular point was discovered !'
      write(7777,*) 'Degree of freedom =',n
      write(7777,*) 'Result cood = ',ind
      stop
 500  format(//5x,75('*'),
     $        /8x,'Message from subroutine skydec !',
     $        /8x,'Matrix has been decompsed in succeed !'
     $        /8x,'Degree of freedom =',i8,
     $         5x,'Result cood = ',i8,
     $        /5x,75('*'))  
      end
*
      subroutine skysol(msizl,n,al,idiag,x)
      implicit real*8(a-h,o-z)
      dimension al(msizl),x(n)
      dimension idiag(0:n)
      if(n.lt.1) goto 420
      do 100 k= 1,n
         ist= k-idiag(k)+idiag(k-1)+1
         iaki0= idiag(k)-k
         do 110 i= ist,k-1
            x(k)= x(k)-al(iaki0+i)*x(i)
 110     continue
 100  continue
      do 120 i= 1,n
         x(i)= al(idiag(i))*x(i)
 120  continue
      do 130 k= n,2,-1
         lst= k-idiag(k)+idiag(k-1)+1
         ilkl0= idiag(k)-k
         do 140 l= lst,k-1
            x(l)= x(l)-x(k)*al(ilkl0+l)
 140     continue
 130  continue
      return
 420  ind= 30000
      write(6,520) n,ind
 520  format(//5x,75('*'),
     $        /8x,'Message from subroutine skysol !',
     $        /8x,'Irregular Input data is expected !',
     $        /8x,'Degree of freedom =',i8,
     $         5x,'Result cood = ',i8,
     $        /5x,75('*'))  
      return
      end