uBLAS 概覽

原理

如果每一種數值處理軟件都可以由C++來實現而且沒有效率損失,那將是一件非常棒的事情,但是除非找到一個實現這個目的但又不會危及C++類型系統的方法,否則依賴於 Fortran,彙編或體系結構特定的擴展可能是更好的選擇 (Bjarne Stroustrup)。

這個C++類庫直接面向基於一些使用矩陣和向量來構造的基本的線性代數及其抽像操作的科學計算。基本的設計目標是:

另一個目的是評估歷為使用這些矩陣和向量類型所產生的抽像懲罰(abstraction penalty)是否是可接受的。

參考資源

這個庫的開發由一系列的類似的庫來指導:

BLAS 似乎是基本的線性代數使用最廣泛的庫,所以它可以被稱作de-facto standard。它的接口是面向過程的(procedural),單獨的函數是對一些基本的線性代數運算的抽像。由於它是使用Fortran語言實現以及它所進行的優化,BLAS似乎也是現在最快的庫之一。因為我無決定以面向對象的方式來設計和實現我們的類庫,所以技術方式上有顯著的不同。然而,每個人都可以使用我們的庫中的運算符表達所有的BLAS抽像並對效率進行比較。

Blitz++ 是一個使用C++實現的讓人映像深刻的(impressive)庫。它的主要設計似乎是面向多維數組以及包括張量(tensor)的相關的運算符。Blitz++庫的作者聲稱由於它的實現技術使用表達式模板(expression template)和模板元編程技術(template metaprogram),他的庫的性能與對應的Fortan代碼的性能相當或更好一些。然而我們有一些理由來開發一個我們自己設計和實現的方法。我們不知道是否有人使用Blitz++庫來實現傳統的線性代數運算。我們同樣也假定由於Blitz++庫所使用的實現風格(idioms),即使在今天也需要最高級的C++編譯器技術。另一方面,Blitz++庫也使用我們相信,使用表達式模板技術(expression templates)是將抽像懲罰降低到一個可接受限度的必需的技術方法。

POOMA 的設計目標似乎是在許多部分對Blitz++庫的大量部分進行並行計算。它從偏微分方程和理論物理領域中提取類來擴展Blitz++的概念。這種實現支持並行體系結構(parallel architectures)。

MTL 是另一個使用C++來支持基本的線性代數運算的類庫。我們共同的觀點是線性代數庫應當提供與BLAS庫相應的功能。另一方面,我們認為C++標準庫的概念並沒有支持所需要的數值計算的概念要求。另一個區別是MTL庫當前並沒有使用表達式模板技術(expression templates)。這可能導致兩個問題中的一個:可能存在表現能力缺失或可能的性能損失。

概念

數學記法(Mathematical Notation)

數學記法的使用可以簡化科學算法的開發。因此一個實現基本的線性代數概念的C++庫應當在矩陣和向量類上仔細地重載可選的C++運算符。

我們決定對下面的原始概念(primitive)使用運算符重載:

說明 運算符
向量和矩陣的索引 vector::operator(size_t i);
matrix::operator(size_t i, size_t j);
向量和矩陣的賦值 vector::operator = (const vector_expression &);
vector::operator += (const vector_expression &);
vector::operator -= (const vector_expression &);
vector::operator *= (const scalar_expression &);
matrix::operator = (const matrix_expression &);
matrix::operator += (const matrix_expression &);
matrix::operator -= (const matrix_expression &);
matrix::operator *= (const scalar_expression &);
向量和矩陣上的一元運算 vector_expression operator - (const vector_expression &);
matrix_expression operator - (const matrix_expression &);
向量和矩陣上的二元運算 vector_expression operator + (const vector_expression &, const vector_expression &);
vector_expression operator - (const vector_expression &, const vector_expression &);
matrix_expression operator + (const matrix_expression &, const matrix_expression &);
matrix_expression operator - (const matrix_expression &, const matrix_expression &);
向量和矩陣與一個標量的乘法 vector_expression operator * (const scalar_expression &, const vector_expression &);
vector_expression operator * (const vector_expression &, const scalar_expression &);
matrix_expression operator * (const scalar_expression &, const matrix_expression &);
matrix_expression operator * (const matrix_expression &, const scalar_expression &);

