! 
!  W.Quapp, Feb 2010
!
! energy + gradient + Hessian from function mod nfk
! Programm bestimmt GE 
!
!     nn = Dimension
!     le = Chainlength
!   pstl = pred-Steplength
!    eps = Convergence in RGF corrector
!   rsel in corrector = tang
       INTEGER nn,le,iterNTmax
       DOUBLE PRECISION pstl
            PARAMETER (nn=2,le=24,iterNTmax=50,pstl=0.15d0)
       DOUBLE PRECISION kette(le+1,nn),eps,tgt,tt,energy,x,y,
     1  rsel(nn),grad(nn),H(nn,nn),Proj(nn,nn),
     2  xstart(nn),sw(nn),norm,dstep(nn),detk,t2(nn),t3(nn),
     4  projNm1(nn-1,nn),dk(nn,nn+1),rgeq(nn-1,nn),
     5  rgrad(nn-1),tang(nn),told(nn),unit(nn,nn),vec_ddot,det,
     7  t1old(nn),t1(nn),rgnorm,point(nn),eival(nn),eivec(nn,nn)
      INTEGER k,j,i,jnr,istatus,icalc,iterNT
      logical omat_inv
ccccccccccccccccccccccccccccccccccccccccc      
      OPEN(8,FILE='dim2kette.dat')
      OPEN(7,FILE='chain.dat')
      OPEN(11,FILE='dim2dir.dat')
      OPEN(53,FILE='dim2Eall.dat')
      OPEN(25,FILE='dim2told.dat')
      OPEN(26,FILE='dim2t1old.dat')
c protocol file: all interesting output of the vri-search
      OPEN(44,FILE='dim2protocol.txt')
cccccccccccccccccccccccccccccccccccccccc
2     FORMAT(2F20.14)
50    FORMAT(1X,3F20.10)
      N=nn
      jnr=1
      iterNT=0
      istatus=2
      Write(44,*)
      write(44,*)' Search TASC solution'
      write(*,*) ' Search TASC solution'
      eps=1.0d-10
c  eps is criterion for rgf convergence to given rsel 
c  jnr is the currend node number
       write(*,*) 'eps = ',eps
       write(44,*)'eps = ',eps
       xstart(1)= 2.85d0
       xstart(2)= -0.1125d0
       do i=1,N
        point(i)=xstart(i)
        kette(1,i)=xstart(i)
       enddo
c other examples: (0.10, -0.10),    (1.0,  -1.30),    (1.80, -2.20)       
       WRITE(44,50)(xstart(I),I=1,N)
       WRITE(*,50) (xstart(I),I=1,N)
       WRITE(44,*) 'ketten-point start is '
       WRITE(44,*) '{', xstart(1), ',', xstart(2) ,' },'
       WRITE(7,*)  '{'
       WRITE(7,*)  '{', xstart(1), ',', xstart(2) ,' },'
       call mat_diag(unit,nn,1.0d0)
       rsel(1)=-1.0d0
       rsel(2)=-1.25d0
       tgt=norm(rsel,nn)
       call vec_dmult_constant(rsel,nn,1.0d0/tgt,rsel)
       write(44,*) 'r_selected '
       write(44,2)(rsel(i),i=1,nn)
       call vec_dmult_constant(rsel,nn,1.0d0,told)
       call vec_dmult_constant(rsel,nn,1.0d0,t1old)
       call vec_dmult_constant(rsel,nn,1.0d0,t2) 
       call vec_dmult_constant(rsel,nn,1.0d0,t3)
       goto 11 
c cccccccccccccccccccccccccccccccccccccccccccccccc
1     continue
c cccccccccccccccccccccccccccccccccccccccccccccccc
        do i=1,nn
         point(i)=kette(jnr,i)
        enddo
        write(44,*) ' point '
        WRITE(44,2) (point(I),I=1,nn)
11        x=point(1)
          y=point(2)
c  0.03- modified NFK test surface, see CPL 2004 
        energy=-9.0d0*Dexp(-(-3.0d0+x)**2 - y**2)
     1         -9.0d0*Dexp(-( 3.0d0+x)**2 - y**2)
     2          +  x*y +0.03d0*(x**2 + y**2)**2  
        write(44,*)'Energy ',energy
        write(53,*) energy
c        rewind 14
c        Read(14,*) (grad(i),i=1,nn)
        grad(1)=18.0d0*(Dexp(-(-3.0d0+x)**2-y**2)*(-3.+x) +
     1                  Dexp(-( 3.0d0+x)**2-y**2)*( 3.+x))+ y+
     2                  0.12d0*x*(x**2 + y**2)
        grad(2)= 
     1   18.0d0*y*(Dexp(-(-3.0d0+x)**2-y**2)+Dexp(-(3.0d0+x)**2-y**2))
     1   + x +  0.12d0*(x**2 + y**2)*y
