Svenska ElektronikForumet

Svenskt forum för elektroniksnack.
https://elektronikforumet.com/forum/

Fungerar va_arg med float* ?

https://elektronikforumet.com/forum/viewtopic.php?t=112045

Sida 2 av 3

Re: Fungerar va_arg med float* ?

Postat: 10 februari 2025, 12:52:37
av DanielM
Det börjar ta sig.
Jag vet att det är VLA som jag använder nu.

row_a = 48. column_a = 10.

Detta är den enda funktionen som jag har kvar nu. Sedan är problemet löst.

Kod: Markera allt

static bool opti(float Q[], float c[], float A[], float b[], float x[], size_t row_a, size_t column_a){		
	/* Declare */
	size_t i, j, k;
	
	/* Use Cholesky factorization to solve x from Qx = c because Q is square and symmetric */
	linsolve_chol(Q, x, c, column_a);

	/* Save address */
	float* Ai = A;

	/* Turn x negative */
	for (i = 0; i < column_a; i++) {
		x[i] = -x[i];

		/* If x was nan - return false - Use higher value on lambda constant! */
#ifdef _MSC_VER
		if (_isnanf(x[i])) {
			return false;
		}
#else
		if (isnanf(x[i])) {
			return false;
		}
#endif
	}

	/* Count how many constraints A*x > b */
	//float *K = (float*)malloc(row_a * sizeof(float));
	float K[row_a];
	size_t violations = 0;
	float value;
	for(i = 0; i < row_a; i++){

		/* Check how many rows are A*x > b */
		value = 0.0f;
		for(j = 0; j < column_a; j++){
			value += Ai[j] * x[j];
			/* value += A[i*column_a + j]*x[j]; */
		}
		Ai += column_a;

		/* Constraints difference */
		K[i] = b[i] - value;

		/* Check constraint violation */
		if (K[i] < 0.0f) {
			violations++;
		}
	}

	/* No violation */
	if (violations == 0) {
		//free(K);
		return true;
	}

	/* Solve QP = A' (Notice that we are using a little trick here so we can avoid A') */
	//float *P = (float*)malloc(row_a * column_a * sizeof(float));
	float p[row_a * column_a];
	float* P = p;
	float *P0 = P;
	Ai = A;
	for (i = 0; i < row_a; i++) {
		linsolve_chol(Q, P, Ai, column_a);
		/* linsolve_gauss(Q, P, Ai, column_a, column_a, 0.0f); - Old */
		P += column_a;
		Ai += column_a;
	}
	P = P0;
	tran(P, row_a, column_a);

	/* Multiply H = A*Q*A' */
	//float *H = (float*)malloc(row_a * row_a * sizeof(float));
	float h[row_a * row_a];
	float* H = h;
	float *Hj;
	mul(A, P, H, row_a, column_a, row_a);

	/* Solve lambda from H*lambda = -K, where lambda >= 0 */
	float *lambda = (float*)malloc(row_a * sizeof(float));
	memset(lambda, 0, row_a * sizeof(float));
	float *lambda_p = (float*)malloc(row_a * sizeof(float));
	float w;
	for(i = 0; i < MAX_ITERATIONS; i++){
		/* Update */
		memcpy(lambda_p, lambda, row_a * sizeof(float));

		/* Use Gauss Seidel */
		Hj = H;
		for (j = 0; j < row_a; j++) {
			/* w = H(i, :)*lambda */
			w = 0.0f;
			for (k = 0; k < row_a; k++) {
				w += Hj[k] * lambda[k];
			}
			
			/* Find a solution */
			w = -1.0f / Hj[j] * (K[j] + w - Hj[j]*lambda[j]);
			Hj += row_a;
			lambda[j] = vmax(0.0f, w);
		}

		/* Check if we should break - Found the optimal solution */
		w = 0.0f;
		for (j = 0; j < row_a; j++) {
			value = lambda[j] - lambda_p[j];
			w += value * value;
		}
#ifdef _MSC_VER
		if (w < MIN_VALUE || _isnanf(w)) {
			break;
		}
#else
		if (w < MIN_VALUE || isnanf(w)) {
			break;
		}
#endif
	}

	/* Solve x = x + P*lambda (Notice that x is negative (see above)) */
	float *Plambda = (float*)malloc(column_a * sizeof(float));
	mul(P, lambda, Plambda, column_a, row_a, 1);
	for (j = 0; j < column_a; j++) {
		x[j] -= Plambda[j];
	}

	/* Free */
	//free(K);
	//free(P);
	//free(H);
	free(lambda);
	free(lambda_p);
	free(Plambda);

	/* If w was nan - return false */
#ifdef _MSC_VER
	if (_isnanf(w)) {
		return false;
	}
#else
	if (isnanf(w)) {
		return false;
	}
#endif

	/* If i equal to MAX_ITERATIONS, then it did not find a solution */
	return i < MAX_ITERATIONS;
}

