const Real EPS = 1e-15;
// Additional (not at this place, in R file)
Real* gam2 = (Real*)malloc(k*sizeof(Real));
- Real* nY21 = (Real*)malloc(n*m*k*sizeof(Real));
Real* sqNorm2 = (Real*)malloc(k*sizeof(Real));
gsl_matrix* matrix = gsl_matrix_alloc(m, m);
gsl_permutation* permutation = gsl_permutation_alloc(m);
dotProduct += X2[ai(i,u,r,n,p,k)] * phi[ai(u,mm,r,p,m,k)];
//ps1[i,mm,r] = Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r])
ps1[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * dotProduct;
- nY21[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * Y2[ai(i,mm,r,n,m,k)];
}
//ps[mm,r] = sum(ps1[,mm,r])
Real sumPs1 = 0.;
for (int u=0; u<n; u++)
sumPs1 += ps1[ai(u,mm,r,n,m,k)];
ps[mi(mm,r,m,k)] = sumPs1;
- //nY2[mm,r] = sum(nY21[,mm,r])
- Real sumNy21 = 0.;
+ //nY2[mm,r] = sum(nY2[,mm,r])
+ Real sumNy2 = 0.;
for (int u=0; u<n; u++)
- sumNy21 += nY21[ai(u,mm,r,n,m,k)];
- nY2[mi(mm,r,m,k)] = sumNy21;
+ sumNy2 += nY2[ai(u,mm,r,n,m,k)];
+ nY2[mi(mm,r,m,k)] = sumNy2;
//rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r]))) / (2*nY2[mm,r])
rho[ai(mm,mm,r,m,m,k)] = ( ps[mi(mm,r,m,k)] + sqrt( ps[mi(mm,r,m,k)]*ps[mi(mm,r,m,k)]
+ 4*nY2[mi(mm,r,m,k)] * gam2[r] ) ) / (2*nY2[mi(mm,r,m,k)]);
ite++;
}
+ //TODO: affec = apply(gam, 1,which.max) à traduire, et retourner affec aussi.
+
//free memory
free(b);
free(gam);
free(ps);
free(nY2);
free(ps1);
- free(nY21);
free(Gram2);
free(ps2);
gsl_matrix_free(matrix);