c        rewind 15
c        Do 15  k=1,nn
c15      Read(15,*) (H(k,i),i=1,nn)
        h(1,1)=
      1  18.0d0*(Dexp(-(-3.0d0+x)**2-y**2)+Dexp(-(3.0d0+x)**2-y**2))-
      2  36.0d0*(Dexp(-(-3.0d0+x)**2-y**2)*(-3.0d0+ x)+
      3          Dexp(-( 3.0d0+x)**2 -y**2)*(3.0d0+ x)) +
      4          0.12d0*(3.0d0*x**2 + y**2)
        h(1,2)=
      1   - 36.0d0*Dexp(-(-3.0d0+x)**2-y**2)*(-3.0d0+x)*y+
      2     1.d0 + 0.24d0*x*y
      3   - 36.0d0*Dexp(-( 3.0d0+x)**2-y**2)*( 3.0d0+x)*y
        h(2,2)=
      1  18.0d0*(Dexp(-(-3.0d0+x)**2 -y**2)+Dexp(-(3.0d0+x)**2-y**2))-
      2  36.d0*y**2*(Dexp(-(-3.d0+x)**2-y**2)+Dexp(-(3.d0+x)**2-y**2))+
      3   0.12d0*( x**2 + 3.0d0*y**2)
        h(2,1)= h(1,2)                                                               
       write(44,*) ' Hessian  '
       write(44,2) ((H(I,K),I=1,nn),K=1,nn)
       detHe=det(H,nn)
       write(44,*) 'det of  Hessian', detHe
       write(*,*)  'det of  Hessian', detHe
c                                 use Proj for a help matrix
       call eigen(H,nn,nn,eival,eivec,Proj,0)
       write(44,*)'eigenvalues  of H'
       write(44,*)(eival(i),i=1,nn)
       write(44,*)'eigenvectors of H'
       do k=1,nn
          write(44,*) 'vec(',k,') =   ',(eivec(k,i),i=1,nn)
       enddo
       write(44,*) ' grad '
       WRITE(44,2) (grad(I),I=1,nn)
       rgnorm=norm(grad,nn)
       WRITE(44,*) ' full gradient norm ',rgnorm
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc  new:  projNm1: n-1 linearly independent lines from Proj
        call ortcompl(rsel,projNm1,unit,nn)
        write(44,*) ' projNm1 '
        Do 13, k=1,nn-1
13      WRITE(44,2) (projNm1(k,I),I=1,nn)
cccccccccccccccccccccccccccccccccccccccccccccccccccc                     
c test print :  is o.k. 
      if(jnr.eq.1) then
        write(44,*) '  rsel '
        WRITE(44,2) (rsel(I),I=1,nn)
        call vec_dinit(sw,nn,0.0d0)
        call matmult(projNm1,nn-1,rsel,1,sw,nn)
        write(44,*) ' Proj * rsel '
        WRITE(44,*) (sw(I),I=1,nn)
        call vec_dinit(sw,nn,0.0d0)
        call matmult(projNm1,nn-1,grad,1,sw,nn)
        write(44,*) ' Proj * grad '
        WRITE(44,*) (sw(I),I=1,nn)
      endif
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc                 
        call matmult(projNm1,nn-1,H,nn,rgeq,nn)
        write(44,*) 'print  rgeq  '
        Do 14, k=1,nn-1
14       WRITE(44,2) (rgeq(k,I),I=1,nn)
        call vec_dinit(rgrad,nn-1,0.0d0)
        call matmult(projNm1,nn-1,grad,1,rgrad,nn)
        write(44,*) ' rgrad '
        WRITE(44,2) (rgrad(I),I=1,nn-1)
        rgnorm=norm(rgrad,nn-1)
        WRITE(44,*) ' PES reduced grad norm  ',rgnorm
c Tangent ccccccccccccccccccccccccccccccccccccccccccccc
         call gettang(rgeq,tang,nn)
         write(44,*) 'tang from rgeq is '
         WRITE(44,2) (tang(I),I=1,nn)
c         ccccccccccccccccccccccccccccccccccccccccccccc
        call matmult(rgeq,nn-1,tang,1,sw,nn)
        write(44,*) ' test: rgeq*tang gives '
        write(44,2) (sw(I),I=1,nn)
        tgt=vec_ddot(tang,told,nn)
        write(44,*)' tang * told     ',tgt
        if (tgt.lt.0.0d0)
     1   call vec_dmult_constant(tang,nn,-1.0d0,tang)
cccccccccccccccccccccccccccccccccccccccccccccc
        call buildk(dk,rgeq,tang,rgrad,nn)
cccccccccccccccccccccccccccccccccccccccccccccc
        write(44,*) ' K mat is '
        Do 71, k=1,nn
71        WRITE(44,*) (dk(k,I),I=1,nn+1)
        detk=det(dk,nn)
        write(*,*)  ' detK is       ', detk
        write(44,*) ' detK is       ', detk
        if (detk.lt.0.0d0) then
           call vec_dmult_constant(tang,nn,-1.0d0,t1)
        else
           call vec_dcopy(tang,t1,nn)
        endif
          tg=vec_ddot(t1,t1old,nn)
          if (tg.lt.0.0d0 .and. jnr.gt. le/2 ) then 
             write(*,*)  ' bif of branches'
             write(44,*) ' bif of branches'
