Dieser Beitrag wurde ursprünglich im Blog von Quansight Labs ver?ffentlicht und wurde hier mit Genehmigung von Quansight ge?ndert und erneut ver?ffentlicht.

Webanwendungen entwickeln sich schnell zu einer neuen Grenze für leistungsstarke wissenschaftliche Berechnungen und KI-gestützte Endbenutzererlebnisse. Grundlage der ML/KI-Revolution ist die lineare Algebra, ein Zweig der Mathematik, der sich mit linearen Gleichungen und deren Darstellungen in Vektorr?umen und über Matrizen besch?ftigt. LAPACK (?Linear Algebra Package“) ist eine grundlegende Softwarebibliothek für numerische lineare Algebra, die robuste, kampferprobte Implementierungen g?ngiger Matrixoperationen bereitstellt . Obwohl LAPACK ein grundlegender Bestandteil der meisten Programmiersprachen und Bibliotheken für numerische Computer ist, muss eine umfassende, qualitativ hochwertige LAPACK-Implementierung, die auf die besonderen Einschr?nkungen des Webs zugeschnitten ist, noch verwirklicht werden. Das ist...bis jetzt.

Anfang dieses Jahres hatte ich das gro?e Glück, ein Sommerpraktikant bei Quansight Labs zu sein, der gemeinnützigen Abteilung von Quansight und führend im wissenschaftlichen Python-?kosystem. W?hrend meines Praktikums habe ich daran gearbeitet, anf?ngliche LAPACK-Unterstützung zu stdlib hinzuzufügen, einer grundlegenden Bibliothek für wissenschaftliche Berechnungen, die in C und JavaScript geschrieben und für die Verwendung in Webbrowsern und anderen webnativen Umgebungen wie Node.js und Deno optimiert ist. In diesem Blogbeitrag bespreche ich meine Reise, einige erwartete und unerwartete (!) Herausforderungen und den weiteren Weg. Ich hoffe, dass diese Arbeit mit etwas Glück einen entscheidenden Baustein dafür liefert, Webbrowser zu einer erstklassigen Umgebung für numerische Berechnungen und maschinelles Lernen zu machen und eine Zukunft mit leistungsf?higeren KI-gestützten Webanwendungen vorwegzunehmen.

Klingt interessant? Auf geht's!

Was ist stdlib?

Leser dieses Blogs, die mit LAPACK vertraut sind, sind wahrscheinlich nicht so genau mit der wilden Welt der Webtechnologien vertraut. Für diejenigen, die aus der Welt der numerischen und wissenschaftlichen Berechnungen kommen und mit dem wissenschaftlichen Python-?kosystem vertraut sind, kann man sich stdlib am einfachsten als eine Open-Source-Bibliothek für wissenschaftliche Berechnungen nach dem Vorbild von NumPy und SciPy vorstellen. Es bietet mehrdimensionale Array-Datenstrukturen und zugeh?rige Routinen für Mathematik, Statistik und lineare Algebra, verwendet jedoch JavaScript anstelle von Python als prim?re Skriptsprache. Daher ist stdlib auf das Web-?kosystem und seine Anwendungsentwicklungsparadigmen ausgerichtet. Dieser Fokus erfordert einige interessante Design- und Projektarchitekturentscheidungen, die stdlib im Vergleich zu traditionelleren Bibliotheken, die für numerische Berechnungen entwickelt wurden, ziemlich einzigartig machen.

Um NumPy als Beispiel zu nehmen: NumPy ist eine einzelne monolithische Bibliothek, in der alle ihre Komponenten, abgesehen von optionalen Abh?ngigkeiten von Drittanbietern wie OpenBLAS, eine einzige, unteilbare Einheit bilden. Man kann nicht einfach NumPy-Routinen zur Array-Manipulation installieren, ohne NumPy vollst?ndig zu installieren. Wenn Sie eine Anwendung bereitstellen, die nur das ndarray-Objekt von NumPy und einige seiner Manipulationsroutinen ben?tigt, bedeutet die Installation und Bündelung von NumPy insgesamt, dass eine betr?chtliche Menge ?toter Code“ einbezogen wird. Im Webentwicklungsjargon würden wir sagen, dass NumPy nicht ?baumschüttelbar“ ist. Für eine normale NumPy-Installation bedeutet dies mindestens 30 MB Festplattenspeicher und mindestens 15 MB Festplattenspeicher für einen benutzerdefinierten Build, der alle Debug-Anweisungen ausschlie?t. Für SciPy k?nnen diese Zahlen auf 130 MB bzw. 50 MB ansteigen. Es erübrigt sich zu erw?hnen, dass die Auslieferung einer 15-MB-Bibliothek in einer Webanwendung für nur wenige Funktionen ein Kinderspiel ist, insbesondere für Entwickler, die Webanwendungen auf Ger?ten mit schlechter Netzwerkkonnektivit?t oder Speicherbeschr?nkungen bereitstellen müssen.

Angesichts der besonderen Einschr?nkungen bei der Entwicklung von Webanwendungen verfolgt stdlib beim Design einen Bottom-up-Ansatz, bei dem jede Funktionseinheit unabh?ngig von nicht verwandten und ungenutzten Teilen der Codebasis installiert und genutzt werden kann. Durch den Einsatz einer zerlegbaren Softwarearchitektur und radikaler Modularit?t bietet stdlib Benutzern die M?glichkeit, genau das zu installieren und zu verwenden, was sie ben?tigen, mit wenig bis gar keinem überschüssigen Code über einen gewünschten Satz von APIs und deren expliziten Abh?ngigkeiten hinaus, wodurch ein geringerer Speicherbedarf gew?hrleistet wird Gr??en und schnellere Bereitstellung.

Angenommen, Sie arbeiten mit zwei Matrizenstapeln (d. h. zweidimensionalen Scheiben dreidimensionaler Würfel) und m?chten jede zweite Scheibe ausw?hlen und die übliche BLAS-Operation y = a * x ausführen. Dabei sind x und y ndarrays und a eine Skalarkonstante. Um dies mit NumPy zu tun, müssen Sie zun?chst NumPy vollst?ndig installieren

pip install numpy

und führen Sie dann die verschiedenen Vorg?nge aus

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

Mit stdlib k?nnen Sie nicht nur das Projekt als monolithische Bibliothek installieren, sondern auch die verschiedenen Funktionseinheiten als separate Pakete installieren

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

und führen Sie dann die verschiedenen Vorg?nge aus

// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';

// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;

// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');

// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);

Wichtig ist, dass Sie nicht nur jedes der über 4.000 Pakete von stdlib unabh?ngig installieren k?nnen, sondern Sie k?nnen auch jedes dieser Pakete reparieren, verbessern und neu mischen, indem Sie ein zugeh?riges GitHub-Repository forken (siehe z. B. @stdlib/ndarray-fancy ). Durch die Definition expliziter Abstraktionsebenen und Abh?ngigkeitsb?ume bietet Ihnen stdlib die Freiheit, die richtige Abstraktionsebene für Ihre Anwendung auszuw?hlen. In mancher Hinsicht ist es eine einfache – und, wenn Sie mit dem herk?mmlichen Design wissenschaftlicher Softwarebibliotheken vertraut sind, vielleicht unorthodoxe – Idee, aber wenn sie eng in die Webplattform integriert ist, hat sie wirkungsvolle Konsequenzen und schafft aufregende neue M?glichkeiten!

Was ist mit WebAssembly?

Okay, vielleicht ist Ihr Interesse geweckt; stdlib scheint faszinierend. Aber was hat das mit LAPACK in Webbrowsern zu tun? Nun, eines unserer Ziele im vergangenen Sommer war es, das stdlib-Ethos – kleine, eng begrenzte Pakete, die eine Sache tun und eine Sache gut machen – anzuwenden, um LAPACK ins Web zu bringen.

