Vad är det som gör så BLAS biblioteket för matemaik är så optimerat?
Jag har alltid undrat varför BLAS är så optimerat? Det finns olika nivåer av BLAS. Nivå 1 som handlar om vektorer * vektorer. Nivå 2 som handlar om matriser * vektorer. Nivå 3 som handlar om matriser * matriser. Som jag har tolkat det som.
BLAS är alltså Basic Linear Algebra Subprograms, dvs multiplikation, subbstraktion, addition, division, invertering, transponat osv. Alltså absolut lägsta nivån inom den linjära algebran. LAPACK finns också, men då handlar det mer om LU, SVD, EIG, QR osv.
När jag läser koden så verkar det inget annat vara än C-kod tidgare än C89-standard, dsv K&R C-standard från forna tider som C-folket körde innan det ens fanns en officiel standard. Koden ser ut som man har rullat en apelsin på ett tangentbord och inget "make sense". Men koden fungerar än idag och körs aktivt inom stora bibliotek. Senaste updatering hos BLAS för C-kod var 2017. Alltså 2 år sedan. Inte illa, med tanke på att BLAS härstammar från 70-talet.
Men då är frågorna:
1. Vad är det som gör BLAS så effektivt?
2. Varför är inte denna kod effektiv? För beräkning av A*x = b, där A och x är kända men inte b.
// Enligt C99 så fungerar det att ha 2D arrayer som argument, så länge dimensionerna nämns före som argument
void mat(int n, int m, A[n][m], x[n], b[n]){
for(int i = 0; i < n; i++){
for(int j = 0; j < m; j++){
b[i] = A[i][j]*x[i]; // Nu är detta indexeringsfel, men hoppas ni förstår vad jag är ute efter.
}
}
}
Här är ett utdrag på hur kod i BLAS ser ut
#include "cblas.h"
#include "cblas_f77.h"
void cblas_cher2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
const void *alpha, const void *A, const int lda,
const void *B, const int ldb, const float beta,
void *C, const int ldc)
{
char UL, TR;
#ifdef F77_CHAR
F77_CHAR F77_TR, F77_UL;
#else
#define F77_TR &TR
#define F77_UL &UL
#endif
#ifdef F77_INT
F77_INT F77_N=N, F77_K=K, F77_lda=lda, F77_ldb=ldb;
F77_INT F77_ldc=ldc;
#else
#define F77_N N
#define F77_K K
#define F77_lda lda
#define F77_ldb ldb
#define F77_ldc ldc
#endif
extern int CBLAS_CallFromC;
extern int RowMajorStrg;
float ALPHA[2];
const float *alp=(float *)alpha;
CBLAS_CallFromC = 1;
RowMajorStrg = 0;
if( Order == CblasColMajor )
{
if( Uplo == CblasUpper) UL='U';
else if ( Uplo == CblasLower ) UL='L';
else
{
cblas_xerbla(2, "cblas_cher2k", "Illegal Uplo setting, %d\n", Uplo);
CBLAS_CallFromC = 0;
RowMajorStrg = 0;
return;
}
if( Trans == CblasTrans) TR ='T';
else if ( Trans == CblasConjTrans ) TR='C';
else if ( Trans == CblasNoTrans ) TR='N';
else
{
cblas_xerbla(3, "cblas_cher2k", "Illegal Trans setting, %d\n", Trans);
CBLAS_CallFromC = 0;
RowMajorStrg = 0;
return;
}
#ifdef F77_CHAR
F77_UL = C2F_CHAR(&UL);
F77_TR = C2F_CHAR(&TR);
#endif
F77_cher2k(F77_UL, F77_TR, &F77_N, &F77_K, alpha, A, &F77_lda, B, &F77_ldb, &beta, C, &F77_ldc);
} else if (Order == CblasRowMajor)
{
RowMajorStrg = 1;
if( Uplo == CblasUpper) UL='L';
else if ( Uplo == CblasLower ) UL='U';
else
{
cblas_xerbla(2, "cblas_cher2k", "Illegal Uplo setting, %d\n", Uplo);
CBLAS_CallFromC = 0;
RowMajorStrg = 0;
return;
}
if( Trans == CblasTrans) TR ='N';
else if ( Trans == CblasConjTrans ) TR='N';
else if ( Trans == CblasNoTrans ) TR='C';
else
{
cblas_xerbla(3, "cblas_cher2k", "Illegal Trans setting, %d\n", Trans);
CBLAS_CallFromC = 0;
RowMajorStrg = 0;
return;
}
#ifdef F77_CHAR
F77_UL = C2F_CHAR(&UL);
F77_TR = C2F_CHAR(&TR);
#endif
ALPHA[0]= *alp;
ALPHA[1]= -alp[1];
F77_cher2k(F77_UL,F77_TR, &F77_N, &F77_K, ALPHA, A, &F77_lda, B, &F77_ldb, &beta, C, &F77_ldc);
}
else cblas_xerbla(1, "cblas_cher2k", "Illegal Order setting, %d\n", Order);
CBLAS_CallFromC = 0;
RowMajorStrg = 0;
return;
}
