-
Notifications
You must be signed in to change notification settings - Fork 0
/
sparse.c
75 lines (69 loc) · 1.94 KB
/
sparse.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
#ifdef COMPCERT
typedef float _Float16; /* _Float16 is a MacOS thing that CompCert doesn't support */
#endif
#include <stdlib.h>
#include <math.h>
#include "sparse.h"
void *surely_malloc(size_t n) {
void *p = malloc(n);
if (!p) exit(1);
return p;
}
double crs_row_vector_multiply(struct crs_matrix *m, double *v, unsigned i) {
double *val = m->val;
unsigned *col_ind = m->col_ind;
unsigned *row_ptr = m->row_ptr;
unsigned h, hi=row_ptr[i+1];
double s=0.0;
for (h=row_ptr[i]; h<hi; h++) {
double x = val[h];
unsigned j = col_ind[h];
double y = v[j];
s = fma(x,y,s);
}
return s;
}
void crs_matrix_vector_multiply_byrows (struct crs_matrix *m, double *v, double *p) {
unsigned i, rows=crs_matrix_rows(m);
for (i=0; i<rows; i++)
p[i]=crs_row_vector_multiply(m,v,i);
}
/* crs_matrix_vector_multiply(m,v,p)
multiplies a sparse matrix m by a dense vector v,
putting the result into the (already allocated) dense vector p
*/
void crs_matrix_vector_multiply (struct crs_matrix *m, double *v, double *p) {
unsigned i, rows=m->rows;
double *val = m->val;
unsigned *col_ind = m->col_ind;
unsigned *row_ptr = m->row_ptr;
unsigned next=row_ptr[0];
for (i=0; i<rows; i++) {
double s=0.0;
unsigned h=next;
next=row_ptr[i+1];
for (; h<next; h++) {
double x = val[h];
unsigned j = col_ind[h];
double y = v[j];
s = fma(x,y,s);
}
p[i]=s;
}
}
/* Let D be a diagonal matrix, whose diagonal is represented
as the vector diag. Let A be a matrix with number of rows equal
to dimension of D. let m represent A.
Then diag_mult(diag,m) sets m to represent D*A */
void diag_mult(double *diag, struct crs_matrix *m) {
unsigned i, h, rows=m->rows;
for (i=0; i<rows; i++) {
unsigned k=m->row_ptr[i+1];
unsigned x = diag[i];
for (h=m->row_ptr[i]; h<k; h++)
m->val[h] *= x;
}
}
unsigned crs_matrix_rows (struct crs_matrix *m) {
return m->rows;
}