Matrix Invertierung

[W2K2005]

Hi,

ich bin total frustiert. Ich versuche diesen ver"...." Code für eine Matrix-Invertierung zu stabilisieren, aber das Programm stürtzt ständig ab.

Vielleicht kann mir jemand sagen, was zu machen ist bzw. eine Alternative nennen.

Vielen Dank.

int INV(double *A, double *Ainv, int dim) {
	int RetVal = 0;
	int ij, i, j, n, m, k;

	double temp;
	double *op_matrix = (double*) malloc(sizeof(double)*2*dim*dim);

	// Zurücksetzen des ZielVectors Ainv
	for (i=0; i<dim;i++) {
		for (j=0; j<dim; j++) {
			ij = i*dim + j;
			Ainv[ij] = 0;
		}
	}

	// Einheitsmatrix erzeugen
	// a11 a12 a13 || 1 0 0
	// a21 a22 a23 || 0 1 0
	// a31 a32 a33 || 0 0 1

	for (i=0; i<dim; i++) {
		for (j=0; j<dim; j++) {
			ij = 2*dim*i + j;
			op_matrix[ij] = A[i*dim+j];
		}
		for (j=0; j<dim; j++) {
			ij = (1+2*i)*dim + j;
			if (i == j) {
				op_matrix[ij] = 1;
			} else {
				op_matrix[ij] = 0;
			}
		}
	}

	for (k=0;k<dim; k++) {
		if ((op_matrix[k*(2*dim+1)] == 0) && (k < dim-1)) {
			// Suche nächste Zeile, die kein Null-Element enthält
			n=k;
			printf("n:%d\n", n);
			do {
				n++;
			} while (op_matrix[n*(2*dim)+k] == 0);
			// Austauschen der Zeilen k und n vornehmen
			for (m=k; m<2*dim-k; m++) {
				temp = op_matrix[k*(2*dim)+m];
				op_matrix[k*(2*dim)+m] = op_matrix[n*(2*dim)+m];

				printf("%d -- %d\n", m, n*(2*dim)+m);
				op_matrix[n*(2*dim)+m] = temp;
			}
			printf("----\n");
		}

		// Zeile mit dem jeweiligen Diagonalelement normieren
		temp = op_matrix[k*(2*dim)+k];
		if (temp != 0) {
			for (n=k; n<2*dim; n++) {
				op_matrix[k*(2*dim)+n] = op_matrix[k*(2*dim)+n] / temp;
			}
		}

		// Null-Element erzeugen
		for (n=0;n<dim;n++) {
			if ((n == k) && (n < 2*dim)) {
				n = n + 1;
			}
			if (op_matrix[n*(2*dim)+k] != 0) {
				temp = op_matrix[n*(2*dim)+n] / op_matrix[k*(2*dim)+k];
				for (m=k;m<2*dim;m++) {
					op_matrix[n*(2*dim)+m] = op_matrix[n*(2*dim)+m] - op_matrix[k*(2*dim)+m]*temp;
				}
			} else {
				break;
			}
		}
	}

	if (RetVal == 0) {
		for (i=0;i<dim;i++) {
			for (j=0;j<dim;j++) {
				ij = i*dim + j;
				Ainv[ij] = op_matrix[(i*2+1)*dim+j];
			}
		}
	}

	free(op_matrix);

	return RetVal;
}

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

void main(void) {
	const int dim = 3;

	int i;
	int RetVal;
	double Res;

	double *matrixVec = (double*) malloc(sizeof(double)*dim*dim);
	double *matrixVecInv = (double*) malloc(sizeof(double)*dim*dim);

	// Schreibe die einzelnen Elemente in das Datenfeld matrixVec
	matrixVec[0] = 0;
	matrixVec[1] = 0;
	matrixVec[2] = 0;
	matrixVec[3] = 1;
	matrixVec[4] = 2;
	matrixVec[5] = 3;
	matrixVec[6] = 4;
	matrixVec[7] = 5;
	matrixVec[8] = 6;

	printf("\n> Invertierung von A:\n-------------------------\n");

	RetVal = INV(matrixVec, matrixVecInv, dim);

	if (RetVal == 0) {
		for (i=0; i<dim*dim; i++) {
			printf("%.2lf\t", matrixVecInv[i]);
			if ((i+1)%dim == 0) {
				printf("\n");
			}
		}
	} else {
		printf("> ERROR: Die Matrix ist nicht invertierbar...\n\n");
	}

	free(matrixVec);
	free(matrixVecInv);
}

[c.rackwitz]

probier mal structs, koennte manches uebersichtlicher machen.

[W2K2005]

Wozu?

Ich sehe den Sinn nicht. Außerdem sollte es auch so funktionieren.

[c.rackwitz]

W2K2005 schrieb:

Außerdem sollte es auch so funktionieren.

tut es das?