1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 /******************************************************************************
9 *
10 * Member functions for hypre_CSRMatrix class.
11 *
12 *****************************************************************************/
13
14 #include "seq_mv.h"
15
16 #ifdef HYPRE_PROFILE
17 HYPRE_Real hypre_profile_times[HYPRE_TIMER_ID_COUNT] = { 0 };
18 #endif
19
20 /*--------------------------------------------------------------------------
21 * hypre_CSRMatrixCreate
22 *--------------------------------------------------------------------------*/
23
24 hypre_CSRMatrix *
hypre_CSRMatrixCreate(HYPRE_Int num_rows,HYPRE_Int num_cols,HYPRE_Int num_nonzeros)25 hypre_CSRMatrixCreate( HYPRE_Int num_rows,
26 HYPRE_Int num_cols,
27 HYPRE_Int num_nonzeros )
28 {
29 hypre_CSRMatrix *matrix;
30
31 matrix = hypre_CTAlloc(hypre_CSRMatrix, 1, HYPRE_MEMORY_HOST);
32
33 hypre_CSRMatrixData(matrix) = NULL;
34 hypre_CSRMatrixI(matrix) = NULL;
35 hypre_CSRMatrixJ(matrix) = NULL;
36 hypre_CSRMatrixBigJ(matrix) = NULL;
37 hypre_CSRMatrixRownnz(matrix) = NULL;
38 hypre_CSRMatrixNumRows(matrix) = num_rows;
39 hypre_CSRMatrixNumRownnz(matrix) = num_rows;
40 hypre_CSRMatrixNumCols(matrix) = num_cols;
41 hypre_CSRMatrixNumNonzeros(matrix) = num_nonzeros;
42 hypre_CSRMatrixMemoryLocation(matrix) = hypre_HandleMemoryLocation(hypre_handle());
43
44 /* set defaults */
45 hypre_CSRMatrixOwnsData(matrix) = 1;
46
47 #if defined(HYPRE_USING_CUSPARSE) || defined(HYPRE_USING_ROCSPARSE)
48 hypre_CSRMatrixSortedJ(matrix) = NULL;
49 hypre_CSRMatrixSortedData(matrix) = NULL;
50 hypre_CSRMatrixCsrsvData(matrix) = NULL;
51 #endif
52
53 return matrix;
54 }
55
56 /*--------------------------------------------------------------------------
57 * hypre_CSRMatrixDestroy
58 *--------------------------------------------------------------------------*/
59
60 HYPRE_Int
hypre_CSRMatrixDestroy(hypre_CSRMatrix * matrix)61 hypre_CSRMatrixDestroy( hypre_CSRMatrix *matrix )
62 {
63 HYPRE_Int ierr = 0;
64
65 if (matrix)
66 {
67 HYPRE_MemoryLocation memory_location = hypre_CSRMatrixMemoryLocation(matrix);
68
69 hypre_TFree(hypre_CSRMatrixI(matrix), memory_location);
70 hypre_TFree(hypre_CSRMatrixRownnz(matrix), HYPRE_MEMORY_HOST);
71
72 if ( hypre_CSRMatrixOwnsData(matrix) )
73 {
74 hypre_TFree(hypre_CSRMatrixData(matrix), memory_location);
75 hypre_TFree(hypre_CSRMatrixJ(matrix), memory_location);
76 /* RL: TODO There might be cases BigJ cannot be freed FIXME
77 * Not so clear how to do it */
78 hypre_TFree(hypre_CSRMatrixBigJ(matrix), memory_location);
79 }
80
81 #if defined(HYPRE_USING_CUSPARSE) || defined(HYPRE_USING_ROCSPARSE)
82 hypre_TFree(hypre_CSRMatrixSortedData(matrix), memory_location);
83 hypre_TFree(hypre_CSRMatrixSortedJ(matrix), memory_location);
84 hypre_CsrsvDataDestroy(hypre_CSRMatrixCsrsvData(matrix));
85 hypre_GpuMatDataDestroy(hypre_CSRMatrixGPUMatData(matrix));
86 #endif
87
88 hypre_TFree(matrix, HYPRE_MEMORY_HOST);
89 }
90
91 return ierr;
92 }
93
94 /*--------------------------------------------------------------------------
95 * hypre_CSRMatrixInitialize
96 *--------------------------------------------------------------------------*/
97
98 HYPRE_Int
hypre_CSRMatrixInitialize_v2(hypre_CSRMatrix * matrix,HYPRE_Int bigInit,HYPRE_MemoryLocation memory_location)99 hypre_CSRMatrixInitialize_v2( hypre_CSRMatrix *matrix, HYPRE_Int bigInit, HYPRE_MemoryLocation memory_location )
100 {
101 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(matrix);
102 HYPRE_Int num_nonzeros = hypre_CSRMatrixNumNonzeros(matrix);
103 /* HYPRE_Int num_rownnz = hypre_CSRMatrixNumRownnz(matrix); */
104
105 HYPRE_Int ierr = 0;
106
107 hypre_CSRMatrixMemoryLocation(matrix) = memory_location;
108
109 /* Caveat: for pre-existing i, j, data, their memory location must be guaranteed to be consistent with `memory_location'
110 * Otherwise, mismatches will exist and problems will be encountered when being used, and freed */
111
112 if ( !hypre_CSRMatrixData(matrix) && num_nonzeros )
113 {
114 hypre_CSRMatrixData(matrix) = hypre_CTAlloc(HYPRE_Complex, num_nonzeros, memory_location);
115 }
116 /*
117 else
118 {
119 //if (PointerAttributes(hypre_CSRMatrixData(matrix))==HYPRE_HOST_POINTER) printf("MATREIX INITIAL WITH JHOST DATA\n");
120 }
121 */
122
123 if ( !hypre_CSRMatrixI(matrix) )
124 {
125 hypre_CSRMatrixI(matrix) = hypre_CTAlloc(HYPRE_Int, num_rows + 1, memory_location);
126 }
127
128 /*
129 if (!hypre_CSRMatrixRownnz(matrix))
130 {
131 hypre_CSRMatrixRownnz(matrix) = hypre_CTAlloc(HYPRE_Int, num_rownnz, memory_location);
132 }
133 */
134
135 if (bigInit)
136 {
137 if ( !hypre_CSRMatrixBigJ(matrix) && num_nonzeros )
138 {
139 hypre_CSRMatrixBigJ(matrix) = hypre_CTAlloc(HYPRE_BigInt, num_nonzeros, memory_location);
140 }
141 }
142 else
143 {
144 if ( !hypre_CSRMatrixJ(matrix) && num_nonzeros )
145 {
146 hypre_CSRMatrixJ(matrix) = hypre_CTAlloc(HYPRE_Int, num_nonzeros, memory_location);
147 }
148 }
149
150 return ierr;
151 }
152
153 HYPRE_Int
hypre_CSRMatrixInitialize(hypre_CSRMatrix * matrix)154 hypre_CSRMatrixInitialize( hypre_CSRMatrix *matrix )
155 {
156 HYPRE_Int ierr;
157
158 ierr = hypre_CSRMatrixInitialize_v2( matrix, 0, hypre_CSRMatrixMemoryLocation(matrix) );
159
160 return ierr;
161 }
162
163 /*--------------------------------------------------------------------------
164 * hypre_CSRMatrixResize
165 *--------------------------------------------------------------------------*/
166
167 HYPRE_Int
hypre_CSRMatrixResize(hypre_CSRMatrix * matrix,HYPRE_Int new_num_rows,HYPRE_Int new_num_cols,HYPRE_Int new_num_nonzeros)168 hypre_CSRMatrixResize( hypre_CSRMatrix *matrix, HYPRE_Int new_num_rows, HYPRE_Int new_num_cols, HYPRE_Int new_num_nonzeros )
169 {
170 HYPRE_MemoryLocation memory_location = hypre_CSRMatrixMemoryLocation(matrix);
171 HYPRE_Int old_num_nonzeros = hypre_CSRMatrixNumNonzeros(matrix);
172 HYPRE_Int old_num_rows = hypre_CSRMatrixNumRows(matrix);
173
174 if (!hypre_CSRMatrixOwnsData(matrix))
175 {
176 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"Error: called hypre_CSRMatrixResize on a matrix that doesn't own the data\n");
177 return 1;
178 }
179
180 hypre_CSRMatrixNumCols(matrix) = new_num_cols;
181
182 if (new_num_nonzeros != hypre_CSRMatrixNumNonzeros(matrix))
183 {
184 hypre_CSRMatrixNumNonzeros(matrix) = new_num_nonzeros;
185
186 if (!hypre_CSRMatrixData(matrix))
187 {
188 hypre_CSRMatrixData(matrix) = hypre_CTAlloc(HYPRE_Complex, new_num_nonzeros, memory_location);
189 }
190 else
191 {
192 hypre_CSRMatrixData(matrix) = hypre_TReAlloc_v2(hypre_CSRMatrixData(matrix), HYPRE_Complex, old_num_nonzeros, HYPRE_Complex, new_num_nonzeros, memory_location);
193 }
194
195 if (!hypre_CSRMatrixJ(matrix))
196 {
197 hypre_CSRMatrixJ(matrix) = hypre_CTAlloc(HYPRE_Int, new_num_nonzeros, memory_location);
198 }
199 else
200 {
201 hypre_CSRMatrixJ(matrix) = hypre_TReAlloc_v2(hypre_CSRMatrixJ(matrix), HYPRE_Int, old_num_nonzeros, HYPRE_Int, new_num_nonzeros, memory_location);
202 }
203 }
204
205 if (new_num_rows != hypre_CSRMatrixNumRows(matrix))
206 {
207 hypre_CSRMatrixNumRows(matrix) = new_num_rows;
208
209 if (!hypre_CSRMatrixI(matrix))
210 {
211 hypre_CSRMatrixI(matrix) = hypre_CTAlloc(HYPRE_Int, new_num_rows + 1, memory_location);
212 }
213 else
214 {
215 hypre_CSRMatrixI(matrix) = hypre_TReAlloc_v2(hypre_CSRMatrixI(matrix), HYPRE_Int, old_num_rows + 1, HYPRE_Int, new_num_rows + 1, memory_location);
216 }
217 }
218
219 return 0;
220 }
221
222 /*--------------------------------------------------------------------------
223 * hypre_CSRMatrixBigInitialize
224 *--------------------------------------------------------------------------*/
225
226 HYPRE_Int
hypre_CSRMatrixBigInitialize(hypre_CSRMatrix * matrix)227 hypre_CSRMatrixBigInitialize( hypre_CSRMatrix *matrix )
228 {
229 HYPRE_Int ierr;
230
231 ierr = hypre_CSRMatrixInitialize_v2( matrix, 1, hypre_CSRMatrixMemoryLocation(matrix) );
232
233 return ierr;
234 }
235
236 /*--------------------------------------------------------------------------
237 * hypre_CSRMatrixBigJtoJ
238 * RL: TODO GPU impl.
239 *--------------------------------------------------------------------------*/
240
241 HYPRE_Int
hypre_CSRMatrixBigJtoJ(hypre_CSRMatrix * matrix)242 hypre_CSRMatrixBigJtoJ( hypre_CSRMatrix *matrix )
243 {
244 HYPRE_Int num_nonzeros = hypre_CSRMatrixNumNonzeros(matrix);
245 HYPRE_BigInt *matrix_big_j = hypre_CSRMatrixBigJ(matrix);
246 HYPRE_Int *matrix_j = NULL;
247
248 if (num_nonzeros && matrix_big_j)
249 {
250 #if defined(HYPRE_MIXEDINT) || defined(HYPRE_BIGINT)
251 HYPRE_Int i;
252 matrix_j = hypre_TAlloc(HYPRE_Int, num_nonzeros, hypre_CSRMatrixMemoryLocation(matrix));
253 for (i = 0; i < num_nonzeros; i++)
254 {
255 matrix_j[i] = (HYPRE_Int) matrix_big_j[i];
256 }
257 hypre_TFree(matrix_big_j, hypre_CSRMatrixMemoryLocation(matrix));
258 #else
259 hypre_assert(sizeof(HYPRE_Int) == sizeof(HYPRE_BigInt));
260 matrix_j = (HYPRE_Int *) matrix_big_j;
261 #endif
262 hypre_CSRMatrixJ(matrix) = matrix_j;
263 hypre_CSRMatrixBigJ(matrix) = NULL;
264 }
265
266 return hypre_error_flag;
267 }
268
269 /*--------------------------------------------------------------------------
270 * hypre_CSRMatrixJtoBigJ
271 * RL: TODO GPU impl.
272 *--------------------------------------------------------------------------*/
273
274 HYPRE_Int
hypre_CSRMatrixJtoBigJ(hypre_CSRMatrix * matrix)275 hypre_CSRMatrixJtoBigJ( hypre_CSRMatrix *matrix )
276 {
277 HYPRE_Int num_nonzeros = hypre_CSRMatrixNumNonzeros(matrix);
278 HYPRE_Int *matrix_j = hypre_CSRMatrixJ(matrix);
279 HYPRE_BigInt *matrix_big_j = NULL;
280
281 if (num_nonzeros && matrix_j)
282 {
283 #if defined(HYPRE_MIXEDINT) || defined(HYPRE_BIGINT)
284 HYPRE_Int i;
285 matrix_big_j = hypre_TAlloc(HYPRE_BigInt, num_nonzeros, hypre_CSRMatrixMemoryLocation(matrix));
286 for (i = 0; i < num_nonzeros; i++)
287 {
288 matrix_big_j[i] = (HYPRE_BigInt) matrix_j[i];
289 }
290 hypre_TFree(matrix_j, hypre_CSRMatrixMemoryLocation(matrix));
291 #else
292 hypre_assert(sizeof(HYPRE_Int) == sizeof(HYPRE_BigInt));
293 matrix_big_j = (HYPRE_BigInt *) matrix_j;
294 #endif
295 hypre_CSRMatrixBigJ(matrix) = matrix_big_j;
296 hypre_CSRMatrixJ(matrix) = NULL;
297 }
298
299 return hypre_error_flag;
300 }
301
302 /*--------------------------------------------------------------------------
303 * hypre_CSRMatrixSetDataOwner
304 *--------------------------------------------------------------------------*/
305
306 HYPRE_Int
hypre_CSRMatrixSetDataOwner(hypre_CSRMatrix * matrix,HYPRE_Int owns_data)307 hypre_CSRMatrixSetDataOwner( hypre_CSRMatrix *matrix,
308 HYPRE_Int owns_data )
309 {
310 HYPRE_Int ierr = 0;
311
312 hypre_CSRMatrixOwnsData(matrix) = owns_data;
313
314 return ierr;
315 }
316
317 /*--------------------------------------------------------------------------
318 * hypre_CSRMatrixSetRownnz
319 *
320 * function to set the substructure rownnz and num_rowsnnz inside the CSRMatrix
321 * it needs the A_i substructure of CSRMatrix to find the nonzero rows.
322 * It runs after the create CSR and when A_i is known..It does not check for
323 * the existence of A_i or of the CSR matrix.
324 *--------------------------------------------------------------------------*/
325
326 HYPRE_Int
hypre_CSRMatrixSetRownnzHost(hypre_CSRMatrix * matrix)327 hypre_CSRMatrixSetRownnzHost( hypre_CSRMatrix *matrix )
328 {
329 HYPRE_Int ierr = 0;
330 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(matrix);
331 HYPRE_Int *A_i = hypre_CSRMatrixI(matrix);
332 HYPRE_Int *Arownnz;
333
334 HYPRE_Int i, adiag;
335 HYPRE_Int irownnz = 0;
336
337 for (i = 0; i < num_rows; i++)
338 {
339 adiag = A_i[i+1] - A_i[i];
340 if (adiag > 0)
341 {
342 irownnz++;
343 }
344 }
345
346 hypre_CSRMatrixNumRownnz(matrix) = irownnz;
347
348 if (irownnz == 0 || irownnz == num_rows)
349 {
350 hypre_CSRMatrixRownnz(matrix) = NULL;
351 }
352 else
353 {
354 Arownnz = hypre_CTAlloc(HYPRE_Int, irownnz, HYPRE_MEMORY_HOST);
355 irownnz = 0;
356 for (i = 0; i < num_rows; i++)
357 {
358 adiag = A_i[i+1] - A_i[i];
359 if (adiag > 0)
360 {
361 Arownnz[irownnz++] = i;
362 }
363 }
364 hypre_CSRMatrixRownnz(matrix) = Arownnz;
365 }
366
367 return ierr;
368 }
369
370 HYPRE_Int
hypre_CSRMatrixSetRownnz(hypre_CSRMatrix * matrix)371 hypre_CSRMatrixSetRownnz( hypre_CSRMatrix *matrix )
372 {
373 HYPRE_Int ierr = 0;
374
375 #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP)
376 HYPRE_ExecutionPolicy exec = hypre_GetExecPolicy1( hypre_CSRMatrixMemoryLocation(matrix) );
377
378 if (exec == HYPRE_EXEC_DEVICE)
379 {
380 // TODO RL: there's no need currently for having rownnz on GPUs
381 }
382 else
383 #endif
384 {
385 ierr = hypre_CSRMatrixSetRownnzHost(matrix);
386 }
387
388 return ierr;
389 }
390
391 /* check if numnonzeros was properly set to be ia[nrow] */
392 HYPRE_Int
hypre_CSRMatrixCheckSetNumNonzeros(hypre_CSRMatrix * matrix)393 hypre_CSRMatrixCheckSetNumNonzeros( hypre_CSRMatrix *matrix )
394 {
395 if (!matrix)
396 {
397 return 0;
398 }
399
400 HYPRE_Int nnz, ierr = 0;
401
402 hypre_TMemcpy(&nnz, hypre_CSRMatrixI(matrix) + hypre_CSRMatrixNumRows(matrix),
403 HYPRE_Int, 1, HYPRE_MEMORY_HOST, hypre_CSRMatrixMemoryLocation(matrix));
404
405 if (hypre_CSRMatrixNumNonzeros(matrix) != nnz)
406 {
407 ierr = 1;
408 hypre_printf("warning: CSR matrix nnz was not set properly (!= ia[nrow], %d %d)\n", hypre_CSRMatrixNumNonzeros(matrix), nnz );
409 hypre_assert(0);
410 hypre_CSRMatrixNumNonzeros(matrix) = nnz;
411 }
412
413 return ierr;
414 }
415
416 /*--------------------------------------------------------------------------
417 * hypre_CSRMatrixRead
418 *--------------------------------------------------------------------------*/
419
420 hypre_CSRMatrix *
hypre_CSRMatrixRead(char * file_name)421 hypre_CSRMatrixRead( char *file_name )
422 {
423 hypre_CSRMatrix *matrix;
424
425 FILE *fp;
426
427 HYPRE_Complex *matrix_data;
428 HYPRE_Int *matrix_i;
429 HYPRE_Int *matrix_j;
430 HYPRE_Int num_rows;
431 HYPRE_Int num_nonzeros;
432 HYPRE_Int max_col = 0;
433
434 HYPRE_Int file_base = 1;
435
436 HYPRE_Int j;
437
438 /*----------------------------------------------------------
439 * Read in the data
440 *----------------------------------------------------------*/
441 fp = fopen(file_name, "r");
442
443 hypre_fscanf(fp, "%d", &num_rows);
444
445 matrix_i = hypre_CTAlloc(HYPRE_Int, num_rows + 1, HYPRE_MEMORY_HOST);
446 for (j = 0; j < num_rows+1; j++)
447 {
448 hypre_fscanf(fp, "%d", &matrix_i[j]);
449 matrix_i[j] -= file_base;
450 }
451
452 num_nonzeros = matrix_i[num_rows];
453
454 matrix = hypre_CSRMatrixCreate(num_rows, num_rows, matrix_i[num_rows]);
455 hypre_CSRMatrixI(matrix) = matrix_i;
456 hypre_CSRMatrixInitialize_v2(matrix, 0, HYPRE_MEMORY_HOST);
457 matrix_j = hypre_CSRMatrixJ(matrix);
458
459 for (j = 0; j < num_nonzeros; j++)
460 {
461 hypre_fscanf(fp, "%d", &matrix_j[j]);
462 matrix_j[j] -= file_base;
463
464 if (matrix_j[j] > max_col)
465 {
466 max_col = matrix_j[j];
467 }
468 }
469
470 matrix_data = hypre_CSRMatrixData(matrix);
471 for (j = 0; j < matrix_i[num_rows]; j++)
472 {
473 hypre_fscanf(fp, "%le", &matrix_data[j]);
474 }
475
476 fclose(fp);
477
478 hypre_CSRMatrixNumNonzeros(matrix) = num_nonzeros;
479 hypre_CSRMatrixNumCols(matrix) = ++max_col;
480
481 return matrix;
482 }
483
484 /*--------------------------------------------------------------------------
485 * hypre_CSRMatrixPrint
486 *--------------------------------------------------------------------------*/
487
488 HYPRE_Int
hypre_CSRMatrixPrint(hypre_CSRMatrix * matrix,const char * file_name)489 hypre_CSRMatrixPrint( hypre_CSRMatrix *matrix,
490 const char *file_name )
491 {
492 FILE *fp;
493
494 HYPRE_Complex *matrix_data;
495 HYPRE_Int *matrix_i;
496 HYPRE_Int *matrix_j;
497 HYPRE_Int num_rows;
498
499 HYPRE_Int file_base = 1;
500
501 HYPRE_Int j;
502
503 HYPRE_Int ierr = 0;
504
505 /*----------------------------------------------------------
506 * Print the matrix data
507 *----------------------------------------------------------*/
508
509 matrix_data = hypre_CSRMatrixData(matrix);
510 matrix_i = hypre_CSRMatrixI(matrix);
511 matrix_j = hypre_CSRMatrixJ(matrix);
512 num_rows = hypre_CSRMatrixNumRows(matrix);
513
514 fp = fopen(file_name, "w");
515
516 hypre_fprintf(fp, "%d\n", num_rows);
517
518 for (j = 0; j <= num_rows; j++)
519 {
520 hypre_fprintf(fp, "%d\n", matrix_i[j] + file_base);
521 }
522
523 for (j = 0; j < matrix_i[num_rows]; j++)
524 {
525 hypre_fprintf(fp, "%d\n", matrix_j[j] + file_base);
526 }
527
528 if (matrix_data)
529 {
530 for (j = 0; j < matrix_i[num_rows]; j++)
531 {
532 #ifdef HYPRE_COMPLEX
533 hypre_fprintf(fp, "%.14e , %.14e\n",
534 hypre_creal(matrix_data[j]), hypre_cimag(matrix_data[j]));
535 #else
536 hypre_fprintf(fp, "%.14e\n", matrix_data[j]);
537 #endif
538 }
539 }
540 else
541 {
542 hypre_fprintf(fp, "Warning: No matrix data!\n");
543 }
544
545 fclose(fp);
546
547 return ierr;
548 }
549
550 HYPRE_Int
hypre_CSRMatrixPrintMM(hypre_CSRMatrix * matrix,HYPRE_Int basei,HYPRE_Int basej,HYPRE_Int trans,const char * file_name)551 hypre_CSRMatrixPrintMM( hypre_CSRMatrix *matrix,
552 HYPRE_Int basei,
553 HYPRE_Int basej,
554 HYPRE_Int trans,
555 const char *file_name )
556 {
557 FILE *fp;
558
559 HYPRE_Complex *matrix_data;
560 HYPRE_Int *matrix_i;
561 HYPRE_Int *matrix_j;
562 HYPRE_Int num_rows, num_cols;
563
564 /* HYPRE_Int file_base = 1; */
565
566 HYPRE_Int j,k;
567
568 HYPRE_Int ierr = 0;
569
570 /*----------------------------------------------------------
571 * Print the matrix data
572 *----------------------------------------------------------*/
573
574 matrix_data = hypre_CSRMatrixData(matrix);
575 matrix_i = hypre_CSRMatrixI(matrix);
576 matrix_j = hypre_CSRMatrixJ(matrix);
577 num_rows = hypre_CSRMatrixNumRows(matrix);
578 num_cols = hypre_CSRMatrixNumCols(matrix);
579
580 if (file_name)
581 {
582 fp = fopen(file_name, "w");
583 }
584 else
585 {
586 fp = stdout;
587 }
588
589 hypre_fprintf(fp, "%%%%MatrixMarket matrix coordinate real general\n");
590
591 hypre_assert(matrix_i[num_rows] == hypre_CSRMatrixNumNonzeros(matrix));
592
593 if (!trans)
594 {
595 hypre_fprintf(fp, "%d %d %d\n", num_rows, num_cols, hypre_CSRMatrixNumNonzeros(matrix));
596 }
597 else
598 {
599 hypre_fprintf(fp, "%d %d %d\n", num_cols, num_rows, hypre_CSRMatrixNumNonzeros(matrix));
600 }
601
602 for (j = 0; j < num_rows; j++)
603 {
604 for (k = matrix_i[j]; k < matrix_i[j+1]; k++)
605 {
606 if (!trans)
607 {
608 hypre_fprintf(fp, "%d %d %.15e\n", j + basei, matrix_j[k] + basej, matrix_data[k]);
609 }
610 else
611 {
612 hypre_fprintf(fp, "%d %d %.15e\n", matrix_j[k] + basej, j + basei, matrix_data[k]);
613 }
614 }
615 }
616
617 if (file_name)
618 {
619 fclose(fp);
620 }
621
622 return ierr;
623 }
624
625 HYPRE_Int
hypre_CSRMatrixPrint2(hypre_CSRMatrix * matrix,const char * file_name)626 hypre_CSRMatrixPrint2( hypre_CSRMatrix *matrix,
627 const char *file_name )
628 {
629 return hypre_CSRMatrixPrintMM(matrix, 0, 0, 0, file_name);
630 }
631
632 /*--------------------------------------------------------------------------
633 * hypre_CSRMatrixPrintHB: print a CSRMatrix in Harwell-Boeing format
634 *--------------------------------------------------------------------------*/
635
636 HYPRE_Int
hypre_CSRMatrixPrintHB(hypre_CSRMatrix * matrix_input,char * file_name)637 hypre_CSRMatrixPrintHB( hypre_CSRMatrix *matrix_input,
638 char *file_name )
639 {
640 FILE *fp;
641 hypre_CSRMatrix *matrix;
642 HYPRE_Complex *matrix_data;
643 HYPRE_Int *matrix_i;
644 HYPRE_Int *matrix_j;
645 HYPRE_Int num_rows;
646 HYPRE_Int file_base = 1;
647 HYPRE_Int j, totcrd, ptrcrd, indcrd, valcrd, rhscrd;
648 HYPRE_Int ierr = 0;
649
650 /*----------------------------------------------------------
651 * Print the matrix data
652 *----------------------------------------------------------*/
653
654 /* First transpose the input matrix, since HB is in CSC format */
655 hypre_CSRMatrixTranspose(matrix_input, &matrix, 1);
656
657 matrix_data = hypre_CSRMatrixData(matrix);
658 matrix_i = hypre_CSRMatrixI(matrix);
659 matrix_j = hypre_CSRMatrixJ(matrix);
660 num_rows = hypre_CSRMatrixNumRows(matrix);
661
662 fp = fopen(file_name, "w");
663
664 hypre_fprintf(fp, "%-70s Key \n", "Title");
665 ptrcrd = num_rows;
666 indcrd = matrix_i[num_rows];
667 valcrd = matrix_i[num_rows];
668 rhscrd = 0;
669 totcrd = ptrcrd + indcrd + valcrd + rhscrd;
670 hypre_fprintf (fp, "%14d%14d%14d%14d%14d\n",
671 totcrd, ptrcrd, indcrd, valcrd, rhscrd);
672 hypre_fprintf (fp, "%-14s%14i%14i%14i%14i\n", "RUA",
673 num_rows, num_rows, valcrd, 0);
674 hypre_fprintf (fp, "%-16s%-16s%-16s%26s\n", "(1I8)", "(1I8)", "(1E16.8)", "");
675
676 for (j = 0; j <= num_rows; j++)
677 {
678 hypre_fprintf(fp, "%8d\n", matrix_i[j] + file_base);
679 }
680
681 for (j = 0; j < matrix_i[num_rows]; j++)
682 {
683 hypre_fprintf(fp, "%8d\n", matrix_j[j] + file_base);
684 }
685
686 if (matrix_data)
687 {
688 for (j = 0; j < matrix_i[num_rows]; j++)
689 {
690 #ifdef HYPRE_COMPLEX
691 hypre_fprintf(fp, "%16.8e , %16.8e\n",
692 hypre_creal(matrix_data[j]), hypre_cimag(matrix_data[j]));
693 #else
694 hypre_fprintf(fp, "%16.8e\n", matrix_data[j]);
695 #endif
696 }
697 }
698 else
699 {
700 hypre_fprintf(fp, "Warning: No matrix data!\n");
701 }
702
703 fclose(fp);
704
705 hypre_CSRMatrixDestroy(matrix);
706
707 return ierr;
708 }
709
710 /*--------------------------------------------------------------------------
711 * hypre_CSRMatrixCopy: copy A to B,
712 * if copy_data = 0 only the structure of A is copied to B.
713 * the routine does not check if the dimensions/sparsity of A and B match !!!
714 *--------------------------------------------------------------------------*/
715
716 HYPRE_Int
hypre_CSRMatrixCopy(hypre_CSRMatrix * A,hypre_CSRMatrix * B,HYPRE_Int copy_data)717 hypre_CSRMatrixCopy( hypre_CSRMatrix *A, hypre_CSRMatrix *B, HYPRE_Int copy_data )
718 {
719 HYPRE_Int ierr = 0;
720 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
721 HYPRE_Int num_nonzeros = hypre_CSRMatrixNumNonzeros(A);
722
723 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
724 HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
725 HYPRE_BigInt *A_bigj = hypre_CSRMatrixBigJ(A);
726 HYPRE_Complex *A_data;
727
728 HYPRE_Int *B_i = hypre_CSRMatrixI(B);
729 HYPRE_Int *B_j = hypre_CSRMatrixJ(B);
730 HYPRE_BigInt *B_bigj = hypre_CSRMatrixBigJ(B);
731 HYPRE_Complex *B_data;
732
733 HYPRE_MemoryLocation memory_location_A = hypre_CSRMatrixMemoryLocation(A);
734 HYPRE_MemoryLocation memory_location_B = hypre_CSRMatrixMemoryLocation(B);
735
736 hypre_TMemcpy(B_i, A_i, HYPRE_Int, num_rows + 1, memory_location_B, memory_location_A);
737
738 if (A_j && B_j)
739 {
740 hypre_TMemcpy(B_j, A_j, HYPRE_Int, num_nonzeros, memory_location_B, memory_location_A);
741 }
742
743 if (A_bigj && B_bigj)
744 {
745 hypre_TMemcpy(B_bigj, A_bigj, HYPRE_BigInt, num_nonzeros, memory_location_B, memory_location_A);
746 }
747
748 if (copy_data)
749 {
750 A_data = hypre_CSRMatrixData(A);
751 B_data = hypre_CSRMatrixData(B);
752 hypre_TMemcpy(B_data, A_data, HYPRE_Complex, num_nonzeros, memory_location_B, memory_location_A);
753 }
754
755 return ierr;
756 }
757
758 /*--------------------------------------------------------------------------
759 * hypre_CSRMatrixClone
760 * Creates and returns a new copy of the argument, A.
761 * Copying is a deep copy in that no pointers are copied; new arrays are
762 * created where necessary.
763 *--------------------------------------------------------------------------*/
764
765 hypre_CSRMatrix*
hypre_CSRMatrixClone_v2(hypre_CSRMatrix * A,HYPRE_Int copy_data,HYPRE_MemoryLocation memory_location)766 hypre_CSRMatrixClone_v2( hypre_CSRMatrix *A, HYPRE_Int copy_data, HYPRE_MemoryLocation memory_location )
767 {
768 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A);
769 HYPRE_Int num_cols = hypre_CSRMatrixNumCols(A);
770 HYPRE_Int num_nonzeros = hypre_CSRMatrixNumNonzeros( A );
771
772 hypre_CSRMatrix *B = hypre_CSRMatrixCreate( num_rows, num_cols, num_nonzeros );
773
774 HYPRE_Int bigInit = hypre_CSRMatrixBigJ(A) != NULL;
775
776 hypre_CSRMatrixInitialize_v2(B, bigInit, memory_location);
777
778 hypre_CSRMatrixCopy(A, B, copy_data);
779
780 return B;
781 }
782
783 hypre_CSRMatrix*
hypre_CSRMatrixClone(hypre_CSRMatrix * A,HYPRE_Int copy_data)784 hypre_CSRMatrixClone( hypre_CSRMatrix *A, HYPRE_Int copy_data )
785 {
786 return hypre_CSRMatrixClone_v2(A, copy_data, hypre_CSRMatrixMemoryLocation(A));
787 }
788
789 /*--------------------------------------------------------------------------
790 * hypre_CSRMatrixUnion
791 * Creates and returns a matrix whose elements are the union of those of A and B.
792 * Data is not computed, only structural information is created.
793 * A and B must have the same numbers of rows.
794 * Nothing is done about Rownnz.
795 *
796 * If col_map_offd_A and col_map_offd_B are zero, A and B are expected to have
797 * the same column indexing. Otherwise, col_map_offd_A, col_map_offd_B should
798 * be the arrays of that name from two ParCSRMatrices of which A and B are the
799 * offd blocks.
800 *
801 * The algorithm can be expected to have reasonable efficiency only for very
802 * sparse matrices (many rows, few nonzeros per row).
803 * The nonzeros of a computed row are NOT necessarily in any particular order.
804 *--------------------------------------------------------------------------*/
805
806 hypre_CSRMatrix*
hypre_CSRMatrixUnion(hypre_CSRMatrix * A,hypre_CSRMatrix * B,HYPRE_BigInt * col_map_offd_A,HYPRE_BigInt * col_map_offd_B,HYPRE_BigInt ** col_map_offd_C)807 hypre_CSRMatrixUnion( hypre_CSRMatrix *A,
808 hypre_CSRMatrix *B,
809 HYPRE_BigInt *col_map_offd_A,
810 HYPRE_BigInt *col_map_offd_B,
811 HYPRE_BigInt **col_map_offd_C )
812 {
813 HYPRE_Int num_rows = hypre_CSRMatrixNumRows( A );
814 HYPRE_Int num_cols_A = hypre_CSRMatrixNumCols( A );
815 HYPRE_Int num_cols_B = hypre_CSRMatrixNumCols( B );
816 HYPRE_Int num_cols;
817 HYPRE_Int num_nonzeros;
818 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
819 HYPRE_Int *A_j = hypre_CSRMatrixJ(A);
820 HYPRE_Int *B_i = hypre_CSRMatrixI(B);
821 HYPRE_Int *B_j = hypre_CSRMatrixJ(B);
822 HYPRE_Int *C_i;
823 HYPRE_Int *C_j;
824 HYPRE_Int *jC = NULL;
825 HYPRE_BigInt jBg, big_jA, big_jB;
826 HYPRE_Int i, jA, jB;
827 HYPRE_Int ma, mb, mc, ma_min, ma_max, match;
828 hypre_CSRMatrix* C;
829
830 HYPRE_MemoryLocation memory_location = hypre_CSRMatrixMemoryLocation(A);
831
832 hypre_assert( num_rows == hypre_CSRMatrixNumRows(B) );
833
834 if ( col_map_offd_B )
835 {
836 hypre_assert( col_map_offd_A );
837 }
838
839 if ( col_map_offd_A )
840 {
841 hypre_assert( col_map_offd_B );
842 }
843
844 /* ==== First, go through the columns of A and B to count the columns of C. */
845 if ( col_map_offd_A==0 )
846 { /* The matrices are diagonal blocks.
847 Normally num_cols_A==num_cols_B, col_starts is the same, etc.
848 */
849 num_cols = hypre_max( num_cols_A, num_cols_B );
850 }
851 else
852 { /* The matrices are offdiagonal blocks. */
853 jC = hypre_CTAlloc(HYPRE_Int, num_cols_B, HYPRE_MEMORY_HOST);
854 num_cols = num_cols_A; /* initialization; we'll compute the actual value */
855 for ( jB=0; jB<num_cols_B; ++jB )
856 {
857 match = 0;
858 jBg = col_map_offd_B[jB];
859 for ( ma=0; ma<num_cols_A; ++ma )
860 {
861 if ( col_map_offd_A[ma]==jBg )
862 {
863 match = 1;
864 }
865 }
866 if ( match==0 )
867 {
868 jC[jB] = num_cols;
869 ++num_cols;
870 }
871 }
872 }
873
874 /* ==== If we're working on a ParCSRMatrix's offd block,
875 make and load col_map_offd_C */
876 if ( col_map_offd_A )
877 {
878 *col_map_offd_C = hypre_CTAlloc( HYPRE_BigInt, num_cols, HYPRE_MEMORY_HOST);
879 for ( jA=0; jA<num_cols_A; ++jA )
880 {
881 (*col_map_offd_C)[jA] = col_map_offd_A[jA];
882 }
883 for ( jB=0; jB<num_cols_B; ++jB )
884 {
885 match = 0;
886 jBg = col_map_offd_B[jB];
887 for ( ma=0; ma<num_cols_A; ++ma )
888 {
889 if ( col_map_offd_A[ma]==jBg )
890 {
891 match = 1;
892 }
893 }
894 if ( match==0 )
895 {
896 (*col_map_offd_C)[ jC[jB] ] = jBg;
897 }
898 }
899 }
900
901
902 /* ==== The first run through A and B is to count the number of nonzero elements,
903 without HYPRE_Complex-counting duplicates. Then we can create C. */
904 num_nonzeros = hypre_CSRMatrixNumNonzeros(A);
905 for ( i=0; i<num_rows; ++i )
906 {
907 ma_min = A_i[i]; ma_max = A_i[i+1];
908 for ( mb=B_i[i]; mb<B_i[i+1]; ++mb )
909 {
910 jB = B_j[mb];
911 if ( col_map_offd_B )
912 {
913 big_jB = col_map_offd_B[jB];
914 }
915 match = 0;
916 for ( ma=ma_min; ma<ma_max; ++ma )
917 {
918 jA = A_j[ma];
919 if ( col_map_offd_A )
920 {
921 big_jA = col_map_offd_A[jA];
922 }
923 if ( big_jB == big_jA )
924 {
925 match = 1;
926 if( ma==ma_min )
927 {
928 ++ma_min;
929 }
930 break;
931 }
932 }
933 if ( match==0 )
934 {
935 ++num_nonzeros;
936 }
937 }
938 }
939
940 C = hypre_CSRMatrixCreate( num_rows, num_cols, num_nonzeros );
941 hypre_CSRMatrixInitialize_v2( C, 0, memory_location );
942
943 /* ==== The second run through A and B is to pick out the column numbers
944 for each row, and put them in C. */
945 C_i = hypre_CSRMatrixI(C);
946 C_i[0] = 0;
947 C_j = hypre_CSRMatrixJ(C);
948 mc = 0;
949 for ( i=0; i<num_rows; ++i )
950 {
951 ma_min = A_i[i];
952 ma_max = A_i[i+1];
953 for ( ma=ma_min; ma<ma_max; ++ma )
954 {
955 C_j[mc] = A_j[ma];
956 ++mc;
957 }
958 for ( mb=B_i[i]; mb<B_i[i+1]; ++mb )
959 {
960 jB = B_j[mb];
961 if ( col_map_offd_B )
962 {
963 big_jB = col_map_offd_B[jB];
964 }
965 match = 0;
966 for ( ma=ma_min; ma<ma_max; ++ma )
967 {
968 jA = A_j[ma];
969 if ( col_map_offd_A )
970 {
971 big_jA = col_map_offd_A[jA];
972 }
973 if ( big_jB == big_jA )
974 {
975 match = 1;
976 if( ma==ma_min )
977 {
978 ++ma_min;
979 }
980 break;
981 }
982 }
983 if ( match==0 )
984 {
985 if ( col_map_offd_A )
986 {
987 C_j[mc] = jC[ B_j[mb] ];
988 }
989 else
990 {
991 C_j[mc] = B_j[mb];
992 }
993 /* ... I don't know whether column indices are required to be in any
994 particular order. If so, we'll need to sort. */
995 ++mc;
996 }
997 }
998 C_i[i+1] = mc;
999 }
1000
1001 hypre_assert( mc == num_nonzeros );
1002
1003 hypre_TFree(jC, HYPRE_MEMORY_HOST);
1004
1005 return C;
1006 }
1007
hypre_CSRMatrixGetLoadBalancedPartitionBoundary(hypre_CSRMatrix * A,HYPRE_Int idx)1008 static HYPRE_Int hypre_CSRMatrixGetLoadBalancedPartitionBoundary(hypre_CSRMatrix *A, HYPRE_Int idx)
1009 {
1010 HYPRE_Int num_nonzerosA = hypre_CSRMatrixNumNonzeros(A);
1011 HYPRE_Int num_rowsA = hypre_CSRMatrixNumRows(A);
1012 HYPRE_Int *A_i = hypre_CSRMatrixI(A);
1013
1014 HYPRE_Int num_threads = hypre_NumActiveThreads();
1015
1016 HYPRE_Int nonzeros_per_thread = (num_nonzerosA + num_threads - 1)/num_threads;
1017
1018 if (idx <= 0)
1019 {
1020 return 0;
1021 }
1022 else if (idx >= num_threads)
1023 {
1024 return num_rowsA;
1025 }
1026 else
1027 {
1028 return (HYPRE_Int)(hypre_LowerBound(A_i, A_i + num_rowsA, nonzeros_per_thread*idx) - A_i);
1029 }
1030 }
1031
hypre_CSRMatrixGetLoadBalancedPartitionBegin(hypre_CSRMatrix * A)1032 HYPRE_Int hypre_CSRMatrixGetLoadBalancedPartitionBegin(hypre_CSRMatrix *A)
1033 {
1034 return hypre_CSRMatrixGetLoadBalancedPartitionBoundary(A, hypre_GetThreadNum());
1035 }
1036
hypre_CSRMatrixGetLoadBalancedPartitionEnd(hypre_CSRMatrix * A)1037 HYPRE_Int hypre_CSRMatrixGetLoadBalancedPartitionEnd(hypre_CSRMatrix *A)
1038 {
1039 return hypre_CSRMatrixGetLoadBalancedPartitionBoundary(A, hypre_GetThreadNum() + 1);
1040 }
1041
1042 HYPRE_Int
hypre_CSRMatrixPrefetch(hypre_CSRMatrix * A,HYPRE_MemoryLocation memory_location)1043 hypre_CSRMatrixPrefetch( hypre_CSRMatrix *A, HYPRE_MemoryLocation memory_location )
1044 {
1045 HYPRE_Int ierr = 0;
1046
1047 #ifdef HYPRE_USING_UNIFIED_MEMORY
1048 if (hypre_CSRMatrixMemoryLocation(A) != HYPRE_MEMORY_DEVICE)
1049 {
1050 return 1;
1051 }
1052
1053 HYPRE_Complex *data = hypre_CSRMatrixData(A);
1054 HYPRE_Int *ia = hypre_CSRMatrixI(A);
1055 HYPRE_Int *ja = hypre_CSRMatrixJ(A);
1056 HYPRE_Int nrow = hypre_CSRMatrixNumRows(A);
1057 HYPRE_Int nnzA = hypre_CSRMatrixNumNonzeros(A);
1058
1059 hypre_MemPrefetch(data, sizeof(HYPRE_Complex)*nnzA, memory_location);
1060 hypre_MemPrefetch(ia, sizeof(HYPRE_Int)*(nrow+1), memory_location);
1061 hypre_MemPrefetch(ja, sizeof(HYPRE_Int)*nnzA, memory_location);
1062 #endif
1063
1064 return ierr;
1065 }
1066
1067 #if defined(HYPRE_USING_CUSPARSE) || defined(HYPRE_USING_ROCSPARSE)
1068 hypre_GpuMatData *
hypre_CSRMatrixGetGPUMatData(hypre_CSRMatrix * matrix)1069 hypre_CSRMatrixGetGPUMatData(hypre_CSRMatrix *matrix)
1070 {
1071 if (!matrix)
1072 {
1073 return NULL;
1074 }
1075
1076 if (!hypre_CSRMatrixGPUMatData(matrix))
1077 {
1078 hypre_CSRMatrixGPUMatData(matrix) = hypre_GpuMatDataCreate();
1079 }
1080
1081 return hypre_CSRMatrixGPUMatData(matrix);
1082 }
1083 #endif
1084
1085