Lapackの複素数連立方程式ルーチンについて、うまく計算できない。

/*計算結果*/
-449.108242 83.321424 -2473.112709 298.807186 
815.629797 3.301997 4438.644075 288.842703 
122.389464 74.621908 640.321241 388.495341 
-30234.634654 21146.882581 -171817.275446 104708.033646 

/*こちらが正しい答え*/
1.000000 1.000000 -1.000000 -2.000000 
2.000000 -3.000000 5.000000 1.000000 
-4.000000 -5.000000 -3.000000 4.000000 
0.000000 6.000000 2.000000 -3.000000 

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <complex.h>
#include <stdlib.h>

int N, M;

// LU分解
void zgetrf_(int* m, int* n, double complex* A, int* LDA, int* IPIV, int* INFO);
void LU_bunkai(int matrix_size, double complex matrix[matrix_size][matrix_size], int IPIV[matrix_size]){ //LU分解
    int m = matrix_size;
    int n = matrix_size;
    double complex A[m * n];
    for(int i = 0; i < m; i++){
        for(int j = 0; j < n; j++){
            A[i * n + j] = matrix[j][i]; // Fixed the assignment
        }
    }
    int LDA = m;
    int INFO;
    zgetrf_(&m, &n, A, &LDA, IPIV, &INFO);
    for(int i = 0; i < m; i++){
        for(int j = 0; j < n; j++){
            matrix[j][i] = A[i * n + j]; // Fixed the assignment
        }
    }
}

//連立方程式
void zgetrs_(char* TRANS, int* n, int* NRHS, double complex* a, int* LDA, int* IPIV, double complex* b, int* LDB, int* INFO);
void simultaneous_eq(int matrix_size_A, int length_B, double complex matrix_A[matrix_size_A][matrix_size_A], double complex matrix_B[matrix_size_A][length_B],
                     double complex matrix_C[matrix_size_A][length_B]){
    int IPIV[matrix_size_A];
    LU_bunkai(matrix_size_A, matrix_A, IPIV);
    // for(int i = 0; i < matrix_size_A; i++){ printf("%d\n", IPIV[i]); }
    char TRANS = 'N';
    int n = matrix_size_A;
    int NRHS = length_B;
    double complex a[matrix_size_A * matrix_size_A];
    for(int i = 0; i < matrix_size_A; i++){
        for(int j = 0; j < matrix_size_A; j++){
            a[i * matrix_size_A + j] = matrix_A[i][j];
            // printf("%1.1lf+%1.1lf ", creal(matrix_A[i][j]), cimag(matrix_A[i][j]));
        }
        // printf("\n");
    }
    double complex b[matrix_size_A * length_B];
    for(int i = 0; i < length_B; i++){
        for(int j = 0; j < matrix_size_A; j++){
            b[i * matrix_size_A + j] = matrix_B[i][j];
        }
    }
    int LDA = matrix_size_A;
    int LDB = matrix_size_A;
    int INFO;
    zgetrs_ (&TRANS, &n, &NRHS, a, &LDA, IPIV, b, &LDB, &INFO);
    for(int i = 0; i < length_B; i++){
        for(int j = 0; j < matrix_size_A; j++){
            matrix_C[j][i] =  b[i * matrix_size_A + j];
            // printf("%lf+%lfI ", creal(matrix_C[i][j]), cimag(matrix_C[i][j]));
        }
        // printf("\n");
    }
}

void make_matrix(double complex A[4][4], double complex B[4][2]){
    A[0][0] = -1.34 + 2.55 * I;   A[0][1] = 0.28 + 3.17 * I;     A[0][2] = -6.39 + -2.20 * I; A[0][3] =  0.72 + -0.92 * I; 
    A[1][0] = -1.70 + -14.10 * I; A[1][1] =  33.10 +  -1.50 * I; A[1][2] = -1.50 + 13.40 * I; A[1][3] = 12.90 + 13.80 * I; 
    A[2][0] = -3.29 +  -2.39 * I; A[2][1] =  -1.91 +   4.42 * I; A[2][2] = -0.14 + -1.35 * I; A[2][3] =  1.72 +  1.35 * I; 
    A[3][0] =  2.41 +   0.39 * I; A[3][1] =  -0.56 +   1.47 * I; A[3][2] = -0.83 + -0.69 * I; A[3][3] = -1.96 +  0.67 * I;  

    B[0][0] = 26.26 +  51.78 * I; B[0][1] = 31.32 +  -6.70 * I;
    B[1][0] = 64.30 + -86.80 * I; B[1][1] = 158.60 + -14.20 * I;
    B[2][0] = -5.75 +  25.31 * I; B[2][1] = -2.15 +  30.19 * I;
    B[3][0] =  1.16 +   2.57 * I; B[3][1] = -2.56 +   7.55 * I; 
}                  

int main(){
    N = 4, M = 2; double complex matrix_A[N][N], matrix_B[N][M], matrix_C[N][M];
    make_matrix(matrix_A, matrix_B);
    simultaneous_eq(N, M, matrix_A, matrix_B, matrix_C);
    for(int i = 0; i < N; i++){
        for(int j = 0; j < M; j++){
            printf("%lf %lf ", creal(matrix_C[i][j]), cimag(matrix_C[i][j]));
        }
        printf("\n");
    }

    return 0;
}