1 /*
2 *+
3 * Name:
4 * palDmat
5
6 * Purpose:
7 * Matrix inversion & solution of simultaneous equations
8
9 * Language:
10 * Starlink ANSI C
11
12 * Type of Module:
13 * Library routine
14
15 * Invocation:
16 * void palDmat( int n, double *a, double *y, double *d, int *jf,
17 * int *iw );
18
19 * Arguments:
20 * n = int (Given)
21 * Number of simultaneous equations and number of unknowns.
22 * a = double[] (Given & Returned)
23 * A non-singular NxN matrix (implemented as a contiguous block
24 * of memory). After calling this routine "a" contains the
25 * inverse of the matrix.
26 * y = double[] (Given & Returned)
27 * On input the vector of N knowns. On exit this vector contains the
28 * N solutions.
29 * d = double * (Returned)
30 * The determinant.
31 * jf = int * (Returned)
32 * The singularity flag. If the matrix is non-singular, jf=0
33 * is returned. If the matrix is singular, jf=-1 & d=0.0 are
34 * returned. In the latter case, the contents of array "a" on
35 * return are undefined.
36 * iw = int[] (Given)
37 * Integer workspace of size N.
38
39 * Description:
40 * Matrix inversion & solution of simultaneous equations
41 * For the set of n simultaneous equations in n unknowns:
42 * A.Y = X
43 * this routine calculates the inverse of A, the determinant
44 * of matrix A and the vector of N unknowns.
45
46 * Authors:
47 * PTW: Pat Wallace (STFC)
48 * TIMJ: Tim Jenness (JAC, Hawaii)
49 * {enter_new_authors_here}
50
51 * History:
52 * 2012-02-11 (TIMJ):
53 * Combination of a port of the Fortran and a comparison
54 * with the obfuscated GPL C routine.
55 * Adapted with permission from the Fortran SLALIB library.
56 * {enter_further_changes_here}
57
58 * Notes:
59 * - Implemented using Gaussian elimination with partial pivoting.
60 * - Optimized for speed rather than accuracy with errors 1 to 4
61 * times those of routines optimized for accuracy.
62
63 * Copyright:
64 * Copyright (C) 2001 Rutherford Appleton Laboratory.
65 * Copyright (C) 2012 Science and Technology Facilities Council.
66 * All Rights Reserved.
67
68 * Licence:
69 * This program is free software: you can redistribute it and/or
70 * modify it under the terms of the GNU Lesser General Public
71 * License as published by the Free Software Foundation, either
72 * version 3 of the License, or (at your option) any later
73 * version.
74 *
75 * This program is distributed in the hope that it will be useful,
76 * but WITHOUT ANY WARRANTY; without even the implied warranty of
77 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
78 * GNU Lesser General Public License for more details.
79 *
80 * You should have received a copy of the GNU Lesser General
81 * License along with this program. If not, see
82 * <http://www.gnu.org/licenses/>.
83
84 * Bugs:
85 * {note_any_bugs_here}
86 *-
87 */
88
89 #include "pal.h"
90
palDmat(int n,double * a,double * y,double * d,int * jf,int * iw)91 void palDmat ( int n, double *a, double *y, double *d, int *jf, int *iw ) {
92
93 const double SFA = 1e-20;
94
95 int k;
96 double*aoff;
97
98 *jf=0;
99 *d=1.0;
100 for(k=0,aoff=a; k<n; k++, aoff+=n){
101 int imx;
102 double * aoff2 = aoff;
103 double amx=fabs(aoff[k]);
104 imx=k;
105 if(k!=n){
106 int i;
107 double *apos2;
108 for(i=k+1,apos2=aoff+n;i<n;i++,apos2+=n){
109 double t=fabs(apos2[k]);
110 if(t>amx){
111 amx=t;
112 imx=i;
113 aoff2=apos2;
114 }
115 }
116 }
117 if(amx<SFA){
118 *jf=-1;
119 } else {
120 if(imx!=k){
121 double t;
122 int j;
123 for(j=0;j<n;j++){
124 t=aoff[j];
125 aoff[j]=aoff2[j];
126 aoff2[j]=t;
127 }
128 t=y[k];
129 y[k]=y[imx];
130 y[imx]=t;*d=-*d;
131 }
132 iw[k]=imx;
133 *d*=aoff[k];
134 if(fabs(*d)<SFA){
135 *jf=-1;
136 } else {
137 double yk;
138 double * apos2;
139 int i, j;
140 aoff[k]=1.0/aoff[k];
141 for(j=0;j<n;j++){
142 if(j!=k){
143 aoff[j]*=aoff[k];
144 }
145 }
146 yk=y[k]*aoff[k];
147 y[k]=yk;
148 for(i=0,apos2=a;i<n;i++,apos2+=n){
149 if(i!=k){
150 for(j=0;j<n;j++){
151 if(j!=k){
152 apos2[j]-=apos2[k]*aoff[j];
153 }
154 }
155 y[i]-=apos2[k]*yk;
156 }
157 }
158 for(i=0,apos2=a;i<n;i++,apos2+=n){
159 if(i!=k){
160 apos2[k]*=-aoff[k];
161 }
162 }
163 }
164 }
165 }
166 if(*jf!=0){
167 *d=0.0;
168 } else {
169 for(k=n;k-->0;){
170 int ki=iw[k];
171 if(k!=ki){
172 int i;
173 double *apos = a;
174 for(i=0;i<n;i++,apos+=n){
175 double t=apos[k];
176 apos[k]=apos[ki];
177 apos[ki]=t;
178 }
179 }
180 }
181 }
182 }
183