Aber warte, sagst du! Das ist ein extremes Unterfangen. LAPACK ist mit etwa 1.700 Routinen umfangreich, und selbst die Umsetzung von 10 % davon innerhalb eines angemessenen Zeitrahmens ist eine gro?e Herausforderung. W?re es nicht besser, LAPACK einfach in WebAssembly zu kompilieren, einem tragbaren Kompilierungsziel für Programmiersprachen wie C, Go und Rust, das die Bereitstellung im Web erm?glicht, und Schluss damit zu machen?

Leider gibt es bei diesem Ansatz mehrere Probleme.

1. Das Kompilieren von Fortran zu WebAssembly ist immer noch ein Bereich der aktiven Entwicklung (siehe 1, 2, 3, 4 und 5). Zum Zeitpunkt dieses Beitrags bestand ein g?ngiger Ansatz darin, Fortran mit f2c in C zu kompilieren und dann einen separaten Kompilierungsschritt durchzuführen, um C in WebAssembly zu konvertieren. Dieser Ansatz ist jedoch problematisch, da f2c nur Fortran 77 vollst?ndig unterstützt und der generierte Code umfangreiche Patches erfordert. Derzeit wird an der Entwicklung eines LLVM-basierten Fortran-Compilers gearbeitet, es bestehen jedoch weiterhin Lücken und komplexe Toolchains.
2. Wie oben in der Diskussion über monolithische Bibliotheken in Webanwendungen angedeutet, ist die Gr??e von LAPACK Teil des Problems. Selbst wenn das Kompilierungsproblem gel?st ist, bedeutet das Einfügen einer einzigen WebAssembly-Bin?rdatei, die das gesamte LAPACK enth?lt, in eine Webanwendung, die nur eine oder zwei LAPACK-Routinen verwenden muss, erheblichen toten Code, was zu langsameren Ladezeiten und einem erh?hten Speicherverbrauch führt.
3. W?hrend man versuchen k?nnte, einzelne LAPACK-Routinen in eigenst?ndige WebAssembly-Bin?rdateien zu kompilieren, k?nnte dies zu einer Aufbl?hung der Bin?rdateien führen, da mehrere eigenst?ndige Bin?rdateien m?glicherweise doppelten Code aus gemeinsamen Abh?ngigkeiten enthalten. Um die bin?re Aufbl?hung zu verringern, k?nnte man versuchen, eine Modulaufteilung durchzuführen. In diesem Szenario gliedert man zun?chst allgemeine Abh?ngigkeiten in eine eigenst?ndige Bin?rdatei mit gemeinsam genutztem Code aus und generiert dann separate Bin?rdateien für einzelne APIs. Obwohl dies in manchen F?llen sinnvoll ist, kann dies schnell unhandlich werden, da dieser Ansatz die Verknüpfung einzelner WebAssembly-Module zur Ladezeit erfordert, indem die Exporte eines oder mehrerer Module mit den Importen eines oder mehrerer anderer Module zusammengefügt werden. Dies kann nicht nur mühsam sein, sondern dieser Ansatz führt auch zu Leistungseinbu?en, da WebAssembly-Routinen beim Aufruf importierter Exporte nun in JavaScript wechseln müssen, anstatt in WebAssembly zu bleiben. Klingt komplex? Das ist es!
4. Abgesehen von WebAssembly-Modulen, die ausschlie?lich mit skalaren Eingabeargumenten arbeiten (z. B. den Sinus einer einzelnen Zahl berechnen), muss jede WebAssembly-Modulinstanz mit WebAssembly-Speicher verknüpft sein, der in festen Schritten von 64 KB (d. h. einer ?Seite“) zugewiesen wird "). Und was noch wichtiger ist: Zum Zeitpunkt dieses Blogbeitrags kann der WebAssembly-Speicher nur wachsen und niemals schrumpfen. Da es derzeit keinen Mechanismus zum Freigeben von Speicher für einen Host gibt, kann der Speicherbedarf einer WebAssembly-Anwendung nur zunehmen. Diese beiden Aspekte zusammen erh?hen die Wahrscheinlichkeit der Zuweisung von Speicher, der nie verwendet wird, und die H?ufigkeit von Speicherlecks.
5. Schlie?lich ist WebAssembly zwar leistungsstark, erfordert aber eine steilere Lernkurve und einen komplexeren Satz sich h?ufig schnell entwickelnder Toolchains. In Endbenutzeranwendungen führt die Schnittstelle zwischen JavaScript – einer webnativen, dynamisch kompilierten Programmiersprache – und WebAssembly zus?tzlich zu einer erh?hten Komplexit?t, insbesondere wenn eine manuelle Speicherverwaltung durchgeführt werden muss.

Um den letzten Punkt zu veranschaulichen, kehren wir zur BLAS-Routine daxpy zurück, die die Operation y = a*x y ausführt und wobei x und y Schrittvektoren und a eine Skalarkonstante sind. Bei einer Implementierung in C k?nnte eine grundlegende Implementierung wie das folgende Codefragment aussehen.

pip install numpy

Nach der Kompilierung in WebAssembly und dem Laden der WebAssembly-Bin?rdatei in unsere Webanwendung müssen wir eine Reihe von Schritten ausführen, bevor wir die c_daxpy-Routine aus JavaScript aufrufen k?nnen. Zuerst müssen wir ein neues WebAssembly-Modul instanziieren.

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

Als n?chstes müssen wir den Modulspeicher definieren und eine neue WebAssembly-Modulinstanz erstellen.

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

Nachdem wir eine Modulinstanz erstellt haben, k?nnen wir nun die exportierte BLAS-Routine aufrufen. Wenn Daten jedoch au?erhalb des Modulspeichers definiert werden, müssen wir diese Daten zun?chst in die Speicherinstanz kopieren und dies immer in Little-Endian-Byte-Reihenfolge tun.

// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';

// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;

// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');

// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);

Da nun die Daten in den Modulspeicher geschrieben wurden, k?nnen wir die c_daxpy-Routine aufrufen.

void c_daxpy(const int N, const double alpha, const double *X, const int strideX, double *Y, const int strideY) {
    int ix;
    int iy;
    int i;
    if (N <= 0) {
        return;
    }
    if (alpha == 0.0) {
        return;
    }
    if (strideX < 0) {
        ix = (1-N) * strideX;
    } else {
        ix = 0;
    }
    if (strideY < 0) {
        iy = (1-N) * strideY;
    } else {
        iy = 0;
    }
    for (i = 0; i < N; i++) {
        Y[iy] += alpha * X[ix];
        ix += strideX;
        iy += strideY;
    }
    return;
}

Und schlie?lich, wenn wir die Ergebnisse zur Visualisierung oder weiteren Analyse an eine Downstream-Bibliothek ohne Unterstützung für WebAssembly-Speicher-?Zeiger“ (d. h. Byte-Offsets) wie D3 übergeben müssen, müssen wir Daten aus dem Modul kopieren Speicher zurück zum ursprünglichen Ausgabearray.

const binary = new UintArray([...]);

const mod = new WebAssembly.Module(binary);

Allein die Berechnung von y = a*x y ist eine Menge Arbeit. Vergleichen Sie es im Gegensatz dazu mit einer einfachen JavaScript-Implementierung, die wie das folgende Code-Snippet aussehen k?nnte.

// Initialize 10 pages of memory and allow growth to 100 pages:
const mem = new WebAssembly.Memory({
    'initial': 10,  // 640KiB, where each page is 64KiB
    'maximum': 100  // 6.4MiB
});

