C @(#)simul.for 17.1.1.1 (ES0-DMD) 01/25/02 17:19:02 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