CUDA ile Matris Çarpımı – 2

Matris Çarpımı CUDA – 2

Bir önceki yazımızda CUDA ile matris çarpımının nasıl yapılabileceğini en basit hali ile görmüştük. Burada geçen sefer kaldığımız yer olan hafıza erişim optimizasyonuna devam ediyoruz.

En son matMul2 isimli kernel içerisinde toplam alınan satırda hafızaya yapılan yazmaları azaltmıştık.

Paylaşık Hafıza (Shared Memory)

Paylaşılmış olmayan hafıza mı vardır dediğinizi duyar gibiyim. Global Hafıza da paylaşılmıştır elbet, herkes yazabilir herkes okuyabilir. Fakat burada bir nVidia terimi olarak geçen “paylaşık hafıza” bir blok içerisindeki threadler arasında paylaşılmıştır. Yani ancak aynı blok içerisindeki threadler aynı yeri okuyabilir ve aynı yere yazabilir.

Bir bloktaki threadler paylaşık hafızaya yardımlaşarak yükledikleri verileri kullanabilirler.

2. Optimizasyon (Paylaşık Hafıza)

Yapılan optimizasyonda paylaşık hafıza kullanılarak global hafızadan okuma sayısı düşürülmüştür.

Zamanlama

512 x 512 boyutlarında iki matrisin çarpım sürelerini karşılaştırabiliriz.

Versiyon Hız (ms) Hızlanma
Host 9619 1x
MatMul1 509 19x
MatMul2 183 53x
MatMul3 72 135x

Referanslar

Proje dosyasını açabilmek için farklı kaydederek uzantısını 7z yapınız.
Kaynakları indir: CUDA matris carpimi 2

CUDA ile Matris Çarpımı – 1

Matris Çarpımı CUDA

Matrislere iki boyutlu vektörler dersek vektörleri tek boyutlu değişkenler gibi görmek, değişkenlere de “boyutsuz” benzetmesini yapmak yerinde olacaktır. Matris nedir sorusunun cevabını aramayacağımız bu yazımızda matris çarpımını C dilinde CPU üzerinde ve GPU üzerinde en basit haliyle nasıl yapabileceğimizi göreceğiz.
Örnek kodlarda kodun basitliği açısından kare matrisler kullanılmıştır.
A: M x N ( M satır, N sütun),
B: N x R
C: M x R boyutlarında seçilmiştir.
Örnek kodlarda matrislerin boyutları #define ile tanımlanacaktır fakat bu yazımızdaki örneklerde M = N = R olarak seçilmelidir.

Tanımlar:

A_{i,j} A matrisinin i. satırının j. sütunundaki eleman
C = AB
C_{i,k}=\sum^N_{j=1}A_{i,j} B_{j,k}
Bu yazımızdaki örnekler:
  • C ile yazılmış naif matris çarpımı
  • Optimize edilmemiş GPU kerneli
  • 1. Optimizasyon
Optimizasyon çok basit olarak yapılmıştır öyleki optimizasyon demeye dilim varmıyor. İleriki yazılarda kernel adım adım optimize edilerek GPU’nun suyu sıkılacaktır.

Naif CPU Matris Çarpımı

C matrisinin her elemanı için A matrisinin N sütunu ve B matrisinin N satırı çarpılarak toplanır. M x N x R kere loop döner ( O(n^3), n kare matrisimizin bir boyutu )

Bu algoritma, matris çarpımının “text-book” çözümüdür.

Optimize Edilmemiş GPU Kerneli

deviceMultiply1 fonksiyonunun çağırdığı matMul1 kerneli, performans hakkında en ufak tedirginlik hissedilmeden yazılmış bir matris çarpım kernelidir. deviceMultiply1 girizgah fonksiyonunda GPU üzerinde hafıza ayrılmakta, matrisler GPU hafızasına kopyalanmakta ve hesaplanan sonuç CPU hafızasına alındıktan sonra GPU üzerinde ayrılan hafızalar bırakılmaktadır. Sonucun hesaplanması matMul1 kernelinde gerçekleşir.

deviceMultiply1 fonksiyonunda matMul1 kerneli çağrılmadan önce grid ve block boyutları hesaplanır. GPU üzerinde çalışacak thread sayısının C matrisindeki eleman sayısı kadar olmasını sağlayacak değerler hesaplanır. Daha sonra matMul1<<< dimGrid, dimBlock >>>( d_C, d_A, d_B, m, n, r); satırıyla kernel çağrılır.

