BLAS - Basic Linear Algebra Subroutines

 

NAME

BLAS - Basic Linear Algebra Subroutines

BESCHREIBUNG

Die Basic Linear Algebra Subprograms (BLAS) bestehen aus drei Gruppen.

Level 1 BLAS

Diese Routinen werden in <1> beschrieben und umfassen Vektoroperationen, z.B. entspricht die Routine SAXPY einer 'linked triad' V1 = V1 + s * V2 und SDOT dient zur Berechnung des Skalarprodukts zweier Vektoren. Die BLAS-1 Routinen sind für den Gebrauch auf Skalarrechnern gedacht, wo z.T. optimierte Versionen existieren, die z.B. im jeweiligen Assembler geschrieben sind. Auf Vektorrechnern können sie i.a. durch einen bzw. einige Vektorbefehle realisiert werden. Um den zusätzlichen Aufwand für einen Unterprogrammaufruf zu vermeiden, sollten die BLAS-1 Routinen deshalb auf
Vektorrechnern möglichst im Quelltext expandiert werden.

Level-2 BLAS

Die BLAS-2 Routinen wurden definiert, um eine Möglichkeit zu schaffen, portable Software zu entwickeln, die auf einer Vielzahl von Vektorrechnern optimal abläuft. Die rechenzeitintensiven Programmteile (jedenfalls für den umfangreichen Bereich der Linearen Algebra) sind in wenigen, optimierten Rou tinen isoliert. BLAS-2 ist in <2> beschrieben. Die BLAS-2 Programme implementieren Matrix-Vektor-Operationen, also typischerweise O(n**2) arithmetische Operationen.

Level-3 BLAS

Für neuere Rechnerarchitekturen (z.B. mehrere Prozessoren, hierarchische Speicherstrukturen) sind die BLAS-2 Routinen nur bedingt geeignet. Deshalb wird in <3> ein Satz von Level-3 BLAS Programmen vorgeschlagen, die Matrix-Matrix-Operationen, also O(n**3) arithmetische Operationen umfassen.
BLAS-3 trägt auch zur besseren Nutzung der Vektorregister bei. Basierend auf BLAS-3 wird ein Programmpaket entwickelt, das den Leistungsumfang der bekannten Pakete EISPACK und LINPACK hat, also den Bereich der Linearen Algebra für dicht gespeicherte Matrizen abdeckt <4>.

BLAS-2 Funktionsumfang

Die folgenden drei Arten von Matrix-Vektor-Operationen werden von Level 2 BLAS durchgeführt.

a) Matrix-Vektor-Multiplikationen der Form:

y <- alpha A x + beta y, t y <- alpha A x + beta y und
_t y <- alpha A x + beta wobei A eine Matrix, alpha und beta Skalare und x und y Vektoren sind, t _t (A bezeichnet die zu A transponierte Matrix und A die transponierte und konjugiert komplexe Matrix zu A.) x <- T x , t x <- T x und _t x <- T x ,

wobei x ein Vektor und T eine obere oder untere Dreiecksmatrix ist.

b) Rank-one und rank-two updates der Form

t A <- alpha x y + A , _t A <- alpha x y + A , _t H <- alpha x x + H und _t _____ _t H <- alpha x y + alpha y x + H ,

wobei H eine Hermitesche Matrix ist.

c) Lösung von Gleichungssystemen mit Dreiecksmatrix als Koeffizientenmatrix

-1 x <- T x , t -1 x <- ( T ) x und _t -1 x <- ( T ) x ,

wobei T eine nicht-singuläre obere oder untere Dreiecksmatrix ist.

Diese Operationen werden soweit sinnvoll auf allgemeine rechteckige Matrizen, Bandmatrizen, Hermitesche Matrizen, Hermitesche Bandmatrizen, Dreiecksmatrizen und Dreiecksbandmatrizen in reeller und komplexer Arithmetik sowie in einfacher und doppelter Genauigkeit angewendet.

Namenskonventionen

Die Namen der Routinen bestehen aus max. 5 Zeichen mit folgender Bedeutung.

Das vierte und fünfte Zeichen beschreiben die Operation:

MV

- Matrix-vector product

R- Rank-one update
R2- Rank-two update
SV- Solve a system of linear equations

Die Zeichen zwei und drei bezeichnen die Art der Matrix:

GEGeneral matrix
GBGeneral band matrix
HEHermitian matrix
SYSymmetric matrix
HPHermitian matrix stored in packed form
SPSymmetric matrix stored in packed form
HBHermitian band matrix
SBSymmetric band matrix
TRTriangular matrix
TPTriangular matrix in packed form
TBTriangular band matrix

Das erste Zeichen beschreibt den FORTRAN Datentyp:

S
REAL
DDOUBLE PRECISION
CCOMPLEX
ZCOMPLEX*16

Alle Kombinationen sind in der folgenden Tabelle angegeben, wobei in der ersten Spalte noch C durch Z und in der zweiten Spalte S durch D ersetzt werden kann.

Datentyp
 
Operation
complexrealMVRR2SV
CGESGE**  
CGBSGB*   
CHESSY*** 
CHPSSP*** 
CHBSSB*   
CTRSTR**  
CTPSTP**  
CTBSTB**  

Für die Operation 'rank-1 update' für allgemeine Matrizen (GER) gibt es zwei komplexe Versionen:

_t CGERC für A <- alpha x y + A und t CGERU für A <- alpha x y + A .

Insgesamt umfaßt BLAS-2 also 66 Routinen, die 16 verschiedene Operationen enthalten, die i.a. sowohl auf die übergebene Matrix als auch auf ihre Transponierte bzw. für komplexe Matrizen auch auf die transponierte und konjugiert komplexe Matrix angewendet werden können.

BLAS-3 Funktionsumfang

Folgende vier Arten von Matrix-Vektoroperationen werden von Level-3 BLAS durchgeführt.


a) Matrix-Matrix-Multiplikationen der Form:

C <- alpha A B + beta C , t C <- alpha A B + beta C , t C <- alpha A B + beta C und t t C <- alpha A B + beta C , wobei A, B und C Matrizen und alpha und beta Skalare sind.

b) Rank-k und rank-2k updates für symmetrische Matrizen: t C <- alpha A A + beta C , t C <- alpha A A + beta C , t t C <- alpha A A + beta C , t t C <- alpha A B + alpha B A + beta C , t t C <- alpha A B + alpha B A + beta Cc) Multiplikation einer Matrix mit einer Dreiecksmatrix B <- alpha T B, t B <- alpha T B, B <- alpha B T und t B <- alpha B Td) Lösung von mehreren Gleichungssystemen mit einer Dreiecksmatrix als Koeffizientenmatrix -1 B <- alpha T B, t -1 B <- alpha ( T ) B, -1 B <- alpha B T und t -1 B <- alpha B ( T )

wobei T eine nicht-singuläre obere oder untere Dreiecksmatrix ist.

Diese Operationen werden soweit sinnvoll auf allgemeine rechteckige Matrizen, Hermitesche Matrizen und Dreiecksmatrizen in reeller und komplexer Arithmetik sowie in einfacher und doppelter Genauigkeit angewendet. Im Falle komplexer Matrizen kann für AT (Transposition) auch AH (transponiert und konjugiert komplex) stehen, so daß in diesem Falle nicht vier sondern neun verschiedene Matrix-Matrix-Multiplikationen vorgesehen sind.

Namenskonventionen

Die Namen der Routinen bestehen aus max. 6 Zeichen mit folgender Bedeutung.

Das vierte bis sechste Zeichen beschreiben die Operation:

MM- Matrix-matrix product
RK- Rank-k update of a symmetric or Hermitian matrix
R2K- Rank-2k update of a symmetric or Hermitian matrix
SM- Solve a system of linear equations for a matrix of right- hand sides

Die Zeichen zwei und drei bezeichnen die Art der Matrix:

GEGeneral matrix
HEHermitian matrix
SYSymmetric matrix
TRTriangular matrix

 Das erste Zeichen beschreibt den FORTRAN Datentyp:

SREAL
DDOUBLE PRECISION
CCOMPLEX
ZCOMPLEX*16

Alle Kombinationen sind in der folgenden Tabelle angegeben, wobei in der ersten Spalte noch C durch Z und in der zweiten Spalte S durch D ersetzt werden kann.

Datentyp
 
Operation
complexrealMMRKR2KSM
CGESGE*   
CSYSSY*** 
CHE *** 
CTRSTR**  

 

BLAS-3 besteht also aus 6 reellen Routinen, die 48 verschiedene Kombinationen von Optionen erlauben und 9 komplexen Routinen mit 77 verschiedenen Kombinationsmöglichkeiten der Optionen.

BLAS Routinen am Rechen- und Kommunikationszentrum der RWTH Aachen

Für die Siemens Vektorrechner wurden optimierte Versionen der BLAS-2 und BLAS-3 Programme erarbeitet, die die Vektorregister des VP sowie die Parallelität der Pipelines optimal nutzen. Diese Programme sollten als Bausteine für die Entwicklung von Anwendungsprogrammen benutzt werden und werden auch bereits in den Unterprogrammbibliotheken IMSL und NAG eingesetzt.

Am Siemens Vektorrechner sind diese Programme in der Datei /usr/local/lib/libblasvp.a gespeichert. Bei der Übersetzung muß diese Bibliothek mit der Option -lblasvp angebunden werden.

Bei der Benutzung der BLAS Routinen und auch anderer Unterprogramme, in denen Matrizen bearbeitet werden, ist die führende Dimension der Felder i.a. von besonderer Bedeutung. Bei den BLAS-Routinen sind dies der Parameter LDA, LDB oder LDC.

