Hur kan MLK göra så snabba matrismultiplikationer?

Citera flera Citera

2023-09-02 20:51

Joppis

Medlem ★

Plats: 127.0.0.1
Registrerad: Apr 2008

●

Är inte detta ett problem som kan köras multitrådat? MKL kanske gör just det i bakgrunden medan du jämför med en egen singeltrådad implementation?

Visa signatur

Citera flera Citera

2023-09-02 20:51

aent

Medlem

Registrerad: Maj 2005

●

Jag är inte insatt i den implementeringen (och det var länge sen jag läste matristeori), men i grund och botten kan du göra smarta omskrivningar där problemet översätts till andra operationer som är enklare att beräkna.

Den här länken kanske är intressant för dig: https://coffeebeforearch.github.io/2020/06/23/mmul.html

Senast redigerat 2023-09-02 20:56 Lade till länk om olika implementeringar av matrismultiplikation.

Citera flera Citera (1)

2023-09-02 21:14

Joppis

Medlem ★

Plats: 127.0.0.1
Registrerad: Apr 2008

●

@heretic16 En till grej kan vara att MKL är implementerat så att matriserna delas upp i block så att CPU kan jobba direkt mot cache-minnet så mycket som möjligt.

Visa signatur

Citera flera Citera (1)

2023-09-02 21:28

heretic16

Medlem

Registrerad: Jun 2008

●

Skrivet av aent:

Jag är inte insatt i den implementeringen (och det var länge sen jag läste matristeori), men i grund och botten kan du göra smarta omskrivningar där problemet översätts till andra operationer som är enklare att beräkna.

Den här länken kanske är intressant för dig: https://coffeebeforearch.github.io/2020/06/23/mmul.html

Så man behöver inte extra trådar för att få det gå snabbare?

Citera flera Citera

2023-09-02 23:23

trudelutt

Medlem ★

Registrerad: Aug 2021

●

Skrivet av heretic16:

Att göra en matrismultiplikation med två matriser i C krävs minst 3 for-satser.
Det tar många minuter att multiplicera en 5000x5000 matris med en annan 5000x5000 matris.
Men med MLK går det på några sekunder.

Vad är det som gör MLK så effektiv när det kommer till matrismultiplikationer?
Notera att det är bara matrismultiplikationer som MLK är effektivast på. Övriga så som FFT så är en vanlig enkel FFT algoritm med en for-sats betydligt snabbare.

Gissningsvis en kombination av många grejer. Effektivare algoritm ( https://en.wikipedia.org/wiki/Computational_complexity_of_mat... ), bättre utnyttjande av cache, vektorisering (SSE), och multitrådning.

Citera flera Citera (1)

2023-09-02 23:34

Elgot

Hedersmedlem

Plats: Linköping
Registrerad: Apr 2004

●

Skrivet av heretic16:

Notera att det är bara matrismultiplikationer som MLK är effektivast på. Övriga så som FFT så är en vanlig enkel FFT algoritm med en for-sats betydligt snabbare.

Vad är det för fall du har tittat på här? FFTW har ju traditionellt varit standardvalet när man behöver hög prestanda, men MKL är ofta ännu bättre (och dessutom kompatibelt).

Citera flera Citera

2023-09-03 09:46

heretic16

Medlem

Registrerad: Jun 2008

●

Skrivet av Elgot:

Vad är det för fall du har tittat på här? FFTW har ju traditionellt varit standardvalet när man behöver hög prestanda, men MKL är ofta ännu bättre (och dessutom kompatibelt).

Körde MLK jämfört med FFTPack....stor skillnad. Det var nog 50% snabbare, men det handlade om dryga 0.01 sekunder.

FFTPack är ju riktig K&R-kod och uppfanns 1982.

Här är ett exempel där x är en 1D vektor av float och n-2 är storleken av 1D vektorn, Typ om vektorn är 12, så ska n vara 10.

	/* Create handle */
	DFTI_DESCRIPTOR_HANDLE handle = NULL;

	/* Prepare handle */
	DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, n);

	/* Commit */
	DftiCommitDescriptor(handle);

	/* Compute forward */
	DftiComputeForward(handle, x);

	/* Delete model*/
	DftiFreeDescriptor(&handle);

Citera flera Citera

2023-09-03 10:01

Elgot

Hedersmedlem

Plats: Linköping
Registrerad: Apr 2004

●

Skrivet av heretic16:

Körde MLK jämfört med FFTPack....stor skillnad. Det var nog 50% snabbare, men det handlade om dryga 0.01 sekunder.

FFTPack är ju riktig K&R-kod och uppfanns 1982.

