ivan@140: /**************************************************************************
ivan@140:  *
ivan@140:  * File name: ompcore.c
ivan@140:  *
ivan@140:  * Ron Rubinstein
ivan@140:  * Computer Science Department
ivan@140:  * Technion, Haifa 32000 Israel
ivan@140:  * ronrubin@cs
ivan@140:  *
ivan@140:  * Last Updated: 25.8.2009
ivan@140:  *
ivan@140:  *************************************************************************/
ivan@140: 
ivan@140: 
ivan@140: #include "ompcore.h"
ivan@140: #include "omputils.h"
ivan@140: #include "ompprof.h"
ivan@140: #include "myblas.h"
ivan@140: #include <math.h>
ivan@140: #include <string.h>
ivan@140: 
ivan@140: 
ivan@140: 
ivan@140: /******************************************************************************
ivan@140:  *                                                                            *
ivan@140:  *                           Batch-OMP Implementation                         *
ivan@140:  *                                                                            *
ivan@140:  ******************************************************************************/  
ivan@140: 
ivan@140: mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L,
ivan@140:                  int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp)
ivan@140: {
ivan@140:   
ivan@140:   profdata pd;
ivan@140:   mxArray *Gamma;
ivan@140:   mwIndex i, j, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count;
ivan@140:   mwSize allocated_coefs, allocated_cols;
ivan@140:   int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms;
ivan@140:   double *alpha, *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; 
ivan@140:   double eps2, resnorm, delta, deltaprev, secs_remain;
ivan@140:   int mins_remain, hrs_remain;
ivan@140:   clock_t lastprint_time, starttime;
ivan@140:  
ivan@140:   
ivan@140:   
ivan@140:   /*** status flags ***/
ivan@140:   
ivan@140:   DtX_specified = (DtX!=0);   /* indicates whether D'*x was provided */
ivan@140:   XtX_specified = (XtX!=0);   /* indicates whether sum(x.*x) was provided */
ivan@140:   
ivan@140:   standardomp = (G==0);       /* batch-omp or standard omp are selected depending on availability of G */
ivan@140:   batchomp = !standardomp;
ivan@140:   
ivan@140:   
ivan@140:   
ivan@140:   /*** allocate output matrix ***/
ivan@140:   
ivan@140:   
ivan@140:   if (gamma_mode == FULL_GAMMA) {
ivan@140:     
ivan@140:     /* allocate full matrix of size m X L */
ivan@140:     
ivan@140:     Gamma = mxCreateDoubleMatrix(m, L, mxREAL);
ivan@140:     gammaPr = mxGetPr(Gamma);
ivan@140:     gammaIr = 0;
ivan@140:     gammaJc = 0;
ivan@140:   }
ivan@140:   else {
ivan@140:     
ivan@140:     /* allocate sparse matrix with room for allocated_coefs nonzeros */
ivan@140:     
ivan@140:     /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */
ivan@140:     allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T;
ivan@140:     Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL);
ivan@140:     gammaPr = mxGetPr(Gamma);
ivan@140:     gammaIr = mxGetIr(Gamma);
ivan@140:     gammaJc = mxGetJc(Gamma);
ivan@140:     gamma_count = 0;
ivan@140:     gammaJc[0] = 0;
ivan@140:   }
ivan@140:   
ivan@140:   
ivan@140:   /*** helper arrays ***/
ivan@140:   
ivan@140:   alpha = (double*)mxMalloc(m*sizeof(double));        /* contains D'*residual */
ivan@140:   ind = (mwIndex*)mxMalloc(n*sizeof(mwIndex));        /* indices of selected atoms */
ivan@140:   selected_atoms = (int*)mxMalloc(m*sizeof(int));     /* binary array with 1's for selected atoms */
ivan@140:   c = (double*)mxMalloc(n*sizeof(double));            /* orthogonal projection result */
ivan@140:   
ivan@140:   /* current number of columns in Dsub / Gsub / Lchol */
ivan@140:   allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T;
ivan@140:   
ivan@140:   /* Cholesky decomposition of D_I'*D_I */
ivan@140:   Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double));
ivan@140: 
ivan@140:   /* temporary vectors for various computations */
ivan@140:   tempvec1 = (double*)mxMalloc(m*sizeof(double));
ivan@140:   tempvec2 = (double*)mxMalloc(m*sizeof(double));
ivan@140:   
ivan@140:   if (batchomp) {
ivan@140:     /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */
ivan@140:     Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double));
ivan@140:   }
ivan@140:   else {
ivan@140:     /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */
ivan@140:     Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double));
ivan@140:     
ivan@140:     /* stores the residual */
ivan@140:     r = (double*)mxMalloc(n*sizeof(double));        
ivan@140:   }
ivan@140:   
ivan@140:   if (!DtX_specified) {
ivan@140:     /* contains D'*x for the current signal */
ivan@140:     DtX = (double*)mxMalloc(m*sizeof(double));  
ivan@140:   }
ivan@140:   
ivan@140:   
ivan@140:   
ivan@140:   /*** initializations for error omp ***/
ivan@140:   
ivan@140:   if (erroromp) {
ivan@140:     eps2 = eps*eps;        /* compute eps^2 */
ivan@140:     if (T<0 || T>n) {      /* unspecified max atom num - set max atoms to n */
ivan@140:       T = n;
ivan@140:     }
ivan@140:   }
ivan@140:   
ivan@140:   
ivan@140:   
ivan@140:   /*** initialize timers ***/
ivan@140:   
ivan@140:   initprofdata(&pd);             /* initialize profiling counters */
ivan@140:   starttime = clock();           /* record starting time for eta computations */
ivan@140:   lastprint_time = starttime;    /* time of last status display */
ivan@140:   
ivan@140:   
ivan@140:   
ivan@140:   /**********************   perform omp for each signal   **********************/
ivan@140:   
ivan@140:   
ivan@140:   
ivan@140:   for (signum=0; signum<L; ++signum) {
ivan@140:     
ivan@140:     
ivan@140:     /* initialize residual norm and deltaprev for error-omp */
ivan@140:     
ivan@140:     if (erroromp) {
ivan@140:       if (XtX_specified) {
ivan@140:         resnorm = XtX[signum];
ivan@140:       }
ivan@140:       else {
ivan@140:         resnorm = dotprod(x+n*signum, x+n*signum, n);
ivan@140:         addproftime(&pd, XtX_TIME);
ivan@140:       }
ivan@140:       deltaprev = 0;     /* delta tracks the value of gamma'*G*gamma */
ivan@140:     }
ivan@140:     else {
ivan@140:       /* ignore residual norm stopping criterion */
ivan@140:       eps2 = 0;
ivan@140:       resnorm = 1;
ivan@140:     }
ivan@140:     
ivan@140:     
ivan@140:     if (resnorm>eps2 && T>0) {
ivan@140:       
ivan@140:       /* compute DtX */
ivan@140:       
ivan@140:       if (!DtX_specified) {
ivan@140:         matT_vec(1, D, x+n*signum, DtX, n, m);
ivan@140:         addproftime(&pd, DtX_TIME);
ivan@140:       }
ivan@140:       
ivan@140:       
ivan@140:       /* initialize alpha := DtX */
ivan@140:       
ivan@140:       memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double));
ivan@140:       
ivan@140:       
ivan@140:       /* mark all atoms as unselected */
ivan@140:       
ivan@140:       for (i=0; i<m; ++i) {
ivan@140:         selected_atoms[i] = 0;
ivan@140:       }
ivan@140:       
ivan@140:     }
ivan@140:     
ivan@140: 
ivan@140:     /* main loop */
ivan@140:     
ivan@140:     i=0;
ivan@140:     while (resnorm>eps2 && i<T) {
ivan@140: 
ivan@140:       /* index of next atom */
ivan@140:       
ivan@140:       pos = maxabs(alpha, m);
ivan@140:       addproftime(&pd, MAXABS_TIME);
ivan@140:       
ivan@140:       
ivan@140:       /* stop criterion: selected same atom twice, or inner product too small */
ivan@140:       
ivan@140:       if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) {
ivan@140:         break;
ivan@140:       }
ivan@140:       
ivan@140:       
ivan@140:       /* mark selected atom */
ivan@140:       
ivan@140:       ind[i] = pos;
ivan@140:       selected_atoms[pos] = 1;
ivan@140:       
ivan@140:       
ivan@140:       /* matrix reallocation */
ivan@140:       
ivan@140:       if (erroromp && i>=allocated_cols) {
ivan@140:         
ivan@140:         allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01);
ivan@140:         
ivan@140:         Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double));
ivan@140:         
ivan@140:         batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) :
ivan@140:                    (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ;
ivan@140:       }
ivan@140:       
ivan@140:       
ivan@140:       /* append column to Gsub or Dsub */
ivan@140:       
ivan@140:       if (batchomp) {
ivan@140:         memcpy(Gsub+i*m, G+pos*m, m*sizeof(double));
ivan@140:       }
ivan@140:       else {
ivan@140:         memcpy(Dsub+i*n, D+pos*n, n*sizeof(double));
ivan@140:       }
ivan@140:       
ivan@140:       
ivan@140:       /*** Cholesky update ***/
ivan@140:       
ivan@140:       if (i==0) {
ivan@140:         *Lchol = 1;
ivan@140:       }
ivan@140:       else {
ivan@140:         
ivan@140:         /* incremental Cholesky decomposition: compute next row of Lchol */
ivan@140:         
ivan@140:         if (standardomp) {
ivan@140:           matT_vec(1, Dsub, D+n*pos, tempvec1, n, i);      /* compute tempvec1 := Dsub'*d where d is new atom */
ivan@140:           addproftime(&pd, DtD_TIME);
ivan@140:         }
ivan@140:         else {
ivan@140:           vec_assign(tempvec1, Gsub+i*m, ind, i);          /* extract tempvec1 := Gsub(ind,i) */
ivan@140:         }
ivan@140:         backsubst('L', Lchol, tempvec1, tempvec2, n, i);   /* compute tempvec2 = Lchol \ tempvec1 */
ivan@140:         for (j=0; j<i; ++j) {                              /* write tempvec2 to end of Lchol */
ivan@140:           Lchol[j*n+i] = tempvec2[j];
ivan@140:         }
ivan@140:         
ivan@140:         /* compute Lchol(i,i) */
ivan@140:         sum = 0;
ivan@140:         for (j=0; j<i; ++j) {         /* compute sum of squares of last row without Lchol(i,i) */
ivan@140:           sum += SQR(Lchol[j*n+i]);
ivan@140:         }
ivan@140:         if ( (1-sum) <= 1e-14 ) {     /* Lchol(i,i) is zero => selected atoms are dependent */
ivan@140:           break;
ivan@140:         }
ivan@140:         Lchol[i*n+i] = sqrt(1-sum);
ivan@140:       }
ivan@140:       
ivan@140:       addproftime(&pd, LCHOL_TIME);
ivan@140: 
ivan@140:       i++;
ivan@140:       
ivan@140:       
ivan@140:       /* perform orthogonal projection and compute sparse coefficients */
ivan@140:       
ivan@140:       vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i);   /* extract tempvec1 = DtX(ind) */
ivan@140:       cholsolve('L', Lchol, tempvec1, c, n, i);                     /* solve LL'c = tempvec1 for c */
ivan@140:       addproftime(&pd, COMPCOEF_TIME);
ivan@140:       
ivan@140: 
ivan@140:       /* update alpha = D'*residual */
ivan@140:       
ivan@140:       if (standardomp) {
ivan@140:         mat_vec(-1, Dsub, c, r, n, i);             /* compute r := -Dsub*c */
ivan@140:         vec_sum(1, x+n*signum, r, n);              /* compute r := x+r */
ivan@140:         
ivan@140:         
ivan@140:         /*memcpy(r, x+n*signum, n*sizeof(double));   /* assign r := x */
ivan@140:         /*mat_vec1(-1, Dsub, c, 1, r, n, i);         /* compute r := r-Dsub*c */
ivan@140:         
ivan@140:         addproftime(&pd, COMPRES_TIME);
ivan@140:         matT_vec(1, D, r, alpha, n, m);            /* compute alpha := D'*r */
ivan@140:         addproftime(&pd, DtR_TIME);
ivan@140:         
ivan@140:         /* update residual norm */
ivan@140:         if (erroromp) {
ivan@140:           resnorm = dotprod(r, r, n);
ivan@140:           addproftime(&pd, UPDATE_RESNORM_TIME);
ivan@140:         }
ivan@140:       }
ivan@140:       else {
ivan@140:         mat_vec(1, Gsub, c, tempvec1, m, i);                              /* compute tempvec1 := Gsub*c */
ivan@140:         memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double));    /* set alpha = D'*x */
ivan@140:         vec_sum(-1, tempvec1, alpha, m);                                  /* compute alpha := alpha - tempvec1 */
ivan@140:         addproftime(&pd, UPDATE_DtR_TIME);
ivan@140:         
ivan@140:         /* update residual norm */
ivan@140:         if (erroromp) {
ivan@140:           vec_assign(tempvec2, tempvec1, ind, i);      /* assign tempvec2 := tempvec1(ind) */
ivan@140:           delta = dotprod(c,tempvec2,i);               /* compute c'*tempvec2 */
ivan@140:           resnorm = resnorm - delta + deltaprev;       /* residual norm update */
ivan@140:           deltaprev = delta;
ivan@140:           addproftime(&pd, UPDATE_RESNORM_TIME);
ivan@140:         }
ivan@140:       }
ivan@140:     }
ivan@140:     
ivan@140:     
ivan@140:     /*** generate output vector gamma ***/
ivan@140: 
ivan@140:     if (gamma_mode == FULL_GAMMA) {    /* write the coefs in c to their correct positions in gamma */
ivan@140:       for (j=0; j<i; ++j) {
ivan@140:         gammaPr[m*signum + ind[j]] = c[j];
ivan@140:       }
ivan@140:     }
ivan@140:     else {
ivan@140:       /* sort the coefs by index before writing them to gamma */
ivan@140:       quicksort(ind,c,i);
ivan@140:       addproftime(&pd, INDEXSORT_TIME);
ivan@140:       
ivan@140:       /* gamma is full - reallocate */
ivan@140:       if (gamma_count+i >= allocated_coefs) {
ivan@140:         
ivan@140:         while(gamma_count+i >= allocated_coefs) {
ivan@140:           allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01);
ivan@140:         }
ivan@140:         
ivan@140:         mxSetNzmax(Gamma, allocated_coefs);
ivan@140:         mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double)));
ivan@140:         mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex)));
ivan@140:         
ivan@140:         gammaPr = mxGetPr(Gamma);
ivan@140:         gammaIr = mxGetIr(Gamma);
ivan@140:       }
ivan@140:       
ivan@140:       /* append coefs to gamma and update the indices */
ivan@140:       for (j=0; j<i; ++j) {
ivan@140:         gammaPr[gamma_count] = c[j];
ivan@140:         gammaIr[gamma_count] = ind[j];
ivan@140:         gamma_count++;
ivan@140:       }
ivan@140:       gammaJc[signum+1] = gammaJc[signum] + i;
ivan@140:     }
ivan@140:     
ivan@140:     
ivan@140:     
ivan@140:     /*** display status messages ***/
ivan@140:     
ivan@140:     if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta)
ivan@140:     {
ivan@140:       lastprint_time = clock();
ivan@140:       
ivan@140:       /* estimated remainig time */
ivan@140:       secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) ,
ivan@140:         &hrs_remain, &mins_remain, &secs_remain);
ivan@140:       
ivan@140:       mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n",        
ivan@140:         signum+1, L, hrs_remain, mins_remain, secs_remain);
ivan@140:       mexEvalString("drawnow;");
ivan@140:     }
ivan@140:     
ivan@140:   }
ivan@140:   
ivan@140:   /* end omp */
ivan@140:   
ivan@140:   
ivan@140:   
ivan@140:   /*** print final messages ***/
ivan@140:   
ivan@140:   if (msg_delta>0) {
ivan@140:     mexPrintf("omp: signal %d / %d\n", signum, L);
ivan@140:   }
ivan@140:   
ivan@140:   if (profile) {
ivan@140:     printprofinfo(&pd, erroromp, batchomp, L);
ivan@140:   }
ivan@140:   
ivan@140:   
ivan@140:   
ivan@140:   /* free memory */
ivan@140:   
ivan@140:   if (!DtX_specified) {
ivan@140:     mxFree(DtX);
ivan@140:   }
ivan@140:   if (standardomp) {
ivan@140:     mxFree(r);
ivan@140:     mxFree(Dsub);
ivan@140:   }
ivan@140:   else {
ivan@140:     mxFree(Gsub);
ivan@140:   }  
ivan@140:   mxFree(tempvec2);
ivan@140:   mxFree(tempvec1);
ivan@140:   mxFree(Lchol);
ivan@140:   mxFree(c);
ivan@140:   mxFree(selected_atoms);
ivan@140:   mxFree(ind);
ivan@140:   mxFree(alpha);
ivan@140:   
ivan@140:   return Gamma;
ivan@140: }