Matrisberäkningar med för STM32?

[Al_Bundy]

Är det inte bara ctrl+c och ctrl+v ? Det är ju det som är tanken med projektet. Att man ska bara dumpa en mapp med massa funktioner i olika utvecklingsverktyg t.ex. Atollic TrueSTUDIO eller Kiel, eller IAR osv.

Ladda ned Eclipse

[hummel]

Om du tänker på din Zip-fil så innehåller den dina .c och .h filer men ingen projektfil med inställningar (eller en make-fil).

[Al_Bundy]

Ska inte behövas. Länka endast "-lm" för mattematiken.

[hummel]

För att få samma exakta inställningar som du har behövs projektfilen. Kanske är den hemlig för ditt test, vad vet jag. Varför jag skriver i den här tråden är för att hjälpa till, inte för att vara elak.

[Al_Bundy]

Den är inte hemlig. Hela mitt projekt finns ju fritt ute.

Men varför inte testa att bara kopiera in mappen vid din main.c fil? Det är så jag gjorde på Atollic TrueSTUDIO och sedan länkande jag -lm i IDE:n.

Svårare än så var det inte.

[Al_Bundy]

Nya uppdateringar.

Kod: Markera allt

void copyval(matrix*a, matrix* b){
	memcpy(b->data, a->data, b->column * b->row * sizeof(float));
}

Kod: Markera allt

matrix* zeros(int n, int m) {

	matrix* out = initMatrix(n, m); // Already zero matrix due to memset
	return out;
}

Kod: Markera allt

matrix* cofact(matrix* a) {

	int n = a->row;
	float* ptr_a = a->data;

	// Create our matrix - cofactor
	matrix* out = initMatrix(n, n);
	float* ptr_out = out->data;

	// Create the minor matrix - Allways one size smaller
	matrix* minor = initMatrix(n-1, n-1);
	float* ptr_minor = minor->data;

	// 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
					*(ptr_minor + minorCellposition) = *((ptr_a + i * n) + j);
					minorCellposition++; // Next cell
				}
			}
		}

		// Insert in our output
		*(ptr_out++) = pow(-1, row+column)* det(minor); //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
		resetMatrix(minor);
		minorCellposition = 0;

	}

	// Delete the minor matrix - we don't need it
	freeMatrix(minor);  // <---- HÄR!

	return out;
}

Kod: Markera allt

#ifndef LINEARALGEBRA_DECLAREFUNCTIONS_H_
#define LINEARALGEBRA_DECLAREFUNCTIONS_H_

#endif /* LINEARALGEBRA_DECLAREFUNCTIONS_H_ */

// Matrix structure
typedef struct _matrix {
    int row;
    int column;
    float* data;
} matrix;

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


/* Use these when you using the library */
matrix* initMatrix(int row, int column);
void freeMatrix(matrix* a);
void printMatrix(matrix* a);
matrix* cut(matrix* a, int start_n, int stop_n, int start_m, int stop_m);
matrix* horzcat(matrix* a, matrix* b);
matrix* vertcat(matrix* a, matrix* b);
float det(matrix* a);
matrix* cofact(matrix* a);
matrix* create(float* arr2D, int n, int m);
matrix* eye(int n, int m);
matrix* tran(matrix* a);
matrix* add(matrix* a, matrix* b);
matrix* sub(matrix* a, matrix* b);
void size(matrix* m, int* row, int* column);
matrix* scale(matrix* m, float value);
matrix* mul(matrix* a, matrix* b, bool elementWise);
matrix* ones(int n, int m);
matrix* zeros(int n, int m);
matrix* diag(float* arr, int length);
matrix* inv(matrix* a);
matrix* linsolve(matrix*a, matrix* b);
float norm(matrix* a);
void qr(matrix* a, matrix* q, matrix* r);
float dot(matrix* a, matrix* b);
void svd(matrix* a, matrix* u, matrix* s, matrix* v);
void copyval(matrix*a, matrix* b);
matrix* triu(matrix* a, int shift);
matrix* tril(matrix* a, int shift);
matrix* vec(matrix* a);
matrix* diagm(matrix* a);
matrix* chol(matrix* a);
void lu(matrix* a, matrix* l, matrix* u);
matrix* repmat(matrix* a, int i, int j);
void max(matrix* a, float* val, int* index);
void min(matrix* a, float* val, int* index);
matrix* summ(matrix* a);
matrix* absm(matrix* a);
matrix* sqrtm(matrix* a);
matrix* hankel(float* arr, int length, int step);
matrix* toeplitz(float* arr, int length);
matrix* pinv(matrix* a);
void resetMatrix(matrix* a); // <-- Här!

Kod: Markera allt