c  to jump over turning points of NTs
            call vec_dmult_constant(tang,nn,-1.0d0,tang)
            call buildk(dk,rgeq,tang,rgrad,nn)
          endif
          call linsolve(dk,dstep,nn)
          WRITE(44,*) 'hypo Corr - dstep by linsolve '
          WRITE(44,2)(dstep(i),i=1,nn)
ccccccc 1.order goto:
333      continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccc Abbruch 
        IF(
     1     rgnorm.lt.eps              .or.  ! Corrector o.k.
     3     iterNT.ge.iterNTmax        ) THEN
            write(44,*)
            write(*,*)
            write(44,*)' Predictor step ----------------  pstl=',pstl
            write(*,*) ' Predictor step ----------------  pstl=',pstl
            write(44,*) 'Nr j in vri2 is =',jnr
            write(*,*) '  |rgnorm|           ', rgnorm
            write(44,*) ' |rgnorm|           ', rgnorm
            iterNT=1
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
           if(jnr.lt.le+1) THEN
cc        c                               do predictorstep
c   this line is TASC !!! ccccccccccccccccccc
            call vec_dcopy(tang,rsel,nn)
c   ccccccccccccccccccccccccccccccccccccccccc           
            do i=1,nn 
              point(i)=kette(jnr,i) 
            enddo
cc predictor step by rgf3
c            dk(nn,nn+1)=pstl
c            call linsolve(dk,dstep,nn)
c            write(44,*)'pred-step by rgf3 '
c  or:
           call vec_dmult_constant(tang,nn,pstl,dstep)
           write(44,*)'pred-step by tangent '  
            tgt=norm(dstep,nn)
            write(*,*)
            write(44,*)'pred-step length would be',tgt
            write(44,2)(dstep(i),i=1,nn)
            do i=1,nn
              point(i)=point(i)+dstep(i)*pstl/tgt
            enddo
            kette(jnr+1,1)=point(1)
            kette(jnr+1,2)=point(2)
       WRITE(44,*) 'ketten-point nach predictor is '
       WRITE(44,22) kette(jnr+1,1), kette(jnr+1,2) 
       WRITE(7,22)  kette(jnr+1,1), kette(jnr+1,2) 
22     Format(4H {  ,F20.14,3H , ,F20.14,5H   },)
            jnr=jnr+1
            call vec_dcopy(t2,t3,nn)
            call vec_dcopy(tang,t2,nn)
            GOTO 66
c noch Abbruch                
           else
            istatus=10 
            goto 66
           ENDIF
ccc AbbruchEnde 
        ENDIF
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccCCCCCCCCCCCCCCCCCC
        tt=norm(dstep,nn)
        write(44,*) ' dstep norm whould be ',tt
        write(44,*) ' dstep is '
ccccc Notbremse  --  emergency break : 
c        if(tt.gt.pstl)
c     1     call vec_dmult_constant(dstep,nn,pstl/tt,dstep)
        write(44,2) (dstep(I),I=1,nn)
        DO 25 i=1,nn
            point(i)=point(i)+dstep(i)
25      CONTINUE
        do i=1,nn
             kette(jnr,i)=point(i)
        enddo 
        WRITE(44,*) 'ketten-point nach dstep addition'
        WRITE(44,22) kette(jnr,1),kette(jnr,2) 
        WRITE(7,22)  kette(jnr,1),kette(jnr,2)
        write(44,*) 'Nr j in vri2 is =',jnr
        write(*,*)  'j vri2 = ',jnr,' iterNT  =  ',iterNT
        write(*,*) ' hypo |dstep|            ', tt,'   rg ',rgnorm
        write(44,*) 'hypo |dstep|            ', tt,'   rg ',rgnorm
        iterNT=iterNT+1
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
66      continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
         WRITE(44,*) jnr
         WRITE(44,2) (kette(jnr,i),i=1,nn) 
         rewind 8
         DO 29 k=1,le+1
29        WRITE(8,2) (kette(k,i),i=1,nn)
c
        if(istatus.lt.10) goto 1
c        
        WRITE(7,23)  kette(jnr,1),kette(jnr,2)
23      Format(4H {  ,F20.14,3H , ,F20.14,5H   }})  
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccc                                       cccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine ginvers(in,NAT,out,omat_inv)
c    calculates the inverse matrix of IN
      real*8 in,out
      logical omat_inv
      integer NAT,i,j,k
      dimension in(NAT,NAT),out(NAT,NAT)
      do 1,k=1,NAT
      do 1,j=1,NAT
  1   out(j,k)=in(j,k)