// Create a new module instance:
const instance = new WebAssembly.Instance(mod, {
    'env': {
        'memory': mem
    }
});

Mit der JavaScript-Implementierung k?nnen wir daxpy dann direkt mit unseren extern definierten Daten aufrufen, ohne die im obigen WebAssembly-Beispiel erforderliche Datenverschiebung.

// External data:
const xdata = new Float64Array([...]);
const ydata = new Float64Array([...]);

// Specify a vector length:
const N = 5;

// Specify vector strides (in units of elements):
const strideX = 2;
const strideY = 4;

// Define pointers (i.e., byte offsets) for storing two vectors:
const xptr = 0;
const yptr = N * 8; // 8 bytes per double

// Create a DataView over module memory:
const view = new DataView(mem.buffer);

// Resolve the first indexed elements in both `xdata` and `ydata`:
let offsetX = 0;
if (strideX < 0) {
    offsetX = (1-N) * strideX;
}
let offsetY = 0;
if (strideY < 0) {
    offsetY = (1-N) * strideY;
}

// Write data to the memory instance:
for (let i = 0; i < N; i++) {
    view.setFloat64(xptr+(i*8), xdata[offsetX+(i*strideX)], true);
    view.setFloat64(yptr+(i*8), ydata[offsetY+(i*strideY)], true);
}

Zumindest in diesem Fall ist der WebAssembly-Ansatz nicht nur weniger ergonomisch, sondern es gibt auch, wie angesichts der erforderlichen Datenbewegung zu erwarten ist, negative Auswirkungen auf die Leistung, wie in der folgenden Abbildung dargestellt.

Abbildung 1: Leistungsvergleich der C-, JavaScript- und WebAssembly (Wasm)-Implementierungen von stdlib für die BLAS-Routine daxpy für zunehmende Array-L?ngen (x-Achse). Im Wasm (Kopie)-Benchmark werden Eingabe- und Ausgabedaten in den Wasm-Speicher und aus dem Wasm-Speicher kopiert, was zu einer schlechteren Leistung führt.

In der Abbildung oben zeige ich einen Leistungsvergleich der C-, JavaScript- und WebAssembly (Wasm)-Implementierungen von stdlib für die BLAS-Routine daxpy für zunehmende Array-L?ngen, wie entlang der x-Achse aufgelistet. Die y-Achse zeigt eine normalisierte Rate relativ zu einer Baseline-C-Implementierung. Im Wasm-Benchmark werden Eingabe- und Ausgabedaten direkt im Speicher des WebAssembly-Moduls zugewiesen und bearbeitet, und im Wasm-Benchmark (Kopieren) werden Eingabe- und Ausgabedaten in den und aus dem Speicher des WebAssembly-Moduls kopiert, wie oben erl?utert. Aus der Tabelle k?nnen wir Folgendes erkennen:

1. Im Allgemeinen kann JavaScript-Code, wenn er sorgf?ltig geschrieben wird, dank hochoptimierter Just-in-Time-Compiler (JIT) nur zwei- bis dreimal langsamer ausgeführt werden als nativer Code. Dieses Ergebnis ist beeindruckend für eine lose typisierte, dynamisch kompilierte Programmiersprache und bleibt, zumindest für Daxpy, über verschiedene Array-L?ngen hinweg konsistent.
2. Wenn die Datengr??e und damit die in einem WebAssembly-Modul verbrachte Zeit zunehmen, kann WebAssembly eine nahezu native (~1,5-fache) Geschwindigkeit erreichen. Dieses Ergebnis stimmt allgemeiner mit der erwarteten WebAssembly-Leistung überein.
3. W?hrend WebAssembly nahezu native Geschwindigkeit erreichen kann, k?nnen sich Datenbewegungsanforderungen negativ auf die Leistung auswirken, wie bei Daxpy beobachtet. In solchen F?llen kann eine gut gestaltete JavaScript-Implementierung, die solche Anforderungen vermeidet, die gleiche, wenn nicht sogar bessere Leistung erzielen, wie es bei daxpy der Fall ist.

Insgesamt kann WebAssembly Leistungsverbesserungen bieten; Die Technologie ist jedoch kein Allheilmittel und muss sorgf?ltig eingesetzt werden, um die gewünschten Vorteile zu erzielen. Und selbst wenn eine überlegene Leistung geboten wird, müssen solche Gewinne gegen die Kosten erh?hter Komplexit?t, potenziell gr??erer Paketgr??en und komplexerer Toolchains abgewogen werden. Für viele Anwendungen reicht eine einfache JavaScript-Implementierung vollkommen aus.

Radikale Modularit?t

Nachdem ich nun den Fall angestrengt habe, einfach das gesamte LAPACK zu WebAssembly zu kompilieren und Schluss zu machen, wo bleibt uns das? Nun, wenn wir das stdlib-Ethos annehmen wollen, brauchen wir eine radikale Modularit?t.

Um radikale Modularit?t anzunehmen, muss man erkennen, dass das Beste in hohem Ma?e kontextabh?ngig ist und dass Entwickler je nach den Anforderungen und Einschr?nkungen von Benutzeranwendungen die Flexibilit?t ben?tigen, die richtige Abstraktion auszuw?hlen. Wenn ein Entwickler eine Node.js-Anwendung schreibt, bedeutet dies m?glicherweise die Bindung an hardwareoptimierte Bibliotheken wie OpenBLAS, Intel MKL oder Apple Accelerate, um eine bessere Leistung zu erzielen. Wenn ein Entwickler eine Webanwendung bereitstellt, die einen kleinen Satz numerischer Routinen ben?tigt, ist JavaScript wahrscheinlich das richtige Werkzeug für diese Aufgabe. Und wenn ein Entwickler an einer gro?en, ressourcenintensiven WebAssembly-Anwendung arbeitet (z. B. für die Bildbearbeitung oder eine Gaming-Engine), ist es von gr??ter Bedeutung, einzelne Routinen einfach als Teil der gr??eren Anwendung kompilieren zu k?nnen. Kurz gesagt, wir wollen ein radikal modulares LAPACK.

Meine Mission bestand darin, den Grundstein für ein solches Unterfangen zu legen, die Hürden zu beseitigen und die Lücken zu finden und uns hoffentlich der leistungsstarken linearen Algebra im Web ein paar Schritte n?her zu bringen. Doch wie sieht radikale Modularit?t aus? Alles beginnt mit der grundlegenden Funktionseinheit, dem Paket.

Jedes Paket in stdlib ist ein eigenst?ndiges Ding, das kolokalisierte Tests, Benchmarks, Beispiele, Dokumentation, Build-Dateien und zugeh?rige Metadaten (einschlie?lich der Aufz?hlung etwaiger Abh?ngigkeiten) enth?lt und eine klare API-Oberfl?che mit der Au?enwelt definiert . Um LAPACK-Unterstützung zu stdlib hinzuzufügen, müssen Sie für jede LAPACK-Routine ein separates eigenst?ndiges Paket mit der folgenden Struktur erstellen:

pip install numpy

Kurz,

• Benchmark: ein Ordner mit Mikro-Benchmarks zur Bewertung der Leistung im Vergleich zu einer Referenzimplementierung (d. h. Referenz-LAPACK).
• docs: ein Ordner mit Hilfsdokumentation, einschlie?lich REPL-Hilfetext und TypeScript-Deklarationen, die typisierte API-Signaturen definieren.
• Beispiele: ein Ordner mit ausführbarem Demonstrationscode, der nicht nur als Dokumentation dient, sondern Entwicklern auch bei der überprüfung des Implementierungsverhaltens hilft.
• include: ein Ordner mit C-Header-Dateien.
• lib: ein Ordner mit JavaScript-Quellimplementierungen, wobei index.js als Paketeinstiegspunkt dient und andere *.js-Dateien interne Implementierungsmodule definieren.
• src: ein Ordner mit C- und Fortran-Quellimplementierungen. Jedes modulare LAPACK-Paket sollte eine leicht modifizierte Fortran-Referenzimplementierung enthalten (F77 für Freiform-Fortran). C-Dateien umfassen eine einfache C-Implementierung, die der Fortran-Referenzimplementierung folgt, einen Wrapper zum Aufrufen der Fortran-Referenzimplementierung, einen Wrapper zum Aufrufen hardwareoptimierter Bibliotheken (z. B. OpenBLAS) in serverseitigen Anwendungen und eine native Bindung zum Aufrufen in kompilierte Dateien C von JavaScript in Node.js oder einer kompatiblen serverseitigen JavaScript-Laufzeitumgebung.
• test: ein Ordner mit Komponententests zum Testen des erwarteten Verhaltens sowohl in JavaScript als auch in nativen Implementierungen. Tests für native Implementierungen werden in JavaScript geschrieben und nutzen die native Bindung für die Zusammenarbeit zwischen JavaScript und C/Fortran.
• binding.gyp/include.gypi: Build-Dateien zum Kompilieren der nativen Node.js-Add-ons, die eine Brücke zwischen JavaScript und nativem Code bilden.
• manifest.json: eine Konfigurationsdatei für die interne Verwaltung von C- und Fortran-kompilierten Quelldateipaketen von stdlib.
• package.json: eine Datei mit Paketmetadaten, einschlie?lich der Aufz?hlung externer Paketabh?ngigkeiten und einem Pfad zu einer einfachen JavaScript-Implementierung zur Verwendung in browserbasierten Webanwendungen.
• README.md: eine Datei mit der Hauptdokumentation eines Pakets, die API-Signaturen und Beispiele für JavaScript- und C-Schnittstellen enth?lt.

Angesichts der anspruchsvollen Dokumentations- und Testanforderungen von stdlib ist das Hinzufügen von Unterstützung für jede Routine ein ordentlicher Arbeitsaufwand, aber das Endergebnis ist robuster, qualitativ hochwertiger und vor allem modularer Code, der als Grundlage für wissenschaftliche Berechnungen geeignet ist im modernen Web. Aber genug mit den Vorrunden! Kommen wir zur Sache!

Ein mehrstufiger Ansatz

Aufbauend auf früheren Bemühungen, die BLAS-Unterstützung zu stdlib hinzuzufügen, haben wir uns entschieden, beim Hinzufügen der LAPACK-Unterstützung einen ?hnlichen mehrstufigen Ansatz zu verfolgen, bei dem wir zun?chst JavaScript-Implementierungen und die damit verbundenen Tests und Dokumentation priorisieren und dann, sobald Tests und Dokumentation vorhanden sind , Backfill-C- und Fortran-Implementierungen und alle damit verbundenen nativen Bindungen an hardwareoptimierte Bibliotheken. Dieser Ansatz erm?glicht es uns, sozusagen einige frühe Punkte auf den Tisch zu bringen, APIs schnell den Benutzern zug?nglich zu machen, robuste Testverfahren und Benchmarks zu etablieren und potenzielle M?glichkeiten für Tools und Automatisierung zu untersuchen, bevor wir uns in die Tiefe der Build-Toolchains und der Leistung stürzen Optimierungen. Aber wo soll man überhaupt anfangen?

Um zu bestimmen, auf welche LAPACK-Routinen zuerst abgezielt werden soll, habe ich den Fortran-Quellcode von LAPACK analysiert, um ein Aufrufdiagramm zu erstellen. Dadurch konnte ich den Abh?ngigkeitsbaum für jede LAPACK-Routine ableiten. Mit dem Diagramm in der Hand führte ich dann eine topologische Sortierung durch, die mir dabei half, Routinen ohne Abh?ngigkeiten zu identifizieren, die unweigerlich Bausteine ??für andere Routinen sein werden. W?hrend ein Tiefenansatz, bei dem ich eine bestimmte High-Level-Routine ausw?hle und rückw?rts arbeite, es mir erm?glichen würde, ein bestimmtes Feature zu landen, k?nnte ein solcher Ansatz dazu führen, dass ich beim Versuch, Routinen mit zunehmender Komplexit?t zu implementieren, stecken bleibe. Durch die Konzentration auf die ?Bl?tter“ des Diagramms konnte ich h?ufig verwendete Routinen (d. h. Routinen mit hohen Ingraden) priorisieren und so meine Wirkung maximieren, indem ich die M?glichkeit freischaltete, sp?ter mehrere Routinen auf h?herer Ebene bereitzustellen meine Bemühungen oder von anderen Mitwirkenden.

Mit meinem Plan in der Hand freute ich mich darauf, mich an die Arbeit zu machen. Für meine erste Routine habe ich dlaswp ausgew?hlt, das eine Reihe von Zeilenaustauschvorg?ngen auf einer allgemeinen rechteckigen Matrix gem?? einer bereitgestellten Liste von Pivot-Indizes durchführt und ein wichtiger Baustein für die LU-Zerlegungsroutinen von LAPACK ist. Und da begannen meine Herausforderungen...

Herausforderungen

Legacy-Fortran

Vor meinem Praktikum bei Quansight Labs habe ich regelm??ig an LFortran, einem modernen interaktiven Fortran-Compiler, der auf LLVM aufbaut, mitgewirkt (und bin es immer noch!), und ich war ziemlich sicher in meinen Fortran-Kenntnissen. Eine meiner ersten Herausforderungen bestand jedoch einfach darin, zu verstehen, was heute als ?alter“ Fortran-Code gilt. Im Folgenden hebe ich drei anf?ngliche Hürden hervor.

Formatierung

LAPACK wurde ursprünglich in FORTRAN 77 (F77) geschrieben. W?hrend die Bibliothek in Version 3.2 (2008) auf Fortran 90 verschoben wurde, bestehen in der Referenzimplementierung weiterhin alte Konventionen. Eine der sichtbarsten dieser Konventionen ist die Formatierung.

Entwickler, die F77-Programme schrieben, verwendeten dazu ein festes Formularlayout, das von Lochkarten übernommen wurde. Dieses Layout stellte strenge Anforderungen an die Verwendung von Zeichenspalten:

• Kommentare, die eine ganze Zeile einnehmen, müssen mit einem Sonderzeichen (z. B. *, ! oder C) in der ersten Spalte beginnen.
• Für Nicht-Kommentarzeilen müssen 1) die ersten fünf Spalten leer sein oder eine numerische Beschriftung enthalten, 2) Spalte sechs ist für Fortsetzungszeichen reserviert, 3) ausführbare Anweisungen müssen in Spalte sieben beginnen und 4) jeglicher Code darüber hinaus Spalte 72 wurde ignoriert.

Fortran 90 führte das Freiform-Layout ein, das die Beschr?nkungen der Spalten- und Zeilenl?nge aufhob und sich entschied ! als Kommentarzeichen. Der folgende Codeausschnitt zeigt die Referenzimplementierung für die LAPACK-Routine dlacpy:

pip install numpy

Der n?chste Codeausschnitt zeigt dieselbe Routine, jedoch implementiert mit dem in Fortran 90 eingeführten Freiform-Layout.

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

Wie zu beobachten ist, ist moderner Fortran-Code durch die Entfernung von Spaltenbeschr?nkungen und die Abkehr von der F77-Konvention, Spezifizierer in GROSSBUCHSTABEN zu schreiben, sichtbar konsistenter und damit besser lesbar.

Beschriftete Kontrollstrukturen

Eine weitere g?ngige Praxis in LAPACK-Routinen ist die Verwendung gekennzeichneter Kontrollstrukturen. Betrachten Sie beispielsweise den folgenden Codeausschnitt, in dem die Bezeichnung 10 mit einem entsprechenden CONTINUE übereinstimmen muss.

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

Fortran 90 machte diese Vorgehensweise überflüssig und verbesserte die Lesbarkeit des Codes, indem es erlaubte, end do zu verwenden, um eine do-Schleife zu beenden. Diese ?nderung wird in der oben bereitgestellten Freiformversion von dlacpy angezeigt.

Arrays mit angenommener Gr??e

Um Flexibilit?t bei der Handhabung von Arrays unterschiedlicher Gr??e zu erm?glichen, arbeiten LAPACK-Routinen üblicherweise mit Arrays mit einer angenommenen Gr??e. In der obigen dlacpy-Routine wird die Eingabematrix A als zweidimensionales Array mit einer angenommenen Gr??e gem?? dem Ausdruck A(LDA, *) deklariert. Dieser Ausdruck deklariert, dass A eine LDA-Anzahl von Zeilen hat und verwendet * als Platzhalter, um anzugeben, dass die Gr??e der zweiten Dimension vom aufrufenden Programm bestimmt wird.

Eine Konsequenz der Verwendung von Arrays mit angenommener Gr??e ist, dass Compiler keine Grenzprüfung für die nicht angegebene Dimension durchführen k?nnen. Daher besteht die aktuelle Best Practice darin, explizite Schnittstellen und Arrays mit angenommener Form (z. B. A(LDA,:)) zu verwenden, um Speicherzugriffe au?erhalb der Grenzen zu verhindern. Dies bedeutet, dass die Verwendung von Arrays mit angenommener Form problematisch sein kann, wenn Untermatrizen an andere Funktionen übergeben werden müssen, da hierfür ein Slicing erforderlich ist, was h?ufig dazu führt, dass Compiler interne Kopien von Array-Daten erstellen.

Migration auf Fortran 95

Unn?tig zu erw?hnen, dass es eine Weile gedauert hat, bis ich mich an die LAPACK-Konventionen gew?hnt und eine LAPACK-Denkweise angenommen habe. Da ich jedoch eher ein Purist bin, wollte ich, wenn ich sowieso Routinen portieren wollte, zumindest die Routinen, die ich portieren konnte, in ein moderneres Zeitalter bringen, in der Hoffnung, die Lesbarkeit des Codes und die zukünftige Wartung zu verbessern. Nachdem ich die Dinge mit den stdlib-Betreuern besprochen hatte, entschied ich mich für die Migration von Routinen auf Fortran 95, das zwar nicht die neueste und beste Fortran-Version war, aber die richtige Balance zwischen der Beibehaltung des Erscheinungsbilds der ursprünglichen Implementierungen und der Gew?hrleistung ( gut genug) Abw?rtskompatibilit?t und Nutzung neuerer syntaktischer Funktionen.

Testabdeckung

Eines der Probleme bei der Verfolgung eines Bottom-up-Ansatzes zum Hinzufügen von LAPACK-Unterstützung besteht darin, dass explizite Unit-Tests für Dienstprogramme auf niedrigerer Ebene in LAPACK oft nicht vorhanden sind. Die Testsuite von LAPACK basiert weitgehend auf einer hierarchischen Testphilosophie, bei der davon ausgegangen wird, dass beim Testen von Routinen auf h?herer Ebene sichergestellt wird, dass ihre abh?ngigen Routinen auf niedrigerer Ebene als Teil eines Gesamtworkflows ordnungsgem?? funktionieren. Man kann zwar argumentieren, dass die Konzentration auf Integrationstests gegenüber Unit-Tests für Routinen auf niedrigerer Ebene sinnvoll ist, da das Hinzufügen von Tests für jede Routine m?glicherweise den Wartungsaufwand und die Komplexit?t des Test-Frameworks von LAPACK erh?hen k?nnte, aber das bedeutet, dass wir uns nicht ohne weiteres auf frühere Tests verlassen konnten Kunst für Unit-Tests und müssten selbst umfassende eigenst?ndige Unit-Tests für jede Routine auf niedrigerer Ebene entwickeln.

Dokumentation

?hnlich wie bei der Testabdeckung war es au?erhalb von LAPACK selbst eine Herausforderung, reale, dokumentierte Beispiele zu finden, die die Verwendung von Routinen auf niedrigerer Ebene veranschaulichen. W?hrend LAPACK-Routinen stets ein Dokumentationskommentar mit Beschreibungen der Eingabeargumente und m?glichen Rückgabewerte vorangestellt ist, kann die Visualisierung und Erfassung der erwarteten Eingabe- und Ausgabewerte ohne Codebeispiele eine Herausforderung darstellen, insbesondere beim Umgang mit speziellen Matrizen. Und obwohl das Fehlen von Unit-Tests oder dokumentierten Beispielen weder das Ende der Welt bedeutet, bedeutete es, dass das Hinzufügen von LAPACK-Unterstützung zu stdlib mühsamer sein würde, als ich erwartet hatte. Das Schreiben von Benchmarks, Tests, Beispielen und Dokumentationen würde einfach mehr Zeit und Mühe erfordern, was m?glicherweise die Anzahl der Routinen einschr?nkte, die ich w?hrend des Praktikums implementieren konnte.

Speicherlayouts

Beim Speichern von Matrixelementen im linearen Speicher hat man zwei M?glichkeiten: entweder Spalten zusammenh?ngend oder Zeilen zusammenh?ngend speichern (siehe Abbildung 2). Das erstere Speicherlayout wird als spaltengro?e-Reihenfolge und das letztere als zeilengro?e-Reihenfolge bezeichnet.

Abbildung 2: Schematische Darstellung der Speicherung von Matrixelementen im linearen Speicher entweder in (a) spaltengro?er (Fortran-Stil) oder (b) zeilengro?er (C-Stil) Reihenfolge. Die Wahl des zu verwendenden Layouts ist gr??tenteils eine Frage der Konvention.

Die Wahl des zu verwendenden Layouts ist gr??tenteils eine Frage der Konvention. Beispielsweise speichert Fortran Elemente in der Reihenfolge der Hauptspalten und C speichert Elemente in der Reihenfolge der Hauptzeilen. Bibliotheken h?herer Ebenen wie NumPy und stdlib unterstützen sowohl spalten- als auch zeilenbezogene Ordnungen, sodass Sie das Layout eines mehrdimensionalen Arrays w?hrend der Array-Erstellung konfigurieren k?nnen.

pip install numpy

W?hrend keines der Speicherlayouts von Natur aus besser ist als das andere, ist die Anordnung der Daten zur Gew?hrleistung eines sequentiellen Zugriffs gem?? den Konventionen des zugrunde liegenden Speichermodells von entscheidender Bedeutung für die Gew?hrleistung einer optimalen Leistung. Moderne CPUs sind in der Lage, sequentielle Daten effizienter zu verarbeiten als nicht sequentielle Daten, was vor allem auf das CPU-Caching zurückzuführen ist, das wiederum die r?umliche Referenzlokalit?t ausnutzt.

Um die Leistungsauswirkungen des sequentiellen vs. nicht-sequentiellen Elementzugriffs zu demonstrieren, betrachten Sie die folgende Funktion, die alle Elemente aus einer MxN-Matrix A in eine andere MxN-Matrix B kopiert und dabei davon ausgeht, dass Matrixelemente im Hauptspaltenformat gespeichert sind bestellen.

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

Seien A und B die folgenden 3x2-Matrizen:

Wenn sowohl A als auch B in der Hauptspaltenreihenfolge gespeichert sind, k?nnen wir die Kopierroutine wie folgt aufrufen:

