← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/invent/libsrc/invert.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/invent/libsrc/invert.for
 
1
C @(#)invert.for        19.1 (ES0-DMD) 02/25/03 13:25:35
 
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
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C.IDENTIFICATION
 
31
C  subroutine INVERT         version 1.0       830818
 
32
C  A. Kruszewski             ESO Garching
 
33
C.PURPOSE
 
34
C  inverts a symmetrical positively defined array SMTR
 
35
C  and keeps the result in the same array SMTR
 
36
C  it was written and tested in the case when array SMTR represents
 
37
C  a covariance matrix of a multimensional random variable
 
38
C  it may fail in other applications
 
39
C.INPUT/OUTPUT
 
40
C  input arguments
 
41
C  SMTR        real*4 array      symmetrical array
 
42
C  N           integer*4         dimension of array SMTR
 
43
C  output arguments
 
44
C  SMTR        real*4 array      inverted input array
 
45
C-----------------------------------------------------------------------
 
46
      SUBROUTINE INVERT(SMTR,N)
 
47
C
 
48
      IMPLICIT NONE
 
49
      INTEGER  N
 
50
      REAL     SMTR(N,N)
 
51
C
 
52
      INTEGER  I, I1
 
53
      INTEGER  J, J1
 
54
      INTEGER  K
 
55
      INTEGER  N1
 
56
      REAL     TEMP
 
57
C
 
58
C calculates triangular square root of the input array
 
59
C
 
60
      SMTR(1,1) = SQRT(SMTR(1,1))
 
61
      DO 40 I = 2,N
 
62
          I1     = I - 1
 
63
          DO 20 J = 1,I1
 
64
              TEMP   = SMTR(I,J)
 
65
              J1     = J - 1
 
66
              DO 10 K = 1,J1
 
67
                  TEMP   = TEMP - SMTR(I,K)*SMTR(J,K)
 
68
   10         CONTINUE
 
69
              SMTR(I,J) = TEMP/SMTR(J,J)
 
70
   20     CONTINUE
 
71
          TEMP   = SMTR(I,I)
 
72
          DO 30 K = 1,I1
 
73
              TEMP   = TEMP - SMTR(I,K)*SMTR(I,K)
 
74
   30     CONTINUE
 
75
          SMTR(I,I) = SQRT(TEMP)
 
76
   40 CONTINUE
 
77
C
 
78
C inverts the triangular square root matrix
 
79
C
 
80
      DO 50 I = 1,N
 
81
          SMTR(I,I) = 1.0/SMTR(I,I)
 
82
   50 CONTINUE
 
83
      N1     = N - 1
 
84
      DO 80 I = 1,N1
 
85
          I1     = I + 1
 
86
          DO 70 J = I1,N
 
87
              J1     = J - 1
 
88
              TEMP   = 0.0
 
89
              DO 60 K = I,J1
 
90
                  TEMP   = TEMP + SMTR(I,K)*SMTR(J,K)
 
91
   60         CONTINUE
 
92
              SMTR(I,J) = -TEMP*SMTR(J,J)
 
93
   70     CONTINUE
 
94
   80 CONTINUE
 
95
C
 
96
C squares inverted triangular matrix
 
97
C
 
98
      DO 110 I = 1,N
 
99
          DO 100 J = 1,I
 
100
              TEMP   = 0.0
 
101
              DO 90 K = I,N
 
102
                  TEMP   = TEMP + SMTR(I,K)*SMTR(J,K)
 
103
   90         CONTINUE
 
104
              SMTR(J,I) = TEMP
 
105
              SMTR(I,J) = TEMP
 
106
  100     CONTINUE
 
107
  110 CONTINUE
 
108
      RETURN
 
109
 
 
110
      END