summaryrefslogtreecommitdiffstats
path: root/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/common/toolbox/MultiNcut/sparsifyc.c')
-rwxr-xr-xSD-VBS/common/toolbox/MultiNcut/sparsifyc.c232
1 files changed, 232 insertions, 0 deletions
diff --git a/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c
new file mode 100755
index 0000000..82fca98
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/sparsifyc.c
@@ -0,0 +1,232 @@
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
20void 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 int *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;