Matrisberäkningar med för STM32?

[Rick81]

Det jag försökte förklara är att GNU octave troligen använder den mer avancerad algoritm med bättre noggrannhet.

[hawkan]

Jag har också sagt det

[Al_Bundy]

Okej!

Men jag lägger ned felsökandet nu då jag inte hittar om det ens är ett fel.

Hur som helst. Nu tänker jag skriva om mitt C-bibliotek och då har jag några frågor:
1. Ska jag ha felmeddelanden som leder till exit(1) ?
2. Hur blir det med matriser som deklareras utan malloc? typ matris[20][40] och nästa gång är det matris[20][41] ?

[AndLi]

Nej du ska inte ha några exit(1), returnera en felkod...

Du skulle ju aldrig ha större matriser än 36x36?

Du måste förklara när du tänker att denna växling i matrisstorlek sker.
Sker den vid implementation är det inga problem, sker det runtime får den ju helt enkelt vara så stor som den max måste vara.

[Al_Bundy]

För ETT program ja. Jag kanske ska göra andra program också.

Vad tycker ni om denna implementering?

Kod: Markera allt

#include <stdio.h>
#include <stdlib.h>
#include "LinearAlgebra/declareFunctions.h"
#include <time.h>

int main() {


	clock_t start, end;
	float cpu_time_used;

	start = clock();

	/*
	 * Create an empty initial matrix with 5 rows and 5 columns
	 */
	printf("Init a matrix\n");
	matrix* A = (matrix*) malloc(sizeof(matrix)); // <--- kan man göra något här så man slipper "malloca" för varje matris?
	init(A, 5,5);


	/*
	 * Insert values in to the matrix
	 */
	float values[5][5] = { { 1, 76, 86, 1, 5 },
				{ 1, 5, -6, 0, 3 },
				{ -5, 7, 3, 87, 3 },
				{ -8, 3, 1, 4, -1 },
				{ 7, 9, 1, 28, 4 } };

	printf("Create matrix\n");
	insert(*values, A);
	print(A);



	end = clock();
	cpu_time_used = ((float) (end - start)) / CLOCKS_PER_SEC;
	printf("Total speed  was %f,", cpu_time_used);

	return EXIT_SUCCESS;

}

Kod: Markera allt

void init(matrix* a, int row, int column) {


	a->column = column;
	a->row = row;
	a->data = (float*) malloc(sizeof(float) * column * row);

	// Set all elements to 0
	memset(a->data, 0, column * row * sizeof(float));

}

Kod: Markera allt

void insert(float* arr2D, matrix* a) {

	int n = a->row;
	int m = a->column;
	float* ptr = a->data;

	// Insert all values
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < m; j++) {
			*(ptr++) = *((arr2D + i * m) + j);
		}
	}


}

Kod: Markera allt

void print(matrix* a) {
	float* ptr = a->data;
	for (int i = 0; i < a->row; i++) {
		for (int j = 0; j < a->column; j++) {
			printf(" %9.6f", *(ptr++));
		}
		printf("\n");
	}
	printf("\n"); // Default new row
}

Tror ni den är bra, eller kan man göra bättre?

Eller är det bättre om init såg ut så här. Där man bara anger en dimensionen och sedan får man tillbaka adressen för matrisen.

Kod: Markera allt

matrix* initMatrix(int row, int column) {
	matrix* out = (matrix*) malloc(sizeof(matrix));

	out->column = column;
	out->row = row;
	out->data = (float*) malloc(sizeof(float) * column * row);

	memset(out->data, 0.0, column * row * sizeof(float));

	return out;
}

[mankan]

Skrolla upp och kolla vad jag föreslog med stackallokering.

[AndLi]

Och byter du program byter du ju bara vektorstorlek till det du ska ha.
Och med lite korrekt utformade #defines behöver du bara ändra på ett ställe...

[Al_Bundy]

Jag tror jag har hittat ett bra förslag. Ska presentera detta sedan.

[Al_Bundy]

Vad tycker ni om denna lösning?

Indata är matris a och utdata är matris co. Dessa två matriser är skapade via malloc. Men det finns en 1D array som är inte allokerad med malloc, men den deklareras ändå och kommer ha olika storlekar beroende på vilken matris a och co som stoppas in.

Detta är kofatormatisen om ni undrar.

Om ni en ekvation:
\(Ax = b\)

Där \(A\) är dina ekvationer och \(x\) är alla skalärer och \(b\) är uppmätt data. Då kan du hitta inversen av \(A\), dvs \(A^{-1}\) igenom denna
\(A^{-1} = \frac{1}{|A|} (co)^T\)

Och därmed så kan du hitta \(x\)
\(x = A^{-1}b\)

Kod: Markera allt

float det_cofact_internal();