對於下面的原始概念(primitives)我們決定不使用運算符重載:

說明 函數
向量和一個矩陣的左乘(Left multiplicationi) vector_expression prod<vector_type > (const matrix_expression &, const vector_expression &);
vector_expression prod (const matrix_expression &, const vector_expression &);
向量和一個矩陣的右乘(Right multiplication) vector_expression prod<vector_type > (const vector_expression &, const matrix_expression &);
vector_expression prod (const vector_expression &, const matrix_expression &);
矩陣乘法 matrix_expression prod<matrix_type > (const matrix_expression &, const matrix_expression &);
matrix_expression prod (const matrix_expression &, const matrix_expression &);
向量的內積( Inner product ) scalar_expression inner_prod (const vector_expression &, const vector_expression &);
向量的外積(Outer product) matrix_expression outer_prod (const vector_expression &, const vector_expression &);
矩陣的轉置 matrix_expression trans (const matrix_expression &);

效率

為了達到數值計算的效率目標,在用C++來公式化(formulating)抽像時需要克服兩個困難,也就是臨時變量和虛函數調用。表達式模板(Expression template)解決了這兩個問題,但是趨向於減緩編譯速度。

避免臨時變量

向量和矩陣的抽像方程通常包含大量的一元和二元運算。最簡單的計算方法是先計算一個組合運算的結點操作(leaf operation)並存放在一個臨時變量中,然後使用另一個臨時變量來進行組合運算。這種方法對於較小的向量和矩陣在時間方面開銷非常大,對於大的向量和矩陣在空間方面的開銷尤其大。解決這個問題的方法是在現代的函數式編譯語言中所稱作的延遲計算(lazy evaluation)。這種方法的原理是明智地計算一個複雜表達式的元素並將其直接賦給目標結果。

兩個有趣並且很危險的因素導致:

別名(Aliases)

在向量和矩陣上使用element wise evaluation會產生嚴重的負作用。考慮矩陣和向量的積x = A x。計算A1x 並賦值給 x1 改變了右邊( right hand side),所以計算A2x 會返回一個錯誤的結果。在這種情況下,元素 xn 在左邊和右邊都有 別名(aliases)

我們解決這個問題的方法是計算等式右邊的值並存儲到一個臨時變量中,然後將這個臨時變量賦給等式的左邊。為了允許更進一步的優化,對於每一個賦值運算符我們都提供了一個對應的成員運算符以及一個 noalias syntax. 。通常使用這種形式,程序員可以確定等式的左邊和右邊是相互獨立的。因此element wise evaluation和直接對計算結果進行賦值是安全的。

複雜度

在某些特定的情況下,計算的複雜度可能是大得不可預料。考慮鏈式矩陣和向量的乘積A (B x)。通常計算A (B x) 是二次方複雜度的。B xi 的延遲求值(Deferred evaluation)是線性的。因為每個元素B xi 依賴於大小都是線性的,一個完全的鏈式矩陣和向量乘積 A (B x) 的延遲求值(Deferred evaluation)是三次方的。在這種情況下,在這個表達式中我們需要重新引入臨時變量。

避免虛函數調用

延遲表達式求值(Lazy expression evaluation)通常會導向(terms)一個類的繼承層次的定義。 這就導致使用動多態(dynamic polymorphism )來訪問一個向量或矩陣中的單一元素,在時間方面這種開銷很大。在多年以前,由 David Vandervoorde 分別發現了一種技術稱作表達式模板(expression templates)。表達式模板( Expression templates contain)包含延遲計算( lazy evaluation)並使用靜多態(static polymorphism )來代替動多態(dynamic polymorphism ),也就是編譯時多態(compile time polymorphism)。表達式模板技術強烈依賴於著名的 Barton-Nackman 技巧(trick),同樣也被Jim Coplien定義為遞歸模板。

表達式模板(Expression templates)構成我們的實現的基礎。

編譯時間

表達式模板技術挑戰當前的編譯器技術同樣也是一個知名的事實。使用 Barton-Nackman 技巧我們可以顯著地減少表達式模板的使用數量。

