Luke's Archive

$L_i$, $D_i$, $U_i$를 $n \times n$ matrix라 하자. ($1\leq i \leq N$)

매트랩에서 다음과 같은 $Nn \times Nn$ block tridiagonal matrix $A$를 만들어보자.

$$A = \begin{bmatrix} D_1 & U_1 & & \\ L_1 & D_2 & U_2 & & \\ & \ddots & \ddots & \ddots & \\ & & L_{N-2} & D_{N-1} & U_{N-1} \\ & & & L_{N-1} & D_N \end{bmatrix}$$

1. Block Diagonal Matrix 만들기

우선 block diagonal matrix를 만드는 방법부터 생각해본다.

$B_i$를 $n \times n$ matrix라 하자.

매트랩으로 다음과 같은 $Nn \times Nn$ block diagonal matrix $M$를 만든다고 생각해보자.

$$M=\begin{bmatrix} B_{ 1 } &  &  &  &  \\  & B_{ 2 } &  &  &  \\  &  & \ddots  &  &  \\  &  &  & B_{ N-1 } &  \\  &  &  &  & B_{ N } \end{bmatrix}$$

1-1. $B = B_i = B_j \; \forall i,j$일 경우

kron(kronecker tensor product)이라는 함수를 이용하면 쉽게 할 수 있다.

kron(A,B)에서 A는 B를 하나의 element로 봤을 때의 matrix를 의미하고, B가 block이 된다.

M = kron(speye(N),B);

1-2. 1-1이 아닐 경우

blkdiag(block diagonal)라는 함수를 이용한다.

우선 $B_i$가 Bi로 저장되어 있다고 가정하자.

(cell로 저장하든 일일히 다른 변수명들로 저장하든 이건 알아서 하자)

blkdiag(B1,B2,B3,...)으로 입력하면 되는데 만약 $N = 100$이라 하면 B1부터 B100까지 일일히 쓰는 건 꽤나 효율적이지 못하다.

그래서 for문을 이용해 M = blkdiag(B1,B2,B3,...,BN)이라는 문자형 data를 만든 뒤 eval이란 함수로 실행시키는 방식으로 한다.

mat_input = [];
for i = 1:N
    mat_input = [mat_input sprintf('B%d,',i)];
end
mat_input = ['M=blkdiag(' ,mat_input(1:end-1), ');'];
eval(mat_input);

또한, 1-1의 경우더라도 repmat을 이용해서 이 방식으로 만들면 kron보다 더 빠르다고 한다.

2. Block Tridiagonal Matrix 만들기

우선 1에서의 방법을 통해 $D_i$로 만든 $Nn \times Nn$ block diagonal matrix $A_D$를 만든다.

만약 $L = L_i = L_j \; \forall i,j$일 경우 1-1의 방법을 응용해 $Nn \times Nn$ block lower diagonal matrix $A_L$을 만든다.

같은 논리로 $U = U_i = U_j \; \forall i,j$일 경우에도 $Nn \times Nn$ block upper diagonal matrix $A_U$를 만들 수 있다.

AL = kron(diag(ones(N-1,1),-1), L);
AU = kron(diag(ones(N-1,1),1), U);

아닐 경우 1-2의 방법을 통해 $(N-1)n \times (N-1)n$ block diagonal matrix $A_L$를 만든 뒤 0을 추가해 $Nn \times Nn$ matrix로 확장한다.

역시 같은 논리로 $A_U$를 만들 수 있다.

AL = [zeros(N*n,n); AL zeros(n,(N-1)*n)];
AU = [AU zeros(n,(N-1)*n); zeros(N*n,n)];

그 다음 $A_L$, $A_D$, $A_U$를 다 더하면 block tridiagonal matrix $A$가 완성된다.

A = AL + AD + AU;

Reference)

http://jasonjuang.blogspot.kr/2014/07/how-to-create-block-diagonal-matrix.html

http://stackoverflow.com/questions/16370371/matlab-create-a-block-diagonal-matrix-with-same-repeating-block

C언어 함수를 MEX라는 걸 이용해 변환시켜 매트랩에서 쓸 수 있다.

이때 변환과정에서 C언어 컴파일러가 쓰이는데, 공식적으로 OS X에서 MEX를 사용할 때 Xcode에 있는 clang만 지원한다.

OS 별 MEX 지원 컴파일러는 여기에서 확인할 수 있다.

mex -setup을 입력하면 현재 설정되어 있는 언어와 컴파일러를 확인할 수 있다.

Xcode에 있는 clang은 알다시피 OpenMP를 지원하지 않는다(위의 링크에서도 확인할 수 있다).

OpenMP를 OS X에서 사용하려면 clang-omp를 쓰거나 Xcode용이 아닌 별도의 gcc를 설치해야 하지만(설치 과정은 여기에서 확인),

공식적으로 매트랩에선 지원하지 않는 컴파일러이기 때문에 추가적인 작업이 필요하다.

이번 포스트에서는 그 과정을 설명하고자 한다.

아래의 기준으로 설명한다.

- gcc 6.2.1

- OS X 10.11 El Capitan

- MATLAB R2016b

1. gcc_maci64.xml 만들기

MEX C 컴파일러를 지정해주는 건 매트랩 내에 있는 clang_maci64.xml이라고 하는 파일에 의해서다.

우선 해당 파일부터 열기 위해 파인더에서 다음의 경로로 간다.

mexopts 폴더 안에 3가지가 있는데(C, C++, fortran), 현재 우리는 C언어만 다룬다는 전제이므로 clang_maci64.xml을 복붙하고 이름을 gcc_maci64.xml로 수정한다.

그 후 gcc_maci64.xml을 자신이 쓰고 있는 텍스트 에디터를 통해 연다.

연 뒤 앞쪽에 보면

이라고 되어 있는 부분이 있는데 이는 clang에 관한 내용이므로 gcc에 맞게 알아서 수정한다.

난 아래와 같이 수정했다.

CC의 루트가 clang으로 아래와 같이 지정되어 있는데,

이를 gcc 경로로 수정한다. Homebrew로 설치했다면 아래와 같다.

마지막으로, OpenMP를 사용하기 위해 CFLAGS, LDFLAGS 뒤에 -fopenmp를 추가해준다.

2. 매트랩에서 컴파일러 변경하기

매트랩으로 돌아가서 mex -setup을 다시 입력해본다.

입력하면 아래와 같이 나온다.

이전과 다르게 C 컴파일러에서 선택권이 생겼다.

GNU GCC를 클릭한다.

이제 mex를 사용할 때 C 컴파일러는 gcc로 지정이 되었다.

테스트를 해보기 위해 mex test.c를 입력한다.

test.c 코드는 여기에 있는 코드를 사용했고, 내용은 다음과 같다.

#include "mex.h"
#include <stdio.h>
#include "omp.h"

void mexFunction(int nlhs,mxArray* plhs[], int nrhs, const mxArray* prhs[])
{
    int i;
    double sum; 
    /* Check for proper number of arguments. */
    if (nrhs!=0)
    {
        mexErrMsgTxt("No input");
    }
    else if (nlhs!=0)
    {
        mexErrMsgTxt("No output");
    }
    
    omp_set_num_threads(omp_get_max_threads());
    mexPrintf("Max number of threads %d\n",omp_get_max_threads());

    #pragma omp parallel for
    for (i = 0; i < omp_get_max_threads(); i++)
    {
        mexPrintf("Hello from thread %d, nthreads %d\n", omp_get_thread_num(), omp_get_num_threads());
    }
    return; 
}

MEX가 성공적으로 완료되었다고는 하지만 위에 두 개의 에러가 나타난다.

두 개의 에러는 본질적으로는 같은 에러이고, 메세지는 다음과 같은 방식으로 뜬다.

저 에러를 해결할 수 있는 방법을 아직 찾지 못했지만, 당장에 문제가 있지 않아 이대로 두고 사용하고 있다.

나중에 저 에러를 해결하게 되면 포스트를 수정하겠다.

변환된 파일을 실행해보기 위해 test라고 입력해보자.

이 코드는 input, output이 없는 거라 test라고만 입력해도 된다.

MEX로 변환된 파일이 문제없이 잘 작동함을 알 수 있다.

Reference)

https://kr.mathworks.com/support/compilers.html

https://kr.mathworks.com/matlabcentral/newsreader/view_thread/334250

http://stackoverflow.com/questions/37362414/openmp-with-mex-in-matlab-on-mac

특정 파트에서 글씨 크기를 키우거나, 작게 쓰거나 할 때 사용할 수 명령어들을 모아두었다.

10pt, 11pt, 12pt는 문서 전체에 설정된 기본 폰트 크기를 의미한다.

Reference)

http://tex.stackexchange.com/questions/24599/what-point-pt-font-size-are-large-etc