1 /**@file
2 *
3 * Simple Gauss elimination with total pivot search for non-degenerate
4 * quadratic matrices.
5 *
6 * author: Claus-Justus Heine
7 * Abteilung fuer Angewandte Mathematik
8 * Albert-Ludwigs-Universitaet Freiburg
9 * Hermann-Herder-Str. 10
10 * 79104 Freiburg
11 * Germany
12 * claus@mathematik.uni-freiburg.de
13 *
14 * Copyright (c) by Claus-Justus Heine (2007)
15 */
16
17 #include <alberta_util.h>
18
19 /**@addtogroup MISCELLANEOUS
20 * @{
21 */
22
23 static inline void
24 swap_rows(REAL *M, int r1, int r2, int n, REAL *b, int nbcol);
25 static inline void
26 swap_cols(REAL *M, int c1, int c2, int *perm, int n);
27
28 /** Simple Gauss elimination with total pivot search for
29 * non-degenerate quadratic matrices. BIG FAT NOTE: M and b are
30 * destroyed by this function.
31 *
32 * @param[in,out] M Quadratic NxN matrix in row-major storage order.
33 *
34 * @param[in] n Dimension of @a M.
35 *
36 * @param[in,out] b Right hand side in column major order.
37 *
38 * @param[in] nbcol Number of columns of @a b and @a x.
39 *
40 * @param[in,out] x Storage for the result. @a x is in column major
41 * order.
42 */
square_gauss(REAL * M,REAL * b,REAL * x,int n,int nbcol)43 void square_gauss(REAL *M, REAL *b, REAL *x, int n, int nbcol)
44 {
45 int perm[n];
46 REAL factor, pivot, max;
47 int maxcol, maxrow, col, row, i;
48
49 for (i = 0; i < n; i++) {
50 perm[i] = i;
51 }
52 for (i = 0; i < n; i++) {
53 /* find the pivot element */
54 max = 0.0;
55 maxrow = maxcol = i;
56 for (row = i; row < n; row++) {
57 for (col = i; col < n; col++) {
58 if (fabs(M[row*n+col]) > max) {
59 max = fabs(M[row*n+col]);
60 maxrow = row;
61 maxcol = col;
62 }
63 }
64 }
65 if (maxrow != i) {
66 swap_rows(M, i, maxrow, n, b, nbcol);
67 }
68 if (maxcol != i) {
69 swap_cols(M, i, maxcol, perm, n);
70 }
71 pivot = M[i*n + i];
72
73 /* clear i-th column of M */
74 for (row = i+1; row < n; row++) {
75 factor = M[row*n + i] / pivot;
76 for (col = i+1; col < n; col++) {
77 M[row*n + col] -= M[i*n + col] * factor;
78 }
79 /* also perform that operation on b */
80 for (col = 0; col < nbcol; col++) {
81 b[col*n + row] -= b[col*n + i] * factor;
82 }
83 }
84 }
85
86 /* at this stage the upper triangular part of M holds the upper
87 * triangular factor of the LR decomposition; just perform
88 * back-substitution to get the solution.
89 */
90 for (row = n-1; row >= 0; --row) {
91 for (col = 0; col < nbcol; col++) {
92 x[col*n + perm[row]] = b[col * n + row];
93 for (i = row+1; i < n; i++) {
94 x[col*n + perm[row]] -= M[row*n+i]*x[col*n + perm[i]];
95 }
96 x[col*n + perm[row]] /= M[row*n+row];
97 }
98 }
99 }
100
101 static inline void
swap_rows(REAL * M,int r1,int r2,int n,REAL * b,int nbcol)102 swap_rows(REAL *M, int r1, int r2, int n, REAL *b, int nbcol)
103 {
104 REAL tmp;
105 REAL *row1 = M+r1*n;
106 REAL *row2 = M+r2*n;
107 int col;
108
109 for (col = 0; col < n; col++) {
110 tmp = *row1;
111 *row1++ = *row2;
112 *row2++ = tmp;
113 }
114 for (col = 0; col < nbcol; col++) {
115 tmp = b[col*n + r1];
116 b[col*n + r1] = b[col*n + r2];
117 b[col*n + r2] = tmp;
118 }
119 }
120
121 static inline void
swap_cols(REAL * M,int c1,int c2,int * perm,int n)122 swap_cols(REAL *M, int c1, int c2, int *perm, int n)
123 {
124 REAL tmp;
125 int tmpi, i;
126
127 tmpi = perm[c1];
128 perm[c1] = perm[c2];
129 perm[c2] = tmpi;
130 for (i = 0; i < n; i++) {
131 tmp = M[i*n + c1];
132 M[i*n + c1] = M[i*n + c2];
133 M[i*n + c2] = tmp;
134 }
135 }
136
137 /**@} MISCELLANEOUS */
138
139 #ifdef TESTING
140
141 #include <stdlib.h>
142 #include <string.h>
143 #include <stdio.h>
144
145 #define DIM 5
146
Mv(const REAL * M,int n,const REAL * x,REAL * result,int nres)147 static inline Mv(const REAL *M, int n, const REAL *x, REAL *result, int nres)
148 {
149 int i, j, k;
150
151 for (k = 0; k < nres; k++) {
152 for (i = 0; i < n; i++) {
153 result[k*n+i] = 0.0;
154 for (j = 0; j < n; j++) {
155 result[k*n+i] += M[i*n+j]*x[k*n+j];
156 }
157 }
158 }
159 }
160
main(int argc,const char * argv[])161 int main(int argc, const char *argv[])
162 {
163 int i, j;
164 REAL M[DIM*DIM];
165 REAL MC[DIM*DIM];
166 REAL b[3*DIM];
167 REAL bC[3*DIM];
168 REAL x[3*DIM];
169
170 for (i = 0; i < 3; i++) {
171 for (j = 0; j < DIM; j++) {
172 b[i*DIM+j] = drand48()-0.5;
173 }
174 }
175 #if 1
176 for (i = 0; i < DIM; i++) {
177 for (j = 0; j < DIM; j++) {
178 M[i*DIM+j] = drand48()-0.5;
179 }
180 M[i*DIM+i] += 1.0;
181 }
182 #else
183 memset(M, 0, sizeof(M));
184 for (i = 0; i < DIM; i++) {
185 M[i*DIM+i] = i+1;
186 #if 0
187 for (j = 0; j <= i; j++) {
188 M[i*DIM+j] = j+1;
189 }
190 #endif
191 }
192 #endif
193 memcpy(MC, M, sizeof(M));
194 memcpy(bC, b, sizeof(b));
195
196 gauss(M, b, x, DIM, 3);
197
198 Mv(MC, DIM, x, b, 3);
199
200 printf("Matrix:\n");
201 for (i = 0; i < DIM; i++) {
202 for (j = 0; j < DIM; j++) {
203 printf("%e ", MC[i*DIM+j]);
204 }
205 printf("\n");
206 }
207
208 printf("RHS:\n");
209 for (i = 0; i < 3; i++) {
210 for (j = 0; j < DIM; j++) {
211 printf("%e ", bC[i*DIM+j]);
212 }
213 printf("\n");
214 }
215 printf("\n");
216 printf("Result:\n");
217 for (i = 0; i < 3; i++) {
218 for (j = 0; j < DIM; j++) {
219 printf("%e ", b[i*DIM+j]);
220 }
221 printf("\n");
222 }
223 }
224 #endif
225 /*
226 * Local Variables: ***
227 * mode: c ***
228 * c-basic-offset: 2 ***
229 * compile-command: "gcc -o gauss -DTESTING -O0 -ggdb3 gausselim.c -lm" ***
230 * End: ***
231 */
232
233