c
      do 10,k=1,NAT
         if (out(k,k).eq.0.0d0) then
            omat_inv=.false.
            return
         endif
         out(k,k)=1.0d0/out(k,k)
         do 20,i=1,NAT
             if (i.ne.k) out(i,k)=-(out(i,k)*out(k,k))
   20    continue
         do 30,j=1,NAT
             if (j.ne.k) then
                do 40,i=1,NAT
                      if (i.ne.k) out(i,j)=out(i,j)+out(i,k)*out(k,j)
   40           continue
                out(k,j)=out(k,j)*out(k,k)
             endif
   30    continue
   10 continue
      omat_inv=.true.
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine mat_diag(d,n,v)
c create a NxN diagonal matrix D with entries V
       real*8  d,v
       integer n,i
       dimension d(n,n)
       if(n .le. 0) return
       call vec_dinit(d,n*n,0.0d0)
       do i=1,n
          d(i,i)=v
       end do
       return
       end
cccccccccccccccccccccccccccccccc
!
! Multiplikation Ax=b
!
      SUBROUTINE multiply(b,A,x,nn)
      DOUBLE PRECISION b(nn),A(nn,nn),x(nn)
      DO 20 i=1,nn
        b(i)=0
        DO 10 j=1,nn        
           b(i)=b(i)+A(i,j)*x(j)
10      CONTINUE
20    CONTINUE
      END
!
! Gibt die Norm von v aus
!
      DOUBLE PRECISION FUNCTION norm(v,nn)
      DOUBLE PRECISION v(nn),s
      s=0.d0
      DO 10 j=1,nn
        s=s+v(j)**2
10    CONTINUE
      IF (s.LT.1.D-23) THEN
        norm=0.d0
      ELSE
        norm=DSQRT(s)
      ENDIF
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
       subroutine ortcompl(vec,compl,dmet,ndim)
c   calculate an Orthonormalcomplement COMPL to vector VEC
       real*8 basis,vec,compl,swap,det,dmgs,dmet
       integer i,j,ndim
       dimension vec(ndim),compl(ndim-1,ndim),basis(ndim,ndim),
     *   swap(ndim),dmet(ndim,ndim)
       call vec_dinit(basis,ndim*ndim,0.0d0)
       do 10,i=1,ndim
         basis(i,i)=1.0d0
 10    continue
       do 20,i=1,ndim
         if (vec(i).ne.0.0d0) then
           do 30,j=1,ndim
             swap(j)=basis(j,1)
             basis(j,1)=vec(j)
             if (i.ne.1) basis(j,i)=swap(j)
 30        continue
           goto 40
         endif
 20    continue
 40    continue
c  note: dmgs is the Gram-Schmidt method 
       det=dmgs(basis,dmet,ndim)
       call vec_dinit(compl,(ndim-1)*ndim,0.0d0)
       call mat_trans(basis,ndim,ndim,basis)
       do 60,i=2,ndim
          call vec_mrowcopy(i,basis,ndim,i-1,compl,ndim-1,ndim)
 60    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function dmgs(a,dmet,dim)
c  modified Gram-Schmidt algorithm return det(a)
       real*8 det,rr,a,q,r,vec_mmdot,dmet,dzero
       integer dim,i,k
       dimension a(dim,dim),q(dim,dim),r(dim,dim),dmet(dim,dim)
       data dzero/0.0d0/
       call vec_dinit(q,dim*dim,0.0d0)
       call vec_dinit(r,dim*dim,0.0d0)
       det=1.0d0
       call mat_dcopy(a,q,dim,dim)
       do 10,k=1,dim
         if (k.eq.1) goto 21
         do 20,i=1,k-1
            r(i,k)=vec_mmdot(i,q,dim,k,q,dim,dim,dmet)
            call vec_madd(k,q,dim,1.0d0,i,q,dim,-r(i,k),k,q,dim,dim)
 20      continue
 21      continue
         rr=vec_mmdot(k,q,dim,k,q,dim,dim,dmet)
         r(k,k)=dsqrt(rr)
         det=det*r(k,k)
         if(rr.lt.dzero) goto 33
         call vec_madd(k,q,dim,1.0d0/r(k,k),k,q,dim,0.0d0,k,q,dim,dim)
 10    continue
       call mat_dcopy(q,a,dim,dim)
       dmgs=det
       return
 33    CONTINUE
cc       call cerr(-23)
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function vec_mmdot
     &(nri,matin,din,nra,matout,dout,dim,dmet)
c     multiply the (nri) col of (matin) and the (nra) col of
c     (matout). (dim) is the number of the rows with metric dmet
       real*8 matin,matout,prod,v1,v2,dmet,vec_mdot
       integer nri,din,nra,dout,dim
       dimension matin(dim,din),matout(dim,dout),dmet(dim,dim),
     +   v1(dim),v2(dim)
ccc      if ((nri.gt.din).or.(nra.gt.dout)) call cerr(-100)
        call vec_mat_centry(v1,matin,nri,din,dim,'o')
        call vec_mat_centry(v2,matout,nra,dout,dim,'o')
        prod=vec_mdot(v1,dmet,v2,dim)
         vec_mmdot=prod
        return
        end
cccccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function vec_mdot(vec1,metric,vec2,dim)
c  metric scalar product of VEC1 and VEC2 with respect to METRIC
       integer dim
       real*8 vec1,metric,vec2,cov2,vec_ddot
       dimension vec1(dim),vec2(dim),cov2(dim),metric(dim,dim)
           call matmult(metric,dim,vec2,1,cov2,dim)
           vec_mdot=vec_ddot(vec1,cov2,dim)
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine gettang(curveq,tangent,nn)
c   This routine calculate the tangent on a curve
c   whereas the tangent is unique determined by:
c
c   H*t=0
c   |t|=1
c   det[H,t]>0
c
      real*8 tangent,curveq,q,ss1,det,zero,norm
      integer nn,i,j
      dimension tangent(nn),curveq(nn-1,nn),q(nn*nn),ss1(nn,nn)
      dimension gmat(nn,nn)
      data zero /0.0d0/
      call qr(curveq,q,nn-1)
      call vec_dcopy(q(nn*(nn-1)+1),tangent(1),nn)
c test     do 1,i=1,nn
c  1   tangent(i)=q(nn,i)
c wq
       call vec_dmult_constant
     +   (tangent,nn,1.0d0/norm(tangent,nn),tangent)
cc       call mat_dcopy(curveq,ss1,nn-1,nn)
       do 10,j = 1,nn-1
       do 20,i = 1,nn
          ss1(j,i) = curveq(j,i)
  20    continue
  10    continue
       do 30,i = 1,nn
  30      ss1(nn,i)=tangent(i)
       if (det(ss1,nn).lt.zero) then
        call vec_dmult_constant(tangent,nn,-1.0d0,tangent)
      endif
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine qr(gex,q,n)
C      QR-Zerlegung von gex
       integer i,j,k,n
       real*8 s1,s2,s,gex,q,r,qq,rr
       dimension gex(n,n+1),q(n+1,n+1),r(n+1,n),
     * qq(n+1,n+1),rr(n+1,n)
       do 10,i=1,n
          do 20,j=1,n+1
             r(j,i)=gex(i,j)
 20       continue
 10    continue
       do 30,i=1,n+1
          do 40,j=1,n+1
             q(i,j)=0.0d0
 40       continue
          q(i,i)=1.0d0
 30    continue
       do 50,i=1,n
          do 60,j=i+1,n+1
             s1=r(i,i)
             s2=r(j,i)
             s=dsqrt(s1*s1+s2*s2)
             if (s.gt.0.0d0) then
                s1=s1/s
                s2=s2/s
                else
                s1=0.0d0
                s2=0.0d0
                endif
             call vec_dcopy(r,rr,(n+1)*n)
             do 70,k=1,n
                rr(i,k)=s1*r(i,k)+s2*r(j,k)
                rr(j,k)=s1*r(j,k)-s2*r(i,k)
 70          continue
             call vec_dcopy(rr,r,(n+1)*n)
             call vec_dcopy(q,qq,(n+1)*(n+1))
             do 80,k=1,n+1
                qq(i,k)=s1*q(i,k)+s2*q(j,k)
                qq(j,k)=s1*q(j,k)-s2*q(i,k)
 80          continue
             call vec_dcopy(qq,q,(n+1)*(n+1))
 60       continue
 50    continue
       do 90,i=1,n+1
          do 100,j=1,n+1
             q(i,j)=qq(j,i)
 100      continue
 90    continue
c       write(44,*) 'in QR '
c       do i=1,n+1
c         write(44,*) (q(i,j),j=1,n+1)
c       enddo
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine tri_mat(dm1,nr1,nc1,dm2,nc2,dm3,nc3,dres)
c     dm1[nr1,nc1]*dm2[nc1,nc2]*dm3[nc2,nc3]=dres[nr1,nc3]
      real*8 dm1,dm2,dm3,dres,r1
      integer nr1,nc1,nc2,nc3
      dimension dm1(nr1,nc1),dm2(nc1,nc2),dm3(nc2,nc3),
     +          r1(nr1,nc2),dres(nr1,nc3)
      call matmult(dm1,nr1,dm2,nc2,r1,nc1)
      call matmult(r1,nr1,dm3,nc3,dres,nc2)
      return
      end
ccccccccccccccc
       subroutine mat_dcopy(arra1,arra2,ndim,mdim)
c  Copies matrix ARRA1 to matrix ARRA2
       real*8  arra1,arra2
       integer i,j,ndim,mdim
       dimension arra1(ndim,mdim),arra2(ndim,mdim)
       do 10,i = 1,ndim
       do 20,j = 1,mdim
        arra2(i,j) = arra1(i,j)
 20    continue
 10    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine buildk(dk,req,tang,rgr,nn)
c build the K-Matrix
      real*8 dk,req,tang,rgr
      integer nn,i,j
      dimension dk(nn,nn+1),req(nn-1,nn),tang(nn),rgr(nn-1)
      do 10,i=1,nn-1
        do 11,j=1,nn
        dk(i,j)=req(i,j)
 11   continue
 10   continue
      do 20,j=1,nn
         dk(nn,j)=tang(j)
 20   continue
      do 30,i=1,nn-1
         dk(i,nn+1)=-rgr(i)
 30   continue
      dk(nn,nn+1)=0.0d0
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
       subroutine linsolve(matein,x,n)
