Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <b.luong@fogale.findmycountry>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Can this run any faster?
Date: Sat, 13 Apr 2013 09:44:08 +0000 (UTC)
Organization: FOGALE nanotech
Lines: 59
Message-ID: <kkb9d8$ok5$1@newscl01ah.mathworks.com>
References: <kk24q4$3j5$1@newscl01ah.mathworks.com> <kk9lr7$782$1@newscl01ah.mathworks.com> <kka7p6$smp$1@newscl01ah.mathworks.com> <kkb8f9$m7v$1@newscl01ah.mathworks.com>
Reply-To: "Bruno Luong" <b.luong@fogale.findmycountry>
NNTP-Posting-Host: www-01-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1365846248 25221 172.30.248.46 (13 Apr 2013 09:44:08 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sat, 13 Apr 2013 09:44:08 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 390839
Xref: news.mathworks.com comp.soft-sys.matlab:793430

Sorry I make a wrong assumption in the normalization. Here is the correct one:

/* MATLAB MEX BCJR_exp_loop_mex.c
 
 CALLING: >> s = BCJR_exp_loop_mex(e, t1, t2, c);
 COMPILE: >>mex -O BCJR_exp_loop_mex.c
 ATTENTION: e() will be modified inplace
 */

#include <math.h>
#include "mex.h"

#define E_INOUT prhs[0]
#define T1_UINT32 prhs[1]
#define T2_UINT32 prhs[2]
#define C prhs[3]
#define S_OUT plhs[0]
void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])
{
    int m, n, i, j;
    int *t1_p, *t2_p, *t1, *t2;
    double *e, *ei, *en, *c, *s, r, t;
    
    n = mxGetM(E_INOUT);
    m = mxGetN(E_INOUT);
    t1_p = (int*)mxGetData(T1_UINT32);
    t2_p = (int*)mxGetData(T2_UINT32);
    e = mxGetPr(E_INOUT);
    c = mxGetPr(C);
    
    S_OUT = mxCreateDoubleMatrix(1, m, mxREAL);
    s = mxGetPr(S_OUT);
    ei = e + n;
    en = ei-1;
    e = e-1;
    for (i=m-1; i--; )
    {
        t1 = t1_p;
        t2 = t2_p;

        for (j=32; j--; )
        {
            r = *(c++) * *(e + *(t1++));
            *(en + *(t2++)) += r;
            t += r;
        }
        
        t = 0;
        for (j=n; j--; ) t += ei[j];
        for (j=n; j--; ) *(ei++) /= t;
        
        *(s+1) = *s + log(t);
        
        s++;        
        e += n;
        en += n;

    }
} /* End of BCJR_exp_loop_mex.c */