source/unres/quindet2ok.F90

   1        subroutine quindet2(c,b,d,n,mu,no)
   2         real (kind=8) :: mu ,x
   3         real (kind=8) ,dimension(n) :: c,b,d
   4         real (kind=8) ,dimension(3,5) :: r
   5         integer ::i,j,k,l,ll,rr,w,n,no
   6 ! the arrays c,b i d contain the diagonal subdiagonal and subsubdiagonal  elements of a simmetric quindiagonalmatrix. the value of the output parameter   no is the number of eigenvalues greater than mu;
   7         r(2,1)=c(1)-mu
   8         r(2,2)=b(2)
   9         r(3,1)=b(2)
  10         r(2,3)=d(3)
  11         r(3,2)=c(2)-mu
  12         r(3,3)=b(3)
  13         r(3,4)=d(4)
  14         r(2,4)=0
  15         r(2,5)=0
  16         r(3,5)=0
  17         if (r(2,1).ge.0) then
  18            no=1
  19         else
  20            no=0
  21         endif
  22 ! k count the major stages
  23         do k=2,n
  24           rr=1
  25           do i=2,1,-1
  26            if (k.gt.i) then
  27              w=3-i
  28 !interchange row w and row 3 if necessary
  29              if (dabs(r(3,1)).gt. dabs(r(w,1))) then
  30                do j=1,i+3
  31                  x=r(3,j)
  32                  r(3,j)=r(w,j)
  33                  r(w,j)=x
  34                enddo
  35                if (r(3,1).ge.0.eqv.r(w,1).ge.0) rr=-rr
  36              endif
  37 !elimination of subdiagonal elements
  38              if (r(w,1).eq.0) then
  39                x=0
  40              else
  41                x=r(3,1)/r(w,1)
  42              endif
  43              do j=2,i+3
  44                r(3,j-1)=r(3,j)-x*r(w,j)
  45              enddo
  46              r(3,i+3)=0
  47            endif
  48           enddo ! i
  49
  50           if (r(3,1).lt.0) rr=-rr
  51           if (rr.gt.0) no=no+1
  52 ! update elements of array r
  53           do i = 1,2
  54             do j = 1,5
  55               r(i,j)=r(i+1,j)
  56             enddo
  57           enddo
  58           if ((k+1).le.n) then
  59             r(3,1)=d(k+1)
  60             r(3,2)=b(k+1)
  61             r(3,3)=c(k+1)-mu
  62           endif
  63           if ((k+2).le.n) then
  64             r(3,4)=b(k+2)
  65           else
  66             r(3,4)=0
  67           endif
  68           if ((k+3).le.n) then
  69            r(3,5)=d(k+3)
  70           else
  71            r(3,5)=0
  72           endif
  73         enddo ! k
  74       endsubroutine