c  solves MATEIN * X = 0
       real*8 matein,dre,x,w
       integer n
       dimension matein(n,n+1),dre(n,n+1),x(n),w(n)
       call householder(matein,dre,w,n,n+1)
c       write(44,*) 'in linsolve '
c       do k=1,n
c         write(44,*) (dre(k,i),i=1,n+1)
c       enddo
       call loesung(dre,x,n)
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_dcopy(arra1,arra2,ndim)
c  Copies ARRA1 to ARRA2
      real*8 arra1,arra2
      integer i,ndim
      dimension arra1(ndim),arra2(ndim)
      if(ndim .le. 0) return
      do 10 i = 1,ndim
           arra2(i) = arra1(i)
 10   continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine householder(matein,mataus,w,n,sp)
c    householder transformation
       integer k,l,m,n,sp
       real*8 alpha,rho,sk,s,w,matein,mataus
       dimension matein(n,sp),mataus(n,sp),w(n)
       call vec_dcopy(matein,mataus,n*sp)
       do 10,k=1,n-1
        s=0.0d0
        do 20,l=k,n
         s=s+mataus(l,k)*mataus(l,k)
 20     continue
        if (mataus(k,k).lt.0) then
         alpha=dsqrt(s)
        else
         alpha=-dsqrt(s)
        endif
        rho=dsqrt(s-mataus(k,k)*mataus(k,k)+
     *       (mataus(k,k)-alpha)*(mataus(k,k)-alpha))
        if(rho .lt. 1.0d-63) then
         rho= 1.0d-63
            w(k)= 0.0d0
            goto 22
        endif
        w(k)=(mataus(k,k)-alpha)/rho
 22     continue
        do 30,l=k+1,n
         w(l)=mataus(l,k)/rho
 30     continue
        do 40,l=k,sp
         sk=0.0d0
         do 55,m=k,n
          sk=sk+w(m)*mataus(m,l)
 55      continue
         do 60,m=k,n
           mataus(m,l)=mataus(m,l)-2*sk*w(m)
 60      continue
 40     continue
 10    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccc
c
       subroutine loesung(matein,x,n)
c   intern linear equation routine
       real*8 matein,x
       integer i,j,n,k
       dimension matein(n,n+1),x(n)
       do 10,i=1,n
 10     x(i)=matein(i,n+1)
       x(n)=x(n)/matein(n,n)
       do 20,i=1,n-1
         k=n-i
         do 30,j=k+1,n
 30        x(k)=x(k)-matein(k,j)*x(j)
           x(k)=x(k)/matein(k,k)
 20    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_mat_centry(vec,mat,npos,ndim,ncol,zkind)
c Copies VEC 'i'n or 'o'ut the NPOS column of matrix MAT
        real*8 vec,mat
        integer npos,ndim,ncol,i
        character zkind
        dimension vec(ndim),mat(ndim,ncol)
        if(npos.gt.ncol) return
        if (zkind.eq.'i') then
         do i=1,ndim
           mat(i,npos)=vec(i)
         end do
        endif
        if (zkind.eq.'o') then
         do i=1,ndim
           vec(i)=mat(i,npos)
         enddo
        endif
        return
        end
cccccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function vec_ddot(arra1,arra2,ndim)
c   cartesian scalar product of ARRA1 and ARRA2
       real*8  arra1,arra2
       integer ndim,i
       dimension arra1(ndim),arra2(ndim)
       vec_ddot = 0
       if(ndim .le. 0)return
       do 70 i = 1,ndim
         vec_ddot = vec_ddot + arra1(i) * arra2(i)
  70   continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_dmult_add(arra1,arra2,ndim,zahl,arra4)
c  For each Komponent: ARRA4 = ARRA1 + ZAHL*ARRA2
      real*8  arra1,arra2,arra4,zahl
      integer ndim,i
      dimension arra1(ndim),arra2(ndim),arra4(ndim)
        if(ndim .le. 0)return
           do 145 i = 1,ndim
             arra4(i) = arra1(i)+arra2(i)*zahl
 145       continue
        return
        end       
ccccccccccccccccccccccccccccccccccccccc
       subroutine matmult(m1,d1,m2,d2,mout,dim)
