C ++********************************************************************
C                                                                      *
C                                                                      *
C                                                                      *
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  Health Research Inc.                      *
C=*                                                                    *
C=* HEALTH RESEARCH INCORPORATED (HRI),                                *   
C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455.                   *
C=*                                                                    *
C=* Email:  spider@wadsworth.org                                       *
C=*                                                                    *
C=* This program is free software; you can redistribute it and/or      *
C=* modify it under the terms of the GNU General Public License as     *
C=* published by the Free Software Foundation; either version 2 of the *
C=* License, or (at your option) any later version.                    *
C=*                                                                    *
C=* This program is distributed in the hope that it will be useful,    *
C=* but WITHOUT ANY WARRANTY; without even the implied warranty of     *
C=* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C **********************************************************************
C                                                                      *
C VPROP                                                                     *
C                                                                      *
C  PURPOSE:                                                            *
C                                                                      *
C  PARAMETERS:                                                         *
C                                                                      *
C
C Double precision version 10/21/96 PP
C======================================================================

      SUBROUTINE VPROP( NDIM, N, W, D, S, KODE ) 
                            
      DIMENSION  W(NDIM,N) , D(N) , S(N) 
                                       
      DATA  SEUIL / 1.0E-7 /
	DOUBLE PRECISION Q,H,T,U,V,P,B
                                                  
      KODE  = 0                                                                 
      CALL  TRIDI (NDIM,N,W,D,S)                                                
        IF (N .EQ. 1)        GO TO  140                                         
        DO   I = 2,N                                                          
      S(I-1)= S(I)                                                              
        ENDDO
      S(N)  = 0.0                                                               
        DO 90  K = 1,N                                                          
      M     = 0                                                                 
   20   DO   J = K,N                                                          
        IF (J .EQ. N)        GO TO   40                                         
      ABJ   = ABS(S(J))                                                         
      EPS   = SEUIL*(ABS(D(J)) + ABS(D(J+1)))                                   
        IF (ABJ .LE. EPS)    GO TO   40                                         
	ENDDO
   40 H     = D(K)                                                              
        IF (J .EQ. K)        GO TO   90                                         
        IF (M .EQ. 30)       GO TO  130                                         
      M     = M + 1                                                             
      Q     = (D(K+1) - H) / (2.0*S(K))                                         
      T     = DSQRT (Q*Q + 1.0)                                                  
      Q     = D(J) - H + S(K) / (Q + DSIGN(T,Q))                                 
      U     = 1.0                                                               
      V     = 1.0                                                               
      H     = 0.0                                                               
      JK    = J - K                                                             
        DO  IJK = 1,JK                                                        
      I     = J - IJK                                                           
      P     = U * S(I)                                                          
      B     = V * S(I)                                                          
        IF (ABS(P).LT.ABS(Q))GO TO   50                                         
      V     = Q / P                                                             
      T     = DSQRT (V*V + 1.0)                                                  
      S(I+1)= P * T                                                             
      U     = 1.0 / T                                                           
      V     = V * U                                                             
                             GO TO   60                                         
   50 U     = P / Q                                                             
      T     = DSQRT (U*U + 1.0)                                                  
      S(I+1)= Q * T                                                             
      V     = 1.0 / T                                                           
      U     = U * V                                                             
   60 Q     = D(I+1) - H                                                        
      T     = (D(I) - Q)*U + 2.0*V*B                                            
      H     = U * T                                                             
      D(I+1)= Q + H                                                             
      Q     = V*T - B                                                           
        DO   L = 1,N                                                          
      P     = W(L,I+1)                                                          
      W(L,I+1)= U*W(L,I) + V*P                                                  
      W(L,I)= V*W(L,I) - U*P                                                    
	ENDDO
	ENDDO
      D(K)  = D(K) - H                                                          
      S(K)  = Q                                                                 
      S(J)  = 0.0                                                               
                             GO TO   20                                         
   90   CONTINUE                                                                
        DO 120  IJ = 2,N                                                        
      I     = IJ - 1                                                            
      L     = I                                                                 
      H     = D(I)                                                              
        DO 100  M = IJ,N                                                        
        IF (D(M) .LE. H)     GO TO  100                                         
      L     = M                                                                 
      H     = D(M)                                                              
  100   CONTINUE                                                                
        IF (L .EQ. I)        GO TO  120                                         
      D(L)  = D(I)                                                              
      D(I)  = H                                                                 
        DO   M = 1,N                                                         
      H     = W(M,I)                                                            
      W(M,I)= W(M,L)                                                            
      W(M,L)= H                                                                 
	ENDDO
  120   CONTINUE                                                                
                             GO TO  140                                         
  130 KODE  = K                                                                 
  140     RETURN                                                                
          END