我們同樣也決定支持雙重的常規實現(dual conventional implementation) (例如,不使用表達式模板) ,使用向量和矩陣的擴展邊界和類型檢查來支持開發週期(development cycle)。從調試模式轉換到發行模式是通常<cassert>中的預編譯符號NDEBUG 來控制的。

功能

每一個支持線性代數的C++庫都會相比於長期存在的Fortran 包 BLAS進行度量。現在我們說明如何將BLAS中的函數調用映射到我們的類中。

矩陣和向量運算概覽頁面對於大多數的向量和矩陣運算給出了一個簡短的總結。

Blas Level 1

BLAS Call Mapped Library Expression Mathematical Description Comment
sasum OR dasum norm_1 (x) sum |xi| 計算實數向量的l1 (sum) 基。
scasum OR dzasum real (sum (v)) + imag (sum (v)) sum re(xi) + sum im(xi) C計算一個複數向量的元素的和。
_nrm2 norm_2 (x) sqrt (sum |xi|2 ) 計算一個向量的 linf (euclidean) 基。
i_amax norm_inf (x)
index_norm_inf (x)
max |xi| 計算一個向量的l2 (maximum) 基
BLAS 計算第一個包含這個值的元素的下標。
_dot
_dotu
_dotc
inner_prod (x, y)or
inner_prod (conj (x), y)
xT y or
xH y
計算兩個向量的內積。
BLAS 實現特定的循環展開。
dsdot
sdsdot
a + prec_inner_prod (x, y) a + xT y 以double精度計算內積。
_copy x = y
y.assign (x)
x <- y 將一個向量拷貝到另一個向量中。
BLAS 實現特定的循環展開。
_swap swap (x, y) x <-> y 交換兩個向量的內容。
BLAS 實現特定的循環展開。
_scal
csscal
zdscal
x *= a x <- a x 擴充一個向量。
BLAS 實現特定的循環展開。
_axpy y += a * x y <- a x + y 增加一個成比例的向量。
BLAS 實現特定的循環展開。
_rot
_rotm
csrot
zdrot
t.assign (a * x + b * y),
y.assign (- b * x + a * y),
x.assign (t)
(x, y) <- (a x + b y, -b x + a y) 應用一個平面旋轉。
_rotg
_rotmg
  (a, b) <-
  (? a / sqrt (a
2 + b2),
    ? b / sqrt (a
2 + b2)) or
(1, 0) <- (0, 0)
構造一個平面旋轉。

Blas Level 2

BLAS 調用 映射為庫中的表達式 數學說明 註釋
_t_mv x = prod (A, x) or
x = prod (trans (A), x)
or
x = prod (herm (A), x)
x <- A x or
x <- A
T x or
x <- A
H x
計算一個矩陣與一個向量的乘積。
_t_sv y = solve (A, x, tag) or
inplace_solve (A, x, tag) or
y = solve (trans (A), x, tag) or
inplace_solve (trans (A), x, tag) or
y = solve (herm (A), x, tag)or
inplace_solve (herm (A), x, tag)
y <- A-1 x or
x <- A
-1 x or
y <- A
T-1 x or
x <- A
T-1 x or
y <- A
H-1 x or
x <- A
H-1 x
使用三角形式來求解線性方程組,例如, A 是三角形的。
_g_mv
_s_mv
_h_mv
y = a * prod (A, x) + b * y or
y = a * prod (trans (A), x) + b * y
or
y = a * prod (herm (A), x) + b * y
y <- a A x + b y or
y <- a A
T x + b y
y <- a A
H x + b y
增加一個矩陣與一個向量的scaled product。
_g_r
_g_ru
_g_rc
A += a * outer_prod (x, y) or
A += a * outer_prod (x, conj (y))
A <- a x yT + A or
A <- a x y
H + A
進行一個秩(rank) 1 變換。
_s_r
_h_r
A += a * outer_prod (x, x) or
A += a * outer_prod (x, conj (x))
A <- a x xT + A or
A <- a x x
H + A
P進行一個對稱的或埃爾米特秩(rank) 1 變換。
_s_r2
_h_r2
A += a * outer_prod (x, y) +
 a * outer_prod (y, x))
or
A += a * outer_prod (x, conj (y)) +
 conj (a) * outer_prod (y, conj (x)))
