Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 10 additions & 10 deletions pygrt/C_extension/include/grt/common/RT_matrix.h
Original file line number Diff line number Diff line change
Expand Up @@ -16,18 +16,18 @@

/** 应该把 PSV 和 SH 的小矩阵统一放在一个结构体中管理 */
typedef struct {
MYCOMPLEX RD[2][2]; ///< P-SV 下传反射系数矩阵
MYCOMPLEX RU[2][2]; ///< P-SV 上传反射系数矩阵
MYCOMPLEX TD[2][2]; ///< P-SV 下传透射系数矩阵
MYCOMPLEX TU[2][2]; ///< P-SV 上传透射系数矩阵
MYCOMPLEX RDL; ///< SH 下传反射系数
MYCOMPLEX RUL; ///< SH 上传反射系数
MYCOMPLEX TDL; ///< SH 下传透射系数
MYCOMPLEX TUL; ///< SH 上传透射系数
cplx_t RD[2][2]; ///< P-SV 下传反射系数矩阵
cplx_t RU[2][2]; ///< P-SV 上传反射系数矩阵
cplx_t TD[2][2]; ///< P-SV 下传透射系数矩阵
cplx_t TU[2][2]; ///< P-SV 上传透射系数矩阵
cplx_t RDL; ///< SH 下传反射系数
cplx_t RUL; ///< SH 上传反射系数
cplx_t TDL; ///< SH 下传透射系数
cplx_t TUL; ///< SH 上传透射系数

/* 一些辅助变量 */
MYCOMPLEX invT[2][2];
MYCOMPLEX invTL;
cplx_t invT[2][2];
cplx_t invTL;

/* 状态变量 */
MYINT stats; ///< 是否有除零错误,非0为异常值
Expand Down
4 changes: 2 additions & 2 deletions pygrt/C_extension/include/grt/common/attenuation.h
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,9 @@
*
* @return atncoef 系数因子,作用在 \f$ k=\omega / c(\omega)\f$的计算
*/
MYCOMPLEX grt_attenuation_law(MYREAL Qinv, MYCOMPLEX omega);
cplx_t grt_attenuation_law(real_t Qinv, cplx_t omega);

/**
* attenuation_law函数在python中被调用的版本,长度2的数组分别表示复数的实部和虚部
*/
void grt_py_attenuation_law(MYREAL Qinv, MYREAL omg[2], MYREAL atte[2]);
void grt_py_attenuation_law(real_t Qinv, real_t omg[2], real_t atte[2]);
4 changes: 2 additions & 2 deletions pygrt/C_extension/include/grt/common/bessel.h
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
* @param[out] bj2 \f$ J_2(x) \f$
*
*/
void grt_bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2);
void grt_bessel012(real_t x, real_t *bj0, real_t *bj1, real_t *bj2);


/**
Expand All @@ -31,4 +31,4 @@ void grt_bessel012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2);
* @param[in,out] bj2 传入 \f$ J_2(x) \f$, 返回\f$ J_2^{'}(x) \f$
*
*/
void grt_besselp012(MYREAL x, MYREAL *bj0, MYREAL *bj1, MYREAL *bj2);
void grt_besselp012(real_t x, real_t *bj0, real_t *bj1, real_t *bj2);
16 changes: 7 additions & 9 deletions pygrt/C_extension/include/grt/common/const.h
Original file line number Diff line number Diff line change
Expand Up @@ -32,19 +32,17 @@
#endif


typedef int MYINT; ///< 整数

// #define GRT_USE_FLOAT ///< 是否使用单精度浮点数
typedef double real_t;
typedef double complex cplx_t;

typedef int MYINT; ///< 整数
#define GRT_NC_MYINT NC_INT
#define GRT_NC_FUNC_MYINT(func) func##_int

#define NC_REAL NC_DOUBLE
#define NC_FUNC_REAL(func) func##_double

#ifdef GRT_USE_FLOAT
typedef float _Complex MYCOMPLEX; ///< 复数
typedef float MYREAL; ///< 浮点数
#else
typedef double _Complex MYCOMPLEX;
typedef double MYREAL;
#endif

// 常数
#define RTWOTHIRD 0.6666666666666667 ///< 2/3
Expand Down
12 changes: 6 additions & 6 deletions pygrt/C_extension/include/grt/common/dwm.h
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,11 @@
*
* @return k 积分截至时的波数
*/
MYREAL grt_discrete_integ(
GRT_MODEL1D *mod1d, MYREAL dk, MYREAL kmax, MYREAL keps,
MYINT nr, MYREAL *rs,
MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
real_t grt_discrete_integ(
GRT_MODEL1D *mod1d, real_t dk, real_t kmax, real_t keps,
MYINT nr, real_t *rs,
cplx_t sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
bool calc_upar,
MYCOMPLEX sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
MYCOMPLEX sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
cplx_t sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
cplx_t sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
FILE *fstats, GRT_KernelFunc kerfunc);
12 changes: 6 additions & 6 deletions pygrt/C_extension/include/grt/common/fim.h
Original file line number Diff line number Diff line change
Expand Up @@ -48,13 +48,13 @@
*
* @return k 积分截至时的波数
*/
MYREAL grt_linear_filon_integ(
GRT_MODEL1D *mod1d, MYREAL k0, MYREAL dk0, MYREAL filondk, MYREAL kmax, MYREAL keps,
MYINT nr, MYREAL *rs,
MYCOMPLEX sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
real_t grt_linear_filon_integ(
GRT_MODEL1D *mod1d, real_t k0, real_t dk0, real_t filondk, real_t kmax, real_t keps,
MYINT nr, real_t *rs,
cplx_t sum_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
bool calc_upar,
MYCOMPLEX sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
MYCOMPLEX sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
cplx_t sum_uiz_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
cplx_t sum_uir_J[nr][GRT_SRC_M_NUM][GRT_INTEG_NUM],
FILE *fstats, GRT_KernelFunc kerfunc);


20 changes: 10 additions & 10 deletions pygrt/C_extension/include/grt/common/integral.h
Original file line number Diff line number Diff line change
Expand Up @@ -25,10 +25,10 @@
*
*/
void grt_int_Pk(
MYREAL k, MYREAL r,
const MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
real_t k, real_t r,
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
bool calc_uir,
MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);



Expand All @@ -40,7 +40,7 @@ void grt_int_Pk(
* @param[out] tol Z、R、T分量结果
*/
void grt_merge_Pk(
const MYCOMPLEX sum_J[GRT_SRC_M_NUM][GRT_INTEG_NUM], MYCOMPLEX tol[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]);
const cplx_t sum_J[GRT_SRC_M_NUM][GRT_INTEG_NUM], cplx_t tol[GRT_SRC_M_NUM][GRT_CHANNEL_NUM]);



Expand All @@ -59,10 +59,10 @@ void grt_merge_Pk(
*
*/
void grt_int_Pk_filon(
MYREAL k, MYREAL r, bool iscos,
const MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
real_t k, real_t r, bool iscos,
const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
bool calc_uir,
MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);


/**
Expand All @@ -78,7 +78,7 @@ void grt_int_Pk_filon(
*
*/
void grt_int_Pk_sa_filon(
const MYREAL k3[3], MYREAL r,
const MYCOMPLEX QWV3[3][GRT_SRC_M_NUM][GRT_QWV_NUM],
const real_t k3[3], real_t r,
const cplx_t QWV3[3][GRT_SRC_M_NUM][GRT_QWV_NUM],
bool calc_uir,
MYCOMPLEX SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);
cplx_t SUM[GRT_SRC_M_NUM][GRT_INTEG_NUM]);
6 changes: 3 additions & 3 deletions pygrt/C_extension/include/grt/common/iostats.h
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@
* 记录其值主要用于参考其变化趋势。
*/
void grt_write_stats(
FILE *f0, MYREAL k, const MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]);
FILE *f0, real_t k, const cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM]);


/**
Expand Down Expand Up @@ -55,8 +55,8 @@ MYINT grt_extract_stats(FILE *bf0, FILE *af0);
*/
void grt_write_stats_ptam(
FILE *f0,
MYREAL Kpt[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX],
MYCOMPLEX Fpt[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX]);
real_t Kpt[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX],
cplx_t Fpt[GRT_SRC_M_NUM][GRT_INTEG_NUM][GRT_PTAM_PT_MAX]);


