수퍼컴퓨터를이용한병렬몬 테카를로시뮬레이션 유현곤 (Hyun-GonRyu) Dept. of Mathematics Yonsei Univ. Email : yhgon@yonsei ac kr Email : yhgon@yonsei.ac.kr 2008 년 7 월 11 일금요일
목차 1. 기초금융공학, 몬테카를로시뮬레이션 2. 병렬컴퓨팅이란? 3. 병렬몬테카를로시뮬레이션 4. OpenMP 를이용한병렬몬테카를로시뮬레이션 5. MPI 를이용한병렬몬테카를로시뮬레이션 6. 결론및추가적연구사항 (CUDA)
1. 기초금융공학, 몬테카를로시뮬레이션
Introduction for Option Pricing European Vanilla Call Option 1. asset dyanmics 2. payoff GBM : [ S K ] + T ds = μs dt + σ S dw t t t t Payoff [ S K ] + T K ST 3. option price is Q rτ + T C =E [ e [ S K] ] t
Introduction for Option Pricing European Vanilla Call Option Q r τ + C =E [ e [ S K ] ] t T 4. By Ito s lemma Q T 1 2 ( r 2 σ ) τ+ σw = BS-PDE : S S e τ ST = S e 0 0 1 2 2 ( r σ ) τ+ σ τn(0,1) 5. By Feynmann-Kac Theorem f + rxf + σ x f = rf 1 2 2 t x 2 xx 6. Closed Form Solution r C = Ke τ t S0Φ( d1) Φ ( d 2 ) x 1 x 2 2 Φ ( x ) = e dz 2π d 1 2 ln( S 0 K ) + ( r σ 2) τ = 2 1 + d = d σ τ σ τ
PDE method for ELS 2 2 2 f f f 1 2 f f 1 2 + rx + ry + f 2σx + ρσ 2 xσ y + 2σ y = rf 2 t x y x x y y With Boundary Value& Terminal Payoff Domain for 2-Star n-chance Stepdown ELS Step down Worst Perf.
Monte Carlo Simulation Law of Large Number : 많은시행을통한평균을통해기대값계산 1 N 1 E[ e rτ [ S K] + ] = lim e rτ [ S K] T T N N i i + We can generate any process S T from given dynamics dst = μstdt + σ StdWt ds = μs dt + f ( Y ) S dw t t t t t
Monte Carlo Simulation in Finance 상품설명서 Easy for even junior quant Simple MC coding 몇시간이내 Excel/VBA, Matlab, C/C++ Library, Rand() 함수이용 debugging Hedging ( 매일 ) 가격, Greeks 요청 Computer Hedging ( 매일 ) 가격, Greeks 출력 몬테카를로시뮬레이션의구현상품발행시상품거래시 Pricing Pricing Overview of Greek Greeks Search Hedge Simulation VaR계산, 위험관리조기상환확률계산
Pseudo Code for MC Parameter Inputs For-Loop of N { XT simulation with Brown Motion Path Transformation (Box-Muller, Inversion) RNG (LCG32, RAND48, MWC, MT19937) Compute Payoff } Sum Average (1/N) Pi Print tresults IF, Max function (payoff, Prepayment)
Monte Carlo Simulation Why is MC useful? 적용이쉽다. Easy to imply 복잡한상품의평가가가능하다. Complex structure : path-dependence, prepayment, complex payoff, etc. 다차원 ( 다중자산 ) 문제를풀수있다. High demension : multi-asset problems (n>4) MC is the only solution (or sparse grid method)
Convergence of Random Number (excel) Vanilla Call Option Monte Carlo Simulation N = 50000 번이상, 10 만번정도가적당 6.0000% 01:43.68 6.000% 01:43.7 4.0000% 2.0000% 0.0000% 01:26.40 01:09.12 4.000% 2.000% 0.000% 01:26.4 01:09.1-2.0000% 00:51.84-2.000% 00:51.8-4.0000% -6.0000% -8.0000% 00:34.56 00:17.28-4.000% -6.000% -8.000% 00:34.6 00:17.3-10.0000% 00:00.00-10.000% 00:00.0 500 0 1,000 0 5,000 0 10,000 0 20,000 0 50,000 0 100,000 0 400,000 0 1,000,000 0 500 0 1,000 0 5,000 0 10,000 0 20,000 0 50,000 0 100,000 0 400,000 0 1,000,000 0 Pseudo Random Number LCG 32 Quasi Random Number 1D Halton Sequence
European Call Option 환경 : Intel CPU 3Ghz 14.6 10.00 14.5 9.00 14.4 8.00 14.3 7.00 14.2 6.00 14.1 5.00 MC 14 13.9 4.00 3.00 Closed 시간 13.8 2.00 13.7 1.00 13.6 0.00 1000 13000 25000 37000 49000 61000 73000 85000 97000 109000 121000 133000 145000 157000 169000 181000 193000 205000 217000 229000 241000 253000 265000 277000 289000 Simple MC 의경우 12 만번 Simulation 하는것이적당 2 sec 하지만계산시간의제약으로약 5 만번정도돌림
250 time step Path simulation 14.8 14.6 14.4 환경 : Intel CPU 3Ghz Black-Scholes Clsoed : 14.23124493415722500 Monte Carlo Simulation : 14.26851272608143400 Error : 0.03726779192420970 14.2 14 13.8 13.6 13.4 13.2 13 12.8 1000 11000 21000 31000 41000 51000 61000 71000 81000 91000 101000 111000 121000 131000 141000 151000 161000 171000 181000 191000 Simple MC 의경우 10 만이상이 Simulation 하는것이적당 7.65sec 하지만계산시간의제약으로약 5 만번정도만돌리면...
Pseudo vs. Quasi 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 02 0.2 04 0.4 06 0.6 08 0.8 0 1 0.2 0.4 0.6 0.8 1 본 발표에서는 생략 3 4 2 3 2 1 1 0-4 -3-2 -1 0-1 -2-3 -4 LCG32 1 2 3 0 4-4 -3-2 -1-1 0 1 2 3 4-2 -3 3-4 Halton Sequence
Monte Carlo Simulation Why did MC not used commonly? Most contracts(95%) can be solved without MC. (closed, FDM) Computation Time : Too slow to get accurate results -1 minutes for each pricing. 4 minutes for greeks 100 contracts : 7 hours to comutes (1 day : 6 hours) 시나리오분석용 Greek plot 그리려면하루종일걸림 Unstable sensitivity : non-smooth greeks plot rather than FDM Convergence Ratio Chart Speed chart
Reasonable Solutions to speed up MC 1. New theory in convergence Malliavin Calculus, Operator Technique, Asymptotics 2. Fast pseudo & quasi-rngs, Transformation 3. Control Variate, Variance Reduction 4. Using Powerful Computer 5. Parallel Computing HPC(use many CPUs) - OpenMP(SMP), MPI (Cluster) Alternative Method - IBM Cell BE, ClearSpeed, GPGPU(CUDA)
2. 병렬컴퓨팅이란?
HPC vs. distribution/load balancing 병렬처리는크게 3 가지로나뉘는데 하나의방법은 HPC로, 하나의작업을여러개로쪼개어계산하도록함으로써계산시간을단축시킴 다른하나의방법은여러개의작업을여러개의서버에분산처리시켜전체작업시간을단축시키는기법 웹서버, 게임서버에사용되는로드벨런싱도분산처리기법의하나임 계산의관점에서병렬처리는 HPC 를의마함 분산처리의예하나의자산당 200초걸리는자산 200개가있다. 총 40000초걸림 200개의서버에서각자산정보를보내서돌린결과를받음 병렬처리할필요없이 200초만에모든 VaR 계산 Quant에게유리 ( 다시코딩할필요없음 ), 회사입장비용문제 ( 수십억원 )
병렬처리 (HPC High performance Computing) 하나의작업을여러개의 CPU(Cores) 로작업을실행시켜계산시간을줄이는방법 Wall Clock Time Task 를병렬처리에의해처리하는데걸린전체시간 Wall Clock Time 1 Original Jobs Parallel Overheads Wall Clock Time 2
KISTI 수퍼컴퓨터 4호기 Tachyon 세계수퍼컴퓨터랭킹 130 위 구분 내용 모델명 SUN Blade 6048 블레이드노드 CPU 성능 : 24TFlops(Rpeak) 노드 : 188 개 ( 컴퓨팅 ) AMD Opteron 2.0GHz 4 개노드당 : 16 Core (Quad Core) 총 : 3,008 개 ( 전체 ) 메모리 32GB( 노드 ) 스토리지 207TB( 디스크 ), 422TB( 테이프 ) 노드간네트워크 Infiniband 4X DDR
금융권 HPC 용 16 코어 SMP 머신 Blue Box : HP, IBM, SUN 등도입비용 1700 만원 ( 프로모션견적 ) 가량 HP ProLiant DL580G5 Rack Type 4U case 사용인텔제온 E7340 CPU 4 개 : 16 코어메모리 8GB White Box : 비슷한도입비용으로동일스펙에 Memory 48GB정도로확장가능단, 유지보수등의문제발생 Node당 100Gflops 유지 Tachyon 1 노드와유사한성능 Tachyon 전체시스템대비 1/250 배성능 N 개의노드추가를통해 MPI 머신으로성능향상가능대형시스템은 Tachyon 과유사한 blade 기반이유리
수퍼컴퓨터 4 호기 Tachyon on KISTI 구분 내용 계정 SRU 100시간 : 100만원 ( 기업 ) SRU 10 시간 : 10 만원 ( 학생 ) 접속환경 ssh 접속 30분 실행환경컴파일러병렬라이브러리수치라이브러리 CentOS 4.4 - Linux 환경 GCC, PGI Compiler, Intel Compiler OpenMPI IMKL, IMSL 등 계정사용을통해 HPC 도입효과를미리검토해볼수있음
작업노드 ( 프로그래밍, 디버깅 ) Tachyon Server 는 SSH 를통해접속해야하고 30분까지로접속시간이제한되어있음 따라서, 다음의시스템에서 OpenMP, MPI 코딩을작성, 디버깅하고최종결과물은 Tachyon에서실행하는방법을사용한다. 1. Windows XP 환경 (PC) : GUI 환경에서 programming MS Visual Studio 2005+ MPICH2 2. Linux 환경 (server) : linux에서실행테스트 ICC 10.1 + OpenMPI 3. Tachyon Server : 최종코드를 job_batch 실행
Windows 환경에서 PuTTY 를통해접속 SSH 원격접속환경
병렬환경 SMP 머신 : OpenMP 프로그래밍 16coe core 지원 : 2 천만원미만 (16 배성능향상보장 ) Cluster 머신 : MPI 프로그래밍 64 core 지원 : 1 억원미만 : (60 배성능향상보장 ) OpenMP, MPI : ( 무료 ) IMSL, IMKL : 병렬라이브러리 : ( 유료 )
Parallel FDM, FEM (matrix) Main For(i=1,i<N,i++){ PE1 PE2 PEn-1 PEn............................ 인접데이터의통신에의한 bottleneck
3. 병렬몬테카를로시뮬레이션
Parallel MC simulation Main Divide & Conquer For(i=1,i<N,i++){ PE1 PE2 PEn-1 PEn For(i=1,i<N, M){... } For(i=1,i<N, M){... } For(i=1,i<N, M){ For(i=1,i<N, M){...... } } 인접데이터의통신이필요없음 : MC 병렬 scailability 가좋다
Parallel MC with OpenMP (SMP machine) 10 만번실행 Excel/VBA 1. Start Function C/C++ DLL Host 2. Receive Input Parameters 4. Initialize Dynamic Creation Seeds wallclock time Core 개수에따라실행시간다름약 12 초 Errors < 10 OpenMP 2~8 Multi Threads 5. Parallel MT Moro s Inversion 6. Stock path 4 10 Option Price Greeks 7. Show Results Price, Greeks, time
Parallel MC with MPI (Cluster Machine) 10 만번실행 Excel/VBA 1. Start Function C/C++ DLL Host 2. Receive Input Parameters 4. Initialize Dynamic Creation Seeds wallclock time 병렬 CPU 갯수에따라실행시간소요 0.05 초 ~30초 Errors < <10 MPI 2~1024 Multi Threads 5. Parallel MT Moro s Inversion 6. Stock path 4 10 Option Price Greeks 7. Show Results Price, Greeks, time
Parallel Monte Carlo Simulation in Finance 상품설명서 Easy for even junior quant Simple MC coding Need veteran parallel coder or library Single Code가잘구축되어있고, 병렬라이브러리가구축되어있다면신상품의병렬화는매우쉽다. (openmp 이용시 ) 몇시간이내 Use own Parallel library & similar codes parallel l MC coding 매일반복계산 Hedging ( 매일 ) 가격, Greeks 요청 병렬실행병렬실행병렬실행병렬실행 Hedging ( 매일 ) Tachyon 급수퍼컴퓨터 0.05 sec 이내결과출력 가격, Greeks
LCG32 X ( a X c) mod m n+ 1 = n + M a c Visual C/C++ 2 32 214013 2531011 GNU C 2 32 69069 5 Unix (LCG48) 2 48 25214903917 11 ANSI C 2 32 1103515245 12345 Numerical Recipes 2 32 1664525 1013904223 Period : 2 32-1 = 4294967295 = 4.2 * 10 9 약 42 억개
Split Method for (i=0; i<n-1; i++ ) for (i=istart; i<iend; i++) CPU ID 0 CPU ID 1 CPU ID 2 34
Multiseed Method for (i=0; i<n-1; i++ ) srand48p(123+(1+cpuid)*45); for (i=0; i<chunksize-1; ) CPU ID 0 CPU ID 1 CPU ID 2 35
Leap Frog Method for (i=0; i<n-1; i++ ) for (i=istart; i<n-chunksize-1; CHUNKSIZE ) CPU ID 0 CPU ID 1 36
Problems in Parallel MC LCG : 2^32-1 42 억개 37
Suitability for Parallel Method Multiseed Leapfrog Splitting DC D.C JAh J.Ah. LCG OK OK OK X - RAND48 OK OK OK X - MCG OK OK OK X - Lagged Pibonachi OK OK OK X - MWC OK OK OK X - MT19937 OK X X OK OK Well RNG OK?? OK 38
Dynamic Creation for (i=0; i<n-1; i++ ) Complex dividing for non-overlapping cycle for (i=0; i<chunksize-1; ) CPU ID 0 CPU ID 1 CPU ID 2 39
for (i=0; i<n-1; i++ ) Jump-Ahead Method Divide sequence with GF(2) polynomial for (i=0; i<chunksize-1; ) CPU ID 0 CPU ID 1 CPU ID 2 40
smt architecture of MT19937 Initial ii Seed 0 1 2 i i+1 smt Generator Extractor U1 i+m i+m+1 623
MT19937 병렬화 Shared pmt architecture of MT19937 Initial Seed U2 Generator Extractor 0 Generator 1 Extractor 2 U1 i i+1 Parallel pmts With Large StepSize 1cycle generate 624 RNGs i+m i+m+1 623 Generator Extractor Un Share Memory
Parallel smt MT19937 Each Cores execute smt independently We will use this model 0 1 2 0 1 2 0 1 2 i i+1 i i+1 i i+1 i+m i+m+1 i+m i+m+1 i+m i+m+1 623 623 623
Parallel Quasi RNG For Loop Nsim simulation For Loop for T simulation of path Method 1 base p 에의한병렬화 Split Halton with Nsim /16 -- 16 is # of core Core has own Base 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59 Scramble each Sequence with T times Method 2 : 준주기성을고려한 n 병렬화 select large prime such as 23 Compute Split K= Nsim/16 * 24 n=(k*coreid + i) Parallely Compute Halton(23,n) Scrample each Sequence with T times
실제병렬화 신상품, 상품설명서 Excel, Matlab, C/C++ 코딩 싱글코드알고리즘 디버깅 Profile 계산로드가많이걸리는부분파악 병렬화시스템결정 SMP, MPI, CUDA, CSCN, etc 병렬코딩 OpenMP, MPI, Matlab, etc. 실행및디버깅 결과출력
금융문제의병렬화기법 사용방법1 총 M회시뮬레이션을병렬처리각 N 개의 Core 가 M/N 회시뮬레이션실행각 Core가독립적으로 Boxmuller 및 RNG 생성 1/N으로계산시간단축 장점 : single code의 90% 이상그대로사용 방법2 1. RNG 를병렬화하여미리생성 2. N개의시뮬레이션을병렬처리 3. N개의시뮬레이션평균 조기상환시 RNG 생성시간낭비 방법3 Table Technique RNG 를병렬화하여미리생성 (1 회 ) memory에 load하여필요시사용각상품별로 (K 회 ) M/N회의시뮬레이션을병렬처리 모듈화작업시 single code 를그대로사용가능 많은상품을계산해야하는시스템에서적합 1/K*N 으로계산시간단축충분한메모리가확보및분산처리 Memory access bottleneck 발생
Single Code 병렬몬테카를로시뮬레이션예제 1. Full-path 유럽형 Vanilla Call 옵션 : reference & benchmark 용 S : 100, K : 100, r : 0.0303, v : 0.3, T : 1 2. ELS pricing : 삼성증권 1909 호주식연계증권 3 년만기 2 star : POSCO 일반주, S-Oil 12 chance : 디지털옵션 Double Barrier Down barrier : 장중체크 ( 둘중하나라도하락한계 ) Up barrier : 매일종가체크 ( 두개다상승한계 ) Step-down : 없음, 동일한조기상환조건 (Digital Option 형태 ) 지급 : 조기상환시 +2 영업일
ELS 예제의경우실무경험부족으로 ( 본인 ) - 실무파트의현실과경험부족으로실무에서사용하는코드와다를가능성이있음 - 이자율, 날짜계산, vol 설정등파라미터설정에대한노하우부족 - 실제업무용시스템과연구용코드의괴리가능성 - 실제포트폴리오에서의 Greeks 관리의복잡성등 아직완성된프로젝트가아님 - 체계적노하우, 코드관리능력등부족
OpenMP 를이용한병렬화
OpenMP 개념 메모리를공유한 SMP 머신 (core2 duo 등의형태 ) 에서병렬코딩함. 16 코어머신을이용하면 MC 의경우약 16 배의속도향상기대 C/C++ 언어와 Fortran 언어를지원함 Windows 환경은 Visual Studio 2003 이상에서기본지원 Unix 환경의경우최신버전기본지원 Linux 의경우 GCC 2.4.2 이상에서지원 OpenMP 명령어를통해컴파일러가자동으로병렬화 For Loop 를쉽게병렬화할수있고, 각 Core 가메모리를공유하기때문에병렬코딩이매우편함
SMP 한대의컴퓨터 Shared Memory Core0 Core0 Core0 CoreN-1 CoreN 각코어가모두하나의메인보드위에있기때문에모두같은메모리를사용특별히, 통신코딩을할필요가없음 OpenMP 코딩의편리함
OpenMP 개념 Fork Master Thread 병렬화영역 Join Fork 병렬화영역 Join
OpenMP 병렬화예약어들 omp_set_num_threads(16); 총 16개의 thread 사용설정 totaln=omp_get_num_threads; 전체병렬화개수파악 tid=omp_get_thread_num(); 각병렬프로세서번호인식 #pragma omp 지시어 #ifdef_openmp 순차프로그래밍에서도사용가능하도록프로그래밍 parallel l for critical private() shared() schedule()
OpenMP 예제 1 #include <stdio.h> #include <omp.h> int main (int argc, char *argv[]) { int nthreads, tid; omp_set_num_threads(4); #pragma omp parallel private(nthreads, tid) { tid = omp_get_thread_num(); nthreads = omp_get_num_threads(); printf(" Hello World from Thread %d of %d \n", tid, nthreads); } return 0; } 병렬화영역
OpenMP 병렬화방법 for(i=0; i<m; i++) { a[i] = b[i*n]*c[0]; for( j=1; j<n; j++) a[i] += b[i*n+j]*c[j]; } 자동병렬화 #include <omp.h> #pragma opm parallel for shared(m,n) private(i,j) for(i=0; i<m; i++) { a[i] = b[i*n]*c[0]; [0] for( j=1; j<n; j++) a[i] += b[i*n+j]*c[j]; }
수퍼컴퓨터 Tachyon OpenMP job batch scheduler 282% [e063rhg@tachyond e063rhg]$ vi a.sh #!/bin/bash #$ -V #$ -cwd #$ -N openmp_job #$ -pe openmp 4 #$ -q small #$ -R yes #$ -wd /work01/e063rhg/ #$ -l h_rt=00:01:00 000100 #$ -M myemailaddress #$ -m e export OMP_NUM_THREADS=4 /work02/e063rhg/omp.exe >result.txt
Job scheduler 에서대기중인사람들 예상실행시간 5sec 의 job 을 dir 40 분기다려야함 Scheduler 조정을통해제일적게기다리는편
Tachyon idle core 현황 6 월 30 일부터현재까지계속 job 이 waiting 상태로, 수퍼컴퓨터사용의의미가없음 OpenMP, MPI test 불가
Job Schedule 실행된모습 ================================= Black-Scholes : 14.23124493415723357 1 3-1073749904 566413.67638386972248554 566413.67638386972248554 3-1073749904 1132306.36511257803067565 565892.68872870830819011 3-1073749904 1701503.34236495988443494 569196.97725238185375929 3-1073749904 2254231.80765507556498051 552728.46529011556413025 14.29527750057961200 0.06403256642237842 elapsed time is 44.892820
#pragma omp parallel shared(fsum, N) private(tid,fsum_local,nnn) { NNN = 0; fsum_local = 0; tid = omp_get_thread_num(); printf("%d \t",omp_get_num_threads()); #pragma omp for for ( i=0; i< Nsim ;i++ ) { xt1=s; xt2=s; for(m=0; m<path; m++) { xx1 = rand()/(rand_max+1.0);if(xx1==0.0) xx1=0.0000000000001; xx2 = rand()/(rand_max+1.0);if(xx2==0.0) xx2=0.0000000000001; normal1=sqrt(-2.0*log(xx1 ))*cos(2.0*3.14159265358979323846*xx2 ); xt1= xt1 + r * xt1 * dt + v*xt1*sqrt(dt)*normal1; } oprice=max(xt1-k,0); fsum_local =fsum_local+ oprice; //printf("%f",op); } #pragma omp barrier #pragma omp critical (update_sum) { fsum +=fsum_local; stop=clock(); printf("\n%d \t %20.17f \t %20.17f \n",tid,fsum, fsum_local) ; } } // end of OpenMP results = (double) 1/Nsim * exp(-1* r * tau)* fsum; stop=clock(); htime = 0.001*difftime(stop,start); //windows printf("\n%19.17f \t %19.17f %7.5f \n",results,results-bsp, htime) ; Loop 병렬화부분 결과 Reduction 부분 전체병렬화영역
OpenMP 결과 1. 수퍼컴퓨터에서는 1~16배까지 core 개수를늘리면 scailability가증가하지만유휴 CPU 가부족하여 scailability 를테스트하기어려웠음. 2. CUDAp, MATHp, Windows PC 등에서 1~4 core 기반의 scailability 체크가능
Cluster 전체시스템 Network 하나의노드 하나의노드 Core0 CoreN Core0 CoreN Shared Memory Shared Memory 각노드에서각각프로그램을실행해줘야함 (mpirun) 각노드의메모리는서로공유하지않으므로서로통신해줘야함 MPI 코딩의복잡함
MPI 개념 MPIRUN MPI로병렬프로그래밍된실행명령을각각의 node에복사하여각각의노드가병렬명령을실행할수있도록지원함 MPI setting 에각각의 node 설정을해주어야함 : 관리의복잡성발생 각 Node는네트웍을통해연결되어있으므로, 서로메모리를공유하지않음. 각메모리는서로독립적으로작동하므로 PDE solving, Matrix 병렬화등에서는서로의메모리정보를서로통신을통해업데이트해줘야함 MPI 는병렬화보다는오히려통신프로그래밍을쉽게해준역할을함. OpenMP 와는다르게프로그래머가병렬화작업, 메모리공유를직접해줘야함 MPIEXEC, POE 등의 MPI Launcher 를통해실행해줘야한다. 전용서버의경우 scheduler 대신 plink 를활용하면유용하다.
MPI 예제 1 #include <stdio.h> #include "mpi.h int main(int argc, char *argv[]) { int rank, size; MPI_Init(&argc, &argv); MPI_Comm_rank(MPI_COMM_WORLD,&rank); MPI_Comm_size(MPI_COMM_WORLD,&size); printf("hello, world I am %d rd core of %d core system\n",rank,size); } MPI_Finalize(); return 0;
MPIexec 다음은 mpiexec 를통해 core 를 1 개, 2 개 4 개로확장시켰을때의실행결과를나타내고있다.
Job 분할방법 How do we divide job ( FOR LOOP)? ID=0 ID=1 ID=2 ID=3 ID=0 ID=1 ID=2 ID=3
Loop 분할방법들 단순분할 for (i = 0; i < Nsim/N; i++) {. } 1 3 순환분할 for (i = Tid; i < Nsim; i+= N) {. } Block 분할 블록 - 순환분할 2 i_start = Tid * (Nsim /N); i_end = i_start + (Nsim /N); if (Tid == (N-1)) i_end = N; for (i = i_start; i < i_end; i++) {. } 4 for (i = n1*block*myrank; i < n2; i+=nprocs*block) { } for (j = jid; j < min(ij+block-1); i+=n2n) {.
para_range(n1,n2,n3,n4,n5,n6) 함수 Para_range 함수는크게 3가지의병렬화방법중 Method 2의기법으로 For Loop를균등분할함 void para_range(int lowest, int highest, int nprocs, int myrank, int *start, int *end) { int wk1, wk2; wk1 = (highest - lowest + 1) / nprocs; wk2 = (highest - lowest + 1) % nprocs; *start = myrank * wk1 + lowest + ( (rank<wk2)? myrank : wk2); *end = *start + wk1-1; if(wk2 > rank) *end = *end + 1; } N ai () n1 1 i= 1 i= = ai () + n 2 ai () + + ai () i= n 1 2 + i nk 1 1 N = + Core 0 Core 1 Core k-1
MPI 예제 2 ( MPI 통신 ) #include <mpi.h> #include <stdio.h> #define n 100000 void para_range(int, int, int, int, int*, int*); int min(int, int); void main (int argc, char *argv[]){ int i, nprocs, myrank ; int ista, iend; double a[n], sum, tmp; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); para_range(1, n, nprocs, myrank, &ista, &iend); for(i = ista-1; i<iend; i++) a[i] = i+1; sum = 0.0; 0; for(i = ista-1; i<iend; i++) sum = sum + a[i]; MPI_Reduce(&sum, &tmp, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); sum = tmp; if(myrank == 0) printf("sum = %f \n", sum); MPI_Finalize(); }
MPI 예제 3 Pi 계산 #include <math.h> #define n 100000 main(){ int i,istep,itotal[10],itemp; i i i double r, seed, pi, x, y, angle; pi = 3.1415926; for(i=0;i<10;i++) itotal[i]=0; seed = 0.5; srand(seed); for(i=0; i<n; i++){ x = 0.0; y = 0.0; for(istep=0;istep<10;istep++){ r = (double)rand(); angle = 2.0*pi*r/32768 r/32768.0; x = x + cos(angle); y = y + sin(angle); } itemp = sqrt(x*x + y*y); itotal[itemp]=itotal[itemp]+1; t t ] } } for(i=0; i<10; i++){ printf(" %d :", i); printf("total=%d\n",itotal[i]); } #include <mpi.h> #include <stdio.h> #include <math.h> #define n 100000 void para_range(int, int, int, int, int*, int*); int min(int, int); main (int argc, char *argv[]){ int i, istep, itotal[10], iitotal[10], itemp; int ista, iend, nprocs, myrank; double r, seed, pi, x, y, angle; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, COMM &myrank); para_range(0, n-1, nprocs, myrank, &ista, &iend); pi = 3.1415926; for(i=0; i<10; i++) itotal[i] = 0; seed = 0.5 + myrank; srand(seed); for(i=ista; i<=iend; i++){ x = 0.0; y = 0.0; for(istep=0; istep<10; istep++){ r = (double)rand(); angle = 2.0*pi*r/32768.0; x = x + cos(angle); y = y + sin(angle); } itemp = sqrt(x*x + y*y); itotal[itemp] = itotal[itemp] + 1; } MPI_Reduce(itotal, iitotal, 10, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); for(i=0; i<10; i++){ printf(" %d :", i); printf(" total = %d\n",iitotal[i]); } MPI_Finalize();
MPI 예제 4 ( Full Path 유렵형 Vanilla Call 옵션 ) for ( i=0; i< Nsim ;i++ ) { xt1=s; xt2=s; for(m=0; m<path; m++) { xx1 = rand()/(rand_max+1.0);if(xx1==0.0) xx1=0.0000000000001; xx2 = rand()/(rand_max+1.0);if(xx2==0.0) xx2=0.0000000000001; normal1=sqrt(-2.0*log(xx1 ))*cos(2.0*3.14159265358979323846*xx2 ); xt1= xt1 + r * xt1 * dt + v*xt1*sqrt(dt)*normal1; } oprice=max(xt1-k,0); fsum =fsum+ oprice; } results = (double) 1/Nsim * exp(-1* r * tau)* fsum; stop=clock(); htime = 0.001*difftime(stop,start); //windows printf("\n%19.17f \t %19.17f %7.5f \n",results,results-bsp, htime) ; Loop 분할필요 병렬 RNG 필요 multiseed method MPI reduce 필요
MPI 화코드 #include <mpi.h> #include <stdio.h> #include <math.h> #define n 150000 void para_range(int, int, int, int, int*, int*); int min(int, int); main (int argc, char *argv[]){ int i, m; int ista, iend, nprocs, myrank; double r, seed, pi, x, y, angle; Nsim = n; 금융관련변수설정은모두생략함 (single 코드와동일 ) 병렬화블록 MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); para_range(0, Nsim, nprocs, myrank, &ista, &iend); Multiseed 방법사용 seed = 0.5 + myrank; srand(seed); for ( i=ista; i< iend ;i++ ) { xt1=s; xt2=s; for(m=0; m<path; m++) 상황에따라추가적 RNG 코딩 { xx1 = rand()/(rand_max+1.0);if(xx1==0.0) xx1=0.0000000000001; Single code 그대로사용 xx2 = rand()/(rand_max+1.0);if(xx2==0.0) xx2=0.0000000000001; normal1=sqrt(-2.0*log(xx1 ))*cos(2.0*3.14159265358979323846*xx2 ); xt1= xt1 + r * xt1 * dt + v*xt1*sqrt(dt)*normal1; } oprice=max(xt1-k,0); fsum_local =fsum_local+ oprice; } Core간데이터통신 MPI_Reduce(&fsum_local, &fsum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Finalize(); results = (double) 1/Nsim * exp(-1* r * tau)* fsum; printf(" Option Price : %23.17f " results); }
Tachyon : all node is full
MPI 결과 수퍼컴퓨터센터의 tachyon, nobel 서버모두 busy 연세대학교수학과 PDE 팀 4 node 8 core server kisti Hamel Cluster (5 년이상된 Cluster 로성능이많이떨어지만사용자적음 ) 에서테스트 Linear Scailability : CPU 개수와속도향상의선형성이보장됨 2. Monte Carlo Simulation 은병렬화의 Scailability 가매우좋은편임 N 개의 Core 사용시 1/N 으로계산시간단축효과, 1/sqrt(N) 의정확도향상 현재의정확도를유지하면서속도를향상시킴 : 1/N
수퍼컴퓨터를이용한몬테카를로병렬화 1. 계산전용서버는자체구축필요 (KISTI 의경우 Job schedule : 최소 30 분 ) 2. MPI의경우직접병렬코딩을하고, MPIRUN을실행해줘야하지만잘만들어진 single 코드가존재하는경우몬테카를로병렬화는어렵지않다. 3. 병렬화를통한속도향상 200초걸리는시뮬레이션문제의경우 2노드 32 core system 구축시 6.3초 0.5억원 8노드 128 core system 구축시 1.6초 2억원 + 알파 16노드 256 core system 구축시 0.8초예상됨 4억원 + 알파 32노드 512 core system 구축시 0.4초예상됨 8억원 + 알파 32 노드의사용시 MC 의편리함과 FDM 수준의속도를얻을수있음 Parallle Quasi MC 를사용할경우 8 노드정도로 FDM 수준의속도향상예상됨
추가적인결과 ( 대안적수퍼컴퓨팅 ) 성능과비용을고려한병렬컴퓨팅 ( 그린컴퓨팅 ) 1. 게임용그래픽카드를이용한병렬MC (CUDA) (2008 년 6 월 13 일확률론워크샵발표 ) 2. PS3 클러스터를통한병렬 MC ( IBM Cell BE chip 에대한연구로확대예정 ) 3. 가속보드를통한연산가속 ( Clearspeed CSe620 보드를통한가속 ) 계산요청 PC의 CPU 그래픽카드 계산결과 게임산업의발전괘개발로그래픽칩이 CPU보다성능이높음
CUDAp workstation at Yonsei Math-Finance Lab 계산용서버 : Fedora Linux + Intel Compiler + CUDA Tesla C870 3 개장착 : 이론상 SP 1.5Tflops in SP 보조용워크스테이션 Windows XP + Visual Studio 2003 + CUDA Tesla C870 1개장착 : 이론상 SP 0.5Tflops in SP 비교 Kisti 수퍼컴퓨터 4호기 : 전체 : DP 24Tflops(40억이상 ), 하나의노드 : DP 0.1Tflops(2천만원 )
Algorithm for ELS Parameter input Loop Parallelization XT simulation Boxmuller, RNG(parallel SMT19937) ELS pricing Routine IF flag : prepayment, payoff average results
Pseuco code for CUDA ELS Main(){ cudamalloc((void**) &MTd, 624*sizeof(float)); Dim3 DimGrid(); id() Dim3 DimBlock(); ELS_body <<<DimGrid,DimBlock>>> (parameters); sum( option[k] ) /N; //N is 128 } global ELS_ body(parameters){ Tid= blockidx.x*blockdim.x + threadidx.x; N=blockDim.x *GridDim; Initialize(); ELS_ kernel(parameters); } device ELS_ kernel(parameters){ For( I<0, I< Nsim/N;I++){ For(j<0,j<Totalday){ For(k=0;k<monitor;k++){ Norm1(i)=Box_Muller(); Norm2(I)=Box_Muller(); Xt1(I)= Xt1+MuT*dt + SigmaT*Norm1 Xt2(I)= Xt1+MuT*dt + SigmaT*Norm2 if (Xt1(I) <DB1 and Xt2(I) <=DB2 ) down_flag=1; } if (Xt1(I) >DB1 and Xt2(I) >=DB2 ) up_flag=1; if(j=pre1){} if(j=pre2){} } option = ; sum = sum+option; option = 1/(Nsim/N)*sum; Return option[tid]; } device Initialize(){ Initialize DC of MT19937 } device float Box_Muller(){ U1= MyRand();U2= MyRand(); If( used =1){ return} else {return } } device float MyRand(){ Return MT19937()/4294967296.0; } device float MT19937(){ //Use static variables for each threads Algorithms for MT RNG Return y; }
Single & CUDA ELS 실행속도 삼성증권제 1909 회 ELS MT19937 이용 (10 만회 ) 현재 CPU 병렬화 1 차최적화최적화중 RNG BM 변환 CPU 53 초, GPU1 7.3 초, GPU2 2.1 초 CPU 20 초, GPU1 3.6 초, GPU2 1.1 초 Copy CPU 0 초, GPU1 18.3 초, GPU2 6.2 초 MC CPU 18 초, GPU1 1.2 초, GPU2 0.4 초 총 90 초, 30 초 (2.5 배 ) 9.8 초 (9 배 ) 7.1 초 (12 배 ) C 언어 : AMD 페놈 2.5Ghz, GPU vs. Tesla C870
Futher Research 1. Hybrid bidmethod with ihopenmp, MPI & CUDA MT19937대신 MT607 이용시추가적으로 4배이상추가속도향상가능 ( 작업중 ) pmc with MultiGPU : 16 배성능의 1node 서버비용으로 16 노드정도의 200 배이상속도향상예상 Pseudo RNG 를이용한 MC simulation : 0.45 sec 이내 pricing 가능함을의미 단점 : C870 은현재 single precision 만지원. 차기 H/W 에서 double precision 지원예정 single precision : 소수점 7자리의누적오차 double precisoin : 소수점 15자리누적오차 MC 기법을이용한 1BP 이하수준의에러를확보하기위해서는 5만 ~10만번정도 double Precision RNG 시뮬레이션이적당할것으로생각됨 바로적용불가 ( H/W 출시 : 6 개월이상, S/W 안정화 : 계속연구중 ) S1070 출시 : double precision 사용시, 2 sec 정도로예상됨 ( 테스트예정 )
3. Module 화작업 Parallel Quasi Monte, New RNGs, Fast Convergence 4. Module화작업좀더쉽게병렬화응용 : CUDA 의경우 99% 이상 single ge code 그대로이용할수있도록작업중금융쪽에서 OpenMP 혹은 MPI의직접적용가능성연구 plink, 소켓통신을이용한 excel base 원격 MPIRUN 환경구축 Greek Search, RM 등금융시스템의병렬화연구 5. 다양한병렬화환경연구 CPU, GPU, 가속보드등의병렬처리노하우병렬화자문경험 - LG 전자연구소병렬화작업자문, KISTI 수퍼컴퓨터센터병렬시스템도입자문등병렬화협력 이화여대그래픽연구실, CUDA 포럼등 6. Parallel FDM & FEM 계획중 Parallel LU, parallel CG solver 연구및병렬 PDE solver
The End 경청해주셔서감사합니다. Resume 연세대수학과박사과정유현곤 E-mail : yhgon@yonsei.ac.kr 학력 : 1997.3-2002.2 연세대학교상경계열 97학번 ( 경영학, 경제학이중전공 ) 2006.3-2008.2 연세대학교수학과석사논문 Option Pricing with Stochastic Volatility under Markovian Regimes Switching 2008.3- 현재연세대학교수학과박사과정 1 년차재학중 ( 약 3 년후졸업예정 ^^) 관심분야 : 금융수학, 금융공학 - Structured derivatives(els 등 ), efficient Greek search, Stochastic Volatility, Parallel Computing (MC, FDM, FEM, Meshless) 등