mneig.f

Go to the documentation of this file.

*
* $Id: mneig.F,v 1.1.1.1 1996/03/07 14:31:29 mclareni Exp $
*
* $Log: mneig.F,v $
* Revision 1.1.1.1  1996/03/07 14:31:29  mclareni
* Minuit
*
*
      SUBROUTINE mneig(A,NDIMA,N,MITS,WORK,PRECIS,IFAULT)
*
* $Id: d506dp.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
*
* $Log: d506dp.inc,v $
* Revision 1.1.1.1  1996/03/07 14:31:32  mclareni
* Minuit
*
*
*
*
* d506dp.inc
*
C ************ DOUBLE PRECISION VERSION *************
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
C
      parameter(zero=0.0,  one=1.0,   two=2.0)
      parameter(tol=1.0e-35)
      dimension a(ndima,*),work(*)
C          PRECIS is the machine precision EPSMAC
      ifault = 1
C
      i = n
      DO 70 i1 = 2,n
      l = i-2
      f = a(i,i-1)
      gl = zero
C
      IF(l .LT. 1) GO TO 25
C
      DO 20 k = 1,l
   20 gl = gl+a(i,k)**2
   25 h = gl + f**2
C
      IF(gl .GT. tol) GO TO 30
C
      work(i) = zero
      work(n+i) = f
      GO TO 65
   30 l = l+1
C
      gl = sqrt(h)
C
      IF(f .GE. zero) gl = -gl
C
      work(n+i) = gl
      h = h-f*gl
      a(i,i-1) = f-gl
      f = zero
      DO 50 j = 1,l
      a(j,i) = a(i,j)/h
      gl = zero
      DO 40 k = 1,j
   40 gl = gl+a(j,k)*a(i,k)
C
      IF(j .GE. l) GO TO 47
C
      j1 = j+1
      DO 45 k = j1,l
   45 gl = gl+a(k,j)*a(i,k)
   47 work(n+j) = gl/h
      f = f+gl*a(j,i)
   50 CONTINUE
      hh = f/(h+h)
      DO 60 j = 1,l
      f = a(i,j)
      gl = work(n+j)-hh*f
      work(n+j) = gl
      DO 60 k = 1,j
      a(j,k) = a(j,k)-f*work(n+k)-gl*a(i,k)
   60 CONTINUE
      work(i) = h
   65 i = i-1
   70 CONTINUE
      work(1) = zero
      work(n+1) = zero
      DO 110 i = 1,n
      l = i-1
C
      IF(work(i) .EQ. zero .OR. l .EQ. 0) GO TO 100
C
      DO 90 j = 1,l
      gl = zero
      DO 80 k = 1,l
   80 gl = gl+a(i,k)*a(k,j)
      DO 90 k = 1,l
      a(k,j) = a(k,j)-gl*a(k,i)
   90 CONTINUE
  100 work(i) = a(i,i)
      a(i,i) = one
C
      IF(l .EQ. 0) GO TO 110
C
      DO 105 j = 1,l
      a(i,j) = zero
      a(j,i) = zero
  105 CONTINUE
  110 CONTINUE
C
C
      n1 = n-1
      DO 130 i = 2,n
      i0 = n+i-1
  130 work(i0) = work(i0+1)
      work(n+n) = zero
      b = zero
      f = zero
      DO 210 l = 1,n
      j = 0
      h = precis*(abs(work(l))+abs(work(n+l)))
C
      IF(b .LT. h) b = h
C
      DO 140 m1 = l,n
      m = m1
C
      IF(abs(work(n+m)) .LE. b) GO TO 150
C
  140 CONTINUE
C
  150 IF(m .EQ. l) GO TO 205
C
  160 IF(j .EQ. mits) RETURN
C
      j = j+1
      pt = (work(l+1)-work(l))/(two*work(n+l))
      r = sqrt(pt*pt+one)
      pr = pt+r
C
      IF(pt .LT. zero) pr=pt-r
C
      h = work(l)-work(n+l)/pr
      DO 170 i=l,n
  170 work(i) = work(i)-h
      f = f+h
      pt = work(m)
      c = one
      s = zero
      m1 = m-1
      i = m
      DO 200 i1 = l,m1
      j = i
      i = i-1
      gl = c*work(n+i)
      h = c*pt
C
      IF(abs(pt) .GE. abs(work(n+i))) GO TO 180
C
      c = pt/work(n+i)
      r = sqrt(c*c+one)
      work(n+j) = s*work(n+i)*r
      s = one/r
      c = c/r
      GO TO 190
  180 c = work(n+i)/pt
      r = sqrt(c*c+one)
      work(n+j) = s*pt*r
      s = c/r
      c = one/r
  190 pt = c*work(i)-s*gl
      work(j) = h+s*(c*gl+s*work(i))
      DO 200 k = 1,n
      h = a(k,j)
      a(k,j) = s*a(k,i)+c*h
      a(k,i) = c*a(k,i)-s*h
  200 CONTINUE
      work(n+l) = s*pt
      work(l) = c*pt
C
      IF(abs(work(n+l)) .GT. b) GO TO 160
C
  205 work(l) = work(l)+f
  210 CONTINUE
      DO 240 i=1,n1
      k = i
      pt = work(i)
      i1 = i+1
      DO 220 j = i1,n
C
      IF(work(j) .GE. pt) GO TO 220
C
      k = j
      pt = work(j)
  220 CONTINUE
C
      IF(k .EQ. i) GO TO 240
C
      work(k) = work(i)
      work(i) = pt
      DO 230 j=1,n
      pt = a(j,i)
      a(j,i) = a(j,k)
      a(j,k) = pt
  230 CONTINUE
  240 CONTINUE
      ifault = 0
C
      RETURN
      END