#include <stdio.h>
#include <stdlib.h>
/*
A
1 0 3 0 2
0 1 0 4 6
2 5 3 0 0
0 0 7 8 0
0 1 0 3 9
*/
/*
x
1
1
1
1
1
*/
int const n = 5;
double val[14] = {1,3,2,1,4,6,2,5,3,7,8,1,3,9};
int col[14] = {0,2,4,1,3,4,0,1,2,2,3,1,3,4};
int row[6] = {0, 3, 6, 9, 11, 14}; /* 14: sentinel */
double x[5] = {1,1,1,1,1};
double b[5] = {0,0,0,0,0};
double matrix_a[25] = {1,0,3,0,2,0,1,0,4,6,2,5,3,0,0,0,0,7,8,0,0,1,0,3,9};
double vector_x[5] = {1,1,1,1,1};
double vector_b[5];
void gemv_crs(double *b, double *x, double *val, int *col, int *row, int n);
void gemv(double *vector_y, double *matrix_a, double *vector_x, int n);
int main(void){
int i;
gemv_crs(b, x, val, col, row, n);
gemv(vector_b, matrix_a, vector_x, n);
for(i = 0; i < n; i++)
printf("crs_b[%d] = %f\n",i,b[i]);
putchar('\n');
for(i = 0; i < n; i++)
printf("b[%d] = %f\n",i,vector_b[i]);
return 0;
}
void gemv(double *vector_y, double *matrix_a, double *vector_x, int n) {
int i, j;
double vmx;
for(i = 0; i < n; i++){
vmx = 0.0;
for(j = 0; j < n; j++)
vmx += matrix_a[i * n + j] * vector_x[j];
vector_y[i] = vmx;
}
}
void gemv_crs(double *b,double *x,double *val,int *col,int *row,int n) {
int i, j;
double vmx;
for(i = 0; i < n ; i++) {
vmx = 0.0;
for(j = row[i] ; j < row[i + 1] ; j++)
vmx += val[j] * x[col[j]];
b[i] = vmx;
}
}
/* end */