c  matrix multiplication:  m1[d1,dim]*m2[dim,d2]=mout[d1,d2]
       real*8 m1,m2,mout,swap
       integer d1,d2,dim,i,j,k
       dimension m1(d1,dim),m2(dim,d2),mout(d1,d2),swap(d1,d2)
       call vec_dinit(swap,d1*d2,0.0d0)
       do 10,i=1,d1
        do 20,j=1,d2
         do 30,k=1,dim
           swap(i,j)=swap(i,j)+m1(i,k)*m2(k,j)
 30      continue
 20     continue
 10    continue
       call vec_dcopy(swap,mout,d1*d2)
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_dinit(arra1,ndim,wert)
c Set all elements of ARRA1 to the double value of WERT
       real*8   arra1,wert
       integer  ndim,i
       dimension arra1(ndim)
       if(ndim .le. 0) return
       do i=1,ndim
         arra1(i)=wert
       end do
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mat_trans(matin,nrow,ncol,matout)
c     transpose MATOUT of matrix MATIN
      real*8 matin,matout,swap
      integer nrow,ncol,i,j
      dimension matin(nrow,ncol),matout(ncol,nrow),swap(ncol,nrow)
      do 10,i=1,nrow
        do 20,j=1,ncol
                swap(j,i)=matin(i,j)
   20   continue
   10 continue
      call mat_dcopy(swap,matout,ncol,nrow)
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_mrowcopy(nri,matin,din,nra,matout,dout,dim)
c  copies the (nri) row of (matin) with (din) row and (dim) col
c  to the (nra) row of (matout) with (dout) row and (dim) col
      real*8 matin,matout
        integer nri,din,nra,dout,dim,i
        dimension matin(din,dim),matout(dout,dim)
ccc      if ((nri.gt.din).or.(nra.gt.dout)) call cerr(17)
      do 10,i=1,dim
        matout(nra,i)=matin(nri,i)
   10 continue
        return
        end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_madd
     * (nr1,mat1,d1,fact1,nr2,mat2,d2,fact2,nrx,matx,dx,dim)
c     matx(nrx)=fact1*mat1(nr1)+fact2*mat2(nr2)
c     d1,d2,dx - columns
c     dim      - rows
       real*8 mat1,mat2,matx,fact1,fact2,vec
       integer nr1,nr2,nrx,d1,d2,dx,dim,i
       dimension mat1(dim,d1),mat2(dim,d2),matx(dim,dx),vec(dim,1)
       call vec_dinit(vec,dim,0.0d0)
       do 10,i=1,dim
          vec(i,1)=fact1*mat1(i,nr1)+fact2*mat2(i,nr2)
 10    continue
       call vec_mat_centry(vec,matx,nrx,dim,dx,'i')
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_dmult_constant(arra1,ndim,zahl,arra2)
c   ARRA2 = ZAHL * ARRA1
      real*8  arra1,arra2,zahl
      integer ndim,i
      dimension arra1(ndim),arra2(ndim)
      if(ndim .le. 0)return
      do i = 1,ndim
         arra2(i) = arra1(i)*zahl
      end do
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function det(matein,nn)
c   calculation of the determinant of MATEIN
       integer i,nn
       real*8 matein,mataus,w
       dimension matein(nn,nn),mataus(nn,nn),w(nn)
       call householder(matein,mataus,w,nn,nn)
       det=1.0d0
       do 10,i=1,nn
       if (dabs(mataus(i,i)).lt.1.0d-63) then
          det= 1.0d-63
          goto 20
       endif
       if (dabs(det) .lt.1.0d-63) then
          det= 1.0d-63
          goto 20
       endif
 10    det=det*mataus(i,i)
       if (iand(nn,1).eq.0) det=-det 
 20    continue
       return
       end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine eigen (a,m,n,d,vec,e,iff)
c  eigenvalues and -vectors of A      
      real*8 a,d,vec,e,eps,tol,h,g,s,f,b,p,r,c
      integer m,iff,n,i,j,ni,ii,l,k,j1
      dimension a(m,m), d(m), vec(m,m)
      dimension e(m)
      eps = 0.278d-16
      tol = 0.211758d-21
      if (n.eq.1) then
         d(1) = a(1,1)
         vec(1,1) = 1.0d0
      else
         do 30 i = 1 , n
c     householder's reduction
c     simulation of loop do 150 i=n,2,(-1)
            do 20 j = 1 , i
               vec(i,j) = a(i,j)
 20         continue
 30      continue
         do 120 ni = 2 , n
            ii = n + 2 - ni
            do 110 i = ii , ii
               l = i - 2
               h = 0.0d0
               g = vec(i,i-1)
               if (l.gt.0) then
                  do 40 k = 1 , l
                     h = h + vec(i,k)**2
 40               continue
                  s = h + g*g
                  if (s.lt.tol) then
                     h = 0.0d0
                  else if (h.gt.0) then
                     l = l + 1
                     f = g
                     g = dsqrt(s)
                     if (f.gt.0) then
                        g = -g
                     end if
                     h = s - f*g
                     vec(i,i-1) = f - g
                     f = 0.0d0
                     do 70 j = 1 , l
                        vec(j,i) = vec(i,j)/h
                        s = 0.0d0
                        do 50 k = 1 , j
                           s = s + vec(j,k)*vec(i,k)
 50                     continue
                        j1 = j + 1
                        if (j1.le.l) then
                           do 60 k = j1 , l
                              s = s + vec(k,j)*vec(i,k)
 60                        continue
                        end if
                        e(j) = s/h
                        f = f + s*vec(j,i)
 70                  continue
                     f = f/(h+h)
                     do 80 j = 1 , l
                        e(j) = e(j) - f*vec(i,j)
 80                  continue
                     do 100 j = 1 , l
                        f = vec(i,j)
                        s = e(j)
                        do 90 k = 1 , j
                           vec(j,k) = vec(j,k) - f*e(k) - vec(i,k)*s
 90                     continue
 100                 continue
                  end if
               end if
