ZPclustering/zero_diag.cpp at master · maalaria/ZPclustering · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
/**********
 *
 *
 *  mex interface to set to zero the main diagonal of a matrix A
 *
 *  To mexify:
 *             mex  zero_diag.cpp;
 *
 *
 *  [A0] = zero_diag(A);
 *
 *  Input:
 *    A = input matrix
 *
 *  Output:
 *    A0 = same as A with main diagonal set to 0
 *
 *
 *  Lihi Zelnik (Caltech) May.2006
 *
 *
 ************/

#include "mex.h"
#include "matrix.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>


/* ////////////// macros to address input  */
#define INMATRIX  prhs[0]

#define DEBUG 0    /* set to 1 and mex to see print outs */


/*   //////////////// main   */

void
mexFunction (
int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
    /* Make sure at most two output arguments are expected */
    if (nlhs < 1) {
        mexErrMsgTxt("Too few output arguments.");
        return;
    }
    if (nlhs > 1) {
        mexErrMsgTxt("Too many output arguments.");
        return;
    }
    /* Make sure input number is sufficient */
    if (nrhs < 1) {
        mexErrMsgTxt("Too few input arguments.");
        return;
    }
    if (nrhs > 1) {
        mexErrMsgTxt("Too many input arguments.");
        return;
    }

    if( DEBUG )
        mexPrintf("Just starting\n");


    /* check matrix dimensions */
    const int *idims =  mxGetDimensions(INMATRIX);
    if(idims[0] != idims[1])
        mexErrMsgTxt("Input matrix must be square");
    const int node_num = idims[0];
    int i,j;

    if(mxIsSparse(INMATRIX)) {   /* Sparse matrix handling */

        if( DEBUG )
            mexPrintf("Detected the matrix is sparse with %d nodes\n",node_num);

        /* Pointers to non zero rows, cols and values */
        int *in_rows = mxGetIr(INMATRIX);
        int *in_cols = mxGetJc(INMATRIX);
        double *in_values = mxGetPr(INMATRIX);
        int non_zero = mxGetNzmax(INMATRIX);

        if( DEBUG )
            mexPrintf("There are %d non zeros\n",non_zero);

        /* prepare output array */
        double *out_values;
        plhs[0] = mxDuplicateArray(prhs[0]);
        out_values = mxGetPr(plhs[0]);

        int starting_row_index, stopping_row_index, current_index;
        /* Go over the input matrix and set its main diagonal to 0 */
        for (j=0; j<node_num; j++)  {   /* loop over all columns */
            starting_row_index = in_cols[j];
            stopping_row_index = in_cols[j+1];
            if (starting_row_index == stopping_row_index) { /* empty column */
                continue;
            }
            else {  /* loop over all rows for current column */
                for (current_index = starting_row_index;
                current_index < stopping_row_index;
                current_index++)  {
                    if( j == in_rows[current_index] )
                        out_values[current_index] = 0;
                    else
                        out_values[current_index] = in_values[current_index];

                }
            }
        }
        if( DEBUG )
            mexPrintf("All matrix processed safely, current_index=%d\n",current_index);

    }
    /*******************   FULL MATRICES ******************************/
    else { /* Full matrix handling  */

        if( DEBUG )
            mexPrintf("Detected the matrix is full with %d nodes\n",node_num);

        /* Pointers to values */
        double *in_values = mxGetPr(INMATRIX);

        /* prepare output arrays */
        double *out_values;
        plhs[0] = mxDuplicateArray(prhs[0]);
        out_values = mxGetPr(plhs[0]);

        /* prepare value arrays */
        /* Go over the input matrix and set its main diagonal to 0 */
        int total = 0;
        for (j=0; j<node_num; j++)  {   /* loop over all columns */
            for (i=0; i<node_num; i++ ){
                if( i == j )
                    out_values[total] = 0;
                else
                    out_values[total] = in_values[total];
                total++;
            }
        }
        if( DEBUG )
            mexPrintf("All matrix processed safely, total=%d\n",total);
    }

    return;
}