/**
Expand Down
4 changes: 2 additions & 2 deletions pygrt/C_extension/include/grt/common/kernel.h
Original file line number Diff line number Diff line change
Expand Up @@ -17,5 +17,5 @@
* 计算核函数的函数指针,动态与静态的接口一致
*/
typedef void (*GRT_KernelFunc) (
const GRT_MODEL1D *mod1d, MYCOMPLEX QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
bool calc_uiz, MYCOMPLEX QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]);
const GRT_MODEL1D *mod1d, cplx_t QWV[GRT_SRC_M_NUM][GRT_QWV_NUM],
bool calc_uiz, cplx_t QWV_uiz[GRT_SRC_M_NUM][GRT_QWV_NUM]);
54 changes: 27 additions & 27 deletions pygrt/C_extension/include/grt/common/matrix.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,12 @@
* @param[out] invM 逆矩阵
* @param[out] stats 状态代码,是否有除零错误,非0为异常值
*/
inline GCC_ALWAYS_INLINE void grt_cmat2x2_inv(const MYCOMPLEX M[2][2], MYCOMPLEX invM[2][2], MYINT *stats) {
MYCOMPLEX M00 = M[0][0];
MYCOMPLEX M01 = M[0][1];
MYCOMPLEX M10 = M[1][0];
MYCOMPLEX M11 = M[1][1];
MYCOMPLEX det = M00*M11 - M01*M10;
inline GCC_ALWAYS_INLINE void grt_cmat2x2_inv(const cplx_t M[2][2], cplx_t invM[2][2], MYINT *stats) {
cplx_t M00 = M[0][0];
cplx_t M01 = M[0][1];
cplx_t M10 = M[1][0];
cplx_t M11 = M[1][1];
cplx_t det = M00*M11 - M01*M10;
if ( det == 0.0 ){
// fprintf(stderr, "%.5e+%.5ej %.5e+%.5ej \n", creal(M[0][0]), cimag(M[0][0]), creal(M[0][1]), cimag(M[0][1]));
// fprintf(stderr, "%.5e+%.5ej %.5e+%.5ej \n", creal(M[1][0]), cimag(M[1][0]), creal(M[1][1]), cimag(M[1][1]));
Expand All @@ -50,7 +50,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_inv(const MYCOMPLEX M[2][2], MYCOMPLEX
* @param[in] M2 矩阵2
* @param[out] M 和矩阵
*/
inline GCC_ALWAYS_INLINE void grt_cmat2x2_add(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2][2], MYCOMPLEX M[2][2]){
inline GCC_ALWAYS_INLINE void grt_cmat2x2_add(const cplx_t M1[2][2], const cplx_t M2[2][2], cplx_t M[2][2]){
M[0][0] = M1[0][0] + M2[0][0];
M[0][1] = M1[0][1] + M2[0][1];
M[1][0] = M1[1][0] + M2[1][0];
Expand All @@ -64,7 +64,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_add(const MYCOMPLEX M1[2][2], const MY
* @param[in] M2 矩阵2
* @param[out] M 差矩阵, M1 - M2
*/
inline GCC_ALWAYS_INLINE void grt_cmat2x2_sub(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2][2], MYCOMPLEX M[2][2]){
inline GCC_ALWAYS_INLINE void grt_cmat2x2_sub(const cplx_t M1[2][2], const cplx_t M2[2][2], cplx_t M[2][2]){
M[0][0] = M1[0][0] - M2[0][0];
M[0][1] = M1[0][1] - M2[0][1];
M[1][0] = M1[1][0] - M2[1][0];
Expand All @@ -76,7 +76,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_sub(const MYCOMPLEX M1[2][2], const MY
*
* @param[in,out] M 差矩阵 I-M2
*/
inline GCC_ALWAYS_INLINE void grt_cmat2x2_one_sub(MYCOMPLEX M[2][2]){
inline GCC_ALWAYS_INLINE void grt_cmat2x2_one_sub(cplx_t M[2][2]){
M[0][0] = 1.0 - M[0][0];
M[0][1] = - M[0][1];
M[1][0] = - M[1][0];
Expand All @@ -90,9 +90,9 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_one_sub(MYCOMPLEX M[2][2]){
* @param[in] M2 矩阵2
* @param[out] M 积矩阵, M1 * M2
*/
inline GCC_ALWAYS_INLINE void grt_cmat2x2_mul(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2][2], MYCOMPLEX M[2][2]){
MYCOMPLEX M011, M012, M021, M022;
MYCOMPLEX M111, M112, M121, M122;
inline GCC_ALWAYS_INLINE void grt_cmat2x2_mul(const cplx_t M1[2][2], const cplx_t M2[2][2], cplx_t M[2][2]){
cplx_t M011, M012, M021, M022;
cplx_t M111, M112, M121, M122;
M011 = M1[0][0]; M012 = M1[0][1];
M021 = M1[1][0]; M022 = M1[1][1];
M111 = M2[0][0]; M112 = M2[0][1];
Expand All @@ -110,7 +110,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_mul(const MYCOMPLEX M1[2][2], const MY
* @param[in] k 常数
* @param[out] M 积矩阵, k * M2
*/
inline GCC_ALWAYS_INLINE void grt_cmat2x2_k(const MYCOMPLEX M1[2][2], MYCOMPLEX k0, MYCOMPLEX M[2][2]){
inline GCC_ALWAYS_INLINE void grt_cmat2x2_k(const cplx_t M1[2][2], cplx_t k0, cplx_t M[2][2]){
M[0][0] = M1[0][0] * k0;
M[0][1] = M1[0][1] * k0;
M[1][0] = M1[1][0] * k0;
Expand All @@ -124,8 +124,8 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_k(const MYCOMPLEX M1[2][2], MYCOMPLEX
* @param[in] M2 向量
* @param[out] M 积矩阵, M1 * M2
*/
inline GCC_ALWAYS_INLINE void grt_cmat2x1_mul(const MYCOMPLEX M1[2][2], const MYCOMPLEX M2[2], MYCOMPLEX M[2]){
MYCOMPLEX M00, M10;
inline GCC_ALWAYS_INLINE void grt_cmat2x1_mul(const cplx_t M1[2][2], const cplx_t M2[2], cplx_t M[2]){
cplx_t M00, M10;
M00 = M1[0][0]*M2[0] + M1[0][1]*M2[1];
M10 = M1[1][0]*M2[0] + M1[1][1]*M2[1];
M[0] = M00;
Expand All @@ -138,7 +138,7 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x1_mul(const MYCOMPLEX M1[2][2], const MY
* @param[in] M1 源矩阵
* @param[out] M2 目标矩阵
*/
inline GCC_ALWAYS_INLINE void grt_cmat2x2_assign(const MYCOMPLEX M1[2][2], MYCOMPLEX M2[2][2]){
inline GCC_ALWAYS_INLINE void grt_cmat2x2_assign(const cplx_t M1[2][2], cplx_t M2[2][2]){
M2[0][0] = M1[0][0];
M2[0][1] = M1[0][1];
M2[1][0] = M1[1][0];
Expand All @@ -155,9 +155,9 @@ inline GCC_ALWAYS_INLINE void grt_cmat2x2_assign(const MYCOMPLEX M1[2][2], MYCOM
* @param[in] M2 M2矩阵
* @param[out] M 积矩阵 M1 * M2
*/
inline GCC_ALWAYS_INLINE void grt_cmatmxn_mul(MYINT m1, MYINT n1, MYINT p1, const MYCOMPLEX M1[m1][n1], const MYCOMPLEX M2[n1][p1], MYCOMPLEX M[m1][p1]){
inline GCC_ALWAYS_INLINE void grt_cmatmxn_mul(MYINT m1, MYINT n1, MYINT p1, const cplx_t M1[m1][n1], const cplx_t M2[n1][p1], cplx_t M[m1][p1]){
MYINT m, n, k;
MYCOMPLEX M0[m1][p1];
cplx_t M0[m1][p1];
for(m=0; m<m1; ++m){
for(n=0; n<p1; ++n){
M0[m][n] = 0.0;
Expand All @@ -167,7 +167,7 @@ inline GCC_ALWAYS_INLINE void grt_cmatmxn_mul(MYINT m1, MYINT n1, MYINT p1, cons
}
}

// memcpy(M, M0, sizeof(MYCOMPLEX)*m1*p1);
// memcpy(M, M0, sizeof(cplx_t)*m1*p1);
for(m=0; m<m1; ++m){
for(n=0; n<p1; ++n){
M[m][n] = M0[m][n];
Expand All @@ -183,16 +183,16 @@ inline GCC_ALWAYS_INLINE void grt_cmatmxn_mul(MYINT m1, MYINT n1, MYINT p1, cons
* @param[in] M1 M1矩阵
* @param[out] M2 M2矩阵 (M1^T)
*/
inline GCC_ALWAYS_INLINE void grt_cmatmxn_transpose(MYINT m1, MYINT n1, const MYCOMPLEX M1[m1][n1], MYCOMPLEX M2[n1][m1]){
inline GCC_ALWAYS_INLINE void grt_cmatmxn_transpose(MYINT m1, MYINT n1, const cplx_t M1[m1][n1], cplx_t M2[n1][m1]){
MYINT m, n;
MYCOMPLEX M0[n1][m1];
cplx_t M0[n1][m1];
for(m=0; m<m1; ++m){
for(n=0; n<n1; ++n){
M0[n][m] = M1[m][n];
}
}

// memcpy(M, M0, sizeof(MYCOMPLEX)*m1*p1);
// memcpy(M, M0, sizeof(cplx_t)*m1*p1);
for(m=0; m<m1; ++m){
for(n=0; n<n1; ++n){
M2[n][m] = M0[n][m];
Expand All @@ -212,7 +212,7 @@ inline GCC_ALWAYS_INLINE void grt_cmatmxn_transpose(MYINT m1, MYINT n1, const MY
* @param[in] ln 子矩阵列数
* @param[out] Q 子矩阵
*/
inline GCC_ALWAYS_INLINE void grt_cmatmxn_block(MYINT m1, MYINT n1, const MYCOMPLEX M[m1][n1], MYINT im, MYINT in, MYINT lm, MYINT ln, MYCOMPLEX Q[lm][ln]){
inline GCC_ALWAYS_INLINE void grt_cmatmxn_block(MYINT m1, MYINT n1, const cplx_t M[m1][n1], MYINT im, MYINT in, MYINT lm, MYINT ln, cplx_t Q[lm][ln]){
for(MYINT m=0; m<lm; ++m){
for(MYINT n=0; n<ln; ++n){
Q[m][n] = M[im+m][in+n];
Expand All @@ -233,7 +233,7 @@ inline GCC_ALWAYS_INLINE void grt_cmatmxn_block(MYINT m1, MYINT n1, const MYCOMP
* @param[in] ln 子矩阵列数
* @param[out] Q 子矩阵
*/
inline GCC_ALWAYS_INLINE void grt_cmatmxn_block_assign(MYINT m1, MYINT n1, MYCOMPLEX M[m1][n1], MYINT im, MYINT in, MYINT lm, MYINT ln, const MYCOMPLEX Q[lm][ln]){
inline GCC_ALWAYS_INLINE void grt_cmatmxn_block_assign(MYINT m1, MYINT n1, cplx_t M[m1][n1], MYINT im, MYINT in, MYINT lm, MYINT ln, const cplx_t Q[lm][ln]){
for(MYINT m=0; m<lm; ++m){
for(MYINT n=0; n<ln; ++n){
M[im+m][in+n] = Q[m][n];
Expand All @@ -250,7 +250,7 @@ inline GCC_ALWAYS_INLINE void grt_cmatmxn_block_assign(MYINT m1, MYINT n1, MYCOM
* @param[in] M1 M1矩阵
*
*/
inline GCC_ALWAYS_INLINE void grt_cmatmxn_print(MYINT m1, MYINT n1, const MYCOMPLEX M1[m1][n1]){
inline GCC_ALWAYS_INLINE void grt_cmatmxn_print(MYINT m1, MYINT n1, const cplx_t M1[m1][n1]){
for(MYINT i=0; i<m1; ++i){
for(MYINT j=0; j<n1; ++j){
fprintf(stderr, " %15.5e + J%-15.5e ", creal(M1[i][j]), cimag(M1[i][j]));
Expand All @@ -267,8 +267,8 @@ inline GCC_ALWAYS_INLINE void grt_cmatmxn_print(MYINT m1, MYINT n1, const MYCOMP
* @param[in] M2 向量
* @param[out] M 积矩阵, M1 * M2
*/
inline GCC_ALWAYS_INLINE void grt_cmatmx1_mul(MYINT m1, MYINT n1, const MYCOMPLEX M1[m1][n1], const MYCOMPLEX M2[n1], MYCOMPLEX M[n1]){
MYCOMPLEX M0[n1];
inline GCC_ALWAYS_INLINE void grt_cmatmx1_mul(MYINT m1, MYINT n1, const cplx_t M1[m1][n1], const cplx_t M2[n1], cplx_t M[n1]){
cplx_t M0[n1];
for(MYINT i=0; i<m1; ++i){
M0[i] = 0.0;
for(MYINT j=0; j<n1; ++j){
Expand Down
Loading