Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
210 views
in Technique[技术] by (71.8m points)

Matlab Mex C of Complex Cholesky Decomposition with malloc

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

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)
Waitting for answers

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...