Re: Fungerar va_arg med float* ?

Postat: 10 februari 2025, 20:10:08
av DanielM
Är det någon som har varit med om att en const size_t har ändrat sig efter en beräkning?
När jag kör denna kod.

Kod: Markera allt

#include <string.h>
#include <float.h>
#include <math.h>
#include <stdbool.h>
#define MIN_VALUE 1e-11f

bool chol(const float A[], float L[], const size_t row) {
	/* Save address */
	float* Li = L;
	float* Lj;
	const float* Ai = A;

	memset(L, 0, row * row * sizeof(float));
	float s;
	size_t i, j, k;
	for (i = 0; i < row; i++) {
		Lj = L;
		for (j = 0; j <= i; j++) {
			s = 0.0f;
			for (k = 0; k < j; k++) {
				s += Li[k] * Lj[k]; /* s += L[row * i + k] * L[row * j + k]; */
			}

			/* We cannot divide with zero L[row * j + j] */
			if (fabsf(Lj[j]) < MIN_VALUE) {
				Lj[j] = FLT_EPSILON;
				/* L[row * j + j] = FLT_EPSILON; // Same as eps command in MATLAB */
			}
			Li[j] = (i == j) ? sqrtf(Ai[i] - s) : (1.0f / Lj[j] * (Ai[j] - s));
			/* L[row * i + j] = (i == j) ? sqrtf(A[row * i + i] - s) : (1.0f / L[row * j + j] * (A[row * i + j] - s)); */
			Lj += row;
		}
		Li += row;
		Ai += row;
	}
	
	return true;
}
Direkt efter memset(L, 0, row * row * sizeof(float)); så ändras row == 10 till row == 0.

Vad kallar vi detta? Minnesläcka? Arrayen L är dock inte NULL :humm:

Re: Fungerar va_arg med float* ?

Postat: 10 februari 2025, 20:49:19
av mankan
Hur är L allokerad? Innehåller den row*row floats?

Jag gissar du skriver sönder stacken för dig själv.

Re: Fungerar va_arg med float* ?

Postat: 10 februari 2025, 20:59:04
av DanielM
Så här är den allokerad.
row är const size_t row = 10

Kod: Markera allt

	float* L = (float*)malloc(row * row * sizeof(float));
	/* Test */
	if(L == NULL){
		x[0] = 1;
	}
	const bool status = chol(A, L, row);
Edit:
Allt verkar fungera bra när jag kör statisk allokering, dvs float L[10*10].
Så den där NULL-kollen verkar lite lurig.
Är det någon som kan ge en fingervisning på vad man ska ändra på för att den ska kunna klara utav att allokera 10*10*4 = 400 bytes på heap?

Eventuellt måste jag hitta en snålare algoritm som kan lösa cholesky faktorisering för linjära system A*x = b där A är symmetrisk och positiv definitiv.

Re: Fungerar va_arg med float* ?

Postat: 10 februari 2025, 23:20:15
av hawkan
Vad du ska ändra på är väl att ta en processor med mycket mera minne. Ram.

Re: Fungerar va_arg med float* ?

Postat: 10 februari 2025, 23:26:06
av DanielM
Jag vet, men jag är så nära nu målet. :) Håller på skriva om en kvadratisk programmeringsalgoritm så den inte använder/använder lite dynamisk allokering.
Det är där problemet sitter i.

Vi kan börja med Matlab koden.