c     accumulation of transformation matrices
               d(i) = h
               e(i-1) = g
 110        continue
 120     continue
         d(1) = vec(1,1)
         vec(1,1) = 1.0d0
         do 170 i = 2 , n
            l = i - 1
            if (d(i).gt.0) then
               do 150 j = 1 , l
                  s = 0.0d0
                  do 130 k = 1 , l
                     s = s + vec(i,k)*vec(k,j)
 130              continue
                  do 140 k = 1 , l
                     vec(k,j) = vec(k,j) - s*vec(k,i)
 140              continue
 150           continue
            end if
            d(i) = vec(i,i)
            vec(i,i) = 1.0d0
            do 160 j = 1 , l
c     diagonalization of the tridiagonal matrix
               vec(i,j) = 0.0d0
               vec(j,i) = 0.0d0
 160        continue
 170     continue
         b = 0.0d0
         f = 0.0d0
         e(n) = 0.0d0
         do 250 l = 1 , n
c     test for splitting
            h = eps*(dabs(d(l))+dabs(e(l)))
            if (h.gt.b) b = h
            do 180 j = l , n
c     test for convergence
               if (dabs(e(j)).le.b) go to 190
 180        continue
 190        if (j.eq.l) then
               d(l) = d(l) + f
            else
 200           p = (d(l+1)-d(l))*0.5d0/e(l)
               r = dsqrt(p*p+1.0d0)
               if (p.lt.0) then
                  p = p - r
               else
                  p = p + r
               end if
               h = d(l) - e(l)/p
               do 210 i = l , n
c     qr transformation
                  d(i) = d(i) - h
 210           continue
               f = f + h
               p = d(j)
c     simulation of loop do 330 i=j-1,l,(-1)
               c = 1.0d0
               s = 0.0d0
               j1 = j - 1
               do 240 ni = l , j1
                  ii = l + j1 - ni
                  do 230 i = ii , ii
c     protection against underflow of exponents
                     g = c*e(i)
                     h = c*p
                     if (dabs(p).lt.dabs(e(i))) then
                        c = p/e(i)
                        r = dsqrt(c*c+1.0d0)
                        e(i+1) = s*e(i)*r
                        s = 1.0d0/r
                        c = c/r
                     else
                        c = e(i)/p
                        r = dsqrt(c*c+1.0d0)
                        e(i+1) = s*p*r
                        s = c/r
                        c = 1.0d0/r
                     end if
                     p = c*d(i) - s*g
                     d(i+1) = h + s*(c*g+s*d(i))
                     do 220 k = 1 , n
                        h = vec(k,i+1)
                        vec(k,i+1) = vec(k,i)*s + h*c
                        vec(k,i) = vec(k,i)*c - h*s
 220                 continue
 230              continue
 240           continue
               e(l) = s*p
c     convergence
               d(l) = c*p
               if (dabs(e(l)).gt.b) go to 200
c     ordering of eigenvalues
               d(l) = d(l) + f
            end if
 250     continue
         if (iff.lt.1) then
            ni = n - 1
            do 280 i = 1 , ni
               k = i
               p = d(i)
               j1 = i + 1
               do 260 j = j1 , n
                  if (d(j).lt.p) then
                     k = j
                     p = d(j)
                  end if
 260           continue
               if (k.ne.i) then
                  d(k) = d(i)
                  d(i) = p
                  do 270 j = 1 , n
                     p = vec(j,i)
                     vec(j,i) = vec(j,k)
                     vec(j,k) = p
 270              continue
               end if
c     special treatment of case n = 1
 280        continue
         end if
      end if
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
      subroutine mat_adj(matin,ndim,matout)
c adjoint matrix of matin
      real*8 matin,matout,matscr,det
      integer ndim,i,j,k,l,ii,jj
      dimension matin(ndim,ndim),matout(ndim,ndim),
     1          matscr(ndim-1,ndim-1)
      call vec_dinit(matscr,(ndim-1)*(ndim-1),0.0d0)
      do 10, i=1 , ndim
        do 20, j=1 , ndim
             do 30, ii=1 , ndim-1
                do 40, jj=1 , ndim-1
                        if (ii.lt.i) then
                         k=ii
                        else
                         k=ii+1
                        endif
                        if (jj.lt.j) then
                         l=jj
                        else
                         l=jj+1
                        endif
                        matscr(ii,jj)=matin(k,l)
                       if(dabs(matscr(ii,jj)).lt.1.0d-63) then
                         matscr(ii,jj)=1.0d-63
                       endif
 40               continue
 30             continue
                matout(i,j)=det(matscr,ndim-1)
                 if(dabs(matout(i,j)).lt.1.0d-70) goto 20
                matout(i,j)= (-1.0d0)**(i+j)*matout(i,j)
 20     continue
 10   continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccc
