I am a new guy learning Matrix, C and Matlab. I try to complete the Complex Cholesky Decomposition with Matlab mexFunction. When I use malloc to create the array , the output is wrong. Where is my problem ?
If I don't use malloc, I can get the right answer.
I would really appreciate if someone could help me .
#include "mex.h"
#include <stdlib.h>
#include <math.h>
typedef struct{
double *real ;
double *imag ;
}Complex;
void cholesky (Complex A ,Complex Ainv, int n){
Complex L;
L.real=malloc(n*n*sizeof(double));
L.imag=malloc(n*n*sizeof(double));
int i,j,k ;
double sumR ;
double sumI ;
for (i=0; i<n; i++){
for (j=0; j<(i+1); j++) {
sumR = 0 ;
sumI = 0 ;
if (i == j){
for (k=0; k<j; k++ ){
sumR += ( L.real[i+k*n] * L.real[j+k*n]
+ L.imag[i+k*n] * L.imag[j+k*n] );
}
L.real[i+j*n]=sqrt(A.real[i+j*n] - sumR) ;
L.imag[i+j*n]=0 ;
}
else{
for (k=0; k<j; k++ ){
sumR += ( L.real[i+k*n] * L.real[j+k*n]
+ L.imag[i+k*n] * L.imag[j+k*n] );
sumI += (L.imag[i+k*n] * L.real[j+k*n]
- L.real[i+k*n] * L.imag[j+k*n] );
}
L.real[i+j*n] = ((A.imag[i+j*n] - sumR)/L.real[i+j*n]);
L.imag[i+j*n] = ((A.imag[i+j*n] - sumI)/L.real[i+j*n]);
}
}
}
for (i=0; i < n; i++){
for (j=0; j<(i+1) ; j++) {
Ainv.real[i+j*n] = L.real [ i+n*j ] ;
Ainv.imag[i+j*n] = L.imag [ i+n*j ] ;
}
}
free(L.real);
free(L.imag);
}
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
Complex A;
Complex Ainv;
int N ;
A.real = mxGetPr(prhs[0]);
A.imag = mxGetPi(prhs[0]);
N = mxGetN(prhs[0]);
plhs[0] = mxCreateDoubleMatrix(N,N,mxCOMPLEX);
Ainv.real = mxGetPr(plhs[0]);
Ainv.imag = mxGetPi(plhs[0]);
cholesky ( A, Ainv, N );
}
question from:
https://stackoverflow.com/questions/65896329/matlab-mex-c-of-complex-cholesky-decomposition-with-malloc 与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…