Innan så hade inv() kommandot olika anrop med dynamiska matriser som sedan frigjordes. Det var detta som många var kritiska till, att man använder malloc och free i t.ex. en for-loop.
Kod: Markera allt
void inv(matrix* a, matrix* b) {
/*
* This uses LU decomposation to find the inverse of matrix a
*/
int n = a->row;
float vectorI[n];
float x[n];
memset(x, 0, sizeof(x)); // reset
float* ptr_out = b->data;
for (int i = 0; i < n; i++) {
// Create one column identity "matrix"
vectorI[i] = 1;
// Solve
linsolve_internal_inv(a, vectorI, x);
// Insert
for (int j = 0; j < n; j++)
*((ptr_out + j * n) + i) = x[j];
// Reset
memset(vectorI, 0, sizeof(vectorI));
memset(x, 0, sizeof(x));
}
}
/*
* This is internal linear solving procedure for just this inv-command.
*/
void linsolve_internal_inv(matrix*a, float* b, float* c) {
/*
* This uses LU-factorization to solve Ax = b equation
*/
// Get data for matrix a
int n = a->row;
int m_a = a->column;
// Do LU-factorization
float l[n*m_a];
float u[n*m_a];
lu_internal_inv(a, l, u);
// Initial sum
float sum = 0;
float y[n*m_a]; // Vector
///float* ptr_x = c->data; // our output
// Solve Ly = b with gaussian elimination
for (int i = 0; i < n; i++) {
sum = 0;
for (int j = 0; j < i; j++)
sum += *((l + i * n) + j) * (*(y + j));
*(y + i) = (*(b + i) - sum) / (*((l + i * n) + i));
}
// Solve Ux = y with gaussian elimination
for (int i = n-1; i >= 0; i--) {
sum = 0;
for (int j = i; j < n; j++)
sum += *((u + i * n) + j) * (*(c + j));
*(c + i) = (*(y + i) - sum) / (*((u + i * n) + i));
}
}
/*
* This is LU-factorization for internal for just this inv-command.
*/
void lu_internal_inv(matrix* a, float* l, float* u) {
// Get info about matrix 'a'
int n = a->row;
int m = a->column;
float* ptr_a = a->data;
float sum = 0;
// Initial values for matrix u and l
memset(u, 1, m * n * sizeof(float));
memset(l, 0, m * n * sizeof(float));
for (int j = 0; j < n; j++) {
for (int i = j; i < n; i++) {
sum = 0;
for (int k = 0; k < j; k++) {
sum = sum + *((l + i*n) + k) * (*((u + k*n) + j));
}
*((l + i*n) + j) = *((ptr_a + i*n) + j) - sum;
}
for (int i = j; i < n; i++) {
sum = 0;
for (int k = 0; k < j; k++) {
sum = sum + *((l + j*n) + k) * (*((u + k*n) + i));
}
if (*((l + j*n) + j) == 0) {
printf("det(L) close to 0!\n Can't divide by 0...\n");
*((l + j*n) + j) = pow(2.2204, -16); // Same as eps command in MATLAB
}
*((u + j*n) + i) = (*((ptr_a + j*n) + i) - sum) / (*((l + j*n) + j));
}
}
}
