Matrisberäkningar med för STM32?

[Al_Bundy]

Hur var det man gjorde profilering i C hos Eclipse?

Kompilera med -pg
Kör programmet
Kör gprof

Nu.
Jag har fått ned körningen till 0.035 sekunder nu.

[hummel]

Det som är avgörande antalet instruktioner (för att vara exakt är det antalet klockcykler din funktion tar) din funktion generar. Titta på assemblerkoden och jämför när du gör dessa ändringar.

[bearing]

Tänker mig något liknande detta.

Kod: Markera allt

//Innan loopen:
data = data + start_n * n + start_m;

...

//Och sedan i loopen:
memcpy(ptr, data, out_m);
data = data + in_m;
ptr = ptr + out_m;

Men då är det ju loop iallafall. Jag förstår inte riktigt hur denna kan vara effektivare än någon enkel for-loop som itererar mindre nu.

Du hade ju en nästlad loop. Ett varv per element. Med ovanstående bara ett varv per rad. Prova och jämför och rapportera tillbaka hur mycket skillnad det gav.

[Al_Bundy]

Förstår fortfarande inte hur detta ska fungera.

[bearing]

Det kan hända att jag vänt på n och m någonstans.

Hur som helst - tidigare:
Nästlad loop
För varje rad:
Och för varje element i denna rad:
Kopiera element 2, 3, 4, 5 osv

Med memcpy:
Ej nästlad loop
För varje rad:
Kopiera element 2-5 med memcpy

Men enligt din profilering tar ju cut inte längre någon tid? Från 95% till 0%?

[Al_Bundy]

Fast 2D matriserna fungerar inte riktigt som du har tänkt det. 2D matriserna är raka som en 1D array. Så två for-loopar behövs.

Men jag ska göra ett försök!

Edit:

Kod: Markera allt

matrix* cut(matrix* a, int start_n, int stop_n, int start_m, int stop_m) {

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

	// Create the output
	matrix* out = initMatrix(stop_n - start_n + 1, stop_m - start_m + 1);
	float* ptr = out->data;

	for (int i = start_n; i < stop_n + 1; i++) {
		memcpy(ptr++, data + i*n + start_m, sizeof(stop_m+1));
	}

	return out;
}

[Al_Bundy]

Resultat: Total speed of SVD was 0.036362 sekunder. Det blev inte sådan stor skillnad.

Nu är det mul.c och cut.c som drar mest minne. Cut.c kan jag inte optimera längre. Den är bra som den är. Men mul.c behövs göras något åt.

Kod: Markera allt

matrix* mul(matrix* a, matrix* b, bool elementWise) {
	// Get the dimensions - they both are the same
	int row_a = a->row;
	int column_a = a->column;
	//int row_b = b->row;
	int column_b = b->column;

	// Data
	float* data_a = a->data;
	float* data_b = b->data;

	if (elementWise == true) {
		matrix* out = initMatrix(row_a, column_a); // Same dimensions for a and b matrix
		float* ptr = out->data;
		// Dot multiply all values
		if (column_b > 1) { // If matrix b is a matrix
			for (int i = 0; i < row_a; i++) {
				for (int j = 0; j < column_a; j++) {
					// Do element wise mutiplication. In MATLAB it is A.*A
					*(ptr++) = *(data_a++) * *(data_b++);
				}
			}
		} else {
			// If matrix b is a vector
			for (int i = 0; i < row_a; i++) {
				for (int j = 0; j < column_a; j++) {
					// Do element wise mutiplication. In MATLAB it is A.*b
					*(ptr++) = *(data_a++) * *(data_b + i);
				}
			}
		}

		return out;
	} else {
		// Do regular mutiplication. In MATLAB it is A*A
		matrix* out = initMatrix(row_a, column_b);
		float* ptr = out->data;

		// Let's take our a matrix
		for (int i = 0; i < a->row; i++) {

			// Then we go through every column of b
			for (int j = 0; j < b->column; j++) {
				data_a = &a->data[i * a->column];
				data_b = &b->data[j];

				*ptr = 0; // Reset
				// And we multiply rows from a with columns of b
				for (int k = 0; k < a->column; k++) {
					*ptr += *data_a * *data_b;
					data_a++;
					data_b += b->column;
				}
				ptr++;
			}
		}

		return out;
	}

}

