SUBROUTINE SMJROT(V,A,N,NDIM1,DIAG,NROT)
*     Eigenvalues and eigenvectors of symmetric N*N matrix V(*), which
*     is destroyed during computation.
*     Returned: the N*N matrix of eigenvectors in array A(NDIM1,N)
*     and the diagonal elements (eigenvalues) in array DIAG(N).  
*     NROT is the performed number of rotations 
*     Limitation to N <= NDD! 
      INTEGER NDD
      PARAMETER (NDD=100)
      REAL V(*),A(NDIM1,N),DIAG(N),DDAG(NDD) 
      REAL AIP,AIQ,VIP,VIQ,AM,SM,DIFF,EPSFAC,TEST
      REAL COSPHI,SINPHI,TANPHI,TANPH2,THETA
      INTEGER I,J,N,II,IJ,P,Q,PQ,NDIM1,NROT,MSY,MDI,ITER
      PARAMETER (EPSFAC=100.0)
*     indices in symmetric matrix
      MSY(I,J)=MIN(I,J)+(MAX(I,J)*MAX(I,J)-MAX(I,J))/2
      MDI(I)  =(I*I+I)/2
      NROT=0 
      DO J=1,N   
       DO I=1,N
        A(I,J)=0.0
       END DO
       A(J,J)=1.0 ! A = unit matrix
      END DO
      II=0
      DO I=1,N
       II=II+I
       DIAG(I)=V(II) ! copy diagonal elements to DIAG
       IF(I.LE.NDD) DDAG(I)=0.0
      END DO 
      DO ITER=-N,50   ! iterations
       SM=0.0
       AM=0.0
       IJ=0
       DO I=1,N
        DO J=1,I-1
         IJ=IJ+1
         SM=SM+ABS(V(IJ))                   ! sum 
         IF(ABS(V(IJ)).GT.AM) AM=ABS(V(IJ)) ! maximum
        END DO 
        IJ=IJ+1 
       END DO  
       IF(SM.EQ.0.0) RETURN  ! test convergence
       PQ=0
       DO Q=1,N
        DO P=1,Q-1       
         PQ=PQ+1
         IF(ITER.GE.0.OR.ABS(V(PQ)).EQ.AM) THEN
            TEST=EPSFAC*ABS(V(PQ))
            IF(ITER.LT.0.OR.
     +         ABS(DIAG(P)).NE.ABS(DIAG(P))+TEST.OR.
     +         ABS(DIAG(Q)).NE.ABS(DIAG(Q))+TEST) THEN
               NROT=NROT+1       ! begin/Jacobi rotation  
               DIFF=DIAG(Q)-DIAG(P) 
               IF(ABS(DIFF).EQ.ABS(DIFF)+TEST) THEN
                  TANPHI=V(PQ)/DIFF
               ELSE 
                  THETA =0.5*DIFF/V(PQ)
                  TANPHI=1.0/(THETA+SIGN(SQRT(1.0+THETA*THETA),THETA))
               END IF
               COSPHI=1.0/SQRT(1.0+TANPHI*TANPHI)
               SINPHI=COSPHI*TANPHI
               TANPH2=SINPHI/(1.0+COSPHI)
               DO I=1,N
                IF(I.NE.Q.AND.I.NE.P) THEN
                   VIP=V(MSY(I,P)) ! transform symmetric
                   VIQ=V(MSY(I,Q)) !  matrix
                   V(MSY(I,P))=VIP-SINPHI*(VIQ+TANPH2*VIP)
                   V(MSY(I,Q))=VIQ+SINPHI*(VIP-TANPH2*VIQ)
                END IF
                AIP=A(I,P) ! update rotation
                AIQ=A(I,Q) !  matrix
                A(I,P)=AIP-SINPHI*(AIQ+TANPH2*AIP)
                A(I,Q)=AIQ+SINPHI*(AIP-TANPH2*AIQ)
               END DO 
               DIAG(P)=DIAG(P)-TANPHI*V(PQ) ! update special elements
               DIAG(Q)=DIAG(Q)+TANPHI*V(PQ) !  of symmetric matrix   
               IF(P.LE.NDD) DDAG(P)=DDAG(P)-TANPHI*V(PQ) 
               IF(Q.LE.NDD) DDAG(Q)=DDAG(Q)+TANPHI*V(PQ) 
            END IF
            V(PQ)=0.0  ! off-diagonal element becomes zero  
         END IF        ! end/Jacobi rotation  
        END DO 
        PQ=PQ+1 
       END DO
       II=0       
       DO I=1,MIN(N,NDD)
        II=II+I
        V(II)=V(II)+DDAG(I) ! improve diagonal elements  
        DIAG(I)=V(II)
        DDAG(I)=0.0
       END DO
      END DO  
      WRITE(*,*) 'SMJROT: no convergence - stop'
      STOP 
      END