matrix* repmat(matrix* a, int i, int j) {

	matrix* horz;
	matrix* vertz;
	matrix* mat;

	// Copy horizontal
	if (i != 0) {
		for (int k = 0; k < i; k++) {
			if (k == 0) {
				horz = horzcat(a, a);
			} else {
				mat = horzcat(horz, a);
				freeMatrix(horz);
				horz = scale(mat, 1); // Copy
				freeMatrix(mat);

			}
		}
	} else {
		horz = scale(a, 1); // Just copy
	}

	// Copy vertical
	if (j != 0) {
		for (int k = 0; k < j; k++) {
			if (k == 0) {
				vertz = vertcat(horz, horz);
			} else {
				mat = vertcat(vertz, horz);
				freeMatrix(vertz);
				vertz = scale(mat, 1); // Copy
				freeMatrix(mat);

			}
		}
	}else{
		vertz = scale(horz, 1); // Just copy
	}

	freeMatrix(horz);  
        // <--- Tog bort printMatrix
	return vertz;
}

Ny fil.

Kod: Markera allt

void resetMatrix(matrix* a){
	memset(a->data, 0, a->column * a->row * sizeof(float));
}

[Al_Bundy]

Nu körde jag SVD på en 36x36 matris....1.867489 sekunder. Verkar som jag har skrivit ett segt bibliotek för SVD.

Kod: Markera allt

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

int main() {

	/*
	 * G(s) = 3.2/(2s^2 + 0.7s + 3.1) - Model
	 */

	float input[72] = {  5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,

			   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,

			   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   5};


	float output[72] = { 0.00000,   1.89649,   3.30777,   3.40331,   3.23615,   3.18151,   3.19049,   3.20013,   3.20132,   3.20031,   3.19990,


			   3.19993,   3.20000,   3.20001,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,


			   3.20000,   3.20000,   3.20000,   3.20000,   3.20000,   3.20000};



	// Create toeplitz matrix
	matrix* toe = toeplitz(input, 72);

	// Create upper triangular matrix of the toeplitz matrix
	matrix* tru = triu(toe, 0);

	// Find the inverse of tru
	matrix* iv = inv(tru);

	// Create vector the horizon - Important! Else, we cannot find the markov-parameters g
    matrix* Y = create(output, 1, 72);

	// Multiply Y with the iv matrix - Find the markov-parameters g
    matrix* G = mul(Y, iv, false);

    // Detete some mats
    freeMatrix(toe);
    freeMatrix(tru);
    freeMatrix(iv);
    freeMatrix(Y);

    /*--------------------------------------*/

    // turn vector G into a normal vector because the function hankel only want 1D vector
    float g[72];
    for(int i = 0; i < 72; i++)
    	g[i] = *(G->data + i);

    // Create hankel matrix H0 and H1
    matrix* H0 = hankel(g, 72, 0);
    matrix* H1 = hankel(g, 72, 1);

    // Cut H1 and H2 to the half - Remember indexing from zero here!
    matrix* H0_half = cut(H0, 0, 35, 0, 35); // 36x36 from 72x72
    matrix* H1_half = cut(H1, 0, 35, 0, 35);

    // Remove H0
    freeMatrix(H0);

    // Measure the size

    /*
    int n;
    int m;
    size(H0_half, &n, &m);
    printf("The H0_half have the size %dx%d\n\n,", n, m);
    */


    // Do SVD on
    matrix* u =  initMatrix(36,36); // We know the size from the size() command above
    matrix* s =  initMatrix(36,36);
    matrix* v =  initMatrix(36,36);
    clock_t start, end;
    float cpu_time_used;

    start = clock();
    svd(H0_half, u, s, v);
    end = clock();
    cpu_time_used = ((float) (end - start)) / CLOCKS_PER_SEC;
    printf("Total speed of SVD was %f,", cpu_time_used);


    //printMatrix(u);

	return EXIT_SUCCESS;
}

GNU Octave är så här snabbt.

Kod: Markera allt

>> tic; [u,s,v] = svd(H0); toc
Elapsed time is 0.000340939 seconds.
>>

Nu får vi inte glömma att GNU Octave gör sina tuffaste beräkningar i Fortran, vilket är snabbare än C.

[sodjan]

> ...i Fortran, vilket är snabbare än C

Det stämmer ju självklart inte. Inte generellt. Det betyder
betydligt mer vilken plattform man kör på, vilket float stöd
plattformen har samt vilket float format som plattformen och
kompilatorerna använder och som applikationen använder.

Om man väljer ett float format som plattformen inte stöder
direkt i hårdvaran så blir det emulering i programvaran och
prestanda går "up through the roof".

Fortran har dock traditionellt haft kompilatorer som har haft
effektiva rutiner för just beräkningar och med bra stöd för
plattformarnas flyttalsstöd i hårdvaran. Men det är inget
som hindrar att man gör samma sak i C.

