Matrisberäkningar med för STM32?

PIC, AVR, Arduino, Raspberry Pi, Basic Stamp, PLC mm.
Användarvisningsbild
Al_Bundy
Inlägg: 2889
Blev medlem: 11 september 2012, 23:59:50
Ort: The U.S - Chicago
Kontakt:

Re: Matrisberäkningar med för STM32?

Inlägg av Al_Bundy »

Någon som vet varför denna kod inte gör rätt lösningen till ekvationen Ax=b?

https://software.intel.com/sites/produc ... v_ex.c.htm

Kod: Markera allt


Intel® Math Kernel Library LAPACK Examples
/*******************************************************************************
*  Copyright (C) 2009-2015 Intel Corporation. All Rights Reserved.
*  The information and material ("Material") provided below is owned by Intel
*  Corporation or its suppliers or licensors, and title to such Material remains
*  with Intel Corporation or its suppliers or licensors. The Material contains
*  proprietary information of Intel or its suppliers and licensors. The Material
*  is protected by worldwide copyright laws and treaty provisions. No part of
*  the Material may be copied, reproduced, published, uploaded, posted,
*  transmitted, or distributed in any way without Intel's prior express written
*  permission. No license under any patent, copyright or other intellectual
*  property rights in the Material is granted to or conferred upon you, either
*  expressly, by implication, inducement, estoppel or otherwise. Any license
*  under such intellectual property rights must be express and approved by Intel
*  in writing.
*
********************************************************************************
*/
/*
   DGESV Example.
   ==============
 
   The program computes the solution to the system of linear
   equations with a square matrix A and multiple
   right-hand sides B, where A is the coefficient matrix:
 
     6.80  -6.05  -0.45   8.32  -9.67
    -2.11  -3.30   2.58   2.71  -5.14
     5.66   5.36  -2.70   4.35  -7.26
     5.97  -4.44   0.27  -7.17   6.08
     8.23   1.08   9.04   2.14  -6.87

   and B is the right-hand side matrix:
 
     4.02  -1.56   9.81
     6.19   4.00  -4.09
    -8.22  -8.67  -4.57
    -7.57   1.75  -8.61
    -3.03   2.86   8.99
 
   Description.
   ============
 
   The routine solves for X the system of linear equations A*X = B,
   where A is an n-by-n matrix, the columns of matrix B are individual
   right-hand sides, and the columns of X are the corresponding
   solutions.

   The LU decomposition with partial pivoting and row interchanges is
   used to factor A as A = P*L*U, where P is a permutation matrix, L
   is unit lower triangular, and U is upper triangular. The factored
   form of A is then used to solve the system of equations A*X = B.

   Example Program Results.
   ========================
 
 DGESV Example Program Results

 Solution
  -0.80  -0.39   0.96
  -0.70  -0.55   0.22
   0.59   0.84   1.90
   1.32  -0.10   5.36
   0.57   0.11   4.04

 Details of LU factorization
   8.23   1.08   9.04   2.14  -6.87
   0.83  -6.94  -7.92   6.55  -3.99
   0.69  -0.67 -14.18   7.24  -5.19
   0.73   0.75   0.02 -13.82  14.19
  -0.26   0.44  -0.59  -0.34  -3.43

 Pivot indices
      5      5      3      4      5
*/
#include <stdlib.h>
#include <stdio.h>

/* DGESV prototype */
extern void dgesv( int* n, int* nrhs, double* a, int* lda, int* ipiv,
                double* b, int* ldb, int* info );
/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, int m, int n, double* a, int lda );
extern void print_int_vector( char* desc, int n, int* a );

/* Parameters */
#define N 5
#define NRHS 3
#define LDA N
#define LDB N