Kod: Markera allt

    /* GNU Octave code:
     *
	% This is quadratic programming with Hildreth's method
	% Min 1/2x^TQx + c^Tx
	% S.t Ax <= b
	%
	% If you want to do maximization, then turn Q and c negative. The constraints are the same
	%
	% Max 1/2x^T(-Q)x + (-c)^Tx
	% S.t Ax <= b
	%
	% Input: Q(Symmetric matrix), c(Objective function), A(Constraint matrix), b(Constraint vector)
	% Output: x(Solution vector), solution(boolean flag)
	% Example 1: [x, solution] = quadprog(Q, c, A, b)

	function [x, solution] = quadprog(Q, c, A, b)
	  % Assume that the solution is true
	  solution = true;

	  % Same as in C code for Functions.h at CControl
	  MIN_VALUE = 1e-11;
	  MAX_ITERATIONS = 10000;

	  % Unconstrained solution
	  x = -linsolve(Q, c);

	  % Constraints difference
	  K = b - A*x;

	  % Check constraint violation
	  if(sum(K < 0) == 0)
		return; % No violation
	  end

	  % Create P
	  P = linsolve(Q, A');

	  % Create H = A*Q*A'
	  H = A*P;

	  % Solve lambda from H*lambda = -K, where lambda >= 0
	  [m, n] = size(K);
	  lambda = zeros(m, n);
	  for km = 1:MAX_ITERATIONS
		lambda_p = lambda;

		% Use Gauss Seidel
		for i = 1:m
		  w = -1.0/H(i,i)*(K(i) + H(i,:)*lambda - H(i,i)*lambda(i));
		  lambda(i) = max(0, w);
		end

		% Check if the minimum convergence has been reached
		w = (lambda - lambda_p)'*(lambda - lambda_p);
		if (w < MIN_VALUE)
		  break;
		end

		% Check if the maximum iteration have been reached
		if(km == MAX_ITERATIONS)
		  solution = false;
		  return;
		end
	  end

	  % Find the solution: x = -inv(Q)*c - inv(Q)*A'*lambda
	  x = x - P*lambda;
	end
 */
Denna kod kan enkelt översättas till detta

Kod: Markera allt

static bool opti(const float Q[], const float c[], const float A[], const float b[], float x[], const size_t row_a, const size_t column_a){
	/* Declare */
	size_t i, j, k;
	
	/* Use Cholesky factorization to solve x from Qx = c because Q is square and symmetric */
	linsolve_chol(Q, x, c, column_a);

	/* Save address */
	const float* Ai = A;

	/* Turn x negative */
	for (i = 0; i < column_a; i++) {
		x[i] = -x[i];

		/* If x was nan - return false - Use higher value on lambda constant! */
#ifdef _MSC_VER
		if (_isnanf(x[i])) {
			return false;
		}
#else
		if (isnanf(x[i])) {
			return false;
		}
#endif
	}

	/* Count how many constraints A*x > b */
	float *K = (float*)malloc(row_a * sizeof(float));
	size_t violations = 0;
	float value;
	for(i = 0; i < row_a; i++){

		/* Check how many rows are A*x > b */
		value = 0.0f;
		for(j = 0; j < column_a; j++){
			value += Ai[j] * x[j];
			/* value += A[i*column_a + j]*x[j]; */
		}
		Ai += column_a;

		/* Constraints difference */
		K[i] = b[i] - value;

		/* Check constraint violation */
		if (K[i] < 0.0f) {
			violations++;
		}
	}

	/* No violation */
	if (violations == 0) {
		free(K);
		return true;
	}

	/* Solve QP = A' (Notice that we are using a little trick here so we can avoid A') */
	float *P = (float*)malloc(row_a * column_a * sizeof(float));
	float *P0 = P;
	Ai = A;
	for (i = 0; i < row_a; i++) {
		linsolve_chol(Q, P, Ai, column_a);
		/* linsolve_gauss(Q, P, Ai, column_a, column_a, 0.0f); - Old */
		P += column_a;
		Ai += column_a;
	}
	P = P0;
	tran(P, row_a, column_a);

	/* Multiply H = A*Q*A' */
	float *H = (float*)malloc(row_a * row_a * sizeof(float));
	float *Hj;
	mul(A, P, H, row_a, column_a, row_a);

	/* Solve lambda from H*lambda = -K, where lambda >= 0 */
	float *lambda = (float*)malloc(row_a * sizeof(float));
	memset(lambda, 0, row_a * sizeof(float));
	float *lambda_p = (float*)malloc(row_a * sizeof(float));
	float w;
	for(i = 0; i < MAX_ITERATIONS; i++){
		/* Update */
		memcpy(lambda_p, lambda, row_a * sizeof(float));

		/* Use Gauss Seidel */
		Hj = H;
		for (j = 0; j < row_a; j++) {
			/* w = H(i, :)*lambda */
			w = 0.0f;
			for (k = 0; k < row_a; k++) {
				w += Hj[k] * lambda[k];
			}
			
			/* Find a solution */
			w = -1.0f / Hj[j] * (K[j] + w - Hj[j]*lambda[j]);
			Hj += row_a;
			lambda[j] = vmax(0.0f, w);
		}

		/* Check if we should break - Found the optimal solution */
		w = 0.0f;
		for (j = 0; j < row_a; j++) {
			value = lambda[j] - lambda_p[j];
			w += value * value;
		}
#ifdef _MSC_VER
		if (w < MIN_VALUE || _isnanf(w)) {
			break;
		}
#else
		if (w < MIN_VALUE || isnanf(w)) {
			break;
		}
#endif
	}

	/* Solve x = x + P*lambda (Notice that x is negative (see above)) */
	float *Plambda = (float*)malloc(column_a * sizeof(float));
	mul(P, lambda, Plambda, column_a, row_a, 1);
	for (j = 0; j < column_a; j++) {
		x[j] -= Plambda[j];
	}

	/* Free */
	free(K);
	free(P);
	free(H);
	free(lambda);
	free(lambda_p);
	free(Plambda);

	/* If w was nan - return false */
#ifdef _MSC_VER
	if (_isnanf(w)) {
		return false;
	}
#else
	if (isnanf(w)) {
		return false;
	}
#endif

	/* If i equal to MAX_ITERATIONS, then it did not find a solution */
	return i < MAX_ITERATIONS;
}
Men för att minimera dynamisk allokering så kan man använda flera for-satser. Som ni ser så måste det fungera säkert väldigt bra att bara ta en rad i taget.

