C @(#)simul.for	19.1 (ES0-DMD) 02/25/03 13:30:48
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 SIMUL(N,A,X)
C++++++++
C   THE GAUSS-JORDAN COMPLETE ELIMINATION METHOD IS EMPLOYED
C   WITH THE MAXIMUM PIVOT STRATEGY.
C--------
      IMPLICIT    NONE
      INTEGER     N
      DOUBLE PRECISION    A(4,4)
      REAL        X(1)
C
      INTEGER     IROW(4)
      INTEGER     JCOL(4)
      INTEGER     I, ISCAN, IROWI, IROWK
      INTEGER     J, JSCAN, JCOLI, JCOLK
      INTEGER     MAX
      INTEGER     K, KM1
      INTEGER     ISTAT
      DOUBLE PRECISION    AIJCK, PIVOT
      DOUBLE PRECISION    EPS
      PARAMETER (EPS=1.D-15)
C
      MAX=N+1   
C   BEGIN ELIMINATION PROCEDURE   
      DO 19 K=1,N                 
      KM1=K-1                     
C   SEARCH FOR THE PIVOT ELEMENT  
      PIVOT=0.D0                    
      DO 10 I=1,N                 
      DO 11 J=1,N                 
C   SEARCH IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCRIPTS
      IF(K.EQ.1)GOTO 9            
      DO 8 ISCAN=1,KM1            
      DO 7 JSCAN=1,KM1            
      IF(I.EQ.IROW(ISCAN))GOTO 11 
      IF(J.EQ.JCOL(JSCAN))GOTO 11 
7     CONTINUE                    
8     CONTINUE                    
9     CONTINUE                    
      IF(DABS(A(I,J)).LE.DABS(PIVOT))GOTO 11
      PIVOT=A(I,J)
      IROW(K)=I   
      JCOL(K)=J   
11    CONTINUE    
10    CONTINUE    
C   INSURE THAT SELECTED PIVOT IS LARGER THAN EPS
      IF(DABS(PIVOT).GT.EPS)GOTO 13
       CALL STTPUT(' WARNING! SINGULAR MATRIX',' ',' ',ISTAT)
      RETURN                      
13    IROWK=IROW(K)               
      JCOLK=JCOL(K)               
C   NORMALIZE PIVOT ROW ELEMENTS  
      DO 14 J=1,MAX               
14    A(IROWK,J)=A(IROWK,J)/PIVOT 
C   CARRY OUT ELIMINATION         
      A(IROWK,JCOLK)=1.D0/PIVOT     
      DO 18 I=1,N                 
      AIJCK=A(I,JCOLK)            
      IF(I.EQ.IROWK)GOTO 18       
      A(I,JCOLK)=-AIJCK/PIVOT     
      DO 17 J=1,MAX               
      IF(J.NE.JCOLK) A(I,J)=A(I,J)-AIJCK*A(IROWK,J) 
17    CONTINUE              
18    CONTINUE
19    CONTINUE   
C   ORDER SOLUTION VALUES   
      DO 20 I=1,N           
      IROWI=IROW(I)         
      JCOLI=JCOL(I)         
      X(JCOLI)= A(IROWI,MAX)
20    CONTINUE
      RETURN
      END

