[C++] Alternativ till for-loop för snabbare uträkningar?

// ...
double F(double x) {

  double y = 10 * sin(x^2) / (x^2);
  return y;

}

double j;
std::vector<double> M;
for (int i = 0; i < 10000; i++) {

  j = F(k[i]);  // k[i] = double value in position i of vector k
  M.push_back(j):

}
// ...

Jag vet inte exakt om detta blir snabbare, har nämligen inte kollat. Men du skulle kunna testa en for each loop istället. Även förändra så du inte mellanlagrar svaren i variabler.
Sen antar jag att du fyller din vector någon annanstans?
Typ såhär:

// ...
double F(double x) {

   return 10 * sin(x^2) / (x^2);
}


std::vector<double> M;
for (int i = 0; i < 10000; i++) {
  M.push_back(F(k[i]));
}
// ...

Det är inte for-loopen i sig som är långsam, utan allt arbete du utför i den. sin är t.ex. en ganska långsam funktion, och push_back kommer leda till en massa minnesallokeringar. Om du måste anropa F för 10000 värden och lagra dem i en vector så kommer det alltså bli långsamt, oavsett hur du gör.

En del knep som du kan använda är att beräkna sin för en massa värden i förväg och lagra dem i en tabell. Det förutsätter dock att du bara behöver kunna använda sin(x) för en begränsad mängd av x. Se t.ex. denna diskussion för lite andra förslag.

Du kan också undvika en massa extra minnesallokeringar genom att reservera plats i vectorn i förväg. En vector i C++ fungerar nämligen så att den allokerar en viss mängd minne. När det minnet tar slut så allokerar den dubbelt så mycket minne som tidigare, och kopierar sedan över den gamla datan från det gamla till det nya minnet. Om man reserverar plats i förväg så allokeras istället bara så mycket minne som du behöver:

std::vector<double> M;
M.reserve(10000);   // Reservera plats för 10000 element.

Om M alltid innehåller 10000 element så kan du också använda en icke-expanderande array. Antingen en vanlig array eller array från standardbiblioteket (kräver C++11):

double M[10000];
std::array<double, 10000> M;

for (int i = 0; i < 10000; ++i) {
  M[i] = F(k[i]);
}

Sen kan du förstås använda flera trådar om du bara vill att det ska gå snabbare utan att du bryr dig om hur mycket processorkraft som används. C++ fick stöd för multitrådning i C++11, och sen finns en massa tredjepartsbibliotek.

På en normal processor vinner du inget på att försöka dig på parallella uträkningar istället för en forloop. Om matlab blir snabbare på det sättet så är matlabs implementation ineffektiv.
Du kan spara några operationer på att reservera kapacitet i vektorn så att du slipper omallokeringar vid pushbacks, och några ytterligare operationer på att använda emplace_back direkt på resultatet av funktionsanropet istället för push_back, om du kompilerar med stöd för c++11.
Men Alltså, du utför 10000 operationer per sekund i en loop, vilket kanske kan översättas till ca 200000 instruktioner per sekund. Den långsammaste icke-jätteembedded-processorn som går att hitta i dagens läge (en raspberry pi zero) klarar ca 600000000 instruktioner per sekund. Det är lugnt.

Mycket annan bra information i tråden redan, men här är ett snabbt exempel på hur du kan parallellisera transform med hjälp av C++11:s async/future utan att behöva underhålla trådar manuellt:

template <typename Iterator, typename T>
vector<T> parallel_transform(Iterator first, Iterator last, const function<T(const T&)>& callback)
{
    auto length = distance(first, last);
    
    // Om listan är kortare än vår (godtyckliga) gräns så anropar vi std::transform() direkt:
    if(length < 20)
    {
        vector<T> res (length);

        transform(first, last, res.begin(), callback);

        return res;
    }
    // ... annars delar vi upp listan i två delar och gör två rekursiva anrop tillbaks till vår
    // parallel_transform(), det ena asynkront och det andra direkt.
    // När de olika grenarnas listor blir korta nog så hanteras de av blocket ovanför.
    else
    {
        Iterator mid = first;
        advance(mid, length / 2);
        
        // Skapa en "future" för att hantera resultatet av den asynkrona halvan av operationen:
        future<vector<T>> future = async(launch::async, &parallel_transform<Iterator, T>, first, mid, callback);
        
        // Medans vi väntar hanterar vi den andra halvan direkt:
        vector<T> first_half  = parallel_transform(first, mid, callback);
        
        // ... och när vi är färdiga med det hämtar vi upp resultatet från vår "future",
        // som förhoppningsvis har utförts på en annan tråd under tiden:
        vector<T> second_half = future.get();

        first_half.insert(first_half.end(), second_half.begin(), second_half.end());

        return first_half;
    }
}

