← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/ccdred/libsrc/mo_gaussj.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2015-06-10 14:20:37 UTC
  • mfrom: (1.2.1) (6.1.9 experimental)
  • Revision ID: package-import@ubuntu.com-20150610142037-6iowpbtyjrpou36o
Tags: 15.02pl1.3-1
* New upstream version
* Add CI tests
* Move back to unstable

Show diffs side-by-side

added added

removed removed

Lines of Context:
stdred/ccdred/libsrc/mo_gaussj.c
1
 
/*===========================================================================
2
 
  Copyright (C) 1995-2009 European Southern Observatory (ESO)
3
 
 
4
 
  This program is free software; you can redistribute it and/or 
5
 
  modify it under the terms of the GNU General Public License as 
6
 
  published by the Free Software Foundation; either version 2 of 
7
 
  the License, or (at your option) any later version.
8
 
 
9
 
  This program is distributed in the hope that it will be useful,
10
 
  but WITHOUT ANY WARRANTY; without even the implied warranty of
11
 
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12
 
  GNU General Public License for more details.
13
 
 
14
 
  You should have received a copy of the GNU General Public 
15
 
  License along with this program; if not, write to the Free 
16
 
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
17
 
  MA 02139, USA.
18
 
 
19
 
  Correspondence concerning ESO-MIDAS should be addressed as follows:
20
 
        Internet e-mail: midas@eso.org
21
 
        Postal address: European Southern Observatory
22
 
                        Data Management Division 
23
 
                        Karl-Schwarzschild-Strasse 2
24
 
                        D 85748 Garching bei Muenchen 
25
 
                        GERMANY
26
 
===========================================================================*/
27
 
 
28
 
/*
29
 
 
30
 
.VERSION
31
 
 090526         last modif
32
 
 
33
 
*/
34
 
 
35
 
 
36
 
#include <midas_def.h>
37
 
#include <math.h>
38
 
 
39
 
#define SWAP(a,b) {float temp=(a);(a)=(b);(b)=temp;}
40
 
 
41
 
void MO_GAUSSJ(a,n,b,m)
42
 
float **a,**b;
43
 
int n,m;
44
 
 
45
 
{
46
 
int   *indxc,*indxr,*ipiv;
47
 
int   i,icol,irow,j,k,l,ll;
48
 
float big,dum,pivinv;
49
 
 
50
 
irow = icol = 0;
51
 
 
52
 
        /* later on, always indexing from 1 to n is used... (KB) */
53
 
 
54
 
        indxc = (int *) osmmget((n+1) * sizeof(int));
55
 
        indxr = (int *) osmmget((n+1) * sizeof(int));
56
 
        ipiv  = (int *) osmmget((n+1) * sizeof(int));
57
 
 
58
 
        for (j = 1; j<= n;j++) ipiv[j]=0;
59
 
        for (i = 1; i <= n; i++) {
60
 
           big=0.0;
61
 
           for (j = 1; j <= n; j++)
62
 
              if (ipiv[j] != 1)
63
 
                 for (k = 1; k <= n; k++) {
64
 
                    if (ipiv[k] == 0) {
65
 
                       if (fabs(a[j][k]) >= big) {
66
 
                          big  = fabs(a[j][k]);
67
 
                          irow = j;
68
 
                          icol = k;
69
 
                       }
70
 
                    } 
71
 
              else if (ipiv[k] > 1) 
72
 
                 SCETER(2, "*** FATAL: GAUSSJ: Singular Matrix-1");
73
 
                                }
74
 
           ++(ipiv[icol]);
75
 
           if (irow != icol) {
76
 
              for (l = 1; l <= n; l++) SWAP(a[irow][l],a[icol][l])
77
 
              for (l = 1; l <= m; l++) SWAP(b[irow][l],b[icol][l])
78
 
            }
79
 
 
80
 
           indxr[i]       = irow;
81
 
           indxc[i]       = icol;
82
 
           if (a[icol][icol] == 0.0) 
83
 
              SCETER(2, "*** FATAL: GAUSSJ - Singular Matrix-2");
84
 
           pivinv         = 1.0/a[icol][icol];
85
 
           a[icol][icol]  = 1.0;
86
 
           for (l = 1; l <= n; l++) a[icol][l] *= pivinv;
87
 
           for (l = 1; l <= m; l++) b[icol][l] *= pivinv;
88
 
           for (ll = 1; ll <= n; ll++)
89
 
           if (ll != icol) {
90
 
              dum         = a[ll][icol];
91
 
              a[ll][icol] = 0.0;
92
 
              for (l = 1; l <= n; l++) a[ll][l] -= a[icol][l]*dum;
93
 
              for (l = 1; l <= m; l++) b[ll][l] -= b[icol][l]*dum;
94
 
           }
95
 
        }
96
 
 
97
 
        for (l = n; l >= 1; l--) {
98
 
           if (indxr[l] != indxc[l])
99
 
              for (k = 1; k <= n; k++)
100
 
                 SWAP(a[k][indxr[l]],a[k][indxc[l]]);
101
 
        }
102
 
 
103
 
        osmmfree((char *) ipiv);
104
 
        osmmfree((char *) indxr);
105
 
        osmmfree((char *) indxc);
106
 
}
107
 
 
108
 
#undef SWAP