void cofact(matrix* a, matrix* co) {

	int n = a->row;
	float* minor_a = a->data;
	float* minor_co = co->data;

	// Create the minor matrix - Always one size smaller
	float minor[(n-1)*(n-1)];

	// Starting position
	int row = 0;
	int column = 0;
	int minorCellposition = 0;

	// What we do is to create sub matrices for every cell. 3x3 means 9 sub matrices
	for (int k = 0; k < n * n; k++) {

		// Loop rows
		for (int i = 0; i < n; i++) {
			// Loop columns
			for (int j = 0; j < n; j++) {
				// Check if we are not on the "forbidden" row and column
				if (i != row && j != column) {
					// Insert
					minor[minorCellposition] = *((minor_a + i * n) + j);
					minorCellposition++; // Next cell
				}
			}
		}

		// Insert in our output
		*(minor_co++) = pow(-1, row+column)* det_cofact_internal(minor, n-1); //pow is needed because turning minor matrix to cofactor matrix

		// Update - Remember n-1 is due to indexing from zero
		if(column < n-1){
			column++; // Next column
		}else{
			column = 0;
			row++; // Next row
		}

		// Reset the minor matrix and jump back to cell position 0
		memset(minor, 0, sizeof(minor));
		minorCellposition = 0;

	}





}

/*
 * This is the same as det function, but the argument is diffrent!
 */

float det_cofact_internal(float* minor, int n) {


	if (n == 1) {
		// 1x1 "matrix" - Just return the value
		return *(minor);
	} else if (n == 2) {
		// 2x2 matrix - Easy
		return *(minor + 0) * *(minor + 3) - *(minor + 1) * *(minor + 2);
	} else if (n == 3) {

		// Don't change this! It works fine for 3x3 - Took me long time to write this
		return (*(minor + 0)) * (*(minor + 4)) * (*(minor + 8))
				+ (*(minor + 1)) * (*(minor + 5)) * (*(minor + 6))
				+ (*(minor + 2)) * (*(minor + 3)) * (*(minor + 7)) -

		(*(minor + 6)) * (*(minor + 4)) * (*(minor + 2))
				- (*(minor + 7)) * (*(minor + 5)) * (*(minor + 0))
				- (*(minor + 8)) * (*(minor + 3)) * (*(minor + 1));

	} else {


		float ratio;
		float determinant = 1; // Initial det value

		// Gauss Elimination Easy and fast algorithm to find the determinant of a nxn matrix
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				if (j > i) {
					// We cannot divide with zero
					if (*((minor + i * n) + i) == 0){
						*((minor + i * n) + i) = pow(2.2204, -16); // Same as MATLAB command eps
					}
					ratio = *((minor + j * n) + i)/ (*((minor + i * n) + i));
					for (int k = 0; k < n; k++) {
						*((minor + j * n) + k) -= ratio * (*((minor + i * n) + k));
					}
				}
			}
		}

		// Find det
		for (int i = 0; i < n; i++) {
			determinant *= *((minor + i * n) + i);
		}

		return determinant;
	}

}

Denna C-kod är ekvivalent med denna MATLAB kod:

Kod: Markera allt

function cof = cof(a);
[r c] = size(a);    %determine size of input           
m = ones(r,c);      %preallocate r x c cofactor matrix        
a_temp=a;           %create temporary matrix equal to input
for i = 1:r
    for k = 1:c
        a_temp([i],:)=[];   %remove ith row
        a_temp(:,[k])=[];   %remove kth row
        m(i,k) = ((-1)^(i+k))*det(a_temp);  %compute cofactor element
        a_temp=a;   %reset elements of temporary matrix to input elements
    end  
end

cof=m;  %return cofactor matrix as output variable

end

[guckrum]

Utan att försöka förstå din kod i detalj kan jag säga att din "gausseliminator" kommer att få problem med numerisk stabilitet i vissa fall om du inte pivoterar eller på annat sätt tar hänsyn till numeriken. Jag är inte den första i tråden att påpeka denna potentiella felkälla.

[Al_Bundy]

Du menar att jag ska använda LU fakrorisering?

[guckrum]

Nej, det har inget att göra med vad du gör för operationer på dina matriser,
utan hur dessa operationer är implementerade "under skalet". Matematiken
är ideal - implementationen av den är det sällan. Begränsad precision i
praktiken hanteras genom "listiga" implementationer. Googla, tex pivotering!

[Al_Bundy]

Så vad händer om jag inte tar hänsyn till pivoting?

[Rick81]

Gissar att du får sämre noggrannhet i beräkningen, kanske det du ser när du jämför.

Jag är fortfarande nyfiken på hur du tar fram systemets ekvation. Du skrev nåt om machine learning, kan du förklara hur denna process går till?

[Icecap]

Som jag fattar det ger dessa matrix-beräkningar en bild av systemets reaktion på input vilket då blir en typ matematisk beskrivning av regleringen.

Alltså inget mashine-learning eller liknande, bara gammal hederlig matematik - som förvisso kan vara mycket effektiv.

Jag antar att man kan ha "kvar" en matrix som summerar under tiden och som är "bilden" av systemet, att dynamisk bild som anpassas eftersom.

Men det är bara vilda gissningar.