Här är ett exempel där x är en 1D vektor av float och n-2 är storleken av 1D vektorn, Typ om vektorn är 12, så ska n vara 10.

	/* Create handle */
	DFTI_DESCRIPTOR_HANDLE handle = NULL;

	/* Prepare handle */
	DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, n);

	/* Commit */
	DftiCommitDescriptor(handle);

	/* Compute forward */
	DftiComputeForward(handle, x);

	/* Delete model*/
	DftiFreeDescriptor(&handle);

Återanvänder du deskriptorerna? Det tar ju lite tid att initiera och så, men det är ju en engångskostnad och själva beräkningen borde vara snabbare.

Citera flera Citera

2023-09-03 10:36

heretic16

Medlem

Registrerad: Jun 2008

●

Skrivet av Elgot:

Återanvänder du deskriptorerna? Det tar ju lite tid att initiera och så, men det är ju en engångskostnad och själva beräkningen borde vara snabbare.

Ja, det gör jag. MLK är känd på internet för att inte vara den snabbaste. FFTW3 är säkert den snabbaste, men FFTPack, trots sin ålder, så är den en av dom snabbaste, och mest lättanvändarvänligaste också.

Enda skillnaden mellan MLK och FFTPack var att FFTPack kunde ge talen 3.99999999999 medan MLK gav 4.0 exakt.

Nu sitter jag med ett problem att jag måste göra en beräkning med 2D med FFTpack. Bara FFTpack 5.1 har stöd för 2D FFT, men denna är skrivet i Fortran 90, och inte ANSI C som jag skriver min kod i.

Då har jag en plan!

1. Om jag har en matris X som är M x N och sedan tar jag transponatet på X för att göra om den till column-major.
2. Jag använder FFT på varje "rad", som är egentligen varje kolumn med tanke på transponantet.

#define row 8
#define column 3

float X[row * column] = { -0.678816,   1.363392, - 0.020819,
							   0.014361,   0.115923,   0.309136,
							  - 0.195106, - 1.376395, - 0.336756,
							   1.678252, - 0.234396, - 0.023950,
							   1.776940,   1.062231, - 0.580292,
							  - 0.723726, - 1.209755,   0.394779,
							  - 0.906230, - 1.120195,   1.207566,
							   2.242360, - 0.311222,   0.433339 };

Efter ha tagit FFT på matrisen X så får jag

3.2080352       -1.7104170      1.3830030
-1.5349649      1.1842327       0.8222664
-0.8341455      -0.7355201      1.9282329
2.1994600       4.9222131       -1.4719210
4.6299772       0.5482139       -0.2945260
-3.3765469      -0.5819106      0.2966796
0.5881023       -1.2479200      -1.1604111
-3.2144592      1.5684826       -0.8436050

I matlab så skulle detta se ut som

>> fft(X)
ans =

   3.2080 +      0i  -1.7104 +      0i   1.3830 +      0i
  -1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i
   2.1995 + 4.6300i   4.9222 + 0.5482i  -1.4719 - 0.2945i
  -3.3765 + 0.5881i  -0.5819 - 1.2479i   0.2967 - 1.1604i
  -3.2145 +      0i   1.5685 +      0i  -0.8436 +      0i
  -3.3765 - 0.5881i  -0.5819 + 1.2479i   0.2967 + 1.1604i
   2.1995 - 4.6300i   4.9222 - 0.5482i  -1.4719 + 0.2945i
  -1.5350 + 0.8341i   1.1842 + 0.7355i   0.8223 - 1.9282i

>>

Här får man vara vaksam! Notera att talet 1.5350 + 0.8341i uppträder två gånger i kolumnen och att datat är speglat. FFTPack ger oss endast halva spegeln.

   3.2080 +      0i  -1.7104 +      0i   1.3830 +      0i
  -1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i
   2.1995 + 4.6300i   4.9222 + 0.5482i  -1.4719 - 0.2945i
  -3.3765 + 0.5881i  -0.5819 - 1.2479i   0.2967 - 1.1604i
  -3.2145 +      0i   1.5685 +      0i  -0.8436 +      0i

Och då är frågan: Hur ska jag använda FFT här på denna data för att göra 2D FFT?
Jo! Det jag måste göra, är att jag måste ta FFT för varje rad, inte kolumn.
Men det ställer till lite problem för mig.

Om jag ska göra FFT på denna rad

-1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i

Så måste jag göra FFT på denna rad

-1.5349649      1.1842327       0.8222664
-0.8341455      -0.7355201      1.9282329

Dessa tal skulle fått byta plats. Andra problemet är att dessa är komplexa tal.
Om jag skulle köra FFT på raden ovan så skulle jag få detta resultat

>> fft([-1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i])
ans =

   0.4715 + 0.3586i  -4.8451 - 1.7439i  -0.2314 - 1.1170i

Om jag kör med FFT Pack

float Y[6] = { -1.5350, -0.8341,  1.1842, 0.7355,   0.8223 , 1.9282 };
fft(Y, 6);
print(Y, 1, 6);

>>Utskrift: 2.3011003       -2.7266998      2.0788074       -2.3498001      2.7056365       -1.3581001

Ni ser att det fungerar inte att köra komplexa tal med FFTPack.

Senast redigerat 2023-09-03 11:00

Citera flera Citera

2023-09-03 11:16

Elgot

Hedersmedlem

Plats: Linköping
Registrerad: Apr 2004

●

Skrivet av heretic16:

Ja, det gör jag. MLK är känd på internet för att inte vara den snabbaste. FFTW3 är säkert den snabbaste, men FFTPack, trots sin ålder, så är den en av dom snabbaste, och mest lättanvändarvänligaste också.

Enda skillnaden mellan MLK och FFTPack var att FFTPack kunde ge talen 3.99999999999 medan MLK gav 4.0 exakt.

Nu sitter jag med ett problem att jag måste göra en beräkning med 2D med FFTpack. Bara FFTpack 5.1 har stöd för 2D FFT, men denna är skrivet i Fortran 90, och inte ANSI C som jag skriver min kod i.

Då har jag en plan!

1. Om jag har en matris X som är M x N och sedan tar jag transponatet på X för att göra om den till column-major.
2. Jag använder FFT på varje "rad", som är egentligen varje kolumn med tanke på transponantet.

#define row 8
#define column 3

float X[row * column] = { -0.678816,   1.363392, - 0.020819,
							   0.014361,   0.115923,   0.309136,
							  - 0.195106, - 1.376395, - 0.336756,
							   1.678252, - 0.234396, - 0.023950,
							   1.776940,   1.062231, - 0.580292,
							  - 0.723726, - 1.209755,   0.394779,
							  - 0.906230, - 1.120195,   1.207566,
							   2.242360, - 0.311222,   0.433339 };

Efter ha tagit FFT på matrisen X så får jag

3.2080352       -1.7104170      1.3830030
-1.5349649      1.1842327       0.8222664
-0.8341455      -0.7355201      1.9282329
2.1994600       4.9222131       -1.4719210
4.6299772       0.5482139       -0.2945260
-3.3765469      -0.5819106      0.2966796
0.5881023       -1.2479200      -1.1604111
-3.2144592      1.5684826       -0.8436050

I matlab så skulle detta se ut som

>> fft(X)
ans =

   3.2080 +      0i  -1.7104 +      0i   1.3830 +      0i
  -1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i
   2.1995 + 4.6300i   4.9222 + 0.5482i  -1.4719 - 0.2945i
  -3.3765 + 0.5881i  -0.5819 - 1.2479i   0.2967 - 1.1604i
  -3.2145 +      0i   1.5685 +      0i  -0.8436 +      0i
  -3.3765 - 0.5881i  -0.5819 + 1.2479i   0.2967 + 1.1604i
   2.1995 - 4.6300i   4.9222 - 0.5482i  -1.4719 + 0.2945i
  -1.5350 + 0.8341i   1.1842 + 0.7355i   0.8223 - 1.9282i

>>

Här får man vara vaksam! Notera att talet 1.5350 + 0.8341i uppträder två gånger i kolumnen och att datat är speglat. FFTPack ger oss endast halva spegeln.

   3.2080 +      0i  -1.7104 +      0i   1.3830 +      0i
  -1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i
   2.1995 + 4.6300i   4.9222 + 0.5482i  -1.4719 - 0.2945i
  -3.3765 + 0.5881i  -0.5819 - 1.2479i   0.2967 - 1.1604i
  -3.2145 +      0i   1.5685 +      0i  -0.8436 +      0i

Och då är frågan: Hur ska jag använda FFT här på denna data för att göra 2D FFT?
Jo! Det jag måste göra, är att jag måste ta FFT för varje rad, inte kolumn.
Men det ställer till lite problem för mig.

Om jag ska göra FFT på denna rad

-1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i

Så måste jag göra FFT på denna rad

-1.5349649      1.1842327       0.8222664
-0.8341455      -0.7355201      1.9282329

Dessa tal skulle fått byta plats. Andra problemet är att dessa är komplexa tal.
Om jag skulle köra FFT på raden ovan så skulle jag få detta resultat

>> fft([-1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i])
ans =

   0.4715 + 0.3586i  -4.8451 - 1.7439i  -0.2314 - 1.1170i

Om jag kör med FFT Pack

float Y[6] = { -1.5350, -0.8341,  1.1842, 0.7355,   0.8223 , 1.9282 };
fft(Y, 6);
print(Y, 1, 6);

>>Utskrift: 2.3011003       -2.7266998      2.0788074       -2.3498001      2.7056365       -1.3581001

Ni ser att det fungerar inte att köra komplexa tal med FFTPack.

Handlar det om väldigt små vektorer (tiotals element)? I så fall lönar det sig troligen att använda något batch-läge med mkl/fftw.

Citera flera Citera

2023-09-03 11:37

heretic16

Medlem

Registrerad: Jun 2008

●

Skrivet av Elgot:

Handlar det om väldigt små vektorer (tiotals element)? I så fall lönar det sig troligen att använda något batch-läge med mkl/fftw.

Jag körde 1000x1000, vilket är 100 tusen array. FFTPack var snabbare än MLK. Inte mycket snabbare, men några hundradelar snabbare.

Om MKL hade 0.02 så hade FFTPack 0.01.

Hur som helst. Tror du det är möjligt att beräkna imaginära och realla tal för sig om man kör FFT? FFT med komplexa tal går inte för FFTPack.

Citera flera Citera

2023-09-03 12:37

Yoshman

Datavetare ★

Plats: Stockholm
Registrerad: Jun 2011

●

Skrivet av heretic16:

Att göra en matrismultiplikation med två matriser i C krävs minst 3 for-satser.
Det tar många minuter att multiplicera en 5000x5000 matris med en annan 5000x5000 matris.
Men med MLK går det på några sekunder.

Vad är det som gör MLK så effektiv när det kommer till matrismultiplikationer?
Notera att det är bara matrismultiplikationer som MLK är effektivast på. Övriga så som FFT så är en vanlig enkel FFT algoritm med en for-sats betydligt snabbare.

Som referens är det bra att veta kör du på för HW.

"Men med MLK går det på några sekunder." låter fruktansvärt långsamt för att multiplicera två 5000x5000 matriser. Sitter med PyTorch på en M1 (så inte MKL, men rätt liknande optimeringar, med Intel Macs hade detta använt MKL), där tar en sådan multiplikation ~8ms på CPU och ~5ms på GPU (är egentligen för små matriser för GPU ska fungera väl, att det ens är snabbare här beror nog helt på att det är en iGPU och inte en dGPU).

De optimeringar som MKL gör verkar vara

Använder Strassen algoritm ("fast matrix multiplication"), tar ned komplexiteten från O(N³) till O(N^log₂(7)) (log₂(7)≈2,8)
Använder hand-optimerad SIMD
Använder alla CPU-kärnor

Har man t.ex. en i7-9900K (8C Skylake) ger detta för N=5000 och antar 32-bit flyttal

Strassen: 5000³/5000^2,8≈x5
SIMD (256-bit AVX med FMA för Skylake) 256/32*2*2 -> "upp till" x32 FLOPS per cykel och kärna, borde vara upp mot x8 snabbare än skalära fallet ("for-loopar")
Alla CPU-kärnor, "upp till" x8 (blir lägre i praktiken då frekvensen minskar av både att använda alla kärnor + använda AVX + är inte nog knappast perfekt skalning till 8 kärnor med så pass små matriser)

Har lite sämre koll på FFT, men "for-loop" fallet där verkar vara O(N²) medan en vettig algoritm (t.ex. Cooley–Tukey algorithm) bör hamna på O(N*log(N)). För en vektor med 5000 element blir det då ett par hundra gånger snabbare med den bättre algoritmen.

Visa signatur

Care About Your Craft: Why spend your life developing software unless you care about doing it well? - The Pragmatic Programmer

Citera flera Citera (7)

2023-09-03 13:19

heretic16

Medlem

Registrerad: Jun 2008

●

Skrivet av Yoshman:

Som referens är det bra att veta kör du på för HW.

"Men med MLK går det på några sekunder." låter fruktansvärt långsamt för att multiplicera två 5000x5000 matriser. Sitter med PyTorch på en M1 (så inte MKL, men rätt liknande optimeringar, med Intel Macs hade detta använt MKL), där tar en sådan multiplikation ~8ms på CPU och ~5ms på GPU (är egentligen för små matriser för GPU ska fungera väl, att det ens är snabbare här beror nog helt på att det är en iGPU och inte en dGPU).

De optimeringar som MKL gör verkar vara

Använder Strassen algoritm ("fast matrix multiplication"), tar ned komplexiteten från O(N³) till O(N^log₂(7)) (log₂(7)≈2,8)
Använder hand-optimerad SIMD
Använder alla CPU-kärnor

Har man t.ex. en i7-9900K (8C Skylake) ger detta för N=5000 och antar 32-bit flyttal

Strassen: 5000³/5000^2,8≈x5
SIMD (256-bit AVX med FMA för Skylake) 256/32*2*2 -> "upp till" x32 FLOPS per cykel och kärna, borde vara upp mot x8 snabbare än skalära fallet ("for-loopar")
Alla CPU-kärnor, "upp till" x8 (blir lägre i praktiken då frekvensen minskar av både att använda alla kärnor + använda AVX + är inte nog knappast perfekt skalning till 8 kärnor med så pass små matriser)

