#include <assert.h>#include <math.h>#include "mex.h"#include <sys/time.h>#include <string.h>/* // export PATH=$PATH:/softs/stow/matlab-14/binsetenv PATH ${PATH}:/softs/stow/matlab-14/binmex CFLAGS=-Wall em_hdda_c.ccp em_hdda_c.mexglx ..*//******************************************************** * prototypes for BLAS and LAPACK functions */#define doublereal double#define integer intextern int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a, integer *lda, doublereal *w, doublereal *work, integer *lwork, integer *info);extern int dsyevx_(char *jobz, char *range, char *uplo, integer *n,  doublereal *a, integer *lda, doublereal *vl, doublereal *vu, integer *  il, integer *iu, doublereal *abstol, integer *m, doublereal *w,  doublereal *z__, integer *ldz, doublereal *work, integer *lwork,  integer *iwork, integer *ifail, integer *info);extern int dgemm_(char* transA, char* transB, integer* M, integer* N, integer* K,       doublereal* alpha,       doublereal* A, integer* lda,       doublereal* B, integer* ldb,       doublereal* beta,       doublereal* C, integer* ldc);extern intdsyrk_(char* uplo, char* trans, integer* N, integer* K,       doublereal* alpha,       doublereal* A, integer* lda,       doublereal* beta,       doublereal* C, integer* ldc);#undef doublereal#undef integer/******************************************************** * some useful definitions *//*#define NEWA(type,n) (type*)mxMalloc((printf("mall %d %d %d\n",n,sizeof(type),__LINE__),sizeof(type)*(n)))*/#define NEWA(type,n) (type*)mxMalloc(sizeof(type)*(n))#define FREE(p) mxFree(p)#define ASSERT(v,comment) mxAssert(v,comment) double getMilliSecs() {  struct timeval tv;  gettimeofday(&tv,NULL);  return tv.tv_sec*1e3+tv.tv_usec*1e-3;}#define TICXX lines[n_t]=__LINE__; times[n_t++]+=getMilliSecs()#define TIC/******************************************************** * Subfunctions */#define MAXD ((p)/2)/* common flags * any combination of the 3 is acceptable */#define Common_a 1#define Common_b 2#define Common_d 4#define Common_p 8#define Free_a 16/* samples k values among 0..n-1 with replacement  * result in cls(k) */void randsample(int n,int k,int *cls,double *rand_input) {  int i;  if(rand_input)     for(i=0;i<k;i++)       cls[i]=ceil(rand_input[i]*n)-1;  else     for(i=0;i<k;i++)       cls[i]=random()%n;}/* access macros for s and t (the same for all the functions) */#define T(i,j) t[(i)+(j)*n]#define S(i,j) s[(i)+(j)*n]#define A(i,j) a[(i)+(j)*MAXD]/* lwork for compute_params*/static int compute_params_lwork(int k,int n,int p) {  int lwork_eig=(32+3)*p;    return p*k    +p*p    +p*p    +lwork_eig    +p*k    +p*p*k    +p    +2*n    +2*p*n    +k    +k    +k    +k*MAXD/*+k*/    +p;}/* computes intrinsic dimensions and estimator parameters  * all matrices are stored by columns (