LAPACK 是 Linear Algebra PACKage 的縮寫，它是一個以 Fortran 語言寫成的函式庫。 本篇的目標，是要以矩陣相乘等常見的函式為例，來介紹 LAPACK 的使用。 由於自己下載原始碼，並編譯出函式庫的步驟比較繁瑣， 所以這裡提供了事先編譯好的 lapack.rar 讓各位下載， 請解壓縮後把它們放在同一個資料夾，例如"C:\lapack"。 (如果對自行編譯有需要的話，可以檢視本網頁的原始碼，裡面用註解的方式記錄了一些比較重要的步驟。)

接著，我們要讓 Code::Blocks 知道函式庫放在哪裡。請先選擇 Settings → Compiler，應該會出現這個畫面：


再來，請找到 linker settings 這個分頁，把三個 *.a 的函示庫加入路徑中

下一步，請找到 Search directories 這個分頁， 加入 C:\lapack 或你放置的其他路徑，才能讓 Code::Blocks 找到相對應的 *.h 檔案。

現在，你可以先試著編譯與執行這個程式：

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <f2c.h>
#include <clapack.h>
 
void matrixDisplay(double mat[], int m, int n){
    int i, j;
    for(i=0;i<m;i++){
        for(j=0;j<n;j++){
            printf("%.0f ",mat[i*n+j]);
        }
        printf("\n");
    }
}
 
int main()
{
    int n;
    double *A,*B,*C;
    double alpha=1.0, beta=0.0;
    char trans='t';
    int i,j;
 
    clock_t tic,toc;
 
    /* space allocation */
    n = 2;
    A = (double*)malloc(n*n*sizeof(double));
    B = (double*)malloc(n*n*sizeof(double));
    C = (double*)malloc(n*n*sizeof(double));
 
    /* value generation/initialization */
    srand(time(NULL));
    j = n*n;
    for(i=0; i<j; i++){
        A[i] = rand()%5+1;
        B[i] = rand()%5+1;
        C[i] = 0.0;
    }
 
    /* computation */
    tic = clock();
    dgemm_(&trans, &trans, &n, &n, &n,&alpha, A, &n, B, &n,&beta, C, &n);
    toc = clock();
    printf("%.4f\n",(double)(toc-tic)/CLOCKS_PER_SEC);
 
    /* display */
    printf("==================\n");
    matrixDisplay(A, n, n);
    printf("==================\n");
    matrixDisplay(B, n, n);
    printf("==================\n");
    matrixDisplay(C, n, n);
    printf("==================\n");
     
    free(A);
    free(B);
    free(C);
    return 0;
}

這個程式會亂數產生兩個 2-by-2 的矩陣，然後利用 LAPACK 計算矩陣乘積。 首先，要使用 LAPACK 的功能，必須含括 f2c.h 和 clapack.h。 time.h 則是計時用，在此範例中，由於矩陣很小，所以時間不明顯。 再來，對於矩陣的部分，我們必須使用一維矩陣來存放資料，然後自行對應至二維。

然後是最重要的，呼叫 LAPACK 的部分，我們使用的是「dgemm_」這個函式。 它會幫你計算這個式子：

C = α*op(A)*op(B) + β*C

其中，op() 代表矩陣要轉置或者不要轉置，由輸入參數決定。 命名的規則稍後來說明，我們先介紹參數如下， 完整的說明可以參考這裡：
所以，如果要進行兩個方陣 A 和 B 的相乘，則矩陣大小皆為 n，α 為 1，β 為 0。 在 trans 的部分，由於在 C 語言裡面，把二維放成一維時，習慣一個 row 接一個 row 放(稱為 row-major)； 而 LAPACK 使用的 Fortran 語言，習慣一個 column 接一個 column 放(稱為 column-major)， 所以 trans 在此處皆使用 't'。 當然，輸出的 C 矩陣是 column-major，而此處在印出前並沒有特別轉換，所以你會覺得看到的矩陣方向不太對。 另外，此函式庫需要傳入記憶體位置，所以若是非指標型態的變數，記得加上「&」。

接下來我們介紹命名規則。由於早期的 Fortran 在函式名稱的字數上有些限制， 因此有規則如下(少數特例本篇不介紹)：

因為 LAPACK 提供的函式很多，本篇就不一一介紹。 以下示範 double 類型的 LU 分解的範例，需要使用的函式是「dgetrf_」，參數說明如下：


完整程式如下：
?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <f2c.h>
#include <clapack.h>
 
void matrixDisplay(double mat[], int m, int n){
    int i, j;
    for(i=0;i<m;i++){
        for(j=0;j<n;j++){
            printf("%.0f ",mat[i*n+j]);
        }
        printf("\n");
    }
}
 
int main()
{
    int n, *ipiv;
    double *A;
    int i, j, info;
     
    /* space allocation */
    n = 3;
    A = (double*)malloc(n*n*sizeof(double));
    ipiv = (int*)malloc(n*sizeof(int));
 
    /* value generation/initialization */
    srand(time(NULL));
    j = n*n;
    for(i=0; i<j; i++){
        A[i] = rand()%5+1;
    }
     
    matrixDisplay(A, n, n);
    dgetrf_(&n, &n, A, &n, ipiv, &info);
    matrixDisplay(A, n, n);
     
    free(A);
 
    return 0;
}