Har lite sämre koll på FFT, men "for-loop" fallet där verkar vara O(N²) medan en vettig algoritm (t.ex. Cooley–Tukey algorithm) bör hamna på O(N*log(N)). För en vektor med 5000 element blir det då ett par hundra gånger snabbare med den bättre algoritmen.

Låter riktigt bra.
Jag gjorde en 100 tusen array med normalfördelande siffror.
0.005000 sekunder tog det att beräkna FFT där.

Med FFTPack. Gör så här. Ladda ned dessa två filer.
https://cs.android.com/android/platform/superproject/+/master...
https://cs.android.com/android/platform/superproject/+/master...

Sedan hittar du på någon array x som har längden n
Kör denna kod

void fft(float x[], size_t n) {
	/* Init */
	float* wsave = (float*)malloc((2 * n + 15) * sizeof(float));
	int* ifac = (int*)malloc((n + 15) * sizeof(float));
	rffti(n, wsave, ifac);

	/* Forward transform */
	rfftf(n, x, wsave, ifac);

	/* Free */
	free(wsave);
	free(ifac);
}

Här har jag hittat Fortran 77 kod som kan beräkna FFT2!
https://github.com/jlokimlin/fftpack5.1

Bara en liten snabb koll, tror ni det går att använda f2c för att konvertera denna Fortran 77 kod till C kod?

Senast redigerat 2023-09-03 13:50

Citera flera Citera

2023-09-03 14:10

Elgot

Hedersmedlem

Plats: Linköping
Registrerad: Apr 2004

●

Skrivet av heretic16:

Låter riktigt bra.
Jag gjorde en 100 tusen array med normalfördelande siffror.
0.005000 sekunder tog det att beräkna FFT där.

Med FFTPack. Gör så här. Ladda ned dessa två filer.
https://cs.android.com/android/platform/superproject/+/master...
https://cs.android.com/android/platform/superproject/+/master...

Sedan hittar du på någon array x som har längden n
Kör denna kod

void fft(float x[], size_t n) {
	/* Init */
	float* wsave = (float*)malloc((2 * n + 15) * sizeof(float));
	int* ifac = (int*)malloc((n + 15) * sizeof(float));
	rffti(n, wsave, ifac);

	/* Forward transform */
	rfftf(n, x, wsave, ifac);

	/* Free */
	free(wsave);
	free(ifac);
}

Den där koden verkar inte riktigt anpassad för de fftpack du länkade till (den klagar på fel antal argument till rfftf/rffti).

Om man kör fftw såhär verkar en 1024-punkters fft-beräkning ta ungefär en mikrosekund (enkeltrådat/halvkass dator).

#include <fftw3.h>
#include <chrono>
#include <iostream>

const int fftSize = 1024;
fftwf_complex* fftInput;
fftwf_complex* fftOutput;
fftwf_plan fftPlan;

int main(int argc, char *argv[])
{

    fftwf_import_wisdom_from_filename("wisdom.dat");
    fftInput = fftwf_alloc_complex(fftSize);
    fftOutput = fftwf_alloc_complex(fftSize);
    fftPlan = fftwf_plan_dft_1d(fftSize, fftInput, fftOutput, FFTW_FORWARD, FFTW_EXHAUSTIVE | FFTW_PRESERVE_INPUT);
    fftwf_export_wisdom_to_filename("wisdom.dat");

    for(int i = 0; i < fftSize; ++i)
    {
        fftInput[i][0] = rand()%1000000;
        fftInput[i][1] = rand()%1000000;
    }

    std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();

    for(int q = 0; q < 1000000; ++q)
        fftwf_execute(fftPlan);

    std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();

    std::cout << "Time: " << std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1).count();
    fftwf_destroy_plan(fftPlan);
    fftwf_free(fftInput);
    fftwf_free(fftOutput);
}

Citera flera Citera

2023-09-03 15:59

heretic16

Medlem

Registrerad: Jun 2008

●

Jag vet att FFT3W är något bättre. Men för mig så är 0.001 riktigt snabbt.
Jag har lyckats konvertera om FFTPack 5.1 nu till C-kod. Men det verkar inte finnas någon seriös manual hur man använder FFTPack när det är konverterat från Fortran 77 till ANSI C.

Exempel:

 SUBROUTINE RFFT1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

int rfftf1_(integer *n, integer *in, real *c__, real *ch, real *wa, real *fac);

Den enda manualen jag hittar är: https://github.com/jlokimlin/fftpack5.1/blob/master/doc/FFTPA...

Citera flera Citera

2023-09-03 16:05

Elgot

Hedersmedlem

Plats: Linköping
Registrerad: Apr 2004

●

Skrivet av heretic16:

Jag vet att FFT3W är något bättre. Men för mig så är 0.001 riktigt snabbt.
Jag har lyckats konvertera om FFTPack 5.1 nu till C-kod. Men det verkar inte finnas någon seriös manual hur man använder FFTPack när det är konverterat från Fortran 77 till ANSI C.

Exempel:

 SUBROUTINE RFFT1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

int rfftf1_(integer *n, integer *in, real *c__, real *ch, real *wa, real *fac);

Den enda manualen jag hittar är: https://github.com/jlokimlin/fftpack5.1/blob/master/doc/FFTPA...

Ja, men mkl är inte sämre så om du får till det får du ju allt på köpet.

Citera flera Citera

2023-09-03 16:23

Yoshman

Datavetare ★

Plats: Stockholm
Registrerad: Jun 2011

●

Skrivet av heretic16:

Låter riktigt bra.
Jag gjorde en 100 tusen array med normalfördelande siffror.
0.005000 sekunder tog det att beräkna FFT där.

Med FFTPack. Gör så här. Ladda ned dessa två filer.
https://cs.android.com/android/platform/superproject/+/master...
https://cs.android.com/android/platform/superproject/+/master...

Sedan hittar du på någon array x som har längden n
Kör denna kod

void fft(float x[], size_t n) {
	/* Init */
	float* wsave = (float*)malloc((2 * n + 15) * sizeof(float));
	int* ifac = (int*)malloc((n + 15) * sizeof(float));
	rffti(n, wsave, ifac);

	/* Forward transform */
	rfftf(n, x, wsave, ifac);

	/* Free */
	free(wsave);
	free(ifac);
}

Här har jag hittat Fortran 77 kod som kan beräkna FFT2!
https://github.com/jlokimlin/fftpack5.1

Bara en liten snabb koll, tror ni det går att använda f2c för att konvertera denna Fortran 77 kod till C kod?

Vad har du för CPU? Bara så vi har en idé om resultaten verkar rimliga eller ej.

FFTW verkar bästa val för x86_64, men det verkar vara något långsammare än MacOS accelerate (som bl.a. PyTorch, NumPy m.fl. använder) om man kör Mac/ARM64.

Använde FFTW biblioteket från Homebrew, med det tar det ~0,90ms att köra FFT på en vektor med 100 000 element. Samma sak tar ~0,61ms med PyTorch (på CPU).

$ brew install fftw
<skapa C++ filen, anta att den heter app.cpp>
$ g++ -std=c++17 -O2 -march=native -flto app.cpp -o app $(pkg-config fftw3f --libs --cflags)
$ ./app

Steg att köra kodsnutten @Elgot postade ovan på MacOS/Linux

Senast redigerat 2023-09-03 16:50 Justerade byggraden så den fungerar även på Linux

Visa signatur

Care About Your Craft: Why spend your life developing software unless you care about doing it well? - The Pragmatic Programmer

Citera flera Citera

2023-09-03 16:34

heretic16

Medlem

Registrerad: Jun 2008

●

Skrivet av Yoshman:

Vad har du för CPU? Bara så vi har en idé om resultaten verkar rimliga eller ej.

FFTW verkar bästa val för x86_64, men det verkar vara något långsammare än MacOS accelerate (som bl.a. PyTorch, NumPy m.fl. använder) om man kör Mac/ARM64.

Använde FFTW biblioteket från Homebrew, med det tar det ~0,90ms att köra FFT på en vektor med 100 000 element. Samma sak tar ~0,61ms med PyTorch (på CPU).

$ brew install fftw
<skapa C++ filen, anta att den heter app.cpp>
$ g++ -std=c++17 $(pkg-config fftw3f --libs --cflags) -O2 -march=native -flto app.cpp -o app
$ ./app

Steg att köra kodsnutten @Elgot postade ovan på MacOS

Jag kör en Intel I5:a

Jag har tänkt att jag ska köra FFT på alla typer av system. Mest ARM.

Jag lyckades få till det där nu med Fortran 77 till C!

int rfft1f_(integer *n, integer *inc, real *r__, integer *
	lenr, real *wsave, integer *lensav, real *work, integer *lenwrk, 
	integer *ier)

Jag kan göra så här. Jag testar hur snabbt FFTpack 5.1 är, som är skrivet i Fortran 77. Sedan återkopplar jag i denna tråd. Då har vi något att jämföra med

Citera flera Citera

2023-09-03 16:41

Elgot

Hedersmedlem

Plats: Linköping
Registrerad: Apr 2004

●

Skrivet av Yoshman:

Använde FFTW biblioteket från Homebrew, med det tar det ~0,90ms att köra FFT på en vektor med 100 000 element. Samma sak tar ~0,61ms med PyTorch (på CPU).

Använde du wisdom ungefär som ovan då (och FFTW_EXHAUSTIVE eller liknande)? Man får inte full hastighet om man inte låter den utvärdera vilka algoritmer som är bäst på aktuell hårdvara.

Citera flera Citera

2023-09-03 16:45

Yoshman

Datavetare ★

Plats: Stockholm
Registrerad: Jun 2011

●

Skrivet av Elgot:

Använde du wisdom ungefär som ovan då (och FFTW_EXHAUSTIVE eller liknande)? Man får inte full hastighet om man inte låter den utvärdera vilka algoritmer som är bäst på aktuell hårdvara.

Körde ditt exempel rakt, men ändrade fftSize till 100000 och körde loopen 1000 gånger (för att få rimlig tid).

Visa signatur

Care About Your Craft: Why spend your life developing software unless you care about doing it well? - The Pragmatic Programmer

Citera flera Citera

2023-09-03 16:45

Elgot

Hedersmedlem

Plats: Linköping
Registrerad: Apr 2004

●

Skrivet av heretic16:

Jag kan göra så här. Jag testar hur snabbt FFTpack 5.1 är, som är skrivet i Fortran 77. Sedan återkopplar jag i denna tråd. Då har vi något att jämföra med

Eventuellt kan du kanske låta biblioteket vara kvar i Fortran? Det brukar väl inte vara omöjligt att länka till sådana bibliotek (särskilt inte om det bara gäller enstaka funktioner)?

Citera flera Citera

2023-09-03 16:59

heretic16

Medlem

Registrerad: Jun 2008

●

Skrivet av Elgot:

Eventuellt kan du kanske låta biblioteket vara kvar i Fortran? Det brukar väl inte vara omöjligt att länka till sådana bibliotek (särskilt inte om det bara gäller enstaka funktioner)?

Jag har ingen aning. Orsaken varför jag väljer FFTPack har med att den är portabel. När jag använder f2c så får jag koden till ANSI C kod. Gammal hederlig C kod som aldrig dör ut.

Citera flera Citera

2023-09-03 17:13

Yoshman

Datavetare ★

Plats: Stockholm
Registrerad: Jun 2011

●

Skrivet av Elgot:

Använde du wisdom ungefär som ovan då (och FFTW_EXHAUSTIVE eller liknande)? Man får inte full hastighet om man inte låter den utvärdera vilka algoritmer som är bäst på aktuell hårdvara.

En FYI: FFTW_EXHAUSTIVE/FFTW_MEASURE verkar hänga sig på ARM64 (i alla fall combon Ubuntu server 22.04 + Orange Pi 4, som har 2 st Cortex A72 (samma CPU som RPi4) och 4 st Cortex A53 (samma CPU som RPi3)).

Men här är FFTW snabbare, PyTorch använder sig av Arms ComputeLibray.

Samma exempel som ovan (FFT på signal med 100 000 samples) tog 4,7ms med FFTW och 5,0ms med PyTorch/ComputeLibray.

Detta fall använder bara en kärna, så är på Cortex A72. Kör man på Cortex A53 (som är en riktig potatis, är ju in-order CPU), tar det 12,2ms med FFTW och 15,1ms med PyTorch/ComputeLibray.

Visa signatur

Care About Your Craft: Why spend your life developing software unless you care about doing it well? - The Pragmatic Programmer

Citera flera Citera

2023-09-03 17:15

Elgot

Hedersmedlem

Plats: Linköping
Registrerad: Apr 2004

●

Skrivet av Yoshman:

En FYI: FFTW_EXHAUSTIVE/FFTW_MEASURE verkar hänga sig på ARM64 (i alla fall combon Ubuntu server 22.04 + Orange Pi 4, som har 2 st Cortex A72 (samma CPU som RPi4) och 4 st Cortex A53 (samma CPU som RPi3)).

Den tar inte bara ohemult lång tid på sig då? Det skalar illa med fft-längd och antal kärnor (om man kör flertrådat).

Citera flera Citera

2023-09-03 17:24

Yoshman

Datavetare ★

Plats: Stockholm
Registrerad: Jun 2011

●

Skrivet av Elgot:

Den tar inte bara ohemult lång tid på sig då? Det skalar illa med fft-längd och antal kärnor (om man kör flertrådat).

Tror inte det, sänkte till 1 000 samples, är detta anrop man fastnar i (med 100 % CPU-last på en kärna)

    fftPlan = fftwf_plan_dft_1d(fftSize, fftInput, fftOutput, FFTW_FORWARD, FFTW_EXHAUSTIVE | FFTW_PRESERVE_INPUT);

Kanske problem med att det är en 2P+4E config...
Edit: verkar inte så, samma problem på en RPi4 med samma OS (Ubuntu 22.04 server)

Verkar som FFT inte skalar med kärnor, även med 100 000 samples används bara en kärna både på Ubuntu 22.04 och MacOS (det är då testat mer 3 olika ramverk, FFTW, Arm ComputeLibrary samt Apple Accelerate).

Senast redigerat 2023-09-03 17:45

Visa signatur

Care About Your Craft: Why spend your life developing software unless you care about doing it well? - The Pragmatic Programmer

Citera flera Citera