/* Main program */
int main() {
        /* Locals */
        int n = N, nrhs = NRHS, lda = LDA, ldb = LDB, info;
        /* Local arrays */
        int ipiv[N];
        double a[LDA*N] = {
            6.80, -2.11,  5.66,  5.97,  8.23,
           -6.05, -3.30,  5.36, -4.44,  1.08,
           -0.45,  2.58, -2.70,  0.27,  9.04,
            8.32,  2.71,  4.35, -7.17,  2.14,
           -9.67, -5.14, -7.26,  6.08, -6.87
        };
        double b[LDB*NRHS] = {
            4.02,  6.19, -8.22, -7.57, -3.03,
           -1.56,  4.00, -8.67,  1.75,  2.86,
            9.81, -4.09, -4.57, -8.61,  8.99
        };
        /* Executable statements */
        printf( " DGESV Example Program Results\n" );
        /* Solve the equations A*X = B */
        dgesv( &n, &nrhs, a, &lda, ipiv, b, &ldb, &info );
        /* Check for the exact singularity */
        if( info > 0 ) {
                printf( "The diagonal element of the triangular factor of A,\n" );
                printf( "U(%i,%i) is zero, so that A is singular;\n", info, info );
                printf( "the solution could not be computed.\n" );
                exit( 1 );
        }
        /* Print solution */
        print_matrix( "Solution", n, nrhs, b, ldb );
        /* Print details of LU factorization */
        print_matrix( "Details of LU factorization", n, n, a, lda );
        /* Print pivot indices */
        print_int_vector( "Pivot indices", n, ipiv );
        exit( 0 );
} /* End of DGESV Example */

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, double* a, int lda ) {
        int i, j;
        printf( "\n %s\n", desc );
        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.2f", a[i+j*lda] );
                printf( "\n" );
        }
}

/* Auxiliary routine: printing a vector of integers */
void print_int_vector( char* desc, int n, int* a ) {
        int j;
        printf( "\n %s\n", desc );
        for( j = 0; j < n; j++ ) printf( " %6i", a[j] );
        printf( "\n" );
}
hawkan
Inlägg: 3445
Blev medlem: 14 augusti 2011, 10:27:40

Re: Matrisberäkningar med för STM32?

Inlägg av hawkan »

Vad är det som blir fel?
Är exemplet i kommentaren fel?
Eller får du inte det resultatet när du kör exemplet?
Användarvisningsbild
Al_Bundy
Inlägg: 2889
Blev medlem: 11 september 2012, 23:59:50
Ort: The U.S - Chicago
Kontakt:

Re: Matrisberäkningar med för STM32?

Inlägg av Al_Bundy »

För det första så har vi två matriser, a och b.

Kod: Markera allt

>> A = [ 6.80, -2.11,  5.66,  5.97,  8.23,
           -6.05, -3.30,  5.36, -4.44,  1.08,
           -0.45,  2.58, -2.70,  0.27,  9.04,
            8.32,  2.71,  4.35, -7.17,  2.14,
           -9.67, -5.14, -7.26,  6.08, -6.87];
>> b = [ 4.02,  6.19, -8.22, -7.57, -3.03,
           -1.56,  4.00, -8.67,  1.75,  2.86,
            9.81, -4.09, -4.57, -8.61,  8.99];
>> x = linsolve(A,b)
error: linsolve: operator \: nonconformant arguments (op1 is 5x5, op2 is 3x5)
error: called from
    linsolve at line 107 column 7
>>
Koden fungerar inte ens.
hawkan skrev:Vad är det som blir fel?
Är exemplet i kommentaren fel?
Eller får du inte det resultatet när du kör exemplet?

Felet är att det blir fel beräkning. Intel verkar ha gjort fel.
hawkan
Inlägg: 3445
Blev medlem: 14 augusti 2011, 10:27:40

Re: Matrisberäkningar med för STM32?

Inlägg av hawkan »

Här funkar det

Kod: Markera allt

>> A
A =

   6.80000  -6.05000  -0.45000   8.32000  -9.67000
  -2.11000  -3.30000   2.58000   2.71000  -5.14000
   5.66000   5.36000  -2.70000   4.35000  -7.26000
   5.97000  -4.44000   0.27000  -7.17000   6.08000
   8.23000   1.08000   9.04000   2.14000  -6.87000

>> B
B =

   4.0200  -1.5600   9.8100
   6.1900   4.0000  -4.0900
  -8.2200  -8.6700  -4.5700
  -7.5700   1.7500  -8.6100
   3.0300   2.8600   8.9900

>> linsolve(A,B)
ans =

  -0.60025  -0.38962   0.95546
  -0.51329  -0.55443   0.22066
   1.17754   0.84223   1.90064
   1.60348  -0.10380   5.35766
   0.80815   0.10571   4.04060
Du transponerar A och B, ska det vara så?
hawkan
Inlägg: 3445
Blev medlem: 14 augusti 2011, 10:27:40

Re: Matrisberäkningar med för STM32?

Inlägg av hawkan »