Memory access sayısının düşürülmesi, memory coalescing, tiling gibi hiçbir teknik kullanılmadan yazılmış olan naif CUDA matris çarpım kernelini birlikte inceleyelim.

1. Optimizasyon (Hafıza Erişimi)

matMul1 kernelinde yapılacak çok ufak bir değişiklik ile performans arttırılabilir. Bu değişikliğin getirdiği hız kazancı, GPU’da hafıza erişim işlemlerinin CPU hafızasına erişimden çok daha yavaş olması ile gün ışığına çıkmaktadır.
matMul1 kernelinde d_C, C matrisinin GPU hafızasındaki yerine işaretçidir(pointer). Bu işaretçinin kullanıldığı her yer GPU hafızasına erişmek demektir. İşaretçi başta sıfırlanmıştır ve toplama işlemi işaretçi üzerinde yapılmıştır. += operatörü, d_C[ Idx ] adresinin önce okunması, toplam işleminden sonra da tekrar yazılması demektir. Yani d_A ve d_B’yi saymazsak for döngüsünde toplam 2*n kere hafıza erişimi yapılmaktadır. Halbuki C matrisinin değerini hesaplayıp tek seferde hafızaya yazabiliriz.
matMul2 kernelinde yaratılan değişken ile global GPU RAM erişimi tek thread başına 1 adet yazmaya (d_A ve d_B hariç) düşmüştür.

Zamanlama

Yapılan bu ufak değişiklikle kazanılan 128 x 128 boyutlarında iki matrisin çarpım sürelerini karşılaştırabiliriz.
matMul1 ve matMul2 zamanlaması

Proje dosyasını açabilmek için farklı kaydederek uzantısını 7z yapınız.
Kaynakları indir: CUDA matris carpimi 1

Visual Studio 2010 ile CUDA SDK 3.2 Projesi Yaratmak

CUDA Projelerine Giriş

For English: CUDA 3.2 on VS2010 in 9 steps

Bu kısa yazımızda CUDA ile programlamaya giriş yapmanın zevkini yaşayacağız. CUDA ile ekran kartındaki GPU üzerinde çalışacak programlar yazabilirsiniz. Bu programlar GPU’daki özel hardware threadler üzerinde paralel olarak çalışırlar. Detaylı bilgi için tıklayınız.

4. adımdaki dosyalar yok ise indiriniz: CUDA 3.2 Build Rules

Aşağıdaki adımları izlemek için gerekenler:

  • Visual Studio 2008 kurulumu
  • Visual Studio 2010 kurulumu
  • nVidia CUDA 3.2 SDK kurulumu
  • Genel Visual Studio menüleri ve kullanımı

9 Adımda CUDA Derleme

    1. Yeni bir Win32 Console Application yaratın.

    1. Projeyi yaratırken empty project seçeneğini işaretleyin. Boş bir projeye kaynak dosyasını kendimiz ekleyeceğiz.

    1. cu uzantılı kaynak dosyası ekleyin. Bu dosya nvidianın compiler driverı olan nvcc.exe dosyası tarafından işlenerek vs2008 ile gelen versiyon 9 C derleyicisine gönderilecektir.

    1. CUDA 3.2 SDK’nın kuruluşuile gelen Build Customizationların yerinde olduğu kontrol edilir.

    1. Platform Toolset v90 olarak değiştirilir. Visual Studio 2010, CUDA dosyalarını derlemeyi destekler fakat altyapı olarak bir önceki versiyon C derleyicisini kullanır. Dolayısıyla Visual Studio 2008’in makinada kurulu olmasına ihtiyaç vardır.

    1. Eklediğimiz CU dosyasının özelliklerinden Item Type olarak CUDA C/C++ seçilir.

    1. Project Build Customizations olarak CUDA 3.2 seçilir. Bende CUDA 3.1 SDK’sı da yüklü olduğu için o da seçilebiliyor.

    1. CUDA’nın kullandığı kütüphaneler eklenir( cuda.lib ve cudart.lib ).

    1. Yaratılan CU dosyasından CUDA fonksiyonlarının çağrılabildiği ve derlenebildiği görülür.

Bézier Curves in Bing Silverlight Maps

Bézier Curves in Bing Silverlight Maps

Connecting geographical locations on Bing Maps with Bézier curves.

Download sample: MapBezier Silverlight Application