2023-09-03 17:25

heretic16

Medlem

Registrerad: Jun 2008

●

Preliminär siffra:
1000000 lång array, normalfördelat brus. 0.17 sekunder. FFTPack 5.1

Senast redigerat 2023-09-03 17:47

Citera flera Citera

2023-09-03 17:39

Yoshman

Datavetare ★

Plats: Stockholm
Registrerad: Jun 2011

●

Skrivet av heretic16:

Preliminär siffra:
1000000 lång array, normalfördelat brus. 0.17 sekunder. FFTPack 5.1

Är det på Arm (vilken Arm) eller i5:an (vilken modell)?

Verkar lite långsamt, Orange Pi 4:an tar det på ~0,099s med FFTW och ~0,078s med PyTorch/ComputeLibrary (så tydligen är Arm ComputeLibrary bättre på större datamängder, FFTW var snabbare med 100 000 samples).

Men i.o.f.s. kommer det vara en rätt stor prestandaskillnad om FFTPack kör ANSI C (så "vanlig" skalär matematik) när FFTW/ComputeLibrary kör med hand-kodad SIMD (NEON i fallet Cortex A72).

Skapar data på samma sätt som i Matlabs FFT exempel, d.v.s. två överlagrade sinus funktioner + pålagt brus (som precis som i ditt fall simuleras med normalfördelade slumptal).

Visa signatur

Care About Your Craft: Why spend your life developing software unless you care about doing it well? - The Pragmatic Programmer

Citera flera Citera

2023-09-03 17:48

heretic16

Medlem

Registrerad: Jun 2008

●

Skrivet av Yoshman:

Är det på Arm (vilken Arm) eller i5:an (vilken modell)?

Verkar lite långsamt, Orange Pi 4:an tar det på ~0,099s med FFTW och ~0,078s med PyTorch/ComputeLibrary (så tydligen är Arm ComputeLibrary bättre på större datamängder, FFTW var snabbare med 100 000 samples).

Men i.o.f.s. kommer det vara en rätt stor prestandaskillnad om FFTPack kör ANSI C (så "vanlig" skalär matematik) när FFTW/ComputeLibrary kör med hand-kodad SIMD (NEON i fallet Cortex A72).

Skapar data på samma sätt som i Matlabs FFT exempel, d.v.s. två överlagrade sinus funktioner + pålagt brus (som precis som i ditt fall simuleras med normalfördelade slumptal).

Det är hos en HP I5:a generation 8.
Ja. FFTPack kör ansi, men FFTpack från android-sidan var något snabbare än FFTPack 5.1. Men nu handlar det bara om tiondelar mellan så. Dessutom kommer jag inte göra FFT på en 1000*1000 array, utan snarare 64*128 array

Slutgiltig siffra:
1000*1000 lång array, normalfördelat brus. FFT -> IFFT tog 0.3680 sekunder. Alltså en fördubbling.

En skillnad mellan FFTpack som jag länkade från Android-sidan, den är snabbare, mer optimerad tydligen.
Så jag skulle rekommendera FFTPack från den där Android-sidan än FFTPack 5.1 från GitHub.

Men, FFTPack-sidan från Github har betydligt flera funktioner så som FFT 2D, vilket inte FFTPack från Android-sidan har.

Nu ska jag testa om FFT 2D verkligen visar resultat.

En sak till!
FFTPack 5.1 gav inte samma resultat för tidsplan -> frekvensplan. Men tidsplan -> frekvensplan -> tidsplan gav samma resultat som i MATLAB. Så jag antar att man får andra typer av imaginära tal med FFTPack 5.1. FFTPack från Android-sidan gav exakt samma resultat som MATLAB, om man normaliserade värderna efter IFFT.

Citera flera Citera

2023-09-03 18:40

heretic16

Medlem

Registrerad: Jun 2008

●

Nu har jag testat FFT2 med FFTPack 5.1. Det fungerar.
Det blir dock inte samma MATLAB-resultat vid tidsplan -> frekvensplan.
Men det blir samma resultat vid tidsplan -> frekvensplan -> tidsplan.

Varför, vet jag inte.
Snabbheten för att göra FFT2 -> IFFT2 för en 500*500 matris var 0.051000 sekunder. Jag tycker dessa siffror är godtagliga.

Jag gillar dock inte MLK's FFT implementering. Men den så måste man antingen välja complex.h, att inkludera, eller så måste man ha +2 längre data array för att andra elementet måste ha en noll-imaginär tal, vilket är rätt värdelöst. Då början och slutet på en FFT-array i frekvensplanet är alltid ett reellt tal.

Det jag har gjort är att för 1D FFT så använder jag det första biblioteket som var snabbare än MLK.
För 2D FFT så använder jag FFTPack 5.1. Nu är det snabbt på båda av dimensionerna.

Senast redigerat 2023-09-03 18:58