Kod: Markera allt

static bool optislim(const float Q[], const float c[], const float A[], const float b[], float x[], const size_t row_a, const size_t column_a){
	/* Declare */
	size_t i, j, k, l;
	
	/* Use Cholesky factorization to solve x from Qx = c because Q is square and symmetric */
	linsolve_chol(Q, x, c, column_a);

	/* Save address */
	const float* Ai = A;

	/* Turn x negative */
	for (i = 0; i < column_a; i++) {
		x[i] = -x[i];

		/* If x was nan - return false - Use higher value on lambda constant! */
#ifdef _MSC_VER
		if (_isnanf(x[i])) {
			return false;
		}
#else
		if (isnanf(x[i])) {
			return false;
		}
#endif
	}

	/* Count how many constraints A*x > b */
	float K, value;
	for (i = 0; i < row_a; i++) {

		/* Check how many rows are A*x > b */
		value = 0.0f;
		for(j = 0; j < column_a; j++){
			value += Ai[j] * x[j];
			/* value += A[i*column_a + j]*x[j]; */
		}

		/* Constraints difference */
		K = b[i] - value;

		/* Check constraint violation */
		if (K < 0.0f) {
			/* Solve QP = A' (Notice that we are using a little trick here so we can avoid A') */
			float* P = (float*)malloc(column_a * sizeof(float));
			linsolve_chol(Q, P, Ai, column_a);

			/* Multiply H = A*Q*A' */
			for (k = 0; k < row_a; k++) {
				value = 0.0f;
				for (l = 0; l < column_a; l++) {
					value += A[k * column_a + l] * P[l];
				}

				for (l = 0; l < MAX_ITERATIONS; l++) {

				}

			}

			/* Free */
			free(P);
		}

		Ai += column_a;
	}

	/* No violation */
	return true;

	
	/* Solve lambda from H*lambda = -K, where lambda >= 0 */
	float *lambda = (float*)malloc(row_a * sizeof(float));
	memset(lambda, 0, row_a * sizeof(float));
	float *lambda_p = (float*)malloc(row_a * sizeof(float));
	float w;
	for(i = 0; i < MAX_ITERATIONS; i++){
		/* Update */
		memcpy(lambda_p, lambda, row_a * sizeof(float));

		/* Use Gauss Seidel */
		Hj = H;
		for (j = 0; j < row_a; j++) {
			/* w = H(i, :)*lambda */
			w = 0.0f;
			for (k = 0; k < row_a; k++) {
				w += Hj[k] * lambda[k];
			}
			
			/* Find a solution */
			w = -1.0f / Hj[j] * (K[j] + w - Hj[j]*lambda[j]);
			Hj += row_a;
			lambda[j] = vmax(0.0f, w);
		}

		/* Check if we should break - Found the optimal solution */
		w = 0.0f;
		for (j = 0; j < row_a; j++) {
			value = lambda[j] - lambda_p[j];
			w += value * value;
		}
#ifdef _MSC_VER
		if (w < MIN_VALUE || _isnanf(w)) {
			break;
		}
#else
		if (w < MIN_VALUE || isnanf(w)) {
			break;
		}
#endif
	}

	/* Solve x = x + P*lambda (Notice that x is negative (see above)) */
	float *Plambda = (float*)malloc(column_a * sizeof(float));
	mul(P, lambda, Plambda, column_a, row_a, 1);
	for (j = 0; j < column_a; j++) {
		x[j] -= Plambda[j];
	}

	/* Free */
	free(K);
	free(P);
	free(H);
	free(lambda);
	free(lambda_p);
	free(Plambda);

	/* If w was nan - return false */
#ifdef _MSC_VER
	if (_isnanf(w)) {
		return false;
	}
#else
	if (isnanf(w)) {
		return false;
	}
#endif

	/* If i equal to MAX_ITERATIONS, then it did not find a solution */
	return i < MAX_ITERATIONS;
}