Introduction

Bing Maps Silverlight Control library has a MapPolyline class for showing connected points on the map. I wanted my points to be smoothly connected but there wasn’t out-of-the-box support so I developed a custom control deriving from MapShapeBase class.

Bézier Curve on map.

Background

Every developer who messed with Expression Blend or Gimp long enough knows by experimentation how a Bézier curve behaves. Basically a cubic Bézier curve has an initial point (P1), two control points (B1, B2) and a final point (P2).

Bézier Curve

Cubic Bézier Curve

The formula that defines a cubic Bézier curve is:

Cubic Bézier Curve

where t is in the interval [0,1].

Terms multiplying P1, B1, B2, P2 are called the basis functions for the cubic Bézier. Our points determine how much of these basis functions does the curve contains.

Basis Functions for Cubic Bézier Curve

Cubic Bézier Basis Functions

The thing is we need to calculate coordinates of the control points such that our points of interest are on the curve.

Polyline, Bézier, Catmull-Rom

Polyline, Bézier, Catmull-Rom

So how can we calculate these control points? After researching(read googling) I have landed on cardinal splines and Catmull-Rom splines. It appears that every control point of a Catmull-Rom spline is on the curve and it is also a Bézier curve which means we can use it as PathGeometry with a Silverlight Path object.

Calculating Control Points

If we rewrite formulas from the cardinal splines page as the following, we can easily calculate control points.

Point Derivative

Control Points:

Control Points

Fitted Bézier with Control Points Visible.

Fitted Bézier with Control Points Visible

Using the code

The method that calculates the Bézier Points is as the following. GetB1 and GetB2 are straight forward implementations of the aforementioned formulas.

You can use the MapBezier class just like MapPolyline and MapPolygon classes in your silverlight XAML file. See the attached sample for an example silverlight application.

Note: Bing Maps Silverlight SDK is needed for compiling the sample application.

Points of Interest

It was fun messing with the Bing Maps Silverlight Control toolkit and I hope MapBezier is what you are looking for.

If you liked this, vote the article at http://www.codeproject.com/KB/silverlight/MapBezier.aspx

Koordinat Dönüşüm Fonksiyonları

Ne?

İlgilendiğimiz koordinat sistemleri, dünya yüzeyindeki herhangi bir noktayı karışıklığa mahal vermeden tarif etmeye yarar. Koordinat sistemi denildiğinde akla gelen ilk sistem x-y koordinatları olarak da hatırlanan bir 2 boyutlu dik koordinat sistemi olan kartezyen koordinat sistemidir.

Peki bahsettiğimiz gibi eğer bir koordinat sistemi karışıklığa mahal vermeden bir noktayı tarif edebiliyorsa farklı koordinat sistemlerine ne gerek var sorusu akla gelebilir, gelmelidir. Farklı koordinat sistemlerine olan ihtiyaç tamamen kullanım yeri ve kullanım kolaylığı ile ilgilidir.

Nedir?

İlgilendiğimiz koordinat sistemlerinden bahsetmeden önce harita yapımının bir aşaması olan izdüşüm (projection) almaktan bahsetmekte fayda var.

İzdüşüm Nedir?

En basit haliyle izdüşüm, bir noktanın bir yüzey üzerine düşen izidir. Biraz detaylandırırsak belli bir merkezden gelip bir noktadan veya bir noktalar kümesinin her bir noktasından geçen doğruların izdüşüm tarafından tanımlanmış olan yüzeyi kestiği noktaların geometrik yeridir diyebiliriz.

Bu biraz karışık olduysa merkezde bir ışık kaynağı olduğunu ve izdüşüm yüzeyinin bizim istediğimiz şekilde bir perde olduğunu düşünün. İzdüşümünü almak istediğimiz noktaları ışık kaynağı ile perde arasına koyarsak, perde üzerinde izdüşümünü elde ederiz.

Bizim ilgilendiğimiz izdüşümü, dünyanın (noktalar) üzerinde geçirilmiş bir silindire (perde), dünyanın merkezinden (ışık kaynağı) olan izdüşümüdür.

Mercator İzdüşümü

Türkiye’de de tapu, ruhsat, askeri ve benzeri alanlarda en çok kullanılan koordinat sistemi UTM şeklinde kısaltılan Universal Transverse Mercator (Evrensel “Uzun eksene dik” Mercator) sistemidir.

Transverse Mercator

Transverse Mercator

KorTransLib