Dieser Parameter gibt an, mit wievielen Zeilen eine Matrix in der DIMENSION- Anweisung vereinbart wurde. Diese Größe sollte am VP für den Datentyp REAL das Doppelte einer ungeraden Zahl sein, d.h. LDA = 2 * (2 * n + 1), und für die übrigen Datentypen sollte LDA eine ungerade Zahl sein. Nur dann kann i.a. ein halbwegs günstiger Zugriff auf eine Zeile einer Matrix erfolgen. Diese Bedingung gilt übrigens auch für alle anderen Programme, in denen mehrdimensionale Felder verarbeitet werden, und auch für andere Vektorrechner.Beschreibungen der einzelnen Routinen findet man an der IBM 3090 in der Datei SYSRZ.BLAS.DOK. In der Datei SYSRZ.BLAS2.FORT sind portable Versionen der BLAS 2 und in der Datei SYSRZ.BLAS3.FORT portable Versionen der BLAS 3 Programme im Quelltext gespeichert. Diese Programme können z. B. auf Institutsrechner übertragen und dort eingesetzt werden.

Beispiel

Berechnung der L-R Zerlegung einer quadratischen Matrix

SUBROUTINE LR(N,EPS,A,LDA,P,IFAIL) C C .. Scalar Arguments .. DOUBLE PRECISION EPS INTEGER IFAIL,LDA,N C .. C .. Array Arguments .. DOUBLE PRECISION A(LDA,*) INTEGER P(*) C .. C .. Local Scalars .. DOUBLE PRECISION S, Y INTEGER I,IMAX,IP,J C .. C .. External Subroutines .. EXTERNAL DGER C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. DO 10 I = 1,N P(I) = I 10 CONTINUE C DO 190 I = 1,N C Pivotzeile bestimmen S = ABS(A(I,I)) IMAX = I DO 90 J = I + 1,N IF (ABS(A(J,I)).GT.S) THEN S = ABS(A(J,I)) IMAX = J END IF 90 CONTINUE IF (S.LT.EPS) GO TO 400 IF (IMAX.NE.I) THEN C Zeilen vertauschen DO 100 J = 1,N Y = A(I,J) A(I,J) = A(IMAX,J) A(IMAX,J) = Y 100 CONTINUE IP = P(I) P(I) = P(IMAX) P(IMAX) = IP END IF C S = 1.D0/A(I,I) DO 111 J = I + 1,N A(I,J) = A(I,J)*S 111 CONTINUE C CALL DGER(N-I,N-I,-1.0D0,A(I+1,I),1,A(I,I+1),LDA,A(I+1,I+1), + LDA) 190 CONTINUE C IFAIL = 0 C RETURN 400 IFAIL = 1 RETURN END

Literatur

  1. C.Lawson, R. Hanson, D, Kincaid, F. Krogh, C.Lawson, R. Hanson, D, Kincaid, F. Krogh, C.Lawson, R. Hanson, D, Kincaid, F. Krogh, C.Lawson, R. Hanson, D, Kincaid, F. Krogh, "Basic Linear Algebra Subprograms for Fortran Usage
  2. J. Dongarra, J. Du Croz, S. Hammarling, R. Hanson, "An Extended Set of Fortran Basic Linear Algebra Subprograms", NAG Technical Report TR3/87, February 1987
  3. J. Dongarra, J. Du Droz, I. Duff, S. Hammarling, "A Set of Level 3 Basic Linear Algebra Subprograms", Argonne National Laboratory, Mathematics and Computer Science Division, Technical Memorandum No. 88 (Revision 1)
  4. J. Demmell, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, D. Sorensen, "Prospectus for the Development of a Linear Algebra Library for High-Performance Computers", ANL/MCS-TM-97

Benutzung der Software

Aufruf:

  •  IBM:
xlf/xlf90 -lessl ... bzw. optimiert für die Power2-Architektur xlf/xlf90 -lesslp2 ...
  •  SGI:
f77/f90 -n32 -mips4 -L/usr/lib32 -lblas ... bzw. f77/f90 -64 -mips4 -L/usr/lib64 -lblas ...
  •  Fujitsu:
frt -L/usr/local_rwth/lib -lblas
  •  HP:
f77/f90 -lblas ...
  •  Sun:
f77/f90 -xlic_lib=sunperf
  •  Linux:
g77 -L/usr/local_rwth/lib -lblas

Hilfe:

man blas bzw. man <blas-routine>
Ein Anwendungsbeispiel inkl. Rechenzeitmessungen finden Sie hier.
Lizensierung

In der Regel wird eine optimierte Blas-Library vom Workstation-Lieferanten mitgeliefert.
Die (nicht optimierten) Programmquellen sind public domain hier verfügbar.
Ansprechstelle

Beratung:

Dieter an Mey

Tel.: +49 241 80 24377

mailto: hpc@rz.rwth-aachen.de

Stand: 19.05.03


Abschlußinformationen

  • Keine Stichwörter