使用 malloc 进行复杂 Cholesky 分解的 Matlab Mex C

问题描述

我是一个学习 Matrix、C 和 Matlab 的新手。我尝试使用 Matlab mexFunction 完成复杂的 Cholesky 分解。当我使用 malloc 创建数组时，输出是错误的。我的问题在哪里？

如果我不使用 malloc，我可以得到正确的答案。

如果有人能帮助我，我将不胜感激。

#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 );
}

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）