1 /*****************************************************************************
2 *
3 * Elmer, A Finite Element Software for Multiphysical Problems
4 *
5 * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library (in file ../LGPL-2.1); if not, write
19 * to the Free Software Foundation, Inc., 51 Franklin Street,
20 * Fifth Floor, Boston, MA 02110-1301 USA
21 *
22 *****************************************************************************/
23
24 /*******************************************************************************
25 *
26 * LU decomposition.
27 *
28 *******************************************************************************
29 *
30 * Author: Juha Ruokolainen
31 *
32 * Address: CSC - IT Center for Science Ltd.
33 * Keilaranta 14, P.O. BOX 405
34 * 02101 Espoo, Finland
35 * Tel. +358 0 457 2723
36 * Telefax: +358 0 457 2302
37 * EMail: Juha.Ruokolainen@csc.fi
38 *
39 * Date: 30 May 1996
40 *
41 * Modified by:
42 *
43 * Date of modification:
44 *
45 ******************************************************************************/
46 /***********************************************************************
47 |
48 | LU.C - Last Edited 8. 8. 1988
49 |
50 ***********************************************************************/
51
52 /*======================================================================
53 |Syntax of the manual pages:
54 |
55 |FUNCTION NAME(...) params ...
56 |
57 $ usage of the function and type of the parameters
58 ? explane the effects of the function
59 = return value and the type of value if not of type int
60 @ globals effected directly by this routine
61 ! current known bugs or limitations
62 & functions called by this function
63 ~ these functions may interest you as an alternative function or
64 | because they control this function somehow
65 ^=====================================================================*/
66
67
68 /*
69 * $Id: lu.c,v 1.1.1.1 2005/04/14 13:29:14 vierinen Exp $
70 *
71 * $Log: lu.c,v $
72 * Revision 1.1.1.1 2005/04/14 13:29:14 vierinen
73 * initial matc automake package
74 *
75 * Revision 1.2 1998/08/01 12:34:45 jpr
76 *
77 * Added Id, started Log.
78 *
79 *
80 */
81
82 #include "elmer/matc.h"
83
84 #define A(i,j) a[n * (i) + (j)]
85
mtr_LUD(var)86 VARIABLE *mtr_LUD(var)
87 VARIABLE *var;
88 {
89 VARIABLE *res;
90 int i, n, *pivot;
91 double *a;
92
93 if (NCOL(var) != NROW(var))
94 {
95 error("LUD: Matrix must be square.\n");
96 }
97
98 res = var_temp_copy(var);
99 a = MATR(res); n = NROW(res);
100
101 pivot = (int *)ALLOCMEM(n * sizeof(int));
102
103 LUDecomp(a, n, pivot);
104
105 FREEMEM((char *)pivot);
106
107 return res;
108 }
109
mtr_det(var)110 VARIABLE *mtr_det(var)
111 VARIABLE *var;
112 {
113 VARIABLE *tmp, *res;
114 int i, n, *pivot;
115 double det, *a;
116
117 if (NCOL(var) != NROW(var))
118 {
119 error("Det: Matrix must be square.\n");
120 }
121
122 tmp = var_temp_copy(var);
123
124 a = MATR(tmp); n = NROW(tmp);
125
126 pivot = (int *)ALLOCMEM(n * sizeof(int));
127
128 LUDecomp(a, n, pivot);
129
130 det = 1.0;
131 for(i = 0; i < n; i++)
132 {
133 det *= A(i,i);
134 if (pivot[i] != i) det = -det;
135 }
136
137 FREEMEM((char *)pivot); var_delete_temp(tmp);
138
139 res = var_temp_new(TYPE_DOUBLE,1,1);
140
141 M(res,0,0) = det;
142
143 return res;
144 }
145
mtr_inv(var)146 VARIABLE *mtr_inv(var)
147 VARIABLE *var;
148 {
149 VARIABLE *ptr;
150
151 int i, j , k, n, *pivot;
152 double s, *a;
153
154 if (NCOL(var) != NROW(var))
155 {
156 error("Inv: Matrix must be square.\n");
157 }
158
159 ptr = var_temp_copy(var);
160
161 a = MATR(ptr); n = NROW(ptr);
162
163 pivot = (int *)ALLOCMEM(n * sizeof(int));
164
165 /*
166 * AP = LU
167 */
168 LUDecomp(a, n, pivot);
169 for(i = 0; i < n; i++)
170 {
171 if (A(i,i) == 0)
172 error("Inv: Matrix is singular.\n");
173 A(i,i) = 1.0 / A(i,i);
174 }
175
176 /*
177 * INV(U)
178 */
179 for(i = n - 2; i >= 0; i--)
180 for(j = n - 1; j > i; j--)
181 {
182 s = 0.0;
183 for(k = i+1; k <= j; k++)
184 if (k != j)
185 s = s - A(i,k) * A(k,j);
186 else
187 s = s - A(i,k);
188 A(i,j) = s;
189 }
190
191 /*
192 * INV(L)
193 */
194 for(i = n - 2; i >= 0; i--)
195 for(j = n - 1; j > i; j--)
196 {
197 s = 0.0;
198 for(k = i + 1; k <= j; k++)
199 s = s - A(j,k) * A(k,i);
200 A(j,i) = A(i,i) * s;
201 }
202
203 /*
204 * A = INV(AP)
205 */
206 for(i = 0; i < n; i++)
207 for(j = 0; j < n; j++)
208 {
209 s = 0.0;
210 for(k = max(i,j); k < n; k++)
211 if (k != i)
212 s += A(i,k) * A(k,j);
213 else
214 s += A(k,j);
215 A(i,j) = s;
216 }
217
218 /*
219 * A = INV(A) (at last)
220 */
221 for(i = 0; i < n; i++)
222 if (pivot[i] != i)
223 for(j = 0; j < n; j++)
224 {
225 s = A(i,j);
226 A(i,j) = A(pivot[i],j);
227 A(pivot[i],j) = s;
228 }
229
230 FREEMEM((char *)pivot);
231
232 return ptr;
233 }
234
235 /*
236 * LU- decomposition by gaussian elimination. Row pivoting is used.
237 *
238 * result : AP = L'U ; L' = LD; pivot[i] is the swapped column
239 * number for column i (for pivoting).
240 *
241 * Result is stored in place of original matrix.
242 *
243 */
LUDecomp(a,n,pivot)244 void LUDecomp(a, n, pivot)
245 double *a;
246 int n, pivot[];
247 {
248 double swap;
249 int i, j, k, l;
250
251 for (i = 0; i < n - 1; i++)
252 {
253 j = i;
254 for(k = i + 1; k < n; k++)
255 if (abs(A(i,k)) > abs(A(j,k))) j = k;
256
257 if (A(i,j) == 0.0)
258 {
259 error("LUDecomp: Matrix is singular.\n");
260 }
261
262 pivot[i] = j;
263
264 if (j != i)
265 {
266 swap = A(i,i);
267 A(i,i) = A(i,j);
268 A(i,j) = swap;
269 }
270
271 for(k = i + 1; k < n; k++)
272 A(i,k) = A(i,k) / A(i,i);
273
274 for(k = i + 1; k < n; k++)
275 {
276 if (j != i)
277 {
278 swap = A(k,i);
279 A(k,i) = A(k,j);
280 A(k,j) = swap;
281 }
282 for(l = i + 1; l < n; l++)
283 A(k,l) = A(k,l) - A(i,l) * A(k,i);
284 }
285 }
286
287 pivot[n - 1] = n - 1;
288 if ( A(n-1,n-1) == 0.0)
289 {
290 error("LUDecomp: Matrix is singular.\n");
291 }
292 }
293