[mankan]

Fast 2D matriserna fungerar inte riktigt som du har tänkt det. 2D matriserna är raka som en 1D array. Så två for-loopar behövs.

Men jag ska göra ett försök!

Edit:

Kod: Markera allt

matrix* cut(matrix* a, int start_n, int stop_n, int start_m, int stop_m) {

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

	// Create the output
	matrix* out = initMatrix(stop_n - start_n + 1, stop_m - start_m + 1);
	float* ptr = out->data;

	for (int i = start_n; i < stop_n + 1; i++) {
		memcpy(ptr++, data + i*n + start_m, sizeof(stop_m+1));
	}

	return out;
}

Funkar verkligen det där? Du skriver ju över dig själv flera gånger om eftersom du bara flyttar ptr 1 steg per varv. sizeof(stop_m+1) ser också skumt ut, borde vara antal floats*sizeof(float).

[Al_Bundy]

Japp. Det fungerar.

Testa vet jag

Jag lägger upp min senaste version som är mer optimerad. Jag tänkte fråga om någon kan kolla hur mycket minne som det läcker här. Trots att jag lägger till fler freeMatrix(matrix* a) funktioner, så ändras inte antalet hur många gånger jag anropar freeMatrix(matrix* a) i Eclipse.

Jag tror att jag nöjer mig med att laga minnesläkaget, om det är något sådant. Jag tror inte jag personligen, kan optimera något mera. Det skulle kräva väldigt mycket teknik att få matriserna att flyta på riktigt bra. Troligtvis så kanske jag måste skriva om allt för mycket så det inte alls blir MATLAB-likt. Jag har ju fått ner min SVD beräkning med 36x36 matris till 0.023663 sekunder.

För GNU Octave tar det.

Kod: Markera allt

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

Alltså 20 gånger snabbare. Nu får vi inte glömma att GNU Octave har över 30 års erfarenhet när det kommer till optimering och Octave har ju väldigt mystiska beräkningssätt i Fortran och C++, för att verkligen "rappa på" beräkningen.

EDIT: Jag har fixat nu repmat.c så den inte ger några varningar.

[Mr Andersson]

Den där profileringen ser konstig ut. Du har 2 samples över hela programmet, och bara 2 funktioner som används?

[Al_Bundy]

Kan inte direkt säga att jag är bra på att profilera koden, trots att det sker endast med ett knapptryck.

Testa om ni kan köra denna kod. Jag får denna vid 0.023731 sekunder när jag har flaggan -O2 i kompileringen.

Jag anropar freeMatrix 130166 gånger.
Jag anropar initMatrix 130174 gånger. Alltså jag initialiserar 8 gånger mer än vad jag tar bort. Det betyder att minnet läker. Men kan inte hitta vad jag har missat.

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);

    // Delete G
    freeMatrix(G);

    // 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

    // 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);


    // Delete H1 and H0_half
    freeMatrix(H1);
    freeMatrix(H0_half);
    freeMatrix(u);
    freeMatrix(s);
    freeMatrix(v);


	return EXIT_SUCCESS;
}

[Mr Andersson]

Din profiler visar fel. Använd valgrind för minnesanalys.
==6189== HEAP SUMMARY:
==6189== in use at exit: 0 bytes in 0 blocks
==6189== total heap usage: 260,355 allocs, 260,355 frees, 25,002,896 bytes allocated
==6189==
==6189== All heap blocks were freed -- no leaks are possible

[Al_Bundy]

Inga minnesläckor?

Hur blir det med den där stora C filen då i själva. ZIP filen ?

[Mr Andersson]

> Använd valgrind för minnesanalys.

Eller förväntar du dig att alla andra ska göra allt jobb och du bara sätter ditt namn på resultatet?

[Al_Bundy]

Har jag sagt det eller ?