THESE/codes/graphe/Ncut_9/sparsifyc.cpp

   1 /*=================================================================
   2 * syntax: SPMX = SPARSIFY(MX, THRES)
   3 *
   4 * SPARSIFY  - sparsify the input matrix, i.e. ignore the values
   5 *                       of the matrix which     are below a threshold
   6 *
   7 *                       Input:  - MX: m-by-n matrix (sparse or full)
   8 *                                       - THRES: threshold value (double)
   9 *
  10 *                       Output: - SPMX: m-by-n sparse matrix only with values
  11 *                                       whose absolut value is above the given threshold
  12 *
  13 *                       Written by Mirko Visontai (10/24/2003)
  14 *=================================================================*/
  15
  16
  17 #include <math.h>
  18 #include "mex.h"
  19
  20 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  21 {
  22     /* Declare variable */
  23     int i,m,n,nzmax,newnnz,col,processed,passed;
  24         int starting_row_index, current_row_index, stopping_row_index;
  25     double *in_pr,*in_pi,*out_pr,*out_pi;
  26     mwIndex *in_ir,*in_jc,*out_ir,*out_jc;
  27         double thres;
  28
  29     /* Check for proper number of input and output arguments */
  30     if ((nlhs != 1) || (nrhs != 2)){
  31                 mexErrMsgTxt("usage: SPMX = SPARSIFY(MX, THRES).");
  32     }
  33         /* if matrix is complex threshold the norm of the numbers */
  34         if (mxIsComplex(prhs[0])){
  35                 /* Check data type of input argument  */
  36                 if (mxIsSparse(prhs[0])){
  37
  38                         /* read input */
  39                         in_pr   = mxGetPr(prhs[0]);
  40                         in_pi   = mxGetPi(prhs[0]);
  41                         in_ir   = mxGetIr(prhs[0]);
  42                         in_jc   = mxGetJc(prhs[0]);
  43                         nzmax   = mxGetNzmax(prhs[0]);
  44                         m               = mxGetM(prhs[0]);
  45                         n               = mxGetN(prhs[0]);
  46                         thres   = mxGetScalar(prhs[1]);
  47
  48                         /* Count new nonzeros */
  49                         newnnz=0;
  50                         for(i=0; i<nzmax; i++){
  51                                 if (sqrt(in_pr[i]*in_pr[i] + in_pi[i]*in_pi[i])>thres) {newnnz++;}
  52                         }
  53
  54                         if (newnnz>0){
  55                                 /* create output */
  56                                 plhs[0]         = mxCreateSparse(m,n,newnnz,mxCOMPLEX);
  57                                 if (plhs[0]==NULL)
  58                                         mexErrMsgTxt("Could not allocate enough memory!\n");
  59                                 out_pr          = mxGetPr(plhs[0]);
  60                                 out_pi          = mxGetPr(plhs[0]);
  61                                 out_ir          = mxGetIr(plhs[0]);
  62                                 out_jc          = mxGetJc(plhs[0]);
  63                                 passed          = 0;
  64                                 out_jc[0]       = 0;
  65                                 for (col=0; col<n; col++){
  66                                         starting_row_index = in_jc[col];
  67                                         stopping_row_index = in_jc[col+1];
  68                                         out_jc[col+1] = out_jc[col];
  69                                         if (starting_row_index == stopping_row_index)
  70                                                 continue;
  71                                         else {
  72                                                 for (current_row_index = starting_row_index;
  73                                                         current_row_index < stopping_row_index;
  74                                                         current_row_index++)  {
  75                                                                 if (sqrt(in_pr[current_row_index]*in_pr[current_row_index] +
  76                                                                              in_pi[current_row_index]*in_pi[current_row_index] ) > thres){
  77
  78                                                                         out_pr[passed]=in_pr[current_row_index];
  79                                                                         out_pi[passed]=in_pi[current_row_index];
  80                                                                         out_ir[passed]=in_ir[current_row_index];
  81                                                                         out_jc[col+1] = out_jc[col+1]+1;
  82                                                                         passed++;
  83                                                                 }
  84                                                 }
  85                                         }
  86                                 }
  87                         }
  88                         else{
  89                                 plhs[0] = mxCreateSparse(m,n,0,mxCOMPLEX);
  90                         }
  91                 }
  92                 else{ /* for full matrices */
  93                         /* read input */
  94                         in_pr   = mxGetPr(prhs[0]);
  95                         in_pi   = mxGetPr(prhs[0]);
  96                         m               = mxGetM(prhs[0]);
  97                         n               = mxGetN(prhs[0]);
  98                         thres   = mxGetScalar(prhs[1]);
  99
 100                         /* Count new nonzeros */
 101                         newnnz=0;
 102                         for(i=0; i<m*n; i++){
 103                                 if (sqrt(in_pr[i]*in_pr[i] + in_pi[i]*in_pi[i])>thres) {newnnz++;}
 104                         }
 105
 106                         if (newnnz>0){
 107                                 /* create output */
 108                                 plhs[0]         = mxCreateSparse(m,n,newnnz,mxCOMPLEX);
 109                                 if (plhs[0]==NULL)
 110                                         mexErrMsgTxt("Could not allocate enough memory!\n");
 111                                 out_pr          = mxGetPr(plhs[0]);
 112                                 out_pi          = mxGetPi(plhs[0]);
 113                                 out_ir          = mxGetIr(plhs[0]);
 114                                 out_jc          = mxGetJc(plhs[0]);
 115                                 passed          = 0;
 116                                 out_jc[0]       = 0;
 117
 118                                 for (col=0; col<n; col++){
 119                                         out_jc[col+1] = out_jc[col];
 120                                         for (current_row_index=0; current_row_index<m; current_row_index++){
 121                                                         if (sqrt(in_pr[current_row_index+m*col]*in_pr[current_row_index+m*col] +
 122                                                                          in_pi[current_row_index+m*col]*in_pi[current_row_index+m*col]) > thres){
 123
 124                                                                 out_pr[passed]=in_pr[current_row_index+m*col];
 125                                                                 out_ir[passed]=current_row_index;
 126                                                                 out_jc[col+1] = out_jc[col+1]+1;
 127                                                                 passed++;
 128                                                         }
 129                                         }
 130                                 }
 131                         }
 132                         else{
 133                                 plhs[0] = mxCreateSparse(m,n,0,mxCOMPLEX);
 134                         }
 135                 }
 136         }
 137         else {
 138         /* Check data type of input argument  */
 139         if (mxIsSparse(prhs[0])){
 140
 141                         /* read input */
 142                         in_pr   = mxGetPr(prhs[0]);
 143                         in_ir   = mxGetIr(prhs[0]);
 144                 in_jc   = mxGetJc(prhs[0]);
 145                         nzmax   = mxGetNzmax(prhs[0]);
 146                         n               = mxGetN(prhs[0]);
 147                         m               = mxGetM(prhs[0]);
 148                         thres   = mxGetScalar(prhs[1]);
 149
 150                         /* Count new nonzeros */
 151                         newnnz=0;
 152                         for(i=0; i<nzmax; i++){
 153                                 if ((fabs(in_pr[i]))>thres) {newnnz++;}
 154                         }
 155
 156                         if (newnnz>0){
 157                                 /* create output */
 158                                 plhs[0]         = mxCreateSparse(m,n,newnnz,mxREAL);
 159                                 if (plhs[0]==NULL)
 160                                         mexErrMsgTxt("Could not allocate enough memory!\n");
 161                                 out_pr          = mxGetPr(plhs[0]);
 162                         out_ir          = mxGetIr(plhs[0]);
 163                         out_jc          = mxGetJc(plhs[0]);
 164                                 passed          = 0;
 165                                 out_jc[0]       = 0;
 166                                 for (col=0; col<n; col++){
 167                                         starting_row_index = in_jc[col];
 168                                         stopping_row_index = in_jc[col+1];
 169                                         out_jc[col+1] = out_jc[col];
 170                                         if (starting_row_index == stopping_row_index)
 171                                                 continue;
 172                                         else {
 173                                                 for (current_row_index = starting_row_index;
 174                                                         current_row_index < stopping_row_index;
 175                                                         current_row_index++)  {
 176                                                                 if (fabs(in_pr[current_row_index])>thres){
 177                                                                         out_pr[passed]=in_pr[current_row_index];
 178                                                                         out_ir[passed]=in_ir[current_row_index];
 179                                                                         out_jc[col+1] = out_jc[col+1]+1;
 180                                                                         passed++;
 181                                                                 }
 182                                                 }
 183                                         }
 184                                 }
 185                         }
 186                         else{
 187                         plhs[0] = mxCreateSparse(m,n,0,mxREAL);
 188                         }
 189                 }
 190                 else{ /* for full matrices */
 191                         /* read input */
 192                         in_pr   = mxGetPr(prhs[0]);
 193                         n               = mxGetN(prhs[0]);
 194                         m               = mxGetM(prhs[0]);
 195                         thres   = mxGetScalar(prhs[1]);
 196
 197                         /* Count new nonzeros */
 198                         newnnz=0;
 199                         for(i=0; i<m*n; i++){
 200                                 if ((fabs(in_pr[i]))>thres) {newnnz++;}
 201                         }
 202
 203                         if (newnnz>0){
 204                                 /* create output */
 205                                 plhs[0]         = mxCreateSparse(m,n,newnnz,mxREAL);
 206                                 if (plhs[0]==NULL)
 207                                         mexErrMsgTxt("Could not allocate enough memory!\n");
 208                                 out_pr          = mxGetPr(plhs[0]);
 209                         out_ir          = mxGetIr(plhs[0]);
 210                         out_jc          = mxGetJc(plhs[0]);
 211                                 passed          = 0;
 212                         out_jc[0]       = 0;
 213
 214                                 for (col=0; col<n; col++){
 215                                         out_jc[col+1] = out_jc[col];
 216                                         for (current_row_index=0; current_row_index<m; current_row_index++){
 217                                                         if (fabs(in_pr[current_row_index+m*col])>thres){
 218                                                                 out_pr[passed]=in_pr[current_row_index+m*col];
 219                                                                 out_ir[passed]=current_row_index;
 220                                                                 out_jc[col+1] = out_jc[col+1]+1;
 221                                                                 passed++;
 222                                                         }
 223                                         }
 224                                 }
 225                         }
 226                         else{
 227                         plhs[0] = mxCreateSparse(m,n,0,mxREAL);
 228                         }
 229                 }
 230         }
 231 }
 232