source: CIVL/examples/amg/csr/seq/csr_add_impl.c@ 1aaefd4

main test-branch
Last change on this file since 1aaefd4 was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 6.8 KB
Line 
1#include <stdlib.h>
2
3/* = = = = = = = = TASS I/O = = = = = = = = */
4#pragma TASS input
5int N;
6#pragma TASS input
7int M;
8
9#pragma TASS input
10double A0[N*M];
11#pragma TASS input
12double B0[N*M];
13#pragma TASS output
14double OUTS[N*M];
15
16/* = = = = = = = = Dense and CSR Matrix Def = = = = = = = = */
17struct CSRM_struct {
18 double *data;
19 int *i;
20 int *j;
21 int num_rows;
22 int num_cols;
23 int num_nonzeros;
24};
25
26typedef struct CSRM_struct hypre_CSRMatrix;
27
28/* = = = = = = = = Hypre_CSR = = = = = = = = */
29hypre_CSRMatrix *
30hypre_CSRMatrixCreate(
31 int num_rows,
32 int num_cols,
33 int num_nonzeros )
34{
35 hypre_CSRMatrix *matrix;
36
37 matrix = (hypre_CSRMatrix *) malloc (sizeof(hypre_CSRMatrix));
38 matrix->data = (double *) NULL;
39 matrix->i = (int *) NULL;
40 matrix->j = (int *) NULL;
41 matrix->num_rows = num_rows;
42 matrix->num_cols = num_cols;
43 matrix->num_nonzeros = num_nonzeros;
44
45 return matrix;
46}
47
48
49int
50hypre_CSRMatrixInitialize(
51 hypre_CSRMatrix *matrix )
52{
53 int num_rows;
54 int num_nonzeros;
55 int ierr;
56
57 ierr = 0;
58 num_rows = matrix->num_rows;
59 num_nonzeros = matrix->num_nonzeros;
60
61 if (matrix->data == (double *)NULL && num_nonzeros != 0 )
62 matrix->data = (double *)malloc(num_nonzeros*sizeof(double));
63 if (matrix->i == (int *)NULL)
64 matrix->i = (int *)malloc((num_rows + 1)*sizeof(int));
65 if (matrix->j == (int *)NULL)
66 matrix->j = (int *)malloc(num_nonzeros*sizeof(int));
67 return ierr;
68}
69
70void free_CSR( hypre_CSRMatrix *matrix ) {
71 if (matrix->data != (double *) NULL)
72 free(matrix->data);
73 if (matrix->i != (int *) NULL)
74 free(matrix->i);
75 if (matrix->j != (int *) NULL)
76 free(matrix->j);
77 free(matrix);
78}
79
80hypre_CSRMatrix *
81hypre_CSRMatrixAdd( hypre_CSRMatrix *A, hypre_CSRMatrix *B) {
82 double * A_data;
83 int * A_i;
84 int * A_j;
85 int nrows_A;
86 int ncols_A;
87 double * B_data;
88 int * B_i;
89 int * B_j;
90 int nrows_B;
91 int ncols_B;
92 hypre_CSRMatrix *C;
93 double *C_data;
94 int *C_i;
95 int *C_j;
96 int ia;
97 int ib;
98 int ic;
99 int jcol;
100 int num_nonzeros;
101 int pos;
102 int * marker;
103
104 A_data = A->data;
105 A_i = A->i;
106 A_j = A->j;
107 nrows_A = A->num_rows;
108 ncols_A = A->num_cols;
109 B_data = B->data;
110 B_i = B->i;
111 B_j = B->j;
112 nrows_B = B->num_rows;
113 ncols_B = B->num_cols;
114
115 marker = (int *) malloc (ncols_A * sizeof(int));
116 C_i = (int *) malloc ((nrows_A+1) * sizeof(int));
117 for (ia = 0; ia < ncols_A; ia++)
118 marker[ia] = -1;
119 num_nonzeros = 0;
120 C_i[0] = 0;
121 for (ic = 0; ic < nrows_A; ic++) {
122 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++) {
123 jcol = A_j[ia];
124 marker[jcol] = ic;
125 num_nonzeros++;
126 }
127 for (ib = B_i[ic]; ib < B_i[ic+1]; ib++) {
128 jcol = B_j[ib];
129 if (marker[jcol] != ic) {
130 marker[jcol] = ic;
131 num_nonzeros++;
132 }
133 }
134 C_i[ic+1] = num_nonzeros;
135 }
136
137 C = hypre_CSRMatrixCreate(nrows_A, ncols_A, num_nonzeros);
138 C->i = C_i;
139 hypre_CSRMatrixInitialize(C);
140 C_j = C->j;
141 C_data = C->data;
142
143 for (ia = 0; ia < ncols_A; ia++)
144 marker[ia] = -1;
145 pos = 0;
146
147#pragma TASS invariant LoopNNZ pos >= A_i[ic] && pos >= B_i[ic] && pos <= A_i[ic] + B_i[ic];
148#pragma TASS joint invariant LoopCondEquiv ic == spec.i && nrows_A == spec.nr && ncols_A == spec.nc;
149#pragma TASS joint invariant LoopMatAEquiv (forall {int xa | A_i[ic]<=xa && xa<A_i[ic+1]} (spec.j != A_j[xa] || spec.A_data[ic*ncols_A + spec.j] == A_data[xa]));
150#pragma TASS joint invariant LoopMatBEquiv (forall {int xb | B_i[ic]<=xb && xb<B_i[ic+1]} (spec.j != B_j[xb] || spec.B_data[ic*ncols_A + spec.j] == B_data[xb]));
151#pragma TASS joint invariant LoopMatCEquiv (forall {int xc | C_i[ic]<=xc && xc<C_i[ic+1]} (spec.j != C_j[xc] || spec.C_data[ic*ncols_A + spec.j] == C_data[xc]));
152#pragma TASS joint invariant LoopCorrect (forall {int r | 0 <= r && r < spec.nr}(forall {int c | 0 <= c && c < spec.nc} (((forall {int nz1 | C_i[r] <= nz1 && nz1 < C_i[r+1]} (c != C_j[nz1] || spec.C_data[r*spec.nc + c] != C_data[nz1])) != true) || (forall {int nz2 | C_i[r] <= nz2 && nz2 < C_i[r+1]} (c != C_j[nz2] || spec.C_data[r*spec.nc + c] != 0.0)))));
153 for (ic = 0; ic < nrows_A; ic++) {
154 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++) {
155 jcol = A_j[ia];
156 C_j[pos] = jcol;
157 C_data[pos] = A_data[ia];
158 marker[jcol] = pos;
159 pos++;
160 }
161 for (ib = B_i[ic]; ib < B_i[ic+1]; ib++) {
162 jcol = B_j[ib];
163 if (marker[jcol] < C_i[ic]) {
164 C_j[pos] = jcol;
165 C_data[pos] = B_data[ib];
166 marker[jcol] = pos;
167 pos++;
168 }else{
169 C_data[marker[jcol]] += B_data[ib];
170 }
171 }
172 }
173 free(marker);
174 return C;
175}
176
177double *
178expand(hypre_CSRMatrix * mat)
179{
180 int i;
181 int j;
182 int k;
183 int nr = mat->num_rows;
184 int nc = mat->num_cols;
185 double * rtn;
186
187 rtn = (double *) malloc ((mat->num_rows) * (mat->num_cols) * sizeof(double));
188 for (i = 0; i < nr;i++)
189 for (j = 0; j < nc; j++)
190 rtn[i*nc + j] = 0.0;
191 k = 0;
192 for (i = 0; i < nr; i++)
193 while(k < mat->i[i+1])
194 for (j = 0; (k < mat->num_nonzeros) && (j < nc); j++)
195 if (j == mat->j[k]){
196 rtn[i*nc + j] = mat->data[k];
197 k++;
198 j = nc; /* break */
199 }
200
201 return rtn;
202}
203
204int main() {
205 int i;
206 int j;
207 int k;
208 int xnz;
209 int ynz;
210 double tmp;
211 double * sum;
212 hypre_CSRMatrix * X;
213 hypre_CSRMatrix * Y;
214 hypre_CSRMatrix * Z;
215
216 xnz = 0;
217 for (i=0; i < N*M; i++)
218 if (A0[i] != 0)
219 xnz = xnz + 1;
220 ynz = 0;
221 for (i=0; i < N*M; i++)
222 if (B0[i] != 0)
223 ynz = ynz + 1;
224 X = hypre_CSRMatrixCreate(N,M, xnz);
225 Y = hypre_CSRMatrixCreate(N,M, ynz);
226 hypre_CSRMatrixInitialize(X);
227 hypre_CSRMatrixInitialize(Y);
228 k = 0;
229 X->i[0] = k;
230 for (i = 0; i < N; i++){
231 for (j = 0; j < M; j++){
232 tmp = A0[i*M + j];
233 if (tmp != 0){
234 X->data[k] = tmp;
235 X->j[k] = j;
236 k++;
237 }
238 }
239 X->i[i+1] = k;
240 }
241 k = 0;
242 Y->i[0] = k;
243 for (i = 0; i < N; i++){
244 for (j = 0; j < M; j++){
245 tmp = B0[i*M + j];
246 if (tmp != 0){
247 Y->data[k] = tmp;
248 Y->j[k] = j;
249 k++;
250 }
251 }
252 Y->i[i+1] = k;
253 }
254 Z = hypre_CSRMatrixAdd(X,Y);
255 sum = expand(Z);
256 for (i=0; i<N*M; i++)
257 OUTS[i] = sum[i];
258 free_CSR(X);
259 free_CSR(Y);
260 free_CSR(Z);
261 free(sum);
262 return 0;
263}
Note: See TracBrowser for help on using the repository browser.