"A designer knows he has achieved perfection not when there is nothing left to add, but when there is nothing left to take away."
- Antoine de Saint-Exupery
Her zamanki gibi, yazı minimal bir Centos6.9 makinayla başlıyor. Kurulumu ve diğer adımları hem daha önceki yazılarda anlattım hem de bu yazının amacıyla bağdaşmıyor. Bir makinanın kurulu olduğunu varsayarak devam ediyorum.
 |
| yum install gcc gcc-gfortran |
Makinayı kurar kurmaz yum update ile güncellemelerini geçip makinayı yeniden başlattım. Bu adım kesinlikle gerekli değil ama olsa çok iyi olur. Kütüphaneyi derlemek için gereken gcc ve gcc-gfortran paketlerini yum'la kurdum.
Ardından LAPACK kütüphanesini indirdim. LAPACK, Linear Algebra Package kelimelerinin kısaltması. Fortran'la geliştirilmiş ve içinde doğrusal denklem çözümleri, en küçük kareler, özdeğer/özvektör hesaplamaları, matris çarpanlara ayırma algoritmaları gibi önemli matris işlemleri için fonksiyonlar var. İşlemler saydıklarımla sınırlı değil, öğrenim hayatımdan aklımda kalanları bunlar. Gündelik hayatta karşılaşılmasa da birçok bilimsel hesaplama programı bu kütüphane üzerine kurulu.
Ardından LAPACK kütüphanesini indirdim. LAPACK, Linear Algebra Package kelimelerinin kısaltması. Fortran'la geliştirilmiş ve içinde doğrusal denklem çözümleri, en küçük kareler, özdeğer/özvektör hesaplamaları, matris çarpanlara ayırma algoritmaları gibi önemli matris işlemleri için fonksiyonlar var. İşlemler saydıklarımla sınırlı değil, öğrenim hayatımdan aklımda kalanları bunlar. Gündelik hayatta karşılaşılmasa da birçok bilimsel hesaplama programı bu kütüphane üzerine kurulu.
İkinci önemli kütüphane BLAS. Basic Linear Algebra Subroutines'in kısaltması olan BLAS vektör ve matris işlemleri sunan daha alt seviye bir kütüphane. LAPACK, derlenebilmek için BLAS'a ihtiyaç duyuyuyor.
İlkin açıklığa kavuşturulması gereken, basit bir matris çarpma veya toplama işlemi için gerçekten bu kütüphanelere ihtiyacımız olup olmadığı.
Bu sorunun yanıtı için matris cebirinin en alt seviye tanım ve teoremlerine girmek istemem. Yapabildiğim kadar çevresinden dolaşmaya çalışacağım. Matrislerde çarpma, birinci matris (m x n) ve ikinci matris (n x p) boyutlu olmak üzere tanımlıdır. Kolaylık için birinci matrise A, ikinci matrise B diyip m = n = p kabulüyle matrislerin kare olduğunu düşünelim. Sonuç matrisine C diyelim. C, (p x p) boyutlu bir matristir. C'nin birinci elemanı, A'nın birinci satır elemanlarıyla
B'nin birinci sütun elemanlarının tek tek çarpımlarının toplamıdır. C matrisinde birinci elemanın sağında bulunan C'nin ikinci elemanı; A'nın ikinci satır elemanlarıyla B'nin birinci sütun elemanlarının tek tek çarpımları toplamıdır. Sözle ifade etmek kolay olmadığından matematiksel ifadelere başvurmak zorundayım. A matrisi yandaki ve B matrisi aşağıdaki gibi olsun:
A ile B'nin çarpımı aşağıdaki gibi tanımlıdır:
Elde edilecek elemanların C matrisindeki yerleri şöyledir:
Eğer işlem yalancı kodla ifade edilmek istenirse aşağıdaki gibi olur:
for i = 1 to p
for j = 1 to p
c(i)(j) = 0
for k = 1 to p
c(i)(j) = c(i)(j) + a(i)(j) * b(i)(j)
for j = 1 to p
c(i)(j) = 0
for k = 1 to p
c(i)(j) = c(i)(j) + a(i)(j) * b(i)(j)
Bu işlemin karmaşıklığı p'nin küpüne bağlı olarak değişir. Büyük O gösterimiyle O(n3) olarak ifade edilir. Başka bir deyişle (p x p) boyutlu iki matris birim zamanda çarpılıyorsa, (2p x 2p) boyutlu iki matris sekiz birim zamanda çarpılır. Literatürde bu işleme çeşitli yaklaşımlar geliştirilmiş ve işlemin karmaşıklığı farklı algoritmalarla O(n2.37)'ye kadar indirilmiştir (ek okuma: Matrix multiplication algorithm).
Konunun kodlama ve optimizasyon boyutu teorisinden oldukça farklıdır. Örn. yukarıda verdiğim yalancı kod C'nin satır yerleşimli dizilerinde cache bellek yerleşiminden ötürü optimal değildir (bilimsel hesaplama veya sayısal algoritmalar dersinin giriş konusudur) bunun yerine indisleri i, k, j olarak sıralamak daha iyi cache hit oranı sağlar. Başka bir deyişle işlemlerin daha fazlası cache bellek içinde yapılır ve algoritma daha hızlı çalışır.
Bu algoritmaların uygulamalarını üniversitedeyken araştırmıştım. Önce yukarıdakine benzer bir kodla başlayıp cache kullanımını iyileştirdim ve değişkenleri işlemciye yakın tutarak kodu hızlandırdım. Sonra algoritmayı değiştirip Strassen Algoritması'nı ekledim. Üçüncü adımda bu algoritmayı özyinelemeli (recursive) duruma getirip teoride elde edilebilecek sınıra getirdim. Kodu derlerken bazı mikroişlemci optimizasyonlarını da devreye alıp elde edebildiğim en iyi sonucu buldum.
Karşılaştırmak için aynı matrisleri BLAS'la çarptığımda aynı derleyici parametrelerini de kullanmama rağmen maalesef BLAS yaklaşık 1/5 oranında daha hızlı çalışıyordu (yani 5 sn yerine 4 sn). Daha da üzücü olan, BLAS'ın yaklaşık 60sn'de çarptığı iki matrisi Matlab'in bir saniyeden biraz daha uzun bir sürede çarpmasıydı. Matlab'in kodu, güncel mikroişlemci mimarisine çok daha uygun yazılmıştı. Araştırmada kullandığım işlemci SSE2 destekliydi ve matris çarpımında yapılan çarpmaların sayısını işlemcinin çekirdek saat hızına böldüğümde neredeyse bire yakın bir oran çıkıyordu. Yani Matlab çarpma yaparken işlemciyi o kadar verimli kullanılıyordu ki başka hiçbir iş yapmıyordu.
Uzun lafın kısası, gerçekten birşeyler yapabilmek için bu kütüphanelere gerek var. BLAS, otuz yıllık bir kütüphane ve bu zamanda yüzlerce insan tarafından geliştirildi. Basit bir araştırma projesinin bir kaç ayda bu kütüphaneyi geçmesini beklemek insafsızlık olurdu. BLAS terminolojisinde birinci seviye işlemler vektörlerle ilgili çarpımlar, ikinci seviye işlemler matrislerin vektörlerle işlemleri ve üçüncü seviye işlemler matris-matris işlemleri olarak adlandırılır. LAPACK, BLAS fonksiyonlarını kullanır ve daha karmaşık işlemleri yapmaya olanak sağlar.
Açık konuşmak gerekirse üniversitedeyken bu kütüphanelerin adını çok duydum. Keza iş hayatında çok kere bu kütüphaneleri ve bunlara bağlı çalışan programları derledim ancak ne derslerde ne de uygulamalarda kütüphanelerin kullanımı verilmedi. Bu kelimeler bilinmeyen diyarlara ait Ebabil kuşu gibi kullanılırdı. BLAS'ı üniversitede bir kaç projemde kullandım. LAPACK'ı hiç kullanmadım. Bugün dışarıdan bir gözle baktığımda bu durumu oldukça olumsuz görüyorum.
Uzun bir teorik aradan sonra tekrar yazının konusuna döneyim. Yazıyı yazdığım zaman itibariyle LAPACK'ın en yenisi 3.8.0 idi. Bunu indirip açtım:
curl http://www.netlib.org/lapack/lapack-3.8.0.tar.gz > lapack-3.8.0.tar.gz
tar xvfz lapack-3.8.0.tar.gz
Not: İndirmek için curl kullandım çünkü minimal CentOS'da wget yok.
tar xvfz lapack-3.8.0.tar.gz
Lapack açıldıktan sonra oluşan dizine girip Makefile'ın içerdiği (include) ayar dosyasını örnek dosyadan kopyalayıp Makefile'daki hedefleri derledim:
cp make.inc.example make.inc
make lib
make blaslib
make cblaslib
make lapackelib
make lib
make blaslib
make cblaslib
make lapackelib
"lib" LAPACK'ın Fortran kütüphanesi. "blaslib", LAPACK ile birlikte gelen BLAS'ın Fortran kütüphanesi. "cblaslib" BLAS'ın C arabirimi ve "lapackelib" de LAPACK'ın C arabirimi. Bu hedeflerin derlenmesi gerek.
Yazdığım kodu Google Drive'da paylaştım. Aşağıdaki bağlantıdan indirilebilir:
En baştaki matrisyaz fonksiyonu, matrisi ekrana yazdırıyor. Çıktıyı Matlab'a uyumlu yazdırıyorum ki sonuçların doğruluğunu kontrol edebileyim. main fonksiyonunun içinde boyut = 3 olarak tanımladım. Bu değişken matrislerin satır/sütun sayılarını tutuyor. Aşağıdaki kod sayesinde programa argüman olarak bir sayı vererek boyutu değiştirmek mümkün.
if(argc > 1) {
boyut = atoi(argv[1]);
}
boyut = atoi(argv[1]);
}
Lapack'ın kendi lapack_int tamsayı tanımlaması olduğundan boyutu bu değişkene atadım. matrisA, üzerinde işlem yaptığım matris. matrisB, ilk başta matrisA'ya eşit. matrisA'nın tersini aldıktan sonra sağlama için matrisB'yle çarpıp sonucun birim matris olduğunu göstereceğim. Matrisleri malloc ile oluşturup matrisA ile matrisB'nin içeriğini rastgele doldurdum. ipiv vektörü işleme girecek olan matrislerin satır sırasını değiştirmekte kullanılacak vektör. Ben boş bırakıyorum, kullanmayacağım.
Kodda üç basit işlem var. İlkin dgetrf fonksiyonuyla matrisi LU çarpanlarına ayırdım. LAPACK'ın fonksiyon adlandırması şöyle: Baştaki 'd' double anlamına geliyor. 's' single, 'c' complex ve 'z' double complex anlamında. İlk harfler işleme girecek matris elemanlarının değerliklerini (precision) belirtiyor. Sonraki iki karakter 'ge' genel matrislerle ilgili fonksiyon anlamında. Matrisin simetik, köşegen, üç-köşegen, üçgen matris vb. türleri farklı iki karakterle kodlanıyor. Ayrıntılı liste LAPACK belgelendirmesinde veya Wikipedia sayfasında bulunabilir. Son üç karakterse fonksiyonun kısaltması. Örn. yanlış hatırlamıyorsam 'trf' triangular factorization'ın kısaltması. LU ayrıştırmasını yapan kod aşağıda:
LAPACKE_dgetrf(LAPACK_ROW_MAJOR, lapack_boyut, lapack_boyut, matrisA, lapack_boyut, ipiv);
C'de bir malloc yapıldığında bellekteki ardışık iki elemanın, matrisin bir satırındaki ardışık iki eleman olduğu düşünülür. Örn. a11 ve a12 gibi. Fortran'daysa bellekteki ardışık iki elemanın, matrisin sütunlarındaki ardışık elemanlar olduğu düşünülür, a11 ve a21 olarak. LAPACK_ROW_MAJOR değeri matrisA'nın satır olarak işlenmesi gerektiğini gösterir. Sonra gelen lapack_boyut değişkenleri matrisin satır ve sütun sayısı; matrisA, matrisin göstericisidir. Bir sonraki lapack_boyut argümanı LDA olarak ifade edilir. Açılımı Leading Dimension of Array olup matrisin tümü değil de bir parçası üzerinde işlem yaparken anlamlıdır. Ben tümünde işlem yapacağım için aynı boyutu verdim (konuyu dağıtmamak adına bunu en son açıklayacağım). Son olarak ipiv adlı boş vektöru fonksiyona verdim.
Ardından, dgetri fonksiyonuyla yukarıda ayrıştırdığım matrisin çarpımsal tersini alıyorum.
LAPACKE_dgetri(LAPACK_ROW_MAJOR, lapack_boyut, matrisA, lapack_boyut, ipiv);
'tri', triangular inverse olması gerekiyor. Cebire girmek istememiştim ama burada kısaca bahsetmeliyim. Sayılardaki çarpmada birim eleman (ilkokuldaki adıyla etkisiz eleman) '1'dir. Bir sayı birim elemanla çarpıldığında sonuç kendisidir. Bir elemanın işlemsel tersi, söz konusu elemanla işleme girdiğinde birim elemanı veren elemandır. Başka bir deyişle sayı, tersiyle çarpıldığında birim elemanı vermelidir. Çarpmada x sayısının tersi 1/x'tir. Sıfır dışında her elemanın, çarpma işlemine göre tersi vardır. Matrislerde ufak istisnalar dışında yukarıdaki kurallar aynen geçerlidir. Matrislerde birim eleman -birim matris- köşegeninde 1'ler diğer elemanları 0'lar olan matristir ve hepsi değilse de birçok matrisin çarpımsal tersi bulunabilir. (Meraklısına: https://proofwiki.org/wiki/Real_Numbers_form_Ring)
Yalnızca kare matrislerin tersi (aslında tek bir tersi demek daha doğru sanırım) olabileceğinden burada lapack_boyut argümanı bir kere veriliyor. Dördüncü argüman yine LDA.
Buraya kadar yapılan, LAPACK fonksiyonlarıyla matrisA'nın tersini bulmaktı. Şimdi üçüncü seviye bir BLAS fonksiyonu kullanarak matrisB'de sakladığım ilk değeri, tersi olduğunu iddia ettiğim matrisA ile çarpacağım ve doğruysa birim matrisi elde edeceğim. LAPACK fonksiyonları gibi BLAS fonksiyonlarının da adlandırma kuralları var. Matris çarpımı için dgemm fonksiyonunu kullandım. Double GEneral Matrix Multiplication ifadesinin kısaltması. Çarpımı yapan kod şöyle:
cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans, lapack_boyut, lapack_boyut, lapack_boyut, 1.0, matrisA, lapack_boyut, matrisB, lapack_boyut, 0.0, matrisC, lapack_boyut);
CBLAS'ın da LAPACK'a benzeyen öntanımlı değerleri var. CblasRowMajor, bunlardan biri ve görevi LAPACK_ROW_MAJOR ile aynı. İki CblasNoTrans argümanı sırasıyla işleme giren matrisA'nın ve matrisB'nin olduğu gibi işleme sokulması gerektiğini belirtiyor. Eğer bu argüman CblasTrans verilirse matrislerin devriği (transpozesi) işleme sokulur. Devrik işlemi, satırların sütun, sütunların satır durumuna gelmesidir. Türkçe wikipedia'da bu Tersçapraz olarak ifade ediliyor. Matris çarpımının tanımını verirken, çarpımın (m x n) boyutlu matrislerle (n x p) boyutlu matrisler arasında tanımlı olduğunu yazmıştım. Bu durumda sonuç (m x p) boyutlu bir matris olur. Bu fonksiyondaki birbirini tekrar eden üç lapack_boyut argümanı (dört, beş ve altıncı argümanlar) sırayla m, n ve p sayılarına karşılık gelir ama basitlik için bu örnekte m = n = p alınmıştır.
dgemm aslında aşağıdaki işlemi yapar:
 |
| dgemm |
Yani A'nın elemanlarını alfa ile çarpıp B matrisiyle matris çarpımı yapar, C'nin elemanlarını beta ile çarpıp sonucu toplar ve C matrisine yazar. Yapmak istediğim yalnız matris çarpımı olduğundan alfa, 1.0 ve beta 0.0 olarak alınmalıdır. İşte matrisA argümanının önündeki 1.0, alfa parametresi olup matrisA'dan sonraki lapack_boyut argümanı da LDA'dir. Sonrasında matrisB'nin göstericisi ve matrisB'ye ait LDB argümanı bulunur. Son bölümde 0.0 değerli beta sabiti, matrisC'nin göstericisi ve matrisC'ye ait LDC yer alır.
Kaynaklar:
Kodun içinde kolaylık olsun diye fonksiyonların prototiplerini de açıklama olarak koydum.
Kodu aşağıdaki komutla derledim:
gcc matris.c -I/root/lapack-3.8.0/CBLAS/include \
-I/root/lapack-3.8.0/LAPACKE/include -L/root/lapack-3.8.0 \
-llapacke -llapack -lrefblas -lcblas -lgfortran -o matris.x
-I/root/lapack-3.8.0/LAPACKE/include -L/root/lapack-3.8.0 \
-llapacke -llapack -lrefblas -lcblas -lgfortran -o matris.x
Kodun çıktısı aşağıdaki gibi:
Sonuçlar GNU Octave'daki sonuçlarla tamamen tutarlı:
Ek: LDA nedir? Ne zaman gereklidir?
Tekrar LU ayrıştırma fonksiyonuna dönelim:
A matrisi = [
0.000000 4.000000 0.000000 ;
3.000000 7.000000 1.000000 ;
5.000000 5.000000 1.000000 ]
LU sonrasi matrisA = [
5.000000 5.000000 1.000000 ;
0.000000 4.000000 0.000000 ;
0.600000 1.000000 0.400000 ]
matrisA'nin tersi = [
0.250000 -0.500000 0.500000 ;
0.250000 -0.000000 0.000000 ;
-2.500000 2.500000 -1.500000 ]
carpim = [
1.000000 -0.000000 0.000000 ;
0.000000 1.000000 0.000000 ;
0.000000 0.000000 1.000000 ]
0.000000 4.000000 0.000000 ;
3.000000 7.000000 1.000000 ;
5.000000 5.000000 1.000000 ]
LU sonrasi matrisA = [
5.000000 5.000000 1.000000 ;
0.000000 4.000000 0.000000 ;
0.600000 1.000000 0.400000 ]
matrisA'nin tersi = [
0.250000 -0.500000 0.500000 ;
0.250000 -0.000000 0.000000 ;
-2.500000 2.500000 -1.500000 ]
carpim = [
1.000000 -0.000000 0.000000 ;
0.000000 1.000000 0.000000 ;
0.000000 0.000000 1.000000 ]
Sonuçlar GNU Octave'daki sonuçlarla tamamen tutarlı:
Son olarak dikkat çekmek istediğim bir nokta var. Çıktıda ikinci eleman -0 olarak görünüyor. Bu ayrı bir yazının konusu, o yüzden gelecek yazılardan birinde değineceğim.
Ek: LDA nedir? Ne zaman gereklidir?
Tekrar LU ayrıştırma fonksiyonuna dönelim:
LAPACKE_dgetrf(LAPACK_ROW_MAJOR, lapack_boyut, lapack_boyut, matrisA, lapack_boyut, ipiv);
Varsayalım ki görece büyük bir matris var. Örn. 20 x 20. Bunun sol üstündeki 5 x 6'lık parçasıyla işlem yapmak istiyorum. Matlab'deki ifadesiyle A(1:5, 1:6). Bu büyük bir matrisin parçası olduğundan, parçanın satır sonundaki elemandan hemen sonra bellekte bir sonraki satırın ilk elemanı bulunmaz. Bir sonraki ilk eleman 20 - 6 = 14 birim (float, complex vb) uzaklıktadır. İşte LDA alt alta olan sütun elemanları arasındaki mesafenin birimidir. Yani yukarıdaki işlem için fonksiyon:
LAPACKE_dgetrf(LAPACK_ROW_MAJOR, 5, 6, matrisA, 20, ipiv );
olmalıdır. Gösterimi basit tutmak için işleme giren matrisi sol üstten seçtim aksi halde parçanın ilk elemanının bellek adresini bulmak için matrisA ile bazı pointer aritmetiği işlemlerine girmem gerekecekti. Dikkat edilmesi gereken bir konu daha var. LDA, LAPACK_ROW_MAJOR seçildiği için bir sonraki satırın aynı sıradaki elemanı arasındaki uzaklığı belirtiyor. Eğer LAPACK_COL_MAJOR seçilseydi o zaman bir sonraki sütunun aynı sıradaki elemanı arasındaki uzaklığı gösterirdi. Özetle LDA büyük bir matrisin parçaları üzerinde işlem yaparken, asıl matrisin satır veya sütun sayısı olmalıdır. Bunlardan hangisinin olması gerektiğini LAPACK_ROW_MAJOR veya LAPACK_COL_MAJOR argümanı belirler.
