← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to contrib/surfphot/libsrc/simul.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/surfphot/libsrc/simul.for
 
1
C @(#)simul.for 19.1 (ES0-DMD) 02/25/03 13:30:48
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Corresponding concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
      SUBROUTINE SIMUL(N,A,X)
 
30
C++++++++
 
31
C   THE GAUSS-JORDAN COMPLETE ELIMINATION METHOD IS EMPLOYED
 
32
C   WITH THE MAXIMUM PIVOT STRATEGY.
 
33
C--------
 
34
      IMPLICIT    NONE
 
35
      INTEGER     N
 
36
      DOUBLE PRECISION    A(4,4)
 
37
      REAL        X(1)
 
38
C
 
39
      INTEGER     IROW(4)
 
40
      INTEGER     JCOL(4)
 
41
      INTEGER     I, ISCAN, IROWI, IROWK
 
42
      INTEGER     J, JSCAN, JCOLI, JCOLK
 
43
      INTEGER     MAX
 
44
      INTEGER     K, KM1
 
45
      INTEGER     ISTAT
 
46
      DOUBLE PRECISION    AIJCK, PIVOT
 
47
      DOUBLE PRECISION    EPS
 
48
      PARAMETER (EPS=1.D-15)
 
49
C
 
50
      MAX=N+1   
 
51
C   BEGIN ELIMINATION PROCEDURE   
 
52
      DO 19 K=1,N                 
 
53
      KM1=K-1                     
 
54
C   SEARCH FOR THE PIVOT ELEMENT  
 
55
      PIVOT=0.D0                    
 
56
      DO 10 I=1,N                 
 
57
      DO 11 J=1,N                 
 
58
C   SEARCH IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCRIPTS
 
59
      IF(K.EQ.1)GOTO 9            
 
60
      DO 8 ISCAN=1,KM1            
 
61
      DO 7 JSCAN=1,KM1            
 
62
      IF(I.EQ.IROW(ISCAN))GOTO 11 
 
63
      IF(J.EQ.JCOL(JSCAN))GOTO 11 
 
64
7     CONTINUE                    
 
65
8     CONTINUE                    
 
66
9     CONTINUE                    
 
67
      IF(DABS(A(I,J)).LE.DABS(PIVOT))GOTO 11
 
68
      PIVOT=A(I,J)
 
69
      IROW(K)=I   
 
70
      JCOL(K)=J   
 
71
11    CONTINUE    
 
72
10    CONTINUE    
 
73
C   INSURE THAT SELECTED PIVOT IS LARGER THAN EPS
 
74
      IF(DABS(PIVOT).GT.EPS)GOTO 13
 
75
       CALL STTPUT(' WARNING! SINGULAR MATRIX',' ',' ',ISTAT)
 
76
      RETURN                      
 
77
13    IROWK=IROW(K)               
 
78
      JCOLK=JCOL(K)               
 
79
C   NORMALIZE PIVOT ROW ELEMENTS  
 
80
      DO 14 J=1,MAX               
 
81
14    A(IROWK,J)=A(IROWK,J)/PIVOT 
 
82
C   CARRY OUT ELIMINATION         
 
83
      A(IROWK,JCOLK)=1.D0/PIVOT     
 
84
      DO 18 I=1,N                 
 
85
      AIJCK=A(I,JCOLK)            
 
86
      IF(I.EQ.IROWK)GOTO 18       
 
87
      A(I,JCOLK)=-AIJCK/PIVOT     
 
88
      DO 17 J=1,MAX               
 
89
      IF(J.NE.JCOLK) A(I,J)=A(I,J)-AIJCK*A(IROWK,J) 
 
90
17    CONTINUE              
 
91
18    CONTINUE
 
92
19    CONTINUE   
 
93
C   ORDER SOLUTION VALUES   
 
94
      DO 20 I=1,N           
 
95
      IROWI=IROW(I)         
 
96
      JCOLI=JCOL(I)         
 
97
      X(JCOLI)= A(IROWI,MAX)
 
98
20    CONTINUE
 
99
      RETURN
 
100
      END
 
101