#include <iostream>
#include <vector>
#include <future>
#include <algorithm>
#include <chrono>

using namespace std;

template <typename Iterator, typename T>
vector<T> parallel_transform(Iterator first, Iterator last, const function<T(const T&)>& callback)
{
    auto length = distance(first, last);

    if(length < 20)
    {
        vector<T> res (length);

        transform(first, last, res.begin(), callback);

        return res;
    }
    else
    {
        Iterator mid = first;
        advance(mid, length / 2);
        
        future<vector<T>> future = async(launch::async, &parallel_transform<Iterator, T>, first, mid, callback);
        
        vector<T> first_half  = parallel_transform(first, mid, callback);
        vector<T> second_half = future.get();

        first_half.insert(first_half.end(), second_half.begin(), second_half.end());

        return first_half;
    }
}

int main()
{
    vector<double> input (10000, 1);

    function<double(const double&)> callback = [](const double& value)
    {
        // Simulera lite arbete
        this_thread::sleep_for(chrono::milliseconds(1));

        return value;
    };

    auto start  = chrono::steady_clock::now();
    auto result = parallel_transform(input.begin(), input.end(), callback);
    auto end    = chrono::steady_clock::now();

    cout
        << chrono::duration<double, milli>(end - start).count()
        << " ms" << endl;
}

Då jag aldrig använt boost.SIMD kändes det lika bra att ändra på det. Tydligen ska boost.SIMD stoppas in i framtida C++ standard-bibliotek och måste säga att det verkligen är det mest genomtänkta SIMD-bibliotek jag sett så här långt.

Enligt jämförelser i presentationer jag hittade ska boost.SIMD prestera väldigt nära motsvarande kod handskriven med s.k. compiler intrinsics för SSE/AVX. Finns massor med problem med att använda intrinsics, det är relativt komplicerat men främst så binder man då koden till specifik instruktioner -> det går överhuvudtaget inte att kompilera koden för en CPU som saknar instruktionerna man använt.

Håller mig till samma exempel som ovan, tar mig också friheten att förutsätta att antal element är jämt delbar med bundle_t::static_size (som är antal flyttal i varje anrop till F())

#include <boost/simd/function/store.hpp>
#include <boost/simd/pack.hpp>
#include <boost/simd/trigonometric.hpp>
...
// Måste vara en potens-av-två stycken flyttal i varje operation
// Får plats 4 st double per register i AVX, men är OK att stoppa in fler än så
// något som ger en liten prestandaboost då man sprider ut kostnaden för att
// anropa F() över lite fler flyttal. Över 8 st gjorde i princip inget på min dator.
#define BUNDLE_ORDER    3
#define BUNDLE_LEN      (1u << BUNDLE_ORDER)

typedef boost::simd::pack<double, BUNDLE_LEN> bundle_t;
...
bundle_t F(bundle_t x) 
{
    return 10 * boost::simd::sin(x * x) / (x * x);
}
...
    double M[10000];
    std::array<double, 10000> M;

    for (int i = 0; i < 10000; i += bundle_t::static_size) {
        // Skapa argument till F av bundle_t::static_size stycken flyttal
        // Lagra svaren (bundle_t::static_size stycken) från M + i 
        boost::simd::store(F(bundle_t{k + i}), M + i);
    }

Notera att ovan utnyttjar SIMD så det kan kombineras med att använda flera CPU-kärnor, t.ex. via OpenMP eller CilkPlus!

Man får inte glömma att skicka med -std=c++11 när man använder boost.SIMD (blir VÄLDIGT många kompileringsfel annars). Vidare måste man endera lägga till -msse2/-msse4.1/-mavx beroende på vad man har. Vill man bara bygga för sin egen CPU är det enklast att köra -march=native så får man högsta möjliga SIMD-stöd (AVX finns i.o.f.s. sedan Sandy Bridge så det är ändå rätt "safe").

Edit: dokumentationen för boost.SIMD säger att man ska använda "vanliga" boost i version 1.61. Det är den absolut senaste versionen av boost, så även om man kör Ubuntu 16.04 så har man inte en tillräckligt ny version av boost (testade i.o.f.s. aldrig med den version jag hade installerad).