Nä det blev visst inte rätt. Tja vad kan det bero på?
Användarvisningsbild
Al_Bundy
Inlägg: 2889
Blev medlem: 11 september 2012, 23:59:50
Ort: The U.S - Chicago
Kontakt:

Re: Matrisberäkningar med för STM32?

Inlägg av Al_Bundy »

Ja. då typ ser det bättre ut.

Dock blir det inte exakt samma resultat, men väldigt nära!

Kod: Markera allt

>> b'
ans =

   4.0200  -1.5600   9.8100
   6.1900   4.0000  -4.0900
  -8.2200  -8.6700  -4.5700
  -7.5700   1.7500  -8.6100
  -3.0300   2.8600   8.9900

>> x = linsolve(A,b')
x =

  -1.46576   0.26841   0.23346
   2.71274  -1.04172  -0.62699
   2.89937   0.22400   0.25025
   1.85694  -0.36079   1.30327
  -0.94597  -0.57073  -0.27915

>> A
A =

   6.80000  -2.11000   5.66000   5.97000   8.23000
  -6.05000  -3.30000   5.36000  -4.44000   1.08000
  -0.45000   2.58000  -2.70000   0.27000   9.04000
   8.32000   2.71000   4.35000  -7.17000   2.14000
  -9.67000  -5.14000  -7.26000   6.08000  -6.87000

>> x = linsolve(A',b')
x =

  -0.80071  -0.38962   0.95546
  -0.69524  -0.55443   0.22066
   0.59391   0.84223   1.90064
   1.32173  -0.10380   5.35766
   0.56576   0.10571   4.04060

>> Losning = [-0.80  -0.39   0.96
  -0.70  -0.55   0.22
   0.59   0.84   1.90
   1.32  -0.10   5.36
   0.57   0.11   4.04]
Losning =

  -0.80000  -0.39000   0.96000
  -0.70000  -0.55000   0.22000
   0.59000   0.84000   1.90000
   1.32000  -0.10000   5.36000
   0.57000   0.11000   4.04000

>> abs(x-Losning)
ans =

   7.1403e-04   3.7861e-04   4.5351e-03
   4.7566e-03   4.4271e-03   6.5963e-04
   3.9150e-03   2.2274e-03   6.3673e-04
   1.7256e-03   3.8019e-03   2.3385e-03
   4.2438e-03   4.2890e-03   6.0266e-04

>>
Men vem säger att man måste ha transponatet av alla matriser? Vad är det för logik? Kan det vara så att Fortran läser matriser radvis, medan C läser matriser kolumnvis?
hawkan
Inlägg: 3445
Blev medlem: 14 augusti 2011, 10:27:40

Re: Matrisberäkningar med för STM32?

Inlägg av hawkan »

B(1,5) har fel tecken. Det ska vara -3.03, då blir det rätt.
Användarvisningsbild
Al_Bundy
Inlägg: 2889
Blev medlem: 11 september 2012, 23:59:50
Ort: The U.S - Chicago
Kontakt:

Re: Matrisberäkningar med för STM32?

Inlägg av Al_Bundy »

Så här blir det för mig.

Kod: Markera allt

 DGESV Example Program Results

 Solution
  -0.80  -0.39   0.96
  -0.70  -0.55   0.22
   0.59   0.84   1.90
   1.32  -0.10   5.36
   0.57   0.11   4.04

 Details of LU factorization
   8.23   1.08   9.04   2.14  -6.87
   0.83  -6.94  -7.92   6.55  -3.99
   0.69  -0.67 -14.18   7.24  -5.19
   0.73   0.75   0.02 -13.82  14.19
  -0.26   0.44  -0.59  -0.34  -3.43

 Pivot indices
      5      5      3      4      5
Mr Andersson
Inlägg: 1409
Blev medlem: 29 januari 2011, 21:06:30
Ort: Lapplandet

Re: Matrisberäkningar med för STM32?

Inlägg av Mr Andersson »

Al_Bundy skrev:Kan det vara så att Fortran läser matriser radvis, medan C läser matriser kolumnvis?
Fortran (och octave) lagrar matriser i column-major-format medans C och dess derivat använder row-major. Det har inget att göra med i vilken ordning datan läses, det är helt upp till programmeraren. Däremot så är det effektivare att läsa kolumnvis i fortran och radvis i C.
Men, detta gäller endast för flerdimensionella vektorer. Intels kod är endimensionell och då är det helt upp till biblioteksdesignen hur datan ska tolkas. (Men radvis är effektivare oavsett språk för endimensionella vektorer så ser du något annat är det konstigt)
Bild
(bildkälla wikipedia)
Notera att detta bara bestämmer hur tvådimensionell data representeras i endimensionellt minne.

Edit: Första delen av den här bloggen tar upp lite av skillnaderna mellan tvådimensionella vektorer i C och fortran.
Användarvisningsbild
Al_Bundy
Inlägg: 2889
Blev medlem: 11 september 2012, 23:59:50
Ort: The U.S - Chicago
Kontakt:

Re: Matrisberäkningar med för STM32?

Inlägg av Al_Bundy »

Misstänkte detta. Så vad jag kan göra för att lösa min kod?

Ska jag ta transponatet av matris A och B?

Kod: Markera allt

/*
 ============================================================================
 Name        : BlasTest.c
 Author      : 
 Version     :
 Copyright   : MIT
 Description : Hello World in C, Ansi-style
 ============================================================================
 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include "blaswrap.h"
#include "f2c.h"
#include "clapack.h"

int main( )
{

	clock_t start, end;
	float cpu_time_used;

	integer N = 5; // Vi har en 5 dimensions matris A
	integer NRHS = 2; // Vi har en 2 kolumns B matris
	integer LDA = 5; // Samma dimension som matris A
	real A[5*5] = {
			3,    4,    2,    5,    5,
		    2,    4,    5,   -5,    5,
		  -13,   -4,    5,    5,    3,
		  -34,   -4 ,  -2,   -5,   -5,
		    2,    4,    5,   -3,    5};

	integer IPIV[5] = {1,1,1,1,1}; // vi sätter denna som en diagonalmatris
	integer LDB = 5; // Raden på matris B
	real B[5*2] = {
			6,    2,
			3,    5,
			1,   -4,
			5,    1,
			5,   44,
	};

	integer INFO; // Inget vi behöver säga nu!


    start = clock();

    sgesv_(&N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO); // Arrayer behöver vi INTE använda &-tecknet på


    end = clock();
    cpu_time_used = ((float) (end - start)) / CLOCKS_PER_SEC;
    printf("\nTotal speed  was %f and INFO = %d\n", cpu_time_used, (int) INFO);

    // Enligt dokumentationen så om info = 0, så betyder det att B = X, dvs våran lösning på AX=B
    for(int i = 0; i < 5; i++){
    	for(int j = 0; j < 2; j++){
    		printf("%f ", (float) B[i*4+j]);
    	}
    	printf("\n");
    }

    return 0;
}

Senast redigerad av Al_Bundy 8 februari 2019, 17:29:53, redigerad totalt 2 gånger.
Mr Andersson
Inlägg: 1409
Blev medlem: 29 januari 2011, 21:06:30
Ort: Lapplandet

Re: Matrisberäkningar med för STM32?

Inlägg av Mr Andersson »

Du måste läsa dokumentationen för att se vilket format de använder.
Användarvisningsbild
Al_Bundy
Inlägg: 2889
Blev medlem: 11 september 2012, 23:59:50
Ort: The U.S - Chicago
Kontakt:

Re: Matrisberäkningar med för STM32?

Inlägg av Al_Bundy »

Med tanke på att denna dokumentation verkar vara dålig så hittar jag den knappt. Men jag ska söka reda på den.
Annars så verkar det som att om jag har matriserna som ett transponat, så fungerar det.

Kod: Markera allt

/*
 ============================================================================
 Name        : BlasTest.c
 Author      : 
 Version     :
 Copyright   : MIT
 Description : Hello World in C, Ansi-style
 ============================================================================
 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include "blaswrap.h"
#include "f2c.h"
#include "clapack.h"

int main( )
{

	clock_t start, end;
	float cpu_time_used;

	integer N = 5; // Vi har en 5 dimensions matris A
	integer NRHS = 2; // Vi har en 2 kolumns B matris
	integer LDA = 5; // Samma dimension som matris A
	real A[5*5] = {
			3,    2,  -13,  -34,    2,
			    4,    4,   -4,   -4,    4,
			    2,    5 ,   5 ,  -2 ,   5,
			    5  , -5  ,  5  , -5 ,  -3,
			    5  ,  5  ,  3  , -5  ,  5};

	integer IPIV[5] = {1,1,1,1,1}; // vi sätter denna som en diagonalmatris
	integer LDB = 5; // Raden på matris B
	real B[5*2] = {
			6,    3,    1,    5,    5,
			    2   , 5  , -4   , 1  , 44,
	};

	integer INFO; // Inget vi behöver säga nu!


    start = clock();

    sgesv_(&N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO); // Arrayer behöver vi INTE använda &-tecknet på


    end = clock();
    cpu_time_used = ((float) (end - start)) / CLOCKS_PER_SEC;
    printf("\nTotal speed  was %f and INFO = %d\n", cpu_time_used, (int) INFO);

    // Enligt dokumentationen så om info = 0, så betyder det att B = X, dvs våran lösning på AX=B
    for(int i = 0; i < 2; i++){
    	for(int j = 0; j < 5; j++){
    		printf("%f ", (float) B[i*5+j]);
    	}
    	printf("\n");
    }

    return 0;
}
Resultat från Lapack

Kod: Markera allt

Total speed  was 0.004028 and INFO = 0
-0.354839 2.854502 2.215054 1.000000 -2.756720 
-0.096775 46.298393 65.967735 19.499998 -82.467735 
Detta ger samma resultat som

Kod: Markera allt

A = [3,    4,    2,    5,    5,
		    2,    4,    5,   -5,    5,
		  -13,   -4,    5,    5,    3,
		  -34,   -4 ,  -2,   -5,   -5,
		    2,    4,    5,   -3,    5];
        
B = [6,    2,
			3,    5,
			1,   -4,
			5,    1,
			5,   44];
      
 X = linsolve(A,B)

>> X = 
   -0.354839   -0.096774
    2.854503   46.298387
    2.215054   65.967742
    1.000000   19.500000
   -2.756720  -82.467742
Tror jag får använda denna funktion istället
http://www.netlib.org/lapack/explore-ht ... 11bd30cdc9

Eller så gör jag en C-funktion som utför transponatet åt mig :)
Användarvisningsbild
Al_Bundy
Inlägg: 2889
Blev medlem: 11 september 2012, 23:59:50
Ort: The U.S - Chicago
Kontakt:

Re: Matrisberäkningar med för STM32?

Inlägg av Al_Bundy »

Jag håller på att lägga in Lapack i mitt C-bibliotek. Det blir mest bara små .c filer som man kan anropa som gör Lapack mer användarvänligt. Dessutom så är Lapack mycket bra när det kommer till just lösa ekvationer som:

\(A = USV^T\) Avklarad
\(A = LU\)
\(A = QR\) Avklarad
\(A = LL^T\) Avklarad
\(Ax = b\) Avklarad
\(Av = \lambda v\)
\(d = |A|\) Avklarad
\(A^{-1}\) Avklarad
\(A^n\)
\(x = A^{\dagger}b\)

Det är ju dessa som har jag haft problem med. Att göra lite transponat, identitetsmatris, klippa och klistra, skriv ut osv är enkla saker.
Då tänkte jag fråga er! Ni har sagt att ni är intresserad utav hjälpa mig. Medan jag håller på skriva C-kod, kan ni kolla om det går att applicera Lapack för STM32? :)
Jag tänkte att vi håller oss till float.

Någon som känner sig manad? :)
bearing
Inlägg: 11672
Blev medlem: 2 mars 2006, 01:01:45
Ort: Ängelholm

Re: Matrisberäkningar med för STM32?

Inlägg av bearing »

Shimonu skrev:Du pratar goja när det kommer till double och float. Att få kompilatorn att tolka double som float, vad betyder det? Antingen kör du double eller så kör du float, det finns inget som kompilatorn ska tolka. Möjligen vill du se om du kan utnyttja FPUn för beräkning med double men annars kommer det användas mjukvarualgoritmer för att göra beräkningar med double.
Det är inte säkert. Jag jobbade med en TI som saknade double-stöd i hårdvaran, men hade stöd för float. Som standard gjorde den om alla double till float vid kompileringen (i TI's egen miljö iaf).
Shimonu
Inlägg: 324
Blev medlem: 21 oktober 2015, 22:44:33

Re: Matrisberäkningar med för STM32?

Inlägg av Shimonu »

Det finns såklart fallet då double == float som du fick, ja. Men det var inte vad Bundy var inne på

.
Skriv svar