Re: Fungerar va_arg med float* ?

Postat: 11 februari 2025, 06:23:26
av hawkan
Nja, alltså finns det tillräckligt med minne kommer malloc att fungera även om man inte borde använda det.
Det avdelas en mängd minne till heap som malloc använder, men både stack och heap tas från ram, stack underifrån, heap från max ner t ex. På något sätt. Ta reda på hur minne avdelas för heap och hur du kan mäta och påverka det. Och se om det alls går att få plats med alla variabler.

Re: Fungerar va_arg med float* ?

Postat: 11 februari 2025, 08:07:00
av DanielM
Jag tänker samma sak. Om det går att allokera statiskt, så måste det gå att allokera dynamiskt.

Re: Fungerar va_arg med float* ?

Postat: 11 februari 2025, 08:09:22
av hawkan
I princip ja, men minnet är uppdelat, malloc tar inte från stacken utan från heap. Kan vara att malloc startar från max minne och går nedåt.
Så du får ha koll på hur uppdelningen görs, om du kan ändra den om det behövs så du får mer minne till heap (malloc).

Re: Fungerar va_arg med float* ?

Postat: 11 februari 2025, 08:59:48
av hummel
DanielM, du får titta i minneskartan och verifiera exakt vad och var ditt minne är allokerat och används. Hur mycket ledigt har du på stack respektive heap under exekvering?

Re: Fungerar va_arg med float* ?

Postat: 11 februari 2025, 18:29:31
av DanielM
Är det någon funktion i IDE:n eller?
Jag använder mig av STM32CubeIDE (Eclipse) och jag har inte sett något där man kan i realtid kolla hur mycket minne förbrukas för varje rad av kod.

Edit:
Jag har kommit vidare nu. Jag återanvände arrayer och minskade storlekar på dom. Nu så kan jag köra chol utan problem. Men ännu är jag inte i mål!

Jag har kommit till "sista bossen" nu kan man väll säga. En array på 48*48 som har 32-bit float. :twisted:

Re: Fungerar va_arg med float* ?

Postat: 12 februari 2025, 09:14:02
av hummel
Vissa utvecklingssystem har det andra inte, kanske finns det tillägg till ditt? Kanske debuggern har den funktionen.

Re: Fungerar va_arg med float* ?

Postat: 13 februari 2025, 10:44:50
av DanielM
Cortex M0 har dock inte detta tyvärr.
Men som tur så har jag rundat av det svåra nu.

Nu allokerar jag som mest (48*10 + 48 + 100)*4 = 2512 bytes under körningen av quadprog.

Jag tror det inte går att optimera mera. Ursprungligen var det 12152 bytes som quadprog tog i bytes.

Ska presentera koden.

P.S Är det någon här som vill lära sig lite om MPC och motiveringen till användningen av MPC?

Re: Fungerar va_arg med float* ?

Postat: 13 februari 2025, 10:45:55
av hawkan
Och processorn har hur mycket ram?

Re: Fungerar va_arg med float* ?

Postat: 13 februari 2025, 12:21:58
av DanielM
8 kB :vissla:
STM32F051 heter processorn.
Alla tidsangivelser är UTC+02:00 Europe/Stockholm
Sida 2 av 3

Drivs av phpBB® Forum Software © phpBB Limited

Swedish translation by phpBB Sweden © 2006-2020