← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/invent/libsrc/lsqsol.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/lsqsol.for
 
1
C @(#)lsqsol.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 LSQSOL         version 1.0       870130
 
32
C  A. Kruszewski             Obs. de Geneve
 
33
C.PURPOSE
 
34
C  Solves a set of normal equations given by array A.
 
35
C  Results are stored in array X. The last element of array X
 
36
C  contains a sum of squared deviations. Mean errors are stored in
 
37
C  array S. The last element of array S contains the mean error
 
38
C  of an observation of unit weight. The dimension N1 is equal to
 
39
C  number of unknowns plus one.
 
40
C.INPUT/OUTPUT
 
41
C  input arguments
 
42
C  A           real*8 array     symmetrical array
 
43
C  N1          integer*4        dimension of array
 
44
C  M           integer*4        number of conditional equations
 
45
C  output arguments
 
46
C  X           real*8 array     unknowns + vv
 
47
C  S           real*8 array     mean errors + sigma0
 
48
C-----------------------------------------------------------------------
 
49
      SUBROUTINE LSQSOL ( N1 , A , M , X , S )
 
50
C
 
51
      IMPLICIT          NONE
 
52
      INTEGER           N1
 
53
      DOUBLE PRECISION  A(N1,N1)
 
54
      INTEGER           M
 
55
      DOUBLE PRECISION  X(N1)
 
56
      DOUBLE PRECISION  S(N1)
 
57
C
 
58
      INTEGER            N, NN
 
59
      INTEGER            I, I1
 
60
      INTEGER            J, J1
 
61
      INTEGER            K
 
62
      DOUBLE PRECISION   TEMP
 
63
C
 
64
C *** N - number of unknowns.
 
65
C
 
66
      N = N1 - 1
 
67
C
 
68
C *** Calculate triangular square root of the input array.
 
69
C
 
70
      IF ( A(1,1) .LE. 0.0 ) THEN
 
71
          S(N1) = -1.0
 
72
          RETURN
 
73
      ENDIF
 
74
      A(1,1) = SQRT( A(1,1) )
 
75
      DO 10 I = 2 , N1
 
76
          I1 = I - 1
 
77
          DO 20 J = 1 , I1
 
78
              IF ( A(J,J) .LE. 0.0 ) THEN
 
79
                  S(N1) = -1.0
 
80
                  RETURN
 
81
              ENDIF 
 
82
              TEMP = A(I,J)
 
83
              J1 = J - 1
 
84
              DO 30 K = 1 , J1
 
85
                  TEMP = TEMP - A(I,K) * A(J,K)
 
86
   30         CONTINUE
 
87
              A(I,J) = TEMP / A(J,J)
 
88
   20     CONTINUE
 
89
          TEMP = A(I,I)
 
90
          DO 40 K = 1 , I1
 
91
              TEMP = TEMP - A(I,K) * A(I,K)
 
92
   40     CONTINUE
 
93
          IF ( TEMP .LE. 0.0 ) THEN
 
94
              S(N1) = -1.0
 
95
              RETURN
 
96
          ENDIF
 
97
          A(I,I) = SQRT( TEMP )
 
98
   10 CONTINUE
 
99
C
 
100
C *** Invert the triangular square root matrix of size N. 
 
101
C
 
102
      DO 50 I = 1 , N
 
103
          A(I,I) = 1.0 / A(I,I)
 
104
   50 CONTINUE
 
105
      NN = N - 1
 
106
      DO 60 I = 1 , NN
 
107
          I1 = I + 1
 
108
          DO 70 J = I1 , N
 
109
              J1 = J - 1
 
110
              TEMP = 0.0
 
111
              DO 80 K = I , J1
 
112
                  TEMP = TEMP + A(I,K) * A(J,K)   
 
113
   80         CONTINUE
 
114
              A(I,J) = -TEMP * A(J,J)
 
115
   70     CONTINUE
 
116
   60 CONTINUE
 
117
C
 
118
C *** Calculate unknowns. 
 
119
C
 
120
      DO 90 I = 1 , N
 
121
          TEMP = 0.0
 
122
          DO 100 J = I , N
 
123
              TEMP=TEMP + A(I,J) * A(N1,J)
 
124
  100     CONTINUE
 
125
          X(I) = SNGL( TEMP )
 
126
   90 CONTINUE
 
127
C
 
128
C *** Calculate vv and sigma0.
 
129
C
 
130
      TEMP = A(N1,N1)
 
131
      X(N1) = SNGL( TEMP * TEMP )
 
132
      S(N1) = SNGL( TEMP / SQRT( DBLE( M-N ) ) )
 
133
C
 
134
C *** Calculate mean errors of unknowns.
 
135
C
 
136
      DO 110 I = 1 , N
 
137
          TEMP = 0.0
 
138
          DO 120 J = I , N
 
139
              TEMP = TEMP + A(I,J) * A(I,J)
 
140
  120     CONTINUE
 
141
          S(I) = SNGL( SQRT(TEMP) * S(N1) )
 
142
  110 CONTINUE
 
143
C
 
144
      RETURN
 
145
C
 
146
      END