6gjdn.c

#include <stdlib.h>
#include <math.h>
#include <stdio.h>

int gjdn(double *a, double *b, int n, int m)
{
	int *js, l, k, i, j, is = 0, p, q;
	double d, t;
	js = malloc(n * sizeof(int));
	l = 1;
	for (k = 0; k <= n - 1; k++) {
		d = 0.0;
		for (i = k; i <= n - 1; i++)
			for (j = k; j <= n - 1; j++) {
				t = fabs(a[i * n + j]);
				if (t > d) {
					d = t;
					js[k] = j;
					is = i;
				}
			}
		if (d + 1.0 == 1.0)
			l = 0;
		else {
			if (js[k] != k)
				for (i = 0; i <= n - 1; i++) {
					p = i * n + k;
					q = i * n + js[k];
					t = a[p];
					a[p] = a[q];
					a[q] = t;
				}
			if (is != k) {
				for (j = k; j <= n - 1; j++) {
					p = k * n + j;
					q = is * n + j;
					t = a[p];
					a[p] = a[q];
					a[q] = t;
				}
				for (j = 0; j <= m - 1; j++) {
					p = k * m + j;
					q = is * m + j;
					t = b[p];
					b[p] = b[q];
					b[q] = t;
				}
			}
		}
		if (l == 0) {
			free(js);
			fprintf(stderr, "ERROR: fail\n");
			return -1;
		}
		d = a[k * n + k];
		for (j = k + 1; j <= n - 1; j++) {
			p = k * n + j;
			a[p] = a[p] / d;
		}
		for (j = 0; j <= m - 1; j++) {
			p = k * m + j;
			b[p] = b[p] / d;
		}
		for (j = k + 1; j <= n - 1; j++)
			for (i = 0; i <= n - 1; i++) {
				p = i * n + j;
				if (i != k)
					a[p] = a[p] - a[i * n + k] * a[k * n + j];
			}
		for (j = 0; j <= m - 1; j++)
			for (i = 0; i <= n - 1; i++) {
				p = i * m + j;
				if (i != k)
					b[p] = b[p] - a[i * n + k] * b[k * m + j];
			}
	}
	for (k = n - 1; k >= 0; k--)
		if (js[k] != k)
			for (j = 0; j <= m - 1; j++) {
				p = k * m + j;
				q = js[k] * m + j;
				t = b[p];
				b[p] = b[q];
				b[q] = t;
			}
	free(js);
	return 0;
}