source: CIVL/mods/dev.civl.com/examples/amg/csr/csr_add_tass.c@ cb4d4f4

main test-branch
Last change on this file since cb4d4f4 was aad342c, checked in by Stephen Siegel <siegel@…>, 3 years ago

Performing huge refactor to incorporate ABC, GMC, and SARL into CIVL repo and use Java modules.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5664 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 3.1 KB
Line 
1
2// Specification: dense matrix
3
4
5
6// from csr_matrix.h:
7typedef struct {
8 double *data;
9 int *i;
10 int *j;
11 int num_rows;
12 int num_cols;
13 int num_nonzeros;
14} hypre_CSRMatrix;
15
16
17hypre_CSRMatrix *
18hypre_CSRMatrixCreate( int num_rows,
19 int num_cols,
20 int num_nonzeros )
21{
22 hypre_CSRMatrix *matrix =
23 (hypre_CSRMatrix *) malloc(sizeof(hypre_CSRMatrix));
24 matrix->data = NULL;
25 matrix->i = NULL;
26 matrix->j = NULL;
27 matrix->num_rows = num_rows;
28 matrix->num_cols = num_cols;
29 matrix->num_nonzeros = num_nonzeros;
30 return matrix;
31}
32
33
34int
35hypre_CSRMatrixInitialize( hypre_CSRMatrix *matrix )
36{
37 int num_rows = matrix->num_rows;
38 int num_nonzeros = matrix->num_nonzeros;
39
40 matrix->data = (double*)malloc(num_nonzeros*sizeof(double));
41 matrix->i = (int*)malloc((num_rows + 1)*sizeof(int));
42 matrix->j = (int*)malloc(num_nonzeros*sizeof(int));
43 return 0;
44}
45
46
47void free_CSR( hypre_CSRMatrix *matrix ) {
48 free(matrix->data);
49 free(matrix->i);
50 free(matrix->j);
51 free(matrix);
52}
53
54hypre_CSRMatrix *
55hypre_CSRMatrixAdd( hypre_CSRMatrix *A, hypre_CSRMatrix *B) {
56 double *A_data = A->data;
57 int *A_i = A->i;
58 int *A_j = A->j;
59 int nrows_A = A->num_rows;
60 int ncols_A = A->num_cols;
61 double *B_data = N->data;
62 int *B_i = B->i;
63 int *B_j = B->j;
64 int nrows_B = B->num_rows;
65 int ncols_B = B->num_cols;
66 hypre_CSRMatrix *C;
67 double *C_data;
68 int *C_i;
69 int *C_j;
70 int ia, ib, ic, jcol, num_nonzeros;
71 int pos;
72 int *marker;
73
74 marker = (int*)malloc(ncols_A*sizeof(int));
75 C_i = (int*)malloc((nrows_A+1)*sizeof(int));
76 for (ia = 0; ia < ncols_A; ia++)
77 marker[ia] = -1;
78 num_nonzeros = 0;
79 C_i[0] = 0;
80 for (ic = 0; ic < nrows_A; ic++) {
81 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++) {
82 jcol = A_j[ia];
83 marker[jcol] = ic;
84 num_nonzeros++;
85 }
86 for (ib = B_i[ic]; ib < B_i[ic+1]; ib++) {
87 jcol = B_j[ib];
88 if (marker[jcol] != ic) {
89 marker[jcol] = ic;
90 num_nonzeros++;
91 }
92 }
93 C_i[ic+1] = num_nonzeros;
94 }
95
96 C = hypre_CSRMatrixCreate(nrows_A, ncols_A, num_nonzeros);
97 C->i = C_i;
98 hypre_CSRMatrixInitialize(C);
99 C_j = C->j;
100 C_data = C->data;
101
102 for (ia = 0; ia < ncols_A; ia++)
103 marker[ia] = -1;
104 pos = 0;
105
106 /*
107#pragma TASS joint invariant ic==spec.i;
108#pragma TASS joint invariant \
109 forall {int i | 0<=i<ic} \
110 forall {int j | 0<=j<ncols_A} \
111 spec.C[i][j] = (exists k.(C_i[i]<=k && k<C_i[i+1] && j=C_j[k])) ? C_data[k] : 0.;
112 */
113
114 for (ic = 0; ic < nrows_A; ic++) {
115 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++) {
116 jcol = A_j[ia];
117 C_j[pos] = jcol;
118 C_data[pos] = A_data[ia];
119 marker[jcol] = pos;
120 pos++;
121 }
122 for (ib = B_i[ic]; ib < B_i[ic+1]; ib++) {
123 jcol = B_j[ib];
124 if (marker[jcol] < C_i[ic]) {
125 C_j[pos] = jcol;
126 C_data[pos] = B_data[ib];
127 marker[jcol] = pos;
128 pos++;
129 }
130 else {
131 C_data[marker[jcol]] += B_data[ib];
132 }
133 }
134 }
135 free(marker);
136 return C;
137}
138
139
Note: See TracBrowser for help on using the repository browser.