-
Notifications
You must be signed in to change notification settings - Fork 6
/
smvp.c
141 lines (120 loc) · 3.88 KB
/
smvp.c
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
#include "mex.h"
#include "matrix.h"
/* compile with -largeArrayDims flag */
/* to use "mwSize", need matrix.h, but doesn't work for R2006a */
/* So, use the following definitions instead: */
#ifndef mwSize
#define mwSize size_t
#endif
#ifndef mwIndex
#define mwIndex size_t /* should make it compatible w/ 64-bit systems */
#endif
void printUsage() {
mexPrintf("usage:\tsmvp(A,b), where A is a sparse m x n matrix\n");
mexPrintf("\tand b is a dense n x 1 vector. Handles real or complex data\n");
}
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[]
)
{ /* Sparse matrix-dense vector product
*
* Inputs:
*
* A - m x n sparse matrix
* b - n x 1 dense vector
*
* Outputs:
*
* x = A*b
*
* Darren Engwirda 2006.
*
* Modified by Stephen Becker, April 2009
* Notes on compiling:
* compile with
* mex -O -v -largeArrayDims smvp.c
*
* Modified by Stephen Becker, May 2009, to allow complex multiplies
*/
double *s, *b, *x, rhs;
double *si, *bi, *xi, rhsi; /* the imaginary parts, if data is complex */
/*int *ir, *jc, i, ncol, k; */
mwIndex *ir, *jc, i, k, stop;
mwSize ncol;
int COMPLEX=0, B_REAL=1;
/* Check I/O number */
if (nlhs>1) {
printUsage();
mexErrMsgTxt("SMVP: Too many outputs");
}
if (nrhs!=2) {
printUsage();
mexErrMsgTxt("SMVP: Incorrect number of inputs");
}
/* Brief error checking */
ncol = mxGetN(prhs[0]);
if ((ncol!=mxGetM(prhs[1])) || (mxGetN(prhs[1])!=1)) {
printUsage();
mexErrMsgTxt("Wrong input dimensions");
}
if (!mxIsSparse(prhs[0])) {
printUsage();
mexErrMsgTxt("Matrix must be sparse");
}
/* handle complex data -- for now, both A and b must be the same type */
if ( mxIsComplex( prhs[0] ) ){
COMPLEX = 1; B_REAL =0;
if (!mxIsComplex( prhs[1] )) {
/*
printUsage();
mexErrMsgTxt("If A is complex, b must be also");
**/
/* this situation arise a lot, so let's handle it */
B_REAL = 1;
}
} else if (mxIsComplex( prhs[1] )) {
printUsage();
mexErrMsgTxt("If b is complex, A must be also");
}
/* Allocate output */
if (COMPLEX)
plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[0]), 1, mxCOMPLEX);
else
plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[0]), 1, mxREAL);
/* I/O pointers */
ir = mxGetIr(prhs[0]); /* Row indexing */
jc = mxGetJc(prhs[0]); /* Column count */
s = mxGetPr(prhs[0]); /* Non-zero elements */
b = mxGetPr(prhs[1]); /* Rhs vector */
x = mxGetPr(plhs[0]); /* Output vector */
if (COMPLEX) {
si = mxGetPi(prhs[0]);
xi = mxGetPi(plhs[0]);
if (!B_REAL) {
bi = mxGetPi( prhs[1] );
}
}
/* Multiplication */
if (COMPLEX) {
for (i=0; i<ncol; i++) { /* Loop through columns */
stop = jc[i+1];
rhs = b[i];
if ( B_REAL )
rhsi = 0.0;
else
rhsi = bi[i];
for (k=jc[i]; k<stop; k++) { /* Loop through non-zeros in ith column */
x[ir[k]] += s[k] * rhs - si[k] * rhsi;
xi[ir[k]]+= si[k]* rhs + s[k] * rhsi;
}
}
} else {
for (i=0; i<ncol; i++) { /* Loop through columns */
stop = jc[i+1];
rhs = b[i];
for (k=jc[i]; k<stop; k++) { /* Loop through non-zeros in ith column */
x[ir[k]] += s[k] * rhs;
}
}
}
}