C @(#)smain.for	17.1.1.1 (ES0-DMD) 01/25/02 17:17:57
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
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 
C the 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
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      SUBROUTINE SMAIN (ARRAY,NP,NORDER,DET)
C+++
C.PURPOSE:
C.INPUT/OUTPUT: ARRAY  : input matrix which is replaced by its inverse
C.              NP     : matrix array dimension
C.              NORDER : degree of matrix (order of determinant)
C.              DET    : determinant of input matrix
C.              DIMENSION STATEMENT VALID FOR NORDER UP TO 500
C.VERSION:      870617
C---
      IMPLICIT  NONE
      INTEGER   NP
      REAL      ARRAY(NP,NP)
      INTEGER   NORDER
      REAL      DET
C
      INTEGER   IK(500),JK(500)
      INTEGER   I, J, K, L
      REAL      AMAX
      REAL      SAVE
C
C ***
      DET = 1.
      K   = 0
   20 CONTINUE
      IF (K.GE.NORDER.OR.DET.EQ.0.) GO TO 10
      K   = K+1
C
C *** find larger element array(i,j) in rest of matrix
      J    = 0
      I    = 0
      AMAX = 0.
C
   30 CONTINUE
      IF (J-K.GE.0) GO TO 40
   50 CONTINUE
      IF (I-K.GE.0) GO TO 60
      DO 901 I=K,NORDER
         DO 911 J=K,NORDER
            IF (ABS(ARRAY(I,J)).GE.ABS(AMAX)) THEN
               AMAX  = ARRAY(I,J)
               IK(K) = I
               JK(K) = J
             END IF
  911    CONTINUE
  901 CONTINUE
C
C *** interchange rows and columns to put amax in array(k,k)
      IF (AMAX.EQ.0.) THEN
         DET=0.
         RETURN
      ELSE
         I = IK(K)
      END IF
      IF (I-K.GT.0) THEN
         DO 902 J=1,NORDER
            SAVE       = ARRAY(K,J)
            ARRAY(K,J) = ARRAY(I,J)
            ARRAY(I,J) = -SAVE
  902    CONTINUE
      END IF
      GO TO 50
C
   60 CONTINUE
      J = JK(K) 
      IF (J-K.GT.0) THEN
         DO 903 I=1,NORDER
            SAVE       = ARRAY(I,K)
            ARRAY(I,K) = ARRAY(I,J) 
            ARRAY(I,J) = -SAVE
  903    CONTINUE
C
C *** accumulate elements of inverse matrix
      END IF
      GO TO 30
C
   40 CONTINUE
      DO 904 I = 1,NORDER
         IF (I-K.NE.0) THEN
            ARRAY(I,K) = -ARRAY(I,K)/AMAX
         END IF
  904 CONTINUE
C
      DO 905 I=1,NORDER
         DO 915 J=1,NORDER
            IF (I-K.NE.0.AND.J-K.NE.0) THEN
               ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
            END IF
  915    CONTINUE
  905 CONTINUE
C
      DO 906 J = 1,NORDER
         IF (J-K.NE.0) THEN
            ARRAY(K,J) = ARRAY(K,J)/AMAX
         END IF
  906 CONTINUE
      ARRAY(K,K) = 1./AMAX
C     DET        = DET*AMAX
      GO TO 20
C
C *** restore ordering of matrix
   10 CONTINUE
      DO 907 L = 1,NORDER
         K = NORDER-L+1
         J = IK(K)
         IF (J-K.GT.0) THEN
            DO 917 I=1,NORDER
               SAVE       = ARRAY(I,K) 
               ARRAY(I,K) = -ARRAY(I,J)
               ARRAY(I,J) = SAVE
  917       CONTINUE
         END IF
C
         I = JK(K)
         IF (I-K.GT.0) THEN
            DO 927 J=1,NORDER
               SAVE       = ARRAY(K,J)
               ARRAY(K,J) = -ARRAY(I,J)
               ARRAY(I,J) = SAVE
  927       CONTINUE
         END IF
  907 CONTINUE
C
      RETURN
      END