pip install numpy

Wenn A und B jedoch beide in Zeilenreihenfolge gespeichert sind, ?ndert sich die Anrufsignatur zu

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

Beachten Sie, dass wir im letzteren Szenario nicht in der innersten Schleife auf Elemente in sequentieller Reihenfolge zugreifen k?nnen, da da0 2 und da1 -5 ist und das Gleiche gilt für db0 und db1. Stattdessen springen die ?Zeiger“ des Array-Index wiederholt weiter, bevor sie zu früheren Elementen im linearen Speicher zurückkehren, wobei ia = {0, 2, 4, 1, 3, 5} und ib gleich sind. In Abbildung 3 zeigen wir die Auswirkungen des nicht sequentiellen Zugriffs auf die Leistung.

Abbildung 3: Leistungsvergleich bei der Bereitstellung quadratischer Spalten-Haupt- und Zeilen-Hauptmatrizen zum Kopieren, wenn Kopieren einen sequentiellen Elementzugriff entsprechend der Spalten-Hauptreihenfolge voraussetzt. Die x-Achse z?hlt zunehmende Matrixgr??en (d. h. Anzahl der Elemente) auf. Alle Raten werden im Verh?ltnis zu den Hauptergebnissen der Spalte für eine entsprechende Matrixgr??e normalisiert.

Aus der Abbildung k?nnen wir ersehen, dass die spalten- und zeilenbezogene Leistung ungef?hr gleich ist, bis wir mit quadratischen Matrizen mit mehr als 1e5 Elementen arbeiten (M = N = ~316). Bei 1e6-Elementen (M = N = ~1000) führt die Bereitstellung einer zu kopierenden Zeilenhauptmatrix zu einem Leistungsabfall von mehr als 25 %. Für 1e7-Elemente (M = N = ~3160) beobachten wir einen Leistungsabfall von mehr als 85 %. Die erhebliche Auswirkung auf die Leistung kann auf eine verringerte Referenzlokalit?t bei der Arbeit mit zeilengro?en Matrizen mit gro?en Zeilengr??en zurückgeführt werden.

Da es in Fortran geschrieben ist, geht LAPACK von einer spaltenorientierten Zugriffsreihenfolge aus und implementiert seine Algorithmen entsprechend. Dies führt zu Problemen für Bibliotheken wie stdlib, die nicht nur die Zeilenreihenfolge unterstützen, sondern diese auch zu ihrem Standardspeicherlayout machen. Würden wir die Fortran-Implementierungen von LAPACK einfach auf JavaScript portieren, würden Benutzer, die zeilenorientierte Matrizen bereitstellen, negative Auswirkungen auf die Leistung aufgrund des nicht sequentiellen Zugriffs erfahren.

Um negative Auswirkungen auf die Leistung abzumildern, haben wir eine Idee von BLIS übernommen, einer BLAS-?hnlichen Bibliothek, die sowohl zeilen- als auch spaltenorientierte Speicherlayouts in BLAS-Routinen unterstützt, und beschlossen, bei der Portierung von Routinen von Fortran nach JavaScript modifizierte LAPACK-Implementierungen zu erstellen C, die explizit sowohl spalten- als auch zeilenorientierte Speicherlayouts durch separate Schrittparameter für jede Dimension berücksichtigen. Bei einigen Implementierungen, wie z. B. dlacpy, das der oben definierten Kopierfunktion ?hnelt, ist die Einbindung separater und unabh?ngiger Schritte unkompliziert und erfordert h?ufig Schritttricks und Schleifenaustausch. Bei anderen erwiesen sich die ?nderungen jedoch aufgrund von als viel weniger einfach spezielle Matrixbehandlung, unterschiedliche Zugriffsmuster und kombinatorische Parametrisierung.

ndarrays

LAPACK-Routinen arbeiten haupts?chlich mit Matrizen, die im linearen Speicher gespeichert sind und auf deren Elemente gem?? den angegebenen Dimensionen und dem Schritt der führenden (d. h. ersten) Dimension zugegriffen wird. Dimensionen geben die Anzahl der Elemente in jeder Zeile bzw. Spalte an. Der Schritt gibt an, wie viele Elemente im linearen Speicher übersprungen werden müssen, um auf das n?chste Element einer Zeile zuzugreifen. LAPACK geht davon aus, dass Elemente, die zur gleichen Spalte geh?ren, immer zusammenh?ngend sind (d. h. benachbart im linearen Speicher). Abbildung 4 bietet eine visuelle Darstellung der LAPACK-Konventionen (insbesondere die Schaltpl?ne (a) und (b)).

Abbildung 4: Schematische Darstellung der Verallgemeinerung der LAPACK-Strided-Array-Konventionen auf nicht zusammenh?ngende Strided-Arrays. a) Eine zusammenh?ngende 5-mal-5-Matrix, die in Spaltenhauptreihenfolge gespeichert ist. b) Eine nicht zusammenh?ngende 3-mal-3-Untermatrix, die in der Hauptspaltenreihenfolge gespeichert ist. Untermatrizen k?nnen in LAPACK bearbeitet werden, indem ein Zeiger auf das erste indizierte Element bereitgestellt und der Schritt der führenden (d. h. ersten) Dimension angegeben wird. In diesem Fall betr?gt der Schritt der führenden Dimension fünf, obwohl nur drei Elemente pro Spalte vorhanden sind, da die Untermatrixelemente im linearen Speicher nicht zusammenh?ngend sind, wenn sie als Teil einer gr??eren Matrix gespeichert werden. In LAPACK wird der Schritt der nachgestellten (d. h. zweiten) Dimension immer als Eins angenommen. c) Eine nicht zusammenh?ngende 3-mal-3-Untermatrix, die in der Hauptspaltenreihenfolge gespeichert ist und keine Einheitsschritte aufweist und die LAPACK-Schrittkonventionen sowohl auf führende als auch auf nachfolgende Dimensionen verallgemeinert. Diese Verallgemeinerung liegt den mehrdimensionalen Arrays von stdlib (auch als ?ndarrays“ bezeichnet) zugrunde.

Bibliotheken wie NumPy und stdlib verallgemeinern die Strided-Array-Konventionen von LAPACK zur Unterstützung

1. Nicht-Einheitsschritte in der letzten Dimension (siehe Abbildung 4 (c)). LAPACK geht davon aus, dass die letzte Dimension einer Matrix immer einen Einheitsschritt hat (d. h. Elemente innerhalb einer Spalte werden zusammenh?ngend im linearen Speicher gespeichert).
2. Negative Fortschritte für jede Dimension. LAPACK erfordert, dass der Schritt einer führenden Matrixdimension positiv ist.
3. Mehrdimensionale Arrays mit mehr als zwei Dimensionen. LAPACK unterstützt explizit nur Strided-Vektoren und (Sub-)Matrizen.

Die Unterstützung von Nicht-Einheitsschritten in der letzten Dimension gew?hrleistet die Unterstützung der O(1)-Erstellung nicht zusammenh?ngender Ansichten des linearen Speichers, ohne dass eine explizite Datenverschiebung erforderlich ist. Diese Ansichten werden oft als ?Slices“ bezeichnet. Betrachten Sie als Beispiel den folgenden Codeausschnitt, der solche Ansichten mithilfe der von stdlib bereitgestellten APIs erstellt.

pip install numpy

Ohne Unterstützung für Nicht-Einheitsschritte in der letzten Dimension w?re die Rückgabe einer Ansicht aus dem Ausdruck x['::2,::2'] nicht m?glich, da ausgew?hlte Elemente in eine neue Linearit?t kopiert werden müssten Speicherpuffer, um die Kontiguit?t sicherzustellen.

Abbildung 5: Schematische Darstellung der Verwendung der Schrittmanipulation zur Erstellung gespiegelter und gedrehter Ansichten von Matrixelementen, die im linearen Speicher gespeichert sind. Für alle Unterschemata werden die Schritte als [Trailing_Dimension, Leading_Dimension] aufgeführt. Implizit für jeden Schaltplan ist ein ?Offset“, der den Index des ersten indizierten Elements angibt, sodass für eine Matrix A das Element A_ij ist aufgel?st gem?? i?strides[1] j?strides[0] Offset. a) Bei einer 3x3-Matrix, die in Spaltenhauptreihenfolge gespeichert ist, kann man die Schritte der führenden und nachfolgenden Dimensionen manipulieren, um Ansichten zu erstellen, in denen auf Matrixelemente entlang einer oder mehrerer Achsen in umgekehrter Reihenfolge zugegriffen wird. b) Mit einer ?hnlichen Schrittmanipulation k?nnen gedrehte Ansichten von Matrixelementen relativ zu ihrer Anordnung im linearen Speicher erstellt werden.

Die Unterstützung negativer Schritte erm?glicht die O(1)-Umkehr und Drehung von Elementen entlang einer oder mehrerer Dimensionen (siehe Abbildung 5). Um beispielsweise eine Matrix von oben nach unten und von links nach rechts umzudrehen, muss man nur die Schritte negieren. Aufbauend auf dem vorherigen Codeausschnitt demonstriert der folgende Codeausschnitt die Umkehrung von Elementen um eine oder mehrere Achsen.

pip install numpy

Die Diskussion über negative Schritte impliziert die Notwendigkeit eines ?Offset“-Parameters, der den Index des ersten indizierten Elements im linearen Speicher angibt. Für ein mehrdimensionales Array mit Schritten A und eine Liste von Schritten s der Index, der dem Element A_ij???n entspricht kann gem?? der Gleichung gel?st werden

wobei N die Anzahl der Array-Dimensionen ist und s_k dem k-ten Schritt entspricht.

In BLAS- und LAPACK-Routinen, die negative Schritte unterstützen – etwas, das nur bei der Arbeit mit Schrittvektoren unterstützt wird (siehe z. B. Daxpy oben) – wird der Index-Offset mithilfe einer Logik berechnet, die dem folgenden Codeausschnitt ?hnelt:

pip install numpy

wobei M die Anzahl der Vektorelemente ist. Dies setzt implizit voraus, dass ein bereitgestellter Datenzeiger auf den Anfang des linearen Speichers für einen Vektor zeigt. In Sprachen, die Zeiger unterstützen, wie z. B. C, passt man zum Bearbeiten eines anderen Bereichs des linearen Speichers typischerweise einen Zeiger mithilfe der Zeigerarithmetik vor dem Funktionsaufruf an, was zumindest für den eindimensionalen Fall relativ kostengünstig und unkompliziert ist.

Wenn wir beispielsweise zu c_daxpy wie oben definiert zurückkehren, k?nnen wir Zeigerarithmetik verwenden, um den Elementzugriff auf fünf Elemente im linearen Speicher zu beschr?nken, beginnend beim elften und sechzehnten Element (Hinweis: nullbasierte Indizierung) eines Eingabe- und Ausgabearrays. bzw. wie im folgenden Codeausschnitt gezeigt.

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

In JavaScript, das keine explizite Zeigerarithmetik für bin?re Puffer unterstützt, muss man jedoch explizit neue typisierte Array-Objekte mit einem gewünschten Byte-Offset instanziieren. Um im folgenden Codeausschnitt die gleichen Ergebnisse wie im obigen C-Beispiel zu erzielen, müssen wir einen typisierten Array-Konstruktor aufl?sen, einen neuen Byte-Offset berechnen, eine neue typisierte Array-L?nge berechnen und eine neue typisierte Array-Instanz erstellen.

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

Bei gro?en Array-Gr??en sind die Kosten für die Instanziierung typisierter Arrays im Vergleich zu der Zeit, die für den Zugriff auf und die Bearbeitung einzelner Array-Elemente aufgewendet wird, vernachl?ssigbar; Bei kleineren Array-Gr??en kann die Objektinstanziierung jedoch erhebliche Auswirkungen auf die Leistung haben.

Um negative Auswirkungen auf die Objektinstanziierung zu vermeiden, entkoppelt stdlib den Datenpuffer eines Ndarrays von der Position des Pufferelements, das dem Anfang einer Ndarray-Ansicht entspricht. Dadurch k?nnen die Slice-Ausdrücke x[2:,3:] und x[3:,1:] neue ndarray-Ansichten zurückgeben, ohne dass neue Pufferinstanzen instanziiert werden müssen, wie im folgenden Codeausschnitt gezeigt.

// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';

// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;

// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');

// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);

Als Konsequenz der Entkopplung eines Datenpuffers vom Anfang einer Ndarray-Ansicht wollten wir ebenfalls vermeiden, dass beim Aufruf von LAPACK-Routinen mit Ndarray-Daten neue typisierte Array-Instanzen instanziiert werden müssen. Dies bedeutete die Erstellung modifizierter LAPACK-API-Signaturen, die explizite Offset-Parameter für alle Schrittvektoren und Matrizen unterstützen.

Der Einfachheit halber kehren wir zur JavaScript-Implementierung von daxpy zurück, die zuvor oben definiert wurde.

pip install numpy

Wie im folgenden Codeausschnitt gezeigt, k?nnen wir die obige Signatur und Implementierung so ?ndern, dass die Verantwortung für die Aufl?sung des ersten indizierten Elements auf den API-Konsumenten verlagert wird.

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

Bei ndarrays erfolgt die Aufl?sung w?hrend der ndarray-Instanziierung, sodass der Aufruf von daxpy_ndarray mit ndarray-Daten eine einfache übergabe der zugeh?rigen ndarray-Metadaten darstellt. Dies wird im folgenden Codeausschnitt demonstriert.

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

?hnlich wie bei BLIS sahen wir sowohl in herk?mmlichen LAPACK-API-Signaturen (z. B. für Abw?rtskompatibilit?t) als auch in modifizierten API-Signaturen (z. B. zur Minimierung nachteiliger Auswirkungen auf die Leistung) einen Wert und entschieden uns daher für einen Plan, sowohl konventionelle als auch andere API-Signaturen bereitzustellen ge?nderte APIs für jede LAPACK-Routine. Um die Codeduplizierung zu minimieren, wollten wir eine gemeinsame ?Basis“-Implementierung auf niedrigerer Ebene implementieren, die dann von APIs auf h?herer Ebene umschlossen werden k?nnte. W?hrend die oben gezeigten ?nderungen für die BLAS-Routine daxpy relativ einfach erscheinen m?gen, war die Transformation einer herk?mmlichen LAPACK-Routine und ihres erwarteten Verhaltens in eine generalisierte Implementierung oft viel weniger einfach.

dlaswp

Genug mit den Herausforderungen! Wie sieht ein Endprodukt aus?!

Lassen Sie uns den Kreis schlie?en und dies zurück zu dlaswp bringen, einer LAPACK-Routine zum Durchführen einer Reihe von Zeilenaustauschen in einer Eingabematrix gem?? einer Liste von Pivot-Indizes. Der folgende Codeausschnitt zeigt die Referenzimplementierung von LAPACK Fortran.

// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';

// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;

// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');

// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);

Um die Schnittstelle mit der Fortran-Implementierung von C aus zu erleichtern, stellt LAPACK eine zweistufige C-Schnittstelle namens LAPACKE bereit, die Fortran-Implementierungen umschlie?t und sowohl zeilen- als auch spaltenbezogene Eingabe- und Ausgabematrizen unterstützt. Die mittlere Schnittstelle für dlaswp wird im folgenden Codeausschnitt gezeigt.

void c_daxpy(const int N, const double alpha, const double *X, const int strideX, double *Y, const int strideY) {
    int ix;
    int iy;
    int i;
    if (N <= 0) {
        return;
    }
    if (alpha == 0.0) {
        return;
    }
    if (strideX < 0) {
        ix = (1-N) * strideX;
    } else {
        ix = 0;
    }
    if (strideY < 0) {
        iy = (1-N) * strideY;
    } else {
        iy = 0;
    }
    for (i = 0; i < N; i++) {
        Y[iy] += alpha * X[ix];
        ix += strideX;
        iy += strideY;
    }
    return;
}

Beim Aufruf mit einer Spaltenhauptmatrix a übergibt der Wrapper LAPACKE_dlaswp_work einfach die bereitgestellten Argumente an die Fortran-Implementierung. Wenn der Wrapper jedoch mit einer Zeilenhauptmatrix a aufgerufen wird, muss er Speicher zuweisen, a explizit in eine tempor?re Matrix a_t transponieren und kopieren, den Schritt der führenden Dimension neu berechnen, dlaswp mit a_t aufrufen, die in a_t gespeicherten Ergebnisse transponieren und kopieren zu einem und schlie?lich frei zugewiesenen Speicher. Das ist ziemlich viel Arbeit und kommt bei den meisten LAPACK-Routinen vor.

Der folgende Codeausschnitt zeigt die auf JavaScript portierte Referenz-LAPACK-Implementierung mit Unterstützung für führende und nachgestellte Dimensionsschritte, Indexoffsets und einen Schrittvektor, der Pivot-Indizes enth?lt.

const binary = new UintArray([...]);

const mod = new WebAssembly.Module(binary);

Um eine API mit konsistentem Verhalten mit herk?mmlichem LAPACK bereitzustellen, habe ich dann die obige Implementierung verpackt und Eingabeargumente an die ?Basis“-Implementierung angepasst, wie im folgenden Codeausschnitt gezeigt.

pip install numpy

Ich habe anschlie?end einen separaten, aber ?hnlichen Wrapper geschrieben, der eine API-Zuordnung direkter zu den mehrdimensionalen Arrays von stdlib bereitstellt und eine spezielle Behandlung durchführt, wenn die Richtung, in die Pivots angewendet werden sollen, negativ ist, wie im folgenden Codeausschnitt gezeigt.

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

Ein paar Punkte, die Sie beachten sollten:

1. Im Gegensatz zur herk?mmlichen LAPACKE-API ist der Parameter ?matrix_layout“ (Reihenfolge) in den dlaswp_ndarray- und Basis-APIs nicht erforderlich, da die Reihenfolge aus den bereitgestellten Schritten abgeleitet werden kann.
2. Im Gegensatz zur herk?mmlichen LAPACKE-API müssen wir bei einer zeilenmajoren Eingabematrix keine Daten in tempor?re Arbeitsbereich-Arrays kopieren, wodurch unn?tige Speicherzuweisung reduziert wird.
3. Im Gegensatz zu Bibliotheken wie NumPy und SciPy, die direkt mit BLAS und LAPACK verbunden sind, müssen wir beim Aufrufen von LAPACK-Routinen in stdlib vorher nicht zusammenh?ngende mehrdimensionale Daten in und aus tempor?ren Arbeitsbereichs-Arrays kopieren bzw. nach dem Aufruf. Au?er bei der Schnittstelle zu hardwareoptimiertem BLAS und LAPACK tr?gt der verfolgte Ansatz dazu bei, Datenbewegungen zu minimieren und die Leistung in ressourcenbeschr?nkten Browseranwendungen sicherzustellen.

Für serverseitige Anwendungen, die hardwareoptimierte Bibliotheken wie OpenBLAS nutzen m?chten, stellen wir separate Wrapper bereit, die verallgemeinerte Signaturargumente an ihre optimierten API-?quivalente anpassen. In diesem Zusammenhang kann sich zumindest bei ausreichend gro?en Arrays die Erstellung tempor?rer Kopien lohnen.

Aktueller Status und n?chste Schritte

Trotz der Herausforderungen, unvorhergesehenen Rückschl?ge und mehreren Design-Iterationen freue ich mich, berichten zu k?nnen, dass ich zus?tzlich zu dlaswp oben 35 PRs ?ffnen konnte, die Unterstützung für verschiedene LAPACK-Routinen und zugeh?rige Dienstprogramme hinzufügen. Natürlich sind es nicht ganz 1.700 übungen, aber ein guter Anfang! :)

Dennoch ist die Zukunft rosig und wir sind sehr gespannt auf diese Arbeit. Es gibt noch viel Raum für Verbesserungen und zus?tzliche Forschung und Entwicklung. Insbesondere liegt uns daran

1. Entdecken Sie Werkzeuge und Automatisierung.
2. Beheben Sie Build-Probleme beim Aufl?sen der Quelldateien von Fortran-Abh?ngigkeiten, die über mehrere stdlib-Pakete verteilt sind.
3. Rollout von C- und Fortran-Implementierungen und nativen Bindungen für die vorhandenen LAPACK-Pakete von stdlib.
4. die Bibliothek modularer LAPACK-Routinen von stdlib weiter ausbauen.
5. Identifizieren Sie zus?tzliche Bereiche zur Leistungsoptimierung.

W?hrend mein Praktikum bei Quansight Labs beendet ist, habe ich vor, weiterhin Pakete hinzuzufügen und diese Bemühungen voranzutreiben. Angesichts des immensen Potenzials und der grundlegenden Bedeutung von LAPACK würden wir uns freuen, wenn diese Initiative, LAPACK ins Internet zu bringen, weiter w?chst. Wenn Sie also daran interessiert sind, dies voranzutreiben, z?gern Sie bitte nicht, uns zu kontaktieren! Und wenn Sie daran interessiert sind, die Entwicklung zu sponsern, stehen Ihnen die Leute von Quansight gerne für ein Gespr?ch zur Verfügung.

Und damit m?chte ich Quansight für die Bereitstellung dieser Praktikumsm?glichkeit danken. Ich bin unglaublich glücklich, so viel gelernt zu haben. Ein Praktikant bei Quansight zu sein war schon lange ein Traum von mir und ich bin sehr dankbar, dass ich ihn erfüllt habe. Mein besonderer Dank gilt Athan Reines und Melissa Mendon?a, die eine gro?artige Mentorin und rundum wundervolle Person ist! Und vielen Dank an alle stdlib-Kernentwickler und alle anderen bei Quansight, die mir auf diesem Weg sowohl im Gro?en als auch im Kleinen geholfen haben.

Prost!

---

stdlib ist ein Open-Source-Softwareprojekt, das sich der Bereitstellung einer umfassenden Suite robuster, leistungsstarker Bibliotheken widmet, um die Entwicklung Ihres Projekts zu beschleunigen und Ihnen die Gewissheit zu geben, dass Sie auf fachm?nnisch erstellte, qualitativ hochwertige Software angewiesen sind.

Wenn Ihnen dieser Beitrag gefallen hat, geben Sie uns einen Stern? auf GitHub und überlegen Sie, das Projekt finanziell zu unterstützen. Ihre Beitr?ge und Ihre kontinuierliche Unterstützung tragen dazu bei, den langfristigen Erfolg des Projekts sicherzustellen und werden sehr gesch?tzt!

Das obige ist der detaillierte Inhalt vonLAPACK in Ihrem Webbrowser. Für weitere Informationen folgen Sie bitte anderen verwandten Artikeln auf der PHP chinesischen Website!