NumPy의 선형대수 모듈에는 LAPACK의 경량화 버전인 lapack_lite가 포함되어 있으며, 그 중 f2c_blas.c 파일은 복소수 연산을 위한 BLAS(Basic Linear Algebra Subprograms) 루틴들을 포함하고 있다. 이 파일은 f2c 도구를 통해 Fortran 소스에서 자동 생성된 코드로, 단정밀도 복소수(single-precision complex) 연산을 수행하는 핵심 수학 루틴들을 제공한다.
CAXPY: 복소수 스케일 벡터 누적
CAXPY는 y = alpha * x + y 연산을 수행하는 루틴으로, 복소수 상수에 벡터를 곱한 뒤 다른 벡터에 더하는 기능을 제공한다. 벡터의 증분값(increment)이 1인 경우와 그렇지 않은 경우를 분기하여 처리하며, 증분이 음수인 경우 시작 인덱스를 역방향으로 조정한다.
int complex_axpy(int vec_len, complex *scale, complex *vec_x,
int step_x, complex *vec_y, int step_y)
{
if (vec_len <= 0) return 0;
if (c_abs1(scale) == 0.0f) return 0;
int pos_x = 1, pos_y = 1;
if (step_x < 0) pos_x = (1 - vec_len) * step_x + 1;
if (step_y < 0) pos_y = (1 - vec_len) * step_y + 1;
for (int idx = 0; idx < vec_len; idx++) {
complex product;
product.re = scale->re * vec_x[pos_x].re - scale->im * vec_x[pos_x].im;
product.im = scale->re * vec_x[pos_x].im + scale->im * vec_x[pos_x].re;
vec_y[pos_y].re += product.re;
vec_y[pos_y].im += product.im;
pos_x += step_x;
pos_y += step_y;
}
return 0;
}
CCOPY: 벡터 복사
CCOPY는 복소수 벡터를 다른 벡터로 복사하는 단순 루틴이다. 증분이 1로 동일한 경우 최적화된 연속 접근 루프를 사용하고, 그렇지 않으면 증분에 따라 인덱스를 계산하며 복사를 수행한다.
int complex_copy(int count, complex *src, int src_step,
complex *dst, int dst_step)
{
if (count <= 0) return 0;
int src_pos = 1, dst_pos = 1;
if (src_step < 0) src_pos = (1 - count) * src_step + 1;
if (dst_step < 0) dst_pos = (1 - count) * dst_step + 1;
for (int idx = 0; idx < count; idx++) {
dst[dst_pos].re = src[src_pos].re;
dst[dst_pos].im = src[src_pos].im;
src_pos += src_step;
dst_pos += dst_step;
}
return 0;
}
CDOTC와 CDOTU: 복소수 내적
CDOTC는 첫 번째 벡터를 켤레 복소수(conjugate)로 변환한 후 내적을 계산하며, CDOTU는 켤레 변환 없이 두 벡터의 내적을 직접 계산한다. 두 루틴 모두 증분이 1인 경우와 그렇지 않은 경우를 분기 처리한다.
void complex_dotc(complex *result, int n, complex *xv, int dx,
complex *yv, int dy)
{
result->re = 0.0f;
result->im = 0.0f;
if (n <= 0) return;
complex accum = {0.0f, 0.0f};
int px = 1, py = 1;
if (dx < 0) px = (1 - n) * dx + 1;
if (dy < 0) py = (1 - n) * dy + 1;
for (int idx = 0; idx < n; idx++) {
complex conj_x = {xv[px].re, -xv[px].im};
complex prod;
prod.re = conj_x.re * yv[py].re - conj_x.im * yv[py].im;
prod.im = conj_x.re * yv[py].im + conj_x.im * yv[py].re;
accum.re += prod.re;
accum.im += prod.im;
px += dx;
py += dy;
}
result->re = accum.re;
result->im = accum.im;
}
CGEMM: 복소수 행렬-행렬 곱셈
CGEMM은 C = alpha * op(A) * op(B) + beta * C 연산을 수행하며, 전치(transpose)와 켤레 전치(conjugate transpose) 옵션을 지원한다. 입력 파라미터 검증 후, alpha가 0인 특수 케이스를 별도로 처리하고, 그렇지 않은 경우 TRANSA와 TRANSB의 조합에 따라 여러 분기로 나뉘어 계산된다.
int complex_gemm(char op_a, char op_b, int rows_m, int cols_n, int inner_k,
complex *alpha, complex *mat_a, int ld_a,
complex *mat_b, int ld_b,
complex *beta, complex *mat_c, int ld_c)
{
int is_normal_a = (op_a == 'N');
int is_normal_b = (op_b == 'N');
int is_conj_a = (op_a == 'C');
int is_conj_b = (op_b == 'C');
int rows_a = is_normal_a ? rows_m : inner_k;
int rows_b = is_normal_b ? inner_k : cols_n;
// 입력 검증
int err = validate_params(op_a, op_b, rows_m, cols_n, inner_k,
ld_a, ld_b, ld_c, rows_a, rows_b);
if (err != 0) { report_error("CGEMM", err); return 0; }
// 빠른 반환 조건 확인
if (rows_m == 0 || cols_n == 0) return 0;
if (is_zero(alpha) && is_one(beta)) return 0;
// alpha가 0인 경우: C = beta * C 또는 C = 0
if (is_zero(alpha)) {
scale_matrix(mat_c, beta, rows_m, cols_n, ld_c);
return 0;
}
// op(A)와 op(B)의 조합에 따른 분기 처리
if (is_normal_b) {
if (is_normal_a) {
multiply_nn(alpha, mat_a, mat_b, beta, mat_c,
rows_m, cols_n, inner_k, ld_a, ld_b, ld_c);
} else if (is_conj_a) {
multiply_cn(alpha, mat_a, mat_b, beta, mat_c,
rows_m, cols_n, inner_k, ld_a, ld_b, ld_c);
} else {
multiply_tn(alpha, mat_a, mat_b, beta, mat_c,
rows_m, cols_n, inner_k, ld_a, ld_b, ld_c);
}
} else {
// TRANSB가 T 또는 C인 경우의 분기들...
}
return 0;
}
CGEMV: 복소수 행렬-벡터 곱셈
CGEMV는 y = alpha * op(A) * x + beta * y 연산을 수행한다. TRANS 파라미터가 'N'인 경우 직접 곱셈을, 'T'인 경우 전치 행렬 곱셈을, 'C'인 경우 켤레 전치 행렬 곱셈을 수행한다. 베타로 y를 먼저 스케일링한 후 알파 연산을 적용하는 구조를 가진다.
CGERC와 CGERU: 랭크-1 갱신
CGERC는 A = alpha * x * conj(y') + A 연산을, CGERU는 A = alpha * x * y' + A 연산을 수행한다. 두 루틴의 유일한 차이는 y 벡터에 켤레 변환을 적용하는지 여부이며, 나머지 구조는 동일하다. 증분이 1인 경우 단순 인덱싱으로, 그렇지 않은 경우 시작 위치를 계산하여 처리한다.
int complex_gerc(int m_dim, int n_dim, complex *alpha,
complex *vec_x, int inc_x, complex *vec_y, int inc_y,
complex *mat_a, int ld_a)
{
if (m_dim == 0 || n_dim == 0 || is_zero(alpha)) return 0;
int y_start = (inc_y > 0) ? 1 : (1 - (n_dim - 1) * inc_y);
if (inc_x == 1) {
for (int col = 1; col <= n_dim; col++) {
if (is_nonzero(&vec_y[y_start])) {
complex conj_y = {vec_y[y_start].re, -vec_y[y_start].im};
complex scaled_alpha;
scaled_alpha.re = alpha->re * conj_y.re - alpha->im * conj_y.im;
scaled_alpha.im = alpha->re * conj_y.im + alpha->im * conj_y.re;
for (int row = 1; row <= m_dim; row++) {
complex outer;
outer.re = vec_x[row].re * scaled_alpha.re
- vec_x[row].im * scaled_alpha.im;
outer.im = vec_x[row].re * scaled_alpha.im
+ vec_x[row].im * scaled_alpha.re;
mat_a[row + col * ld_a].re += outer.re;
mat_a[row + col * ld_a].im += outer.im;
}
}
y_start += inc_y;
}
}
return 0;
}
CHEMV: 에르미트 행렬-벡터 곱셈
CHEMV는 y = alpha * A * x + beta * y 연산을 수행하며, 여기서 A는 n×n 에르미트 행렬이다. UPLO 파라미터에 따라 상삼각 또는 하삼각 부분만 참조하며, 대각선 원소는 실수로 간주한다. 상삼각 저장 시 대각선 위쪽 원소와 그 켤레를 이용하여 하삼각 부분을 암시적으로 처리한다.
CHER2와 CHERK: 에르미트 랭크-k 갱신
CHER2는 랭크-2 갱신 A = alpha * x * conj(y') + conj(alpha) * y * conj(x') + A를 수행하며, CHERK는 랭크-k 갱신 C = alpha * A * conj(A') + beta * C 또는 C = alpha * conj(A') * A + beta * C를 수행한다. CHERK에서 beta는 실수 스칼라이며, 대각선 원소의 허수부는 항상 0으로 설정된다.
int complex_herk(char uplo, char trans, int n, int k,
float alpha, complex *a, int lda,
float beta, complex *c, int ldc)
{
int is_upper = (uplo == 'U');
int is_normal = (trans == 'N');
int rows_a = is_normal ? n : k;
// 검증 생략...
if (n == 0) return 0;
if (alpha == 0.0f && beta == 1.0f) return 0;
if (alpha == 0.0f) {
// C = beta * C (대각선은 실수로 설정)
scale_hermitian(c, beta, n, ldc, is_upper);
return 0;
}
if (is_normal) {
// C = alpha * A * conj(A') + beta * C
for (int j = 1; j <= n; j++) {
// 대각선 및 비대각선 원소 갱신
if (is_upper) {
for (int l = 1; l <= k; l++) {
complex conj_al = conj(a[j + l * lda]);
complex scaled = {alpha * conj_al.re, alpha * conj_al.im};
for (int i = 1; i < j; i++) {
complex term;
term.re = scaled.re * a[i + l * lda].re
- scaled.im * a[i + l * lda].im;
term.im = scaled.re * a[i + l * lda].im
+ scaled.im * a[i + l * lda].re;
c[i + j * ldc].re += term.re;
c[i + j * ldc].im += term.im;
}
// 대각선: 실수부만 누적, 허수부는 0
}
}
}
}
return 0;
}
CSCAL과 CSSCAL: 벡터 스케일링
CSCAL은 복소수 스칼라로 복소수 벡터를 스케일링하고, CSSCAL은 실수 스칼라로 복소수 벡터를 스케일링한다. 두 루틴 모두 증분이 1인 경우와 그렇지 않은 경우를 분기하며, 증분이 음수이거나 0 이하인 경우를 적절히 처리한다.
int complex_scal(int n, complex *factor, complex *vec, int inc)
{
if (n <= 0 || inc <= 0) return 0;
if (inc == 1) {
for (int i = 1; i <= n; i++) {
complex result;
result.re = factor->re * vec[i].re - factor->im * vec[i].im;
result.im = factor->re * vec[i].im + factor->im * vec[i].re;
vec[i] = result;
}
} else {
int total = n * inc;
for (int i = 1; i <= total; i += inc) {
complex result;
result.re = factor->re * vec[i].re - factor->im * vec[i].im;
result.im = factor->re * vec[i].im + factor->im * vec[i].re;
vec[i] = result;
}
}
return 0;
}
CSROT: 복소수 평면 회전
CSROT는 실수 코사인(c)과 사인(s) 값을 사용하여 두 복소수 벡터에 평면 회전을 적용한다. 변환 공식은 x = c*x + s*y, y = c*y - s*x이며, 복소수 벡터에 실수 스칼라를 곱하는 형태로 구현된다.
CSWAP: 벡터 교환
CSWAP은 두 복소수 벡터의 내용을 교환하는 루틴이다. 임시 변수를 사용하여 원소 단위로 스왑을 수행하며, 증분 처리 방식은 CCOPY와 동일한 패턴을 따른다.
CTRMM: 삼각 행렬-행렬 곱셈
CTRMM은 B = alpha * op(A) * B 또는 B = alpha * B * op(A) 연산을 수행한다. SIDE, UPLO, TRANSA, DIAG 파라미터의 조합에 따라 다양한 분기가 발생하며, 단위 삼각(unit triangular) 여부에 따라 대각선 원소를 1로 가정하거나 실제 값을 사용한다. 상삼각 행렬의 경우 전방 순회로, 하삼각 행렬의 경우 역방향 순회로 처리한다.
int complex_trmm(char side, char uplo, char transa, char diag,
int m, int n, complex *alpha,
complex *a, int lda, complex *b, int ldb)
{
int left_side = (side == 'L');
int is_upper = (uplo == 'U');
int no_conj = (transa == 'T');
int non_unit = (diag == 'N');
int a_rows = left_side ? m : n;
// 검증 생략...
if (m == 0 || n == 0) return 0;
if (is_zero(alpha)) { zero_matrix(b, m, n, ldb); return 0; }
if (left_side) {
if (transa == 'N') {
// B := alpha * A * B
if (is_upper) {
for (int j = 1; j <= n; j++) {
for (int k = 1; k <= m; k++) {
if (is_nonzero(&b[k + j * ldb])) {
complex temp;
temp.re = alpha->re * b[k + j * ldb].re
- alpha->im * b[k + j * ldb].im;
temp.im = alpha->re * b[k + j * ldb].im
+ alpha->im * b[k + j * ldb].re;
for (int i = 1; i < k; i++) {
complex prod;
prod.re = temp.re * a[i + k * lda].re
- temp.im * a[i + k * lda].im;
prod.im = temp.re * a[i + k * lda].im
+ temp.im * a[i + k * lda].re;
b[i + j * ldb].re += prod.re;
b[i + j * ldb].im += prod.im;
}
if (non_unit) {
// temp /= A[k][k]
}
b[k + j * ldb] = temp;
}
}
}
}
}
}
return 0;
}
CTRMV와 CTRSV: 삼각 행렬 연산
CTRMV는 x = op(A) * x 행렬-벡터 곱셈을, CTRSV는 op(A) * x = b 삼각 방정식 시스템을 푼다. 두 루틴은 동일한 파라미터 구조를 가지며, UPLO에 따라 상삼각/하삼각, TRANS에 따라 일반/전치/켤레전치, DIAG에 따라 단위/비단위 삼각 행렬을 지정한다.
int complex_trsv(char uplo, char trans, char diag, int n,
complex *a, int lda, complex *x, int incx)
{
int no_conj = (trans == 'T');
int non_unit = (diag == 'N');
int is_upper = (uplo == 'U');
if (n == 0) return 0;
int start_x = (incx <= 0) ? 1 - (n - 1) * incx : 1;
if (trans == 'N') {
// x := inv(A) * x
if (is_upper) {
int jx = start_x + (n - 1) * incx;
for (int j = n; j >= 1; j--) {
if (is_nonzero(&x[jx])) {
if (non_unit) {
complex_div(&x[jx], &x[jx], &a[j + j * lda]);
}
complex temp = x[jx];
int ix = jx;
for (int i = j - 1; i >= 1; i--) {
ix -= incx;
complex prod;
prod.re = temp.re * a[i + j * lda].re
- temp.im * a[i + j * lda].im;
prod.im = temp.re * a[i + j * lda].im
+ temp.im * a[i + j * lda].re;
x[ix].re -= prod.re;
x[ix].im -= prod.im;
}
}
jx -= incx;
}
}
}
return 0;
}
CTRSM: 삼각 행렬 해법
CTRSM은 op(A) * X = alpha * B 또는 X * op(A) = alpha * B를 풀어 해 행렬 X를 구하며, 결과는 B를 덮어쓴다. SIDE가 'L'인 경우 좌측 해법을, 'R'인 경우 우측 해법을 수행하며, TRANSA 옵션에 따라 전치/켤레전치를 지원한다. alpha가 0인 경우 B를 0으로 채우고 즉시 반환한다.
int complex_trsm(char side, char uplo, char transa, char diag,
int m, int n, complex *alpha,
complex *a, int lda, complex *b, int ldb)
{
int left_side = (side == 'L');
int is_upper = (uplo == 'U');
int no_conj = (transa == 'T');
int non_unit = (diag == 'N');
int a_rows = left_side ? m : n;
// 검증 생략...
if (m == 0 || n == 0) return 0;
if (is_zero(alpha)) {
zero_matrix(b, m, n, ldb);
return 0;
}
if (left_side) {
if (transa == 'N') {
// B := alpha * inv(A) * B
if (is_upper) {
for (int j = 1; j <= n; j++) {
// B 열을 alpha로 스케일
scale_vector(&b[1 + j * ldb], alpha, m, ldb);
for (int k = m; k >= 1; k--) {
if (is_nonzero(&b[k + j * ldb])) {
if (non_unit) {
complex_div(&b[k + j * ldb],
&b[k + j * ldb],
&a[k + k * lda]);
}
complex bkj = b[k + j * ldb];
for (int i = 1; i < k; i++) {
complex prod;
prod.re = bkj.re * a[i + k * lda].re
- bkj.im * a[i + k * lda].im;
prod.im = bkj.re * a[i + k * lda].im
+ bkj.im * a[i + k * lda].re;
b[i + j * ldb].re -= prod.re;
b[i + j * ldb].im -= prod.im;
}
}
}
}
}
}
}
return 0;
}
구현 패턴 요약
이 파일의 모든 루틴은 공통적인 구조적 패턴을 따른다:
- 파라미터 검증: 호출 전에 모든 입력 파라미터의 유효성을 검사하고, 오류 시
xerbla_를 호출한다. - 빠른 반환: 차원이 0이거나 스칼라가 0/1인 특수 케이스를 조기에 처리한다.
- 증분 처리: 증분이 1인 경우 최적화된 루프를, 음수/양수 증분인 경우 시작 위치를 계산하여 처리한다.
- 복소수 연산: f2c가 생성한
q__1,q__2등의 임시 변수를 통해 복소수 곱셈과 덧셈을 단계별로 수행한다. - 켤레 처리:
r_cnjg매크로를 사용하여 복소수의 켤레를 구하고, 에르미트 연산에 활용한다. - 삼각 행렬: 상삼각은 전방 순회, 하삼각은 역방향 순회로 처리하여 올바른 연산 순서를 보장한다.
이러한 루틴들은 NumPy의 numpy.linalg 모듈이 LAPACK 기능을 독립적으로 제공할 수 있도록 하는 기반 계층으로 작용하며, f2c 변환을 통해 Fortran BLAS의 정확한 수학적 동작을 C 코드로 보존하고 있다.