Koordinat Dönüşüm Kütüphanesinin yazılış amacı ileride düşünülen projelerde kullanılmak üzere olduğundan sadece 6 derece dilimli UTM, enlem-boylam ve derece-dakika-saniye olmak üzere üç adet koordinat sınıfı olmakla beraber istenilen koordinat sistemleri kolaylıkla eklenebilecek şekilde tasarlanmıştır.

Coordinate Sınıfı ve Türetilmiş Sınıfları

Coordinate Sınıfı ve Türetilmiş Sınıfları

Koordinat Sınıfları

  • DecimalDegreeCoordinate
  • DegreeMinuteSecondCoordinate
  • UTMCoordinate

DECIMALDEGREECOORDINATE

Ondalık dereceli gösterim. Matematiksel hesap yaparken kolaylıkla kullanılır.

DEGREEMINUTESECONDCOORDINATE

Derece, dakika, saniye sisteminde gösterim.

UTMCOORDINATE

6 derecelik dilimli Universal Transverse Sekant Mercator izdüşüm sistemindeki koordinatlar. UTM koordinat sisteminde sağa ve yukarı değerden başka dilim bölgesi (UTM Zone) değeri de vardır. Bu değer sağa ve yukarı değerlerinin anlamını tanımlar. Türkiye’nin UTM dilimleri aşağıda gösterilmiştir.

Türkiye UTM Dilimleri

Türkiye UTM Dilimleri

Koordinatlar Sınıfları Arasında Matematiksel İşlemler

Örnek olarak herhangi iki koordinatın orta noktasının bulunmasını aşağıdaki şekilde yapabilmemiz için kullanılan TypeConverter sınıfı, kodun yazılmasında da kullanılmasında da büyük kolaylık sağlamıştır.

Yukarıda gösterildiği gibi koordinat dönüşüm işlemi sadece bir değişkene atama yapılmadan önce yapılıyor.

İki koordinatın ortalaması, karşılıklı iki köşe koordinatı bilinen bir paftanın merkezini bulurken kullanılabilir.

TypeConverter Kullanımı

  1. System.ComponentModel içerisindeki TypeConverter sınıfı miras alınarak CanConvertFrom, ConvertFrom, CanConvertTo, ConvertTo fonksiyonları yazılır.
  2. Yaratılan TypeConverter mirasçısı sınıfının hangi tipi dönüştürdüğü Nitelik(Attribute) sınıfları kullanılarak tanımlanır.
  3. TypeDescriptor sınıfı yardımı ile kullanılır.

Project Euler Problem 53

Project Euler 53. Problem ve Çözümü

53. Problem’de 1 <= n <= 100 ve C(n,r) > 1000000 koşullarını sağlayan kaç adet C(n,r) şeklinde kombinasyon olduğunu bulmamız isteniyor.

Bu problemin çözümü için akla gelen pek çok yol içerisinden hızlı ve basit olan bir yol ile problemi çözmeye çalışalım. Problemimizde kombinasyonu hesaplayarak 1000000’dan büyük olup olmadığını sınamak yerine kombinasyonun ve kontrol değerinin logaritmasını alarak logaritmaları karşılatırıyoruz.  Öncelikle bazı kombinasyon formüllerini hatırlayalım:

Faktöriyel ve kombinasyonun formülü

ve bazı logaritma formüllerini:

Çarpımın logaritması

Bölümün logaritması

Yukarıdaki formülleri kullanarak faktoriyelin logaritmasını aşağıdaki gibi yazabiliriz:

Faktöriyelin logaritması

Bu fonksiyon project euler ile doğrudan ilgili olmasa da gammaln fonksiyonunun analogudur.

Böylece kombinasyonun logaritmasını da aşağıdaki şekilde yazabiliriz:

Kombinasyonun logaritması

Buraya kadar yazdıklarımızdan çıkarmak istediğimiz sonuç olan problem 53‘ün logaritmalar ile ifadesini aşağıda görebiliriz:Logarithm of one million

Euler Problem 53

Böylece logaritmanın da yardımı ile kendimizi big-integer çarpım, bölüm ve karşılaştırmalarından kurtararak işimizi double toplamlara ve karşılaştırmalara indirgemiş olduk. Aynı zamanda faktoriyelin logaritma fonksiyonunu tasarlarken memoization tekniğini de kullanabiliriz. Problemi çözen örnek C# kodu aşağıda verilmiştir: