SUBROUTINE SPIN(Q,N,MM,DET,IDEXP)                                         
C                                                                               
C     SUBROUTINE SPIN IS A MATRIX INVERSION ROUTINE FOR SYMMETRIC               
C     POSITIVE-DEFINITE MATRICES.  THE MATRIX INVERTED IS THE UPPER             
C     N BY N PORTION OF THE MATRIX Q WHICH IS DIMENSIONED MM BY MM              
C     IN THE CALLING ROUTINE.                                                   
C                                                                               
C  INPUT:                                                                       
C       Q - THE MATRIX DIMENSIONED MM BY MM WHICH CONTAINS THE                  
C           MATRIX TO BE INVERTED.                                              
C                                                                               
C       N - THE DIMENSION OF THE ACTUAL PART( UPPER LEFT CORNER)                
C           OF Q WHICH IS TO BE INVERTED. ( N MAY BE EQUAL BUT MUST             
C           NOT BE LARGER THAN MM) .                                            
C       MM- DIMENSIONED SIZE OF Q IN THE CALLING ROUTINE.                       
C                                                                               
C                                                                               
C  OUTPUT:                                                                      
C                                                                               
C       Q - THE UPPER LEFT N BY N PORTION CONTAINS THE INVERSE OF               
C           THE INPUT UPPER LEFT N BY N PORTION.                                
C                                                                               
C       DET - THE NON-EXPONENT PORTION OF THE DETERMINANT OF THE                
C             INPUT N BY N (UPPER LEFT PORTION OF Q) MATRIX. SEE                
C             IDEXP  BELOW.                                                     
C                                                                               
C       IDEXP - THE EXPONENT (OF 10) PART OF THE DETERMINANT DESCRIBED          
C               UNDER DET ABOVE. THUS THE DETERMINANT IS RETURNED IN            
C               TWO PARTS CORRESPONDING TO                                      
C                     DETERMINANT = DET * 10 ** IDEXP .                         
C               THIS IS DONE TO AVOID UNDER OR OVERFLOW IN THE                  
C               COMPUTATION OF THE DETERMINANT.  TO PRINT THE DETERM-           
C               INANT THE USER SHOULD PRINT BOTH NUMBERS AS FOLLOWS;            
C               (FOR EXAMPLE)                                                   
C                  PRINT 10,DET,IDEXP                                           
C               10 FORMAT(' ','DETERMINANT= ',F7.4,'D',I4)                      
C                                                                               
C                                                R.R. STEEVES                   
C                                                SEPT., 1979                    
C                                                                               
C                                                                               
      REAL*8 Q(MM,MM),DET,SUM,DSQRT,DABS,RPART,APART                            
      DET=0.D0                                                                  
      DO 4 J=1,N                                                                
      DO 4 I=1,J                                                                
      IF(I.EQ.1) GO TO 2                                                      20
      M=I-1                                                                     
      SUM=0.0D0                                                                 
      DO 1 K=1,M                                                                
1     SUM=SUM+Q(K,I)*Q(K,J)                                                     
      Q(I,J)=Q(I,J)-SUM                                                         
2     IF(I.EQ.J) GO TO 3                                                        
      Q(I,J)=Q(I,J)/Q(I,I)                                                      
      GO TO 4                                                                   
3     CONTINUE                                                                  
      DET=DET+DLOG10(Q(I,I))                                                    
      Q(I,I)=DSQRT(Q(I,I))                                                      
4     CONTINUE                                                                  
      IDEXP=DET                                                                 
      RPART=DET-IDEXP                                                           
      APART=DABS(RPART)                                                         
      IF(APART.LT.1.D-20) DET=1.D0                                              
      IF(APART.LT.1.D-20) GO TO 10                                              
      DET=10**RPART                                                             
   10 CONTINUE                                                                  
      DO 7 J=1,N                                                                
      DO 7 I=1,J                                                                
      IF(I.LT.J) GO TO 5                                                        
      Q(J,J)=1.0D0/Q(J,J)                                                       
      GO TO 7                                                                   
5     SUM=0.0D0                                                                 
      M=J-1                                                                     
      DO 6 K=I,M                                                                
6     SUM=SUM-Q(I,K)*Q(K,J)                                                     
      Q(I,J)=SUM/Q(J,J)                                                         
7     CONTINUE                                                                  
      DO 9 J=1,N                                                                
      DO 9 I=1,J                                                                
      SUM=0.0D0                                                                 
      DO 8 K=J,N                                                                
8     SUM=SUM+Q(I,K)*Q(J,K)                                                     
      Q(I,J)=SUM                                                                
      IF(I.EQ.J) GO TO 9                                                        
      Q(J,I)=SUM                                                                
9     CONTINUE                                                                  
      RETURN                                                                    
      END