[bearing]

Använder koden verkligen processorns flyttalsoperationer om det tar så lång tid?

[Al_Bundy]

Men nu körde jag på en vanlig PC.

Det är väll inte for-looparna som segar ned? Om det är det, då vet jag inte alls vad jag ska göra för att optimera. Jag kör ju endast for-loopar i alla funktioner och alla funktioner ser likadana ut.

Det vore tragiskt om jag tvingas gå tillbaka till reglersystem som vilar på ett operativsystem. Risken är ojämn samplingstid är hög, vilket orsakar ostabilitet.

Det är denna kod i svd.c filen som suger kraften. Testa att undvika köra for-loopen och man kommer sjukt lågt ned i tid. Logiskt dock! Men det bekräftar att det är for-looparna som suger kraften.

Kod: Markera allt

for(int i = 0; i < 100; i++){
		// Solve QR from transpose of s
		t = tran(s);
		qr(t, q, r);
		freeMatrix(t);

		// Copy values from r to s
		copyval(r, s); // We do not need to delete r matrix yet

		// Copy u = u*q
		t = mul(u, q, false); // Not element wise
		copyval(t, u);
		freeMatrix(t);

		// Solve QR from transpose of s
		t = tran(s);
		qr(t, q, r);
		freeMatrix(t);

		// Copy values from r to s
		copyval(r, s); // We do not need to delete r matrix yet

		// Copy v = v*q
		t = mul(v, q, false);
		copyval(t, v);
		freeMatrix(t);

	}

[mankan]

Du körde alltså ditt testprogram på samma maskin som du kör Octave på?

Nu får du chansen att lära dig profilerare och då kommer du förhoppningsvis inse vad jag och andra har påpekat när det gäller minneshanteringen.

Kul att du har något som kör.

[Al_Bundy]

Jag körde SVD parallellt via GNU Octave. Supersnabbt.

Det har inte med minneshanteringen att göra, utan for-loopar. Det var JAG som påpekade detta innan, att for-loopar är ej vektoriserad kod, men många här slog mig på fingrarna och sade att det spelar ingen roll.

Så for-loopar måste optimeras på något sätt.

Jag ser också att qr.c är boven i dramat. Utan den så kommer jag ned till några hudradelar av en sekund.

[sodjan]

> Jag körde SVD parallellt via GNU Octave. Supersnabbt.

Alltså, det var *Octave* var "supersnabbt"?

> Det har inte med minneshanteringen att göra, utan for-loopar. Det var JAG som påpekade detta innan,
> att for-loopar är ej vektoriserad kod, men många här slog mig på fingrarna och sade att det spelar ingen roll.

Jag menar att du är helt ute ock cyklar. Bara för att dina C kod *ser* enklare ut så betyder
inte det att plattformen du kör på stöder vektorberäkningar. Någonstans i en lib körs fortfarande
samma loopar. Visst, den som skrev det kan ju ha varit duktigare på optimering än vad du är,
men om din algoritm som sådan är korrekt så måste ju ända samma operationer utföras.

> Så for-loopar måste optimeras på något sätt.

Ja, det finns vissa problem med att cachen kan bli ineffektiv om man läser igenom
en array i "fel ordning", så att säga. Man bör läsa/loopa så som den är lagrad i minnet.

> Jag ser också att qr.c är boven i dramat. Utan den så kommer jag ned till några hundradelar av en sekund.

Jag har inte kollat vad qr.c gör, men får du alltså fortfarande rätt resultat utan den?
Om inte får det så spelar väl körtiden utan den inte så stor roll...

Jag tror att skillnaderna du ser antingen beror på ineffektiv kod rent generellt eller att
du använder en float typ som inte är optimal på din plattform.

[gkar]

Nu har jag inte orkat bygga koden, men:

1) Kolla så att du genererar flyttalskod om du har flyttalshårdvarustöd, annars får du emulerad kod som aldrig blir klar.

2) Din kod har troligen pointer aliasing problem.
https://en.wikipedia.org/wiki/Pointer_aliasing

Tex detta:

Kod: Markera allt

for (int k = 0; k < a->column; k++) {
	*ptr += *data_a * *data_b;
	data_a++;
	data_b += b->column;
}

Undvik att ha skrivpekare i samma flöden som du har läspekare.
ptr är ett problem eftersom kompilatorn troligtvis inte vet om ptr och data_a eller data_b är samma adress.
Använd lokal variabel istället.

Jag avråder från att använda restrict (keyword)
https://en.cppreference.com/w/c/language/restrict

[bearing]

Men nu körde jag på en vanlig PC.

Det innebär inte att kompilatorn automatiskt genererar flyttalsinstruktioner.