| 1 | #include <civlc.cvh>
|
|---|
| 2 | #include <pointer.cvh>
|
|---|
| 3 |
|
|---|
| 4 | #include "HYPRE_utilities.h"
|
|---|
| 5 | #include "csr_matrix.h"
|
|---|
| 6 | #include "hypre_memory.h"
|
|---|
| 7 |
|
|---|
| 8 | $input int BOUND;
|
|---|
| 9 |
|
|---|
| 10 | $input int M;
|
|---|
| 11 | $assume(1 <= M);
|
|---|
| 12 | //$assume(1 <= M && M <= BOUND);
|
|---|
| 13 |
|
|---|
| 14 | $input int N;
|
|---|
| 15 | $assume(1 <= N && N <= BOUND);
|
|---|
| 16 | $assume(1 <= N);
|
|---|
| 17 |
|
|---|
| 18 | $input HYPRE_Real A[M][N], B[M][N];
|
|---|
| 19 |
|
|---|
| 20 | /** Original **/
|
|---|
| 21 |
|
|---|
| 22 | /* hypre_CSRMatrix* */
|
|---|
| 23 | /* hypre_CSRMatrixAddHost ( hypre_CSRMatrix *A, */
|
|---|
| 24 | /* hypre_CSRMatrix *B ) */
|
|---|
| 25 | /* { */
|
|---|
| 26 | /* HYPRE_Complex *A_data = hypre_CSRMatrixData(A); */
|
|---|
| 27 | /* HYPRE_Int *A_i = hypre_CSRMatrixI(A); */
|
|---|
| 28 | /* HYPRE_Int *A_j = hypre_CSRMatrixJ(A); */
|
|---|
| 29 | /* HYPRE_Int nrows_A = hypre_CSRMatrixNumRows(A); */
|
|---|
| 30 | /* HYPRE_Int ncols_A = hypre_CSRMatrixNumCols(A); */
|
|---|
| 31 | /* HYPRE_Complex *B_data = hypre_CSRMatrixData(B); */
|
|---|
| 32 | /* HYPRE_Int *B_i = hypre_CSRMatrixI(B); */
|
|---|
| 33 | /* HYPRE_Int *B_j = hypre_CSRMatrixJ(B); */
|
|---|
| 34 | /* HYPRE_Int nrows_B = hypre_CSRMatrixNumRows(B); */
|
|---|
| 35 | /* HYPRE_Int ncols_B = hypre_CSRMatrixNumCols(B); */
|
|---|
| 36 | /* hypre_CSRMatrix *C; */
|
|---|
| 37 | /* HYPRE_Complex *C_data; */
|
|---|
| 38 | /* HYPRE_Int *C_i; */
|
|---|
| 39 | /* HYPRE_Int *C_j; */
|
|---|
| 40 |
|
|---|
| 41 | /* HYPRE_Int ia, ib, ic, jcol, num_nonzeros; */
|
|---|
| 42 | /* HYPRE_Int pos; */
|
|---|
| 43 | /* HYPRE_Int *marker; */
|
|---|
| 44 |
|
|---|
| 45 | /* HYPRE_MemoryLocation memory_location_A = hypre_CSRMatrixMemoryLocation(A); */
|
|---|
| 46 | /* HYPRE_MemoryLocation memory_location_B = hypre_CSRMatrixMemoryLocation(B); */
|
|---|
| 47 |
|
|---|
| 48 | /* /\* RL: TODO cannot guarantee, maybe should never assert */
|
|---|
| 49 | /* hypre_assert(memory_location_A == memory_location_B); */
|
|---|
| 50 | /* *\/ */
|
|---|
| 51 |
|
|---|
| 52 | /* /\* RL: in the case of A=H, B=D, or A=D, B=H, let C = D, */
|
|---|
| 53 | /* * not sure if this is the right thing to do. */
|
|---|
| 54 | /* * Also, need something like this in other places */
|
|---|
| 55 | /* * TODO *\/ */
|
|---|
| 56 | /* HYPRE_MemoryLocation memory_location_C = hypre_max(memory_location_A, memory_location_B); */
|
|---|
| 57 |
|
|---|
| 58 | /* if (nrows_A != nrows_B || ncols_A != ncols_B) */
|
|---|
| 59 | /* { */
|
|---|
| 60 | /* hypre_error_w_msg(HYPRE_ERROR_GENERIC,"Warning! incompatible matrix dimensions!\n"); */
|
|---|
| 61 | /* return NULL; */
|
|---|
| 62 | /* } */
|
|---|
| 63 |
|
|---|
| 64 | /* marker = hypre_CTAlloc(HYPRE_Int, ncols_A, HYPRE_MEMORY_HOST); */
|
|---|
| 65 | /* C_i = hypre_CTAlloc(HYPRE_Int, nrows_A+1, memory_location_C); */
|
|---|
| 66 |
|
|---|
| 67 | /* for (ia = 0; ia < ncols_A; ia++) */
|
|---|
| 68 | /* { */
|
|---|
| 69 | /* marker[ia] = -1; */
|
|---|
| 70 | /* } */
|
|---|
| 71 |
|
|---|
| 72 | /* num_nonzeros = 0; */
|
|---|
| 73 | /* C_i[0] = 0; */
|
|---|
| 74 | /* for (ic = 0; ic < nrows_A; ic++) */
|
|---|
| 75 | /* { */
|
|---|
| 76 | /* for (ia = A_i[ic]; ia < A_i[ic+1]; ia++) */
|
|---|
| 77 | /* { */
|
|---|
| 78 | /* jcol = A_j[ia]; */
|
|---|
| 79 | /* marker[jcol] = ic; */
|
|---|
| 80 | /* num_nonzeros++; */
|
|---|
| 81 | /* } */
|
|---|
| 82 | /* for (ib = B_i[ic]; ib < B_i[ic+1]; ib++) */
|
|---|
| 83 | /* { */
|
|---|
| 84 | /* jcol = B_j[ib]; */
|
|---|
| 85 | /* if (marker[jcol] != ic) */
|
|---|
| 86 | /* { */
|
|---|
| 87 | /* marker[jcol] = ic; */
|
|---|
| 88 | /* num_nonzeros++; */
|
|---|
| 89 | /* } */
|
|---|
| 90 | /* } */
|
|---|
| 91 | /* C_i[ic+1] = num_nonzeros; */
|
|---|
| 92 | /* } */
|
|---|
| 93 |
|
|---|
| 94 | /* C = hypre_CSRMatrixCreate(nrows_A, ncols_A, num_nonzeros); */
|
|---|
| 95 | /* hypre_CSRMatrixI(C) = C_i; */
|
|---|
| 96 | /* hypre_CSRMatrixInitialize_v2(C, 0, memory_location_C); */
|
|---|
| 97 | /* C_j = hypre_CSRMatrixJ(C); */
|
|---|
| 98 | /* C_data = hypre_CSRMatrixData(C); */
|
|---|
| 99 |
|
|---|
| 100 | /* for (ia = 0; ia < ncols_A; ia++) */
|
|---|
| 101 | /* { */
|
|---|
| 102 | /* marker[ia] = -1; */
|
|---|
| 103 | /* } */
|
|---|
| 104 | /* pos = 0; */
|
|---|
| 105 | /* for (ic = 0; ic < nrows_A; ic++) */
|
|---|
| 106 | /* { */
|
|---|
| 107 | /* for (ia = A_i[ic]; ia < A_i[ic+1]; ia++) */
|
|---|
| 108 | /* { */
|
|---|
| 109 | /* jcol = A_j[ia]; */
|
|---|
| 110 | /* C_j[pos] = jcol; */
|
|---|
| 111 | /* C_data[pos] = A_data[ia]; */
|
|---|
| 112 | /* marker[jcol] = pos; */
|
|---|
| 113 | /* pos++; */
|
|---|
| 114 | /* } */
|
|---|
| 115 | /* for (ib = B_i[ic]; ib < B_i[ic+1]; ib++) */
|
|---|
| 116 | /* { */
|
|---|
| 117 | /* jcol = B_j[ib]; */
|
|---|
| 118 | /* if (marker[jcol] < C_i[ic]) */
|
|---|
| 119 | /* { */
|
|---|
| 120 | /* C_j[pos] = jcol; */
|
|---|
| 121 | /* C_data[pos] = B_data[ib]; */
|
|---|
| 122 | /* marker[jcol] = pos; */
|
|---|
| 123 | /* pos++; */
|
|---|
| 124 | /* } */
|
|---|
| 125 | /* else */
|
|---|
| 126 | /* { */
|
|---|
| 127 | /* C_data[marker[jcol]] += B_data[ib]; */
|
|---|
| 128 | /* } */
|
|---|
| 129 | /* } */
|
|---|
| 130 | /* } */
|
|---|
| 131 |
|
|---|
| 132 | /* hypre_TFree(marker, HYPRE_MEMORY_HOST); */
|
|---|
| 133 |
|
|---|
| 134 | /* return C; */
|
|---|
| 135 | /* } */
|
|---|
| 136 |
|
|---|
| 137 | hypre_CSRMatrix*
|
|---|
| 138 | hypre_CSRMatrixAddHost ( hypre_CSRMatrix *A,
|
|---|
| 139 | hypre_CSRMatrix *B )
|
|---|
| 140 | {
|
|---|
| 141 | HYPRE_Complex *A_data = hypre_CSRMatrixData(A);
|
|---|
| 142 | HYPRE_Int *A_i = hypre_CSRMatrixI(A);
|
|---|
| 143 | HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
|
|---|
| 144 | HYPRE_Int nrows_A = hypre_CSRMatrixNumRows(A);
|
|---|
| 145 | HYPRE_Int ncols_A = hypre_CSRMatrixNumCols(A);
|
|---|
| 146 | HYPRE_Complex *B_data = hypre_CSRMatrixData(B);
|
|---|
| 147 | HYPRE_Int *B_i = hypre_CSRMatrixI(B);
|
|---|
| 148 | HYPRE_Int *B_j = hypre_CSRMatrixJ(B);
|
|---|
| 149 | HYPRE_Int nrows_B = hypre_CSRMatrixNumRows(B);
|
|---|
| 150 | HYPRE_Int ncols_B = hypre_CSRMatrixNumCols(B);
|
|---|
| 151 | hypre_CSRMatrix *C;
|
|---|
| 152 | HYPRE_Complex *C_data;
|
|---|
| 153 | HYPRE_Int *C_i;
|
|---|
| 154 | HYPRE_Int *C_j;
|
|---|
| 155 |
|
|---|
| 156 | HYPRE_Int ia, ib, ic, jcol, num_nonzeros;
|
|---|
| 157 | HYPRE_Int pos;
|
|---|
| 158 | HYPRE_Int *marker;
|
|---|
| 159 |
|
|---|
| 160 | hypre_CTAlloc(HYPRE_Int, marker, ncols_A);
|
|---|
| 161 | hypre_CTAlloc(HYPRE_Int, C_i, nrows_A+1);
|
|---|
| 162 |
|
|---|
| 163 | for (ia = 0; ia < ncols_A; ia++)
|
|---|
| 164 | {
|
|---|
| 165 | marker[ia] = -1;
|
|---|
| 166 | }
|
|---|
| 167 |
|
|---|
| 168 | num_nonzeros = 0;
|
|---|
| 169 | C_i[0] = 0;
|
|---|
| 170 | //@ transform flatten tag1;
|
|---|
| 171 | for (ic = 0; ic < nrows_A; ic++)
|
|---|
| 172 | {
|
|---|
| 173 | for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
|
|---|
| 174 | {
|
|---|
| 175 | jcol = A_j[ia];
|
|---|
| 176 | marker[jcol] = ic;
|
|---|
| 177 | num_nonzeros++;
|
|---|
| 178 | }
|
|---|
| 179 | for (ib = B_i[ic]; ib < B_i[ic+1]; ib++)
|
|---|
| 180 | {
|
|---|
| 181 | jcol = B_j[ib];
|
|---|
| 182 | if (marker[jcol] != ic)
|
|---|
| 183 | {
|
|---|
| 184 | marker[jcol] = ic;
|
|---|
| 185 | num_nonzeros++;
|
|---|
| 186 | }
|
|---|
| 187 | }
|
|---|
| 188 | C_i[ic+1] = num_nonzeros;
|
|---|
| 189 | }
|
|---|
| 190 |
|
|---|
| 191 | C = hypre_CSRMatrixCreate(nrows_A, ncols_A, num_nonzeros);
|
|---|
| 192 | hypre_CSRMatrixI(C) = C_i;
|
|---|
| 193 | hypre_CSRMatrixInitialize_v2(C);
|
|---|
| 194 | C_j = hypre_CSRMatrixJ(C);
|
|---|
| 195 | C_data = hypre_CSRMatrixData(C);
|
|---|
| 196 |
|
|---|
| 197 | for (ia = 0; ia < ncols_A; ia++)
|
|---|
| 198 | {
|
|---|
| 199 | marker[ia] = -1;
|
|---|
| 200 | }
|
|---|
| 201 | pos = 0;
|
|---|
| 202 | for (ic = 0; ic < nrows_A; ic++)
|
|---|
| 203 | {
|
|---|
| 204 | for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
|
|---|
| 205 | {
|
|---|
| 206 | jcol = A_j[ia];
|
|---|
| 207 | C_j[pos] = jcol;
|
|---|
| 208 | C_data[pos] = A_data[ia];
|
|---|
| 209 | marker[jcol] = pos;
|
|---|
| 210 | pos++;
|
|---|
| 211 | }
|
|---|
| 212 | for (ib = B_i[ic]; ib < B_i[ic+1]; ib++)
|
|---|
| 213 | {
|
|---|
| 214 | jcol = B_j[ib];
|
|---|
| 215 | if (marker[jcol] < C_i[ic])
|
|---|
| 216 | {
|
|---|
| 217 | C_j[pos] = jcol;
|
|---|
| 218 | C_data[pos] = B_data[ib];
|
|---|
| 219 | marker[jcol] = pos;
|
|---|
| 220 | pos++;
|
|---|
| 221 | }
|
|---|
| 222 | else
|
|---|
| 223 | {
|
|---|
| 224 | C_data[marker[jcol]] += B_data[ib];
|
|---|
| 225 | }
|
|---|
| 226 | }
|
|---|
| 227 | }
|
|---|
| 228 |
|
|---|
| 229 | hypre_TFree(marker);
|
|---|
| 230 |
|
|---|
| 231 | return C;
|
|---|
| 232 | }
|
|---|
| 233 |
|
|---|
| 234 | hypre_CSRMatrix *
|
|---|
| 235 | DenseToCSR (int m, int n, HYPRE_Real a[][])
|
|---|
| 236 | {
|
|---|
| 237 | int numnonzeros = 0;
|
|---|
| 238 | for (int i = 0; i < m; ++i)
|
|---|
| 239 | {
|
|---|
| 240 | for (int j = 0; j < n; ++j)
|
|---|
| 241 | {
|
|---|
| 242 | if (a[i][j] != 0.0)
|
|---|
| 243 | {
|
|---|
| 244 | ++numnonzeros;
|
|---|
| 245 | }
|
|---|
| 246 | }
|
|---|
| 247 | }
|
|---|
| 248 |
|
|---|
| 249 | hypre_CSRMatrix * csr_a = hypre_CSRMatrixCreate(m, n, numnonzeros);
|
|---|
| 250 | hypre_CSRMatrixInitialize_v2(csr_a);
|
|---|
| 251 |
|
|---|
| 252 | int k = 0;
|
|---|
| 253 |
|
|---|
| 254 | for (int i = 0; i < m; ++i)
|
|---|
| 255 | {
|
|---|
| 256 | hypre_CSRMatrixI(csr_a)[i] = k;
|
|---|
| 257 | for (int j = 0; j < n; ++j)
|
|---|
| 258 | {
|
|---|
| 259 | if (a[i][j] != 0.0)
|
|---|
| 260 | {
|
|---|
| 261 | hypre_CSRMatrixData(csr_a)[k] = a[i][j];
|
|---|
| 262 | hypre_CSRMatrixJ(csr_a)[k] = j;
|
|---|
| 263 | ++k;
|
|---|
| 264 | }
|
|---|
| 265 | }
|
|---|
| 266 | }
|
|---|
| 267 | hypre_CSRMatrixI(csr_a)[m] = numnonzeros;
|
|---|
| 268 |
|
|---|
| 269 | return csr_a;
|
|---|
| 270 | }
|
|---|
| 271 |
|
|---|
| 272 | void
|
|---|
| 273 | CSRToDense(HYPRE_Real dense[][], hypre_CSRMatrix * csr)
|
|---|
| 274 | {
|
|---|
| 275 | int m = hypre_CSRMatrixNumRows(csr), n = hypre_CSRMatrixNumCols(csr);
|
|---|
| 276 |
|
|---|
| 277 | for (int i = 0; i < m; ++i)
|
|---|
| 278 | {
|
|---|
| 279 | for (int j = 0; j < n; ++j)
|
|---|
| 280 | {
|
|---|
| 281 | dense[i][j] = 0.0;
|
|---|
| 282 | }
|
|---|
| 283 | }
|
|---|
| 284 |
|
|---|
| 285 | for (int i = 0; i < m; ++i)
|
|---|
| 286 | {
|
|---|
| 287 | int start = hypre_CSRMatrixI(csr)[i], end = hypre_CSRMatrixI(csr)[i+1];
|
|---|
| 288 | for (int k = start; k < end; ++k)
|
|---|
| 289 | {
|
|---|
| 290 | dense[i][hypre_CSRMatrixJ(csr)[k]] = hypre_CSRMatrixData(csr)[k];
|
|---|
| 291 | }
|
|---|
| 292 | }
|
|---|
| 293 | }
|
|---|
| 294 |
|
|---|
| 295 | void
|
|---|
| 296 | AddSpec(int m, int n, HYPRE_Real c[][], HYPRE_Real a[][], HYPRE_Real b[][])
|
|---|
| 297 | {
|
|---|
| 298 | for (int i = 0; i < m; ++i)
|
|---|
| 299 | {
|
|---|
| 300 | for (int j = 0; j < n; ++j)
|
|---|
| 301 | {
|
|---|
| 302 | c[i][j] = a[i][j] + b[i][j];
|
|---|
| 303 | }
|
|---|
| 304 | }
|
|---|
| 305 | }
|
|---|
| 306 |
|
|---|
| 307 | void free_CSR(hypre_CSRMatrix * matrix) {
|
|---|
| 308 | free(matrix->data);
|
|---|
| 309 | free(matrix->i);
|
|---|
| 310 | free(matrix->j);
|
|---|
| 311 | free(matrix);
|
|---|
| 312 | }
|
|---|
| 313 |
|
|---|
| 314 | int main(int argc, char* argv[])
|
|---|
| 315 | {
|
|---|
| 316 | $atomic {
|
|---|
| 317 | $elaborate(M);
|
|---|
| 318 | $elaborate(N);
|
|---|
| 319 | HYPRE_Real C[M][N];
|
|---|
| 320 | AddSpec(M, N, C, A, B);
|
|---|
| 321 |
|
|---|
| 322 | hypre_CSRMatrix * csr_A = DenseToCSR(M, N, A);
|
|---|
| 323 | hypre_CSRMatrix * csr_B = DenseToCSR(M, N, B);
|
|---|
| 324 |
|
|---|
| 325 | HYPRE_Real impl_A[M][N];
|
|---|
| 326 | CSRToDense(impl_A, csr_A);
|
|---|
| 327 | //$assert($equals(&A, &impl_A));
|
|---|
| 328 |
|
|---|
| 329 | hypre_CSRMatrix * csr_C = hypre_CSRMatrixAddHost(csr_A, csr_B);
|
|---|
| 330 |
|
|---|
| 331 | HYPRE_Real impl_C[M][N];
|
|---|
| 332 | CSRToDense(impl_C, csr_C);
|
|---|
| 333 |
|
|---|
| 334 | $assert($equals(&C, &impl_C));
|
|---|
| 335 |
|
|---|
| 336 | free_CSR (csr_A), free_CSR (csr_B), free_CSR (csr_C);
|
|---|
| 337 | }
|
|---|
| 338 | }
|
|---|
