source: CIVL/examples/mpi-omp/AMG2013/seq_mv/csr_matop.c

main
Last change on this file 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: 12.1 KB
Line 
1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14
15/******************************************************************************
16 *
17 * Matrix operation functions for hypre_CSRMatrix class.
18 *
19 *****************************************************************************/
20
21#include "headers.h"
22
23/*--------------------------------------------------------------------------
24 * hypre_CSRMatrixAdd:
25 * adds two CSR Matrices A and B and returns a CSR Matrix C;
26 * Note: The routine does not check for 0-elements which might be generated
27 * through cancellation of elements in A and B or already contained
28 in A and B. To remove those, use hypre_CSRMatrixDeleteZeros
29 *--------------------------------------------------------------------------*/
30
31hypre_CSRMatrix *
32hypre_CSRMatrixAdd( hypre_CSRMatrix *A,
33 hypre_CSRMatrix *B)
34{
35 double *A_data = hypre_CSRMatrixData(A);
36 int *A_i = hypre_CSRMatrixI(A);
37 int *A_j = hypre_CSRMatrixJ(A);
38 int nrows_A = hypre_CSRMatrixNumRows(A);
39 int ncols_A = hypre_CSRMatrixNumCols(A);
40 double *B_data = hypre_CSRMatrixData(B);
41 int *B_i = hypre_CSRMatrixI(B);
42 int *B_j = hypre_CSRMatrixJ(B);
43 int nrows_B = hypre_CSRMatrixNumRows(B);
44 int ncols_B = hypre_CSRMatrixNumCols(B);
45 hypre_CSRMatrix *C;
46 double *C_data;
47 int *C_i;
48 int *C_j;
49
50 int ia, ib, ic, jcol, num_nonzeros;
51 int pos;
52 int *marker;
53
54 if (nrows_A != nrows_B || ncols_A != ncols_B)
55 {
56 printf("Warning! incompatible matrix dimensions!\n");
57 return NULL;
58 }
59
60
61 marker = hypre_CTAlloc(int, ncols_A);
62 C_i = hypre_CTAlloc(int, nrows_A+1);
63
64 for (ia = 0; ia < ncols_A; ia++)
65 marker[ia] = -1;
66
67 num_nonzeros = 0;
68 C_i[0] = 0;
69 for (ic = 0; ic < nrows_A; ic++)
70 {
71 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
72 {
73 jcol = A_j[ia];
74 marker[jcol] = ic;
75 num_nonzeros++;
76 }
77 for (ib = B_i[ic]; ib < B_i[ic+1]; ib++)
78 {
79 jcol = B_j[ib];
80 if (marker[jcol] != ic)
81 {
82 marker[jcol] = ic;
83 num_nonzeros++;
84 }
85 }
86 C_i[ic+1] = num_nonzeros;
87 }
88
89 C = hypre_CSRMatrixCreate(nrows_A, ncols_A, num_nonzeros);
90 hypre_CSRMatrixI(C) = C_i;
91 hypre_CSRMatrixInitialize(C);
92 C_j = hypre_CSRMatrixJ(C);
93 C_data = hypre_CSRMatrixData(C);
94
95 for (ia = 0; ia < ncols_A; ia++)
96 marker[ia] = -1;
97
98 pos = 0;
99 for (ic = 0; ic < nrows_A; ic++)
100 {
101 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
102 {
103 jcol = A_j[ia];
104 C_j[pos] = jcol;
105 C_data[pos] = A_data[ia];
106 marker[jcol] = pos;
107 pos++;
108 }
109 for (ib = B_i[ic]; ib < B_i[ic+1]; ib++)
110 {
111 jcol = B_j[ib];
112 if (marker[jcol] < C_i[ic])
113 {
114 C_j[pos] = jcol;
115 C_data[pos] = B_data[ib];
116 marker[jcol] = pos;
117 pos++;
118 }
119 else
120 {
121 C_data[marker[jcol]] += B_data[ib];
122 }
123 }
124 }
125
126 hypre_TFree(marker);
127 return C;
128}
129
130/*--------------------------------------------------------------------------
131 * hypre_CSRMatrixMultiply
132 * multiplies two CSR Matrices A and B and returns a CSR Matrix C;
133 * Note: The routine does not check for 0-elements which might be generated
134 * through cancellation of elements in A and B or already contained
135 in A and B. To remove those, use hypre_CSRMatrixDeleteZeros
136 *--------------------------------------------------------------------------*/
137
138hypre_CSRMatrix *
139hypre_CSRMatrixMultiply( hypre_CSRMatrix *A,
140 hypre_CSRMatrix *B)
141{
142 double *A_data = hypre_CSRMatrixData(A);
143 int *A_i = hypre_CSRMatrixI(A);
144 int *A_j = hypre_CSRMatrixJ(A);
145 int nrows_A = hypre_CSRMatrixNumRows(A);
146 int ncols_A = hypre_CSRMatrixNumCols(A);
147 double *B_data = hypre_CSRMatrixData(B);
148 int *B_i = hypre_CSRMatrixI(B);
149 int *B_j = hypre_CSRMatrixJ(B);
150 int nrows_B = hypre_CSRMatrixNumRows(B);
151 int ncols_B = hypre_CSRMatrixNumCols(B);
152 hypre_CSRMatrix *C;
153 double *C_data;
154 int *C_i;
155 int *C_j;
156
157 int ia, ib, ic, ja, jb, num_nonzeros=0;
158 int row_start, counter;
159 double a_entry, b_entry;
160 int *B_marker;
161
162 if (ncols_A != nrows_B)
163 {
164 printf("Warning! incompatible matrix dimensions!\n");
165 return NULL;
166 }
167
168
169 B_marker = hypre_CTAlloc(int, ncols_B);
170 C_i = hypre_CTAlloc(int, nrows_A+1);
171
172 for (ib = 0; ib < ncols_B; ib++)
173 B_marker[ib] = -1;
174
175 for (ic = 0; ic < nrows_A; ic++)
176 {
177 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
178 {
179 ja = A_j[ia];
180 for (ib = B_i[ja]; ib < B_i[ja+1]; ib++)
181 {
182 jb = B_j[ib];
183 if (B_marker[jb] != ic)
184 {
185 B_marker[jb] = ic;
186 num_nonzeros++;
187 }
188 }
189 }
190 C_i[ic+1] = num_nonzeros;
191 }
192
193 C = hypre_CSRMatrixCreate(nrows_A, ncols_B, num_nonzeros);
194 hypre_CSRMatrixI(C) = C_i;
195 hypre_CSRMatrixInitialize(C);
196 C_j = hypre_CSRMatrixJ(C);
197 C_data = hypre_CSRMatrixData(C);
198
199 for (ib = 0; ib < ncols_B; ib++)
200 B_marker[ib] = -1;
201
202 counter = 0;
203 for (ic = 0; ic < nrows_A; ic++)
204 {
205 row_start = C_i[ic];
206 for (ia = A_i[ic]; ia < A_i[ic+1]; ia++)
207 {
208 ja = A_j[ia];
209 a_entry = A_data[ia];
210 for (ib = B_i[ja]; ib < B_i[ja+1]; ib++)
211 {
212 jb = B_j[ib];
213 b_entry = B_data[ib];
214 if (B_marker[jb] < row_start)
215 {
216 B_marker[jb] = counter;
217 C_j[B_marker[jb]] = jb;
218 C_data[B_marker[jb]] = a_entry*b_entry;
219 counter++;
220 }
221 else
222 C_data[B_marker[jb]] += a_entry*b_entry;
223
224 }
225 }
226 }
227 hypre_TFree(B_marker);
228 return C;
229}
230
231hypre_CSRMatrix *
232hypre_CSRMatrixDeleteZeros( hypre_CSRMatrix *A, double tol)
233{
234 double *A_data = hypre_CSRMatrixData(A);
235 int *A_i = hypre_CSRMatrixI(A);
236 int *A_j = hypre_CSRMatrixJ(A);
237 int nrows_A = hypre_CSRMatrixNumRows(A);
238 int ncols_A = hypre_CSRMatrixNumCols(A);
239 int num_nonzeros = hypre_CSRMatrixNumNonzeros(A);
240
241 hypre_CSRMatrix *B;
242 double *B_data;
243 int *B_i;
244 int *B_j;
245
246 int zeros;
247 int i, j;
248 int pos_A, pos_B;
249
250 zeros = 0;
251 for (i=0; i < num_nonzeros; i++)
252 if (fabs(A_data[i]) <= tol)
253 zeros++;
254
255 if (zeros)
256 {
257 B = hypre_CSRMatrixCreate(nrows_A,ncols_A,num_nonzeros-zeros);
258 hypre_CSRMatrixInitialize(B);
259 B_i = hypre_CSRMatrixI(B);
260 B_j = hypre_CSRMatrixJ(B);
261 B_data = hypre_CSRMatrixData(B);
262 B_i[0] = 0;
263 pos_A = 0;
264 pos_B = 0;
265 for (i=0; i < nrows_A; i++)
266 {
267 for (j = A_i[i]; j < A_i[i+1]; j++)
268 {
269 if (fabs(A_data[j]) <= tol)
270 {
271 pos_A++;
272 }
273 else
274 {
275 B_data[pos_B] = A_data[pos_A];
276 B_j[pos_B] = A_j[pos_A];
277 pos_B++;
278 pos_A++;
279 }
280 }
281 B_i[i+1] = pos_B;
282 }
283 return B;
284 }
285 else
286 return NULL;
287}
288
289
290/******************************************************************************
291 *
292 * Finds transpose of a hypre_CSRMatrix
293 *
294 *****************************************************************************/
295
296/*--------------------------------------------------------------------------
297 * hypre_CSRMatrixTranspose
298 *--------------------------------------------------------------------------*/
299
300
301int hypre_CSRMatrixTranspose(hypre_CSRMatrix *A, hypre_CSRMatrix **AT,
302 int data)
303
304{
305 double *A_data = hypre_CSRMatrixData(A);
306 int *A_i = hypre_CSRMatrixI(A);
307 int *A_j = hypre_CSRMatrixJ(A);
308 int num_rowsA = hypre_CSRMatrixNumRows(A);
309 int num_colsA = hypre_CSRMatrixNumCols(A);
310 int num_nonzerosA = hypre_CSRMatrixNumNonzeros(A);
311
312 double *AT_data;
313 int *AT_i;
314 int *AT_j;
315 int num_rowsAT;
316 int num_colsAT;
317 int num_nonzerosAT;
318
319 int max_col;
320 int i, j;
321
322 /*--------------------------------------------------------------
323 * First, ascertain that num_cols and num_nonzeros has been set.
324 * If not, set them.
325 *--------------------------------------------------------------*/
326
327 if (! num_nonzerosA)
328 {
329 num_nonzerosA = A_i[num_rowsA];
330 }
331
332 if (num_rowsA && ! num_colsA)
333 {
334 max_col = -1;
335 for (i = 0; i < num_rowsA; ++i)
336 {
337 for (j = A_i[i]; j < A_i[i+1]; j++)
338 {
339 if (A_j[j] > max_col)
340 max_col = A_j[j];
341 }
342 }
343 num_colsA = max_col+1;
344 }
345
346 num_rowsAT = num_colsA;
347 num_colsAT = num_rowsA;
348 num_nonzerosAT = num_nonzerosA;
349
350 *AT = hypre_CSRMatrixCreate(num_rowsAT, num_colsAT, num_nonzerosAT);
351
352 AT_i = hypre_CTAlloc(int, num_rowsAT+1);
353 AT_j = hypre_CTAlloc(int, num_nonzerosAT);
354 hypre_CSRMatrixI(*AT) = AT_i;
355 hypre_CSRMatrixJ(*AT) = AT_j;
356 if (data)
357 {
358 AT_data = hypre_CTAlloc(double, num_nonzerosAT);
359 hypre_CSRMatrixData(*AT) = AT_data;
360 }
361
362 /*-----------------------------------------------------------------
363 * Count the number of entries in each column of A (row of AT)
364 * and fill the AT_i array.
365 *-----------------------------------------------------------------*/
366
367 for (i = 0; i < num_nonzerosA; i++)
368 {
369 ++AT_i[A_j[i]+1];
370 }
371
372 for (i = 2; i <= num_rowsAT; i++)
373 {
374 AT_i[i] += AT_i[i-1];
375 }
376
377 /*----------------------------------------------------------------
378 * Load the data and column numbers of AT
379 *----------------------------------------------------------------*/
380
381 for (i = 0; i < num_rowsA; i++)
382 {
383 for (j = A_i[i]; j < A_i[i+1]; j++)
384 {
385 hypre_assert( AT_i[A_j[j]] >= 0 );
386 hypre_assert( AT_i[A_j[j]] < num_nonzerosAT );
387 AT_j[AT_i[A_j[j]]] = i;
388 if (data) AT_data[AT_i[A_j[j]]] = A_data[j];
389 AT_i[A_j[j]]++;
390 }
391 }
392
393 /*------------------------------------------------------------
394 * AT_i[j] now points to the *end* of the jth row of entries
395 * instead of the beginning. Restore AT_i to front of row.
396 *------------------------------------------------------------*/
397
398 for (i = num_rowsAT; i > 0; i--)
399 {
400 AT_i[i] = AT_i[i-1];
401 }
402
403 AT_i[0] = 0;
404
405 return(0);
406}
407
408
409/*--------------------------------------------------------------------------
410 * hypre_CSRMatrixReorder:
411 * Reorders the column and data arrays of a square CSR matrix, such that the
412 * first entry in each row is the diagonal one.
413 *--------------------------------------------------------------------------*/
414
415int hypre_CSRMatrixReorder(hypre_CSRMatrix *A)
416{
417 int i, j, tempi, row_size;
418 double tempd;
419
420 double *A_data = hypre_CSRMatrixData(A);
421 int *A_i = hypre_CSRMatrixI(A);
422 int *A_j = hypre_CSRMatrixJ(A);
423 int num_rowsA = hypre_CSRMatrixNumRows(A);
424 int num_colsA = hypre_CSRMatrixNumCols(A);
425
426 /* the matrix should be square */
427 if (num_rowsA != num_colsA)
428 return -1;
429
430 for (i = 0; i < num_rowsA; i++)
431 {
432 row_size = A_i[i+1]-A_i[i];
433
434 for (j = 0; j < row_size; j++)
435 {
436 if (A_j[j] == i)
437 {
438 if (j != 0)
439 {
440 tempi = A_j[0];
441 A_j[0] = A_j[j];
442 A_j[j] = tempi;
443
444 tempd = A_data[0];
445 A_data[0] = A_data[j];
446 A_data[j] = tempd;
447 }
448 break;
449 }
450
451 /* diagonal element is missing */
452 if (j == row_size-1)
453 return -2;
454 }
455
456 A_j += row_size;
457 A_data += row_size;
458 }
459
460 return 0;
461}
462
463/*--------------------------------------------------------------------------
464 * hypre_CSRMatrixSumElts:
465 * Returns the sum of all matrix elements.
466 *--------------------------------------------------------------------------*/
467
468double hypre_CSRMatrixSumElts( hypre_CSRMatrix *A )
469{
470 double sum = 0;
471 double * data = hypre_CSRMatrixData( A );
472 int num_nonzeros = hypre_CSRMatrixNumNonzeros(A);
473 int i;
474
475 for ( i=0; i<num_nonzeros; ++i ) sum += data[i];
476
477 return sum;
478}
Note: See TracBrowser for help on using the repository browser.