A <- a x yT + a y xT + A or
A <- a x y
H + a- y xH + A
進行一個對稱的或埃爾米特秩(rank) 2 變換。

Blas Level 3

BLAS 調用 映射為庫中的表達式 數學說明 註釋
_t_mm B = a * prod (A, B) or
B = a * prod (trans (A), B) or
B = a * prod (A, trans (B)) or
B = a * prod (trans (A), trans (B)) or
B = a * prod (herm (A), B) or
B = a * prod (A, herm (B)) or
B = a * prod (herm (A), trans (B)) or
B = a * prod (trans (A), herm (B)) or
B = a * prod (herm (A), herm (B))
B <- a op (A) op (B) with
  op (X) = X or
  op (X) = XT or
  op (X) = XH
C計算兩個矩陣的scaled product。
_t_sm C = solve (A, B, tag) or
inplace_solve (A, B, tag) or
C = solve (trans (A), B, tag) or
inplace_solve (trans (A), B, tag)
or
C = solve (herm (A), B, tag)
or
inplace_solve (herm (A), B, tag)
C <- A-1 B or
B <- A
-1 B or
C <- A
T-1 B or
B <- A
-1 B or
C <- A
H-1 B or
B <- A
H-1 B
使用三角形式來求解線性方程組,例如, A 是三角形的。
_g_mm
_s_mm
_h_mm
C = a * prod (A, B) + b * C or
C = a * prod (trans (A), B) + b * C or
C = a * prod (A, trans (B)) + b * C or
C = a * prod (trans (A), trans (B)) + b * C or
C = a * prod (herm (A), B) + b * C or
C = a * prod (A, herm (B)) + b * C or
C = a * prod (herm (A), trans (B)) + b * C or
C = a * prod (trans (A), herm (B)) + b * C or
C = a * prod (herm (A), herm (B)) + b * C
C <- a op (A) op (B) + b C with
  op (X) = X or
  op (X) = XT or
  op (X) = XH
增加兩個矩陣的scaled product。
_s_rk
_h_rk
B = a * prod (A, trans (A)) + b * B or
B = a * prod (trans (A), A) + b * B or
B = a * prod (A, herm (A)) + b * B or
B = a * prod (herm (A), A) + b * B
B <- a A AT + b B or
B <- a A
T A + b B or
B <- a A AH + b B or
B <- a A
H A + b B
進行一個對稱的或埃爾米特秩(rank)k變換。
_s_r2k
_h_r2k
C = a * prod (A, trans (B)) +
 a * prod (B, trans (A)) + b * C
or
C = a * prod (trans (A), B) +
 a * prod (trans (B), A) + b * C
or
C = a * prod (A, herm (B)) +
 conj (a) * prod (B, herm (A)) + b * C
or
C = a * prod (herm (A), B) +
 conj (a) * prod (herm (B), A) + b * C
C <- a A BT + a B AT + b C or
C <- a A
T B + a BT A + b C or
C <- a A B
H + a- B AH + b C or
C <- a A
H B + a- BH A + b C
進行一個對稱的或埃爾米特秩(rank) 2 k 變換。

存儲佈局(Storage Layout)

uBLAS 支持許多不同的存儲佈局(storage layouts)。可以在類型概覽中找到詳細信息。大多數類似於vector<double>matrix<double> 都默認與C風格的數組相兼容,但是可以設置為包含與 FORTAN 兼容的數據。

兼容性

出於兼容性的考慮,我們提供了類似於數組的下標索引。對於一些類型 (埃爾米特矩陣,稀疏矩陣等等),由於會使用臨時的代理對象,這會產生很大的開銷。

uBLAS 使用與 STL兼容的分配器(allocators)來分配容器所需要的存儲空間。

基準結果(Benchmark Results)

下面的表中包含了我們的一個基準中的結果。這個基準比較一個C實現(「C數組」)和一些基於庫的實現。 安全的變體(variants)假定使用別名(aliasing),最快的變體(variants)不使用臨時變量並且在功能上等價於C數組實現。除了一般化的向量和矩陣類之外,比較的基準(benchmark)還使用特殊類c_vectorc_matrix, 通過泛化來避免每一個過度開銷(overhead)。

基準程序(benchmark program) bench1 使用 GCC 4.0 編譯並運行在一台 Athlon 64 3000+的機器上。通過運行 bench1 100來將時間調整到一個合理的精度。

首先我們分別註釋一個double類型的3維向量和一個double類型的3X3的矩陣。

註釋
inner_prod C array 0.61 782 少量的抽像懲罰
c_vector 0.86 554
vector<unbounded_array> 1.02 467
vector + vector C array 0.51 1122 抽像懲罰: 因子為 2
c_vector fast 1.17 489
vector<unbounded_array> fast 1.32 433
c_vector safe 2.02 283
vector<unbounded_array> safe 6.95 82
outer_prod C array 0.59 872 少量的抽像懲罰
c_matrix, c_vector fast 0.88 585
matrix<unbounded_array>, vector<unbounded_array> fast 0.90 572
c_matrix, c_vector safe 1.66 310
matrix<unbounded_array>, vector<unbounded_array> safe 2.95 175
prod (matrix, vector) C array 0.64 671 沒有明顯的抽像懲罰
c_matrix, c_vector fast 0.70 613
matrix<unbounded_array>, vector<unbounded_array> fast 0.79 543
c_matrix, c_vector safe 0.95 452
matrix<unbounded_array>, vector<unbounded_array> safe 2.61 164
matrix + matrix C array 0.75 686 沒有明顯的抽像懲罰
c_matrix fast 0.99 520
matrix<unbounded_array> fast 1.29 399
c_matrix safe 1.7 303
matrix<unbounded_array> safe 3.14 164
prod (matrix, matrix) C array 0.94 457 沒有明顯的抽像懲罰
c_matrix fast 1.17 367
matrix<unbounded_array> fast 1.34 320
c_matrix safe 1.56 275
matrix<unbounded_array> safe 2.06 208

對於較小的向量和矩陣我們發現了一個25%(two fold )的性能損失:首先對於使用類而產生的抽像懲罰,然後是使用泛化的向量和矩陣類所產生的較小的性能損失。 w.r.t. 差別假定也是顯著的。

接下來,我們分別註釋一個double類型的100維向量和一個double類型的100X100的矩陣。

運算 實現 用時 [s] MFLOP/s 註釋
inner_prod C array 0.64 889 沒有明顯的抽像懲罰
c_vector 0.66 862
vector<unbounded_array> 0.66 862
vector + vector C array 0.64 894 沒有明顯的抽像懲罰
c_vector fast 0.66 867
vector<unbounded_array> fast 0.66 867
c_vector safe 1.14 501
vector<unbounded_array> safe 1.23 465
outer_prod C array 0.50 1144 沒有明顯的抽像懲罰
c_matrix, c_vector fast 0.71 806
matrix<unbounded_array>, vector<unbounded_array> fast 0.57 1004
c_matrix, c_vector safe 1.91 300
matrix<unbounded_array>, vector<unbounded_array> safe 0.89 643
prod (matrix, vector) C array 0.65 876 沒有明顯的抽像懲罰
c_matrix, c_vector fast 0.65 876
matrix<unbounded_array>, vector<unbounded_array> fast 0.66 863
c_matrix, c_vector safe 0.66 863
matrix<unbounded_array>, vector<unbounded_array> safe 0.66 863
matrix + matrix C array 0.96 596 沒有明顯的抽像懲罰
c_matrix fast 1.21 473
matrix<unbounded_array> fast 1.00 572
c_matrix safe 2.44 235
matrix<unbounded_array> safe 1.30 440
prod (matrix, matrix) C array 0.70 813 沒有明顯的抽像懲罰
c_matrix fast 0.73 780
matrix<unbounded_array> fast 0.76 749
c_matrix safe 0.75 759
matrix<unbounded_array> safe 0.76 749

對於較大的向量和矩陣,使用類的一笛膜的抽像懲罰似乎在降低,但是使用泛化的向量和矩陣而產生的微量的性能損失似乎仍然存在。 w.r.t. 差別假定也是可見的(visible)。


Copyright (©) 2000-2002 Joerg Walter, Mathias Koch
Use, modification and distribution are subject to the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt ).