source: CIVL/examples/mpi-omp/AMG2013/parcsr_mv/par_csr_matrix.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: 71.6 KB
RevLine 
[2aa6644]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 * Member functions for hypre_ParCSRMatrix class.
18 *
19 *****************************************************************************/
20
21#include "headers.h"
22
23#include "../seq_mv/HYPRE_seq_mv.h"
24/* In addition to publically accessible interface in HYPRE_mv.h, the implementation
25 in this file uses accessor macros into the sequential matrix structure, and
26 so includes the .h that defines that structure. Should those accessor functions
27 become proper functions at some later date, this will not be necessary. AJC 4/99
28*/
29#include "../seq_mv/csr_matrix.h"
30
31
32#ifdef HYPRE_NO_GLOBAL_PARTITION
33int hypre_FillResponseParToCSRMatrix(void*, int, int, void*, MPI_Comm, void**, int*);
34#endif
35
36/*--------------------------------------------------------------------------
37 * hypre_ParCSRMatrixCreate
38 *--------------------------------------------------------------------------*/
39
40/* If create is called for HYPRE_NO_GLOBAL_PARTITION and row_starts and col_starts are NOT null,
41 then it is assumed that they are array of length 2 containing the start row of
42 the calling processor followed by the start row of the next processor - AHB 6/05 */
43
44hypre_ParCSRMatrix *
45hypre_ParCSRMatrixCreate( MPI_Comm comm,
46 HYPRE_BigInt global_num_rows,
47 HYPRE_BigInt global_num_cols,
48 HYPRE_BigInt *row_starts,
49 HYPRE_BigInt *col_starts,
50 int num_cols_offd,
51 int num_nonzeros_diag,
52 int num_nonzeros_offd)
53{
54 hypre_ParCSRMatrix *matrix;
55 int num_procs, my_id;
56 int local_num_rows, local_num_cols;
57 HYPRE_BigInt first_row_index, first_col_diag;
58
59 matrix = hypre_CTAlloc(hypre_ParCSRMatrix, 1);
60
61 MPI_Comm_rank(comm,&my_id);
62 MPI_Comm_size(comm,&num_procs);
63
64 if (!row_starts)
65 {
66
67#ifdef HYPRE_NO_GLOBAL_PARTITION
68 hypre_GenerateLocalPartitioning(global_num_rows, num_procs, my_id, &row_starts);
69#else
70 hypre_GeneratePartitioning(global_num_rows,num_procs,&row_starts);
71#endif
72 }
73 if (!col_starts)
74 {
75 if (global_num_rows == global_num_cols)
76 {
77 col_starts = row_starts;
78 }
79 else
80 {
81#ifdef HYPRE_NO_GLOBAL_PARTITION
82 hypre_GenerateLocalPartitioning(global_num_cols, num_procs, my_id, &col_starts);
83#else
84 hypre_GeneratePartitioning(global_num_cols,num_procs,&col_starts);
85#endif
86 }
87 }
88
89#ifdef HYPRE_NO_GLOBAL_PARTITION
90 /* row_starts[0] is start of local rows. row_starts[1] is start of next
91 processor's rows */
92 first_row_index = row_starts[0];
93 local_num_rows = row_starts[1]-first_row_index ;
94 first_col_diag = col_starts[0];
95 local_num_cols = col_starts[1]-first_col_diag;
96#else
97 first_row_index = row_starts[my_id];
98 local_num_rows = row_starts[my_id+1]-first_row_index;
99 first_col_diag = col_starts[my_id];
100 local_num_cols = col_starts[my_id+1]-first_col_diag;
101#endif
102
103
104 hypre_ParCSRMatrixComm(matrix) = comm;
105 hypre_ParCSRMatrixDiag(matrix) = hypre_CSRMatrixCreate(local_num_rows,
106 local_num_cols,num_nonzeros_diag);
107 hypre_ParCSRMatrixOffd(matrix) = hypre_CSRMatrixCreate(local_num_rows,
108 num_cols_offd,num_nonzeros_offd);
109 hypre_ParCSRMatrixGlobalNumRows(matrix) = global_num_rows;
110 hypre_ParCSRMatrixGlobalNumCols(matrix) = global_num_cols;
111 hypre_ParCSRMatrixFirstRowIndex(matrix) = first_row_index;
112 hypre_ParCSRMatrixFirstColDiag(matrix) = first_col_diag;
113
114 hypre_ParCSRMatrixLastRowIndex(matrix) = first_row_index + local_num_rows - 1;
115 hypre_ParCSRMatrixLastColDiag(matrix) = first_col_diag + local_num_cols - 1;
116
117 hypre_ParCSRMatrixColMapOffd(matrix) = NULL;
118
119 hypre_ParCSRMatrixAssumedPartition(matrix) = NULL;
120
121
122 /* When NO_GLOBAL_PARTITION is set we could make these null, instead
123 of leaving the range. If that change is made, then when this create
124 is called from functions like the matrix-matrix multiply, be careful
125 not to generate a new partition */
126
127 hypre_ParCSRMatrixRowStarts(matrix) = row_starts;
128 hypre_ParCSRMatrixColStarts(matrix) = col_starts;
129
130 hypre_ParCSRMatrixCommPkg(matrix) = NULL;
131 hypre_ParCSRMatrixCommPkgT(matrix) = NULL;
132
133 /* set defaults */
134 hypre_ParCSRMatrixOwnsData(matrix) = 1;
135 hypre_ParCSRMatrixOwnsRowStarts(matrix) = 1;
136 hypre_ParCSRMatrixOwnsColStarts(matrix) = 1;
137 if (row_starts == col_starts)
138 hypre_ParCSRMatrixOwnsColStarts(matrix) = 0;
139 hypre_ParCSRMatrixRowindices(matrix) = NULL;
140 hypre_ParCSRMatrixRowvalues(matrix) = NULL;
141 hypre_ParCSRMatrixGetrowactive(matrix) = 0;
142
143 return matrix;
144}
145
146/*--------------------------------------------------------------------------
147 * hypre_ParCSRMatrixDestroy
148 *--------------------------------------------------------------------------*/
149
150int
151hypre_ParCSRMatrixDestroy( hypre_ParCSRMatrix *matrix )
152{
153
154 if (matrix)
155 {
156 if ( hypre_ParCSRMatrixOwnsData(matrix) )
157 {
158 hypre_CSRMatrixDestroy(hypre_ParCSRMatrixDiag(matrix));
159 hypre_CSRMatrixDestroy(hypre_ParCSRMatrixOffd(matrix));
160 if (hypre_ParCSRMatrixColMapOffd(matrix))
161 hypre_TFree(hypre_ParCSRMatrixColMapOffd(matrix));
162 if (hypre_ParCSRMatrixCommPkg(matrix))
163 hypre_MatvecCommPkgDestroy(hypre_ParCSRMatrixCommPkg(matrix));
164 if (hypre_ParCSRMatrixCommPkgT(matrix))
165 hypre_MatvecCommPkgDestroy(hypre_ParCSRMatrixCommPkgT(matrix));
166 }
167 if ( hypre_ParCSRMatrixOwnsRowStarts(matrix) )
168 hypre_TFree(hypre_ParCSRMatrixRowStarts(matrix));
169 if ( hypre_ParCSRMatrixOwnsColStarts(matrix) )
170 hypre_TFree(hypre_ParCSRMatrixColStarts(matrix));
171
172 hypre_TFree(hypre_ParCSRMatrixRowindices(matrix));
173 hypre_TFree(hypre_ParCSRMatrixRowvalues(matrix));
174
175
176 if (hypre_ParCSRMatrixAssumedPartition(matrix))
177 hypre_ParCSRMatrixDestroyAssumedPartition(matrix);
178
179
180 hypre_TFree(matrix);
181 }
182
183 return hypre_error_flag;
184}
185
186/*--------------------------------------------------------------------------
187 * hypre_ParCSRMatrixInitialize
188 *--------------------------------------------------------------------------*/
189
190int
191hypre_ParCSRMatrixInitialize( hypre_ParCSRMatrix *matrix )
192{
193 if (!matrix)
194 {
195 hypre_error_in_arg(1);
196 return hypre_error_flag;
197 }
198
199 hypre_CSRMatrixInitialize(hypre_ParCSRMatrixDiag(matrix));
200 hypre_CSRMatrixInitialize(hypre_ParCSRMatrixOffd(matrix));
201 hypre_ParCSRMatrixColMapOffd(matrix) =
202 hypre_CTAlloc(HYPRE_BigInt,hypre_CSRMatrixNumCols(
203 hypre_ParCSRMatrixOffd(matrix)));
204
205 return hypre_error_flag;
206}
207
208/*--------------------------------------------------------------------------
209 * hypre_ParCSRMatrixSetNumNonzeros
210 *--------------------------------------------------------------------------*/
211
212int
213hypre_ParCSRMatrixSetNumNonzeros( hypre_ParCSRMatrix *matrix)
214{
215 MPI_Comm comm;
216 hypre_CSRMatrix *diag;
217 int *diag_i;
218 hypre_CSRMatrix *offd;
219 int *offd_i;
220 int local_num_rows;
221 HYPRE_BigInt total_num_nonzeros;
222 HYPRE_BigInt local_num_nonzeros;
223 if (!matrix)
224 {
225 hypre_error_in_arg(1);
226 return hypre_error_flag;
227 }
228 comm = hypre_ParCSRMatrixComm(matrix);
229 diag = hypre_ParCSRMatrixDiag(matrix);
230 diag_i = hypre_CSRMatrixI(diag);
231 offd = hypre_ParCSRMatrixOffd(matrix);
232 offd_i = hypre_CSRMatrixI(offd);
233 local_num_rows = hypre_CSRMatrixNumRows(diag);
234
235 local_num_nonzeros = diag_i[local_num_rows] + offd_i[local_num_rows];
236 MPI_Allreduce(&local_num_nonzeros, &total_num_nonzeros, 1, MPI_HYPRE_BIG_INT,
237 MPI_SUM, comm);
238 hypre_ParCSRMatrixNumNonzeros(matrix) = total_num_nonzeros;
239 return hypre_error_flag;
240}
241
242/*--------------------------------------------------------------------------
243 * hypre_ParCSRMatrixSetDNumNonzeros
244 *--------------------------------------------------------------------------*/
245
246int
247hypre_ParCSRMatrixSetDNumNonzeros( hypre_ParCSRMatrix *matrix)
248{
249 MPI_Comm comm;
250 hypre_CSRMatrix *diag;
251 int *diag_i;
252 hypre_CSRMatrix *offd;
253 int *offd_i;
254 int local_num_rows;
255 double total_num_nonzeros;
256 double local_num_nonzeros;
257 if (!matrix)
258 {
259 hypre_error_in_arg(1);
260 return hypre_error_flag;
261 }
262 comm = hypre_ParCSRMatrixComm(matrix);
263 diag = hypre_ParCSRMatrixDiag(matrix);
264 diag_i = hypre_CSRMatrixI(diag);
265 offd = hypre_ParCSRMatrixOffd(matrix);
266 offd_i = hypre_CSRMatrixI(offd);
267 local_num_rows = hypre_CSRMatrixNumRows(diag);
268 local_num_nonzeros = (double) diag_i[local_num_rows]
269 + (double) offd_i[local_num_rows];
270 MPI_Allreduce(&local_num_nonzeros, &total_num_nonzeros, 1, MPI_DOUBLE,
271 MPI_SUM, comm);
272 hypre_ParCSRMatrixDNumNonzeros(matrix) = total_num_nonzeros;
273 return hypre_error_flag;
274}
275
276/*--------------------------------------------------------------------------
277 * hypre_ParCSRMatrixSetDataOwner
278 *--------------------------------------------------------------------------*/
279
280int
281hypre_ParCSRMatrixSetDataOwner( hypre_ParCSRMatrix *matrix,
282 int owns_data )
283{
284 if (!matrix)
285 {
286 hypre_error_in_arg(1);
287 return hypre_error_flag;
288 }
289
290 hypre_ParCSRMatrixOwnsData(matrix) = owns_data;
291
292 return hypre_error_flag;
293}
294
295/*--------------------------------------------------------------------------
296 * hypre_ParCSRMatrixSetRowStartsOwner
297 *--------------------------------------------------------------------------*/
298
299int
300hypre_ParCSRMatrixSetRowStartsOwner( hypre_ParCSRMatrix *matrix,
301 int owns_row_starts )
302{
303 if (!matrix)
304 {
305 hypre_error_in_arg(1);
306 return hypre_error_flag;
307 }
308
309 hypre_ParCSRMatrixOwnsRowStarts(matrix) = owns_row_starts;
310
311 return hypre_error_flag;
312}
313
314/*--------------------------------------------------------------------------
315 * hypre_ParCSRMatrixSetColStartsOwner
316 *--------------------------------------------------------------------------*/
317
318int
319hypre_ParCSRMatrixSetColStartsOwner( hypre_ParCSRMatrix *matrix,
320 int owns_col_starts )
321{
322 if (!matrix)
323 {
324 hypre_error_in_arg(1);
325 return hypre_error_flag;
326 }
327
328 hypre_ParCSRMatrixOwnsColStarts(matrix) = owns_col_starts;
329
330 return hypre_error_flag;
331}
332
333/*--------------------------------------------------------------------------
334 * hypre_ParCSRMatrixRead
335 *--------------------------------------------------------------------------*/
336
337hypre_ParCSRMatrix *
338hypre_ParCSRMatrixRead( MPI_Comm comm,
339 const char *file_name )
340{
341 hypre_ParCSRMatrix *matrix;
342 hypre_CSRMatrix *diag;
343 hypre_CSRMatrix *offd;
344 int my_id, i, num_procs;
345 char new_file_d[80], new_file_o[80], new_file_info[80];
346 HYPRE_BigInt global_num_rows, global_num_cols;
347 int num_cols_offd;
348 int local_num_rows;
349 HYPRE_BigInt *row_starts;
350 HYPRE_BigInt *col_starts;
351 HYPRE_BigInt *col_map_offd;
352 FILE *fp;
353 int equal = 1;
354#ifdef HYPRE_NO_GLOBAL_PARTITION
355 HYPRE_BigInt row_s, row_e, col_s, col_e;
356#endif
357
358
359 MPI_Comm_rank(comm,&my_id);
360 MPI_Comm_size(comm,&num_procs);
361#ifdef HYPRE_NO_GLOBAL_PARTITION
362 row_starts = hypre_CTAlloc(HYPRE_BigInt, 2);
363 col_starts = hypre_CTAlloc(HYPRE_BigInt, 2);
364#else
365 row_starts = hypre_CTAlloc(HYPRE_BigInt, num_procs+1);
366 col_starts = hypre_CTAlloc(HYPRE_BigInt, num_procs+1);
367#endif
368 sprintf(new_file_d,"%s.D.%d",file_name,my_id);
369 sprintf(new_file_o,"%s.O.%d",file_name,my_id);
370 sprintf(new_file_info,"%s.INFO.%d",file_name,my_id);
371 fp = fopen(new_file_info, "r");
372#ifdef HYPRE_LONG_LONG
373 fscanf(fp, "%lld", &global_num_rows);
374 fscanf(fp, "%lld", &global_num_cols);
375#else
376 fscanf(fp, "%d", &global_num_rows);
377 fscanf(fp, "%d", &global_num_cols);
378#endif
379 fscanf(fp, "%d", &num_cols_offd);
380#ifdef HYPRE_NO_GLOBAL_PARTITION
381 /* the bgl input file should only contain the EXACT range for local processor */
382#ifdef HYPRE_LONG_LONG
383 fscanf(fp, "%lld %lld %lld %lld", &row_s, &row_e, &col_s, &col_e);
384#else
385 fscanf(fp, "%d %d %d %d", &row_s, &row_e, &col_s, &col_e);
386#endif
387 row_starts[0] = row_s;
388 row_starts[1] = row_e;
389 col_starts[0] = col_s;
390 col_starts[1] = col_e;
391
392#else
393 for (i=0; i < num_procs; i++)
394#ifdef HYPRE_LONG_LONG
395 fscanf(fp, "%lld %lld", &row_starts[i], &col_starts[i]);
396#else
397 fscanf(fp, "%d %d", &row_starts[i], &col_starts[i]);
398#endif
399 row_starts[num_procs] = global_num_rows;
400 col_starts[num_procs] = global_num_cols;
401#endif
402
403 col_map_offd = hypre_CTAlloc(HYPRE_BigInt, num_cols_offd);
404
405 for (i=0; i < num_cols_offd; i++)
406#ifdef HYPRE_LONG_LONG
407 fscanf(fp, "%lld", &col_map_offd[i]);
408#else
409 fscanf(fp, "%d", &col_map_offd[i]);
410#endif
411
412 fclose(fp);
413
414#ifdef HYPRE_NO_GLOBAL_PARTITION
415 for (i=1; i >= 0; i--)
416 if (row_starts[i] != col_starts[i])
417 {
418 equal = 0;
419 break;
420 }
421#else
422 for (i=num_procs; i >= 0; i--)
423 if (row_starts[i] != col_starts[i])
424 {
425 equal = 0;
426 break;
427 }
428#endif
429 if (equal)
430 {
431 hypre_TFree(col_starts);
432 col_starts = row_starts;
433 }
434
435
436
437
438 diag = hypre_CSRMatrixRead(new_file_d);
439 local_num_rows = hypre_CSRMatrixNumRows(diag);
440
441 if (num_cols_offd)
442 {
443 offd = hypre_CSRMatrixRead(new_file_o);
444 }
445 else
446 {
447 offd = hypre_CSRMatrixCreate(local_num_rows,0,0);
448 hypre_CSRMatrixInitialize(offd);
449 }
450
451
452 matrix = hypre_CTAlloc(hypre_ParCSRMatrix, 1);
453
454 hypre_ParCSRMatrixComm(matrix) = comm;
455 hypre_ParCSRMatrixGlobalNumRows(matrix) = global_num_rows;
456 hypre_ParCSRMatrixGlobalNumCols(matrix) = global_num_cols;
457#ifdef HYPRE_NO_GLOBAL_PARTITION
458 hypre_ParCSRMatrixFirstRowIndex(matrix) = row_s;
459 hypre_ParCSRMatrixFirstColDiag(matrix) = col_s;
460 hypre_ParCSRMatrixLastRowIndex(matrix) = row_e - 1;
461 hypre_ParCSRMatrixLastColDiag(matrix) = col_e - 1;
462#else
463 hypre_ParCSRMatrixFirstRowIndex(matrix) = row_starts[my_id];
464 hypre_ParCSRMatrixFirstColDiag(matrix) = col_starts[my_id];
465 hypre_ParCSRMatrixLastRowIndex(matrix) = row_starts[my_id+1]-1;
466 hypre_ParCSRMatrixLastColDiag(matrix) = col_starts[my_id+1]-1;
467#endif
468
469 hypre_ParCSRMatrixRowStarts(matrix) = row_starts;
470 hypre_ParCSRMatrixColStarts(matrix) = col_starts;
471 hypre_ParCSRMatrixCommPkg(matrix) = NULL;
472
473 /* set defaults */
474 hypre_ParCSRMatrixOwnsData(matrix) = 1;
475 hypre_ParCSRMatrixOwnsRowStarts(matrix) = 1;
476 hypre_ParCSRMatrixOwnsColStarts(matrix) = 1;
477 if (row_starts == col_starts)
478 hypre_ParCSRMatrixOwnsColStarts(matrix) = 0;
479
480 hypre_ParCSRMatrixDiag(matrix) = diag;
481 hypre_ParCSRMatrixOffd(matrix) = offd;
482 if (num_cols_offd)
483 hypre_ParCSRMatrixColMapOffd(matrix) = col_map_offd;
484 else
485 hypre_ParCSRMatrixColMapOffd(matrix) = NULL;
486
487 return matrix;
488}
489
490/*--------------------------------------------------------------------------
491 * hypre_ParCSRMatrixPrint
492 *--------------------------------------------------------------------------*/
493
494int
495hypre_ParCSRMatrixPrint( hypre_ParCSRMatrix *matrix,
496 const char *file_name )
497{
498 MPI_Comm comm;
499 HYPRE_BigInt global_num_rows;
500 HYPRE_BigInt global_num_cols;
501 HYPRE_BigInt *col_map_offd;
502#ifndef HYPRE_NO_GLOBAL_PARTITION
503 HYPRE_BigInt *row_starts;
504 HYPRE_BigInt *col_starts;
505#endif
506 int my_id, i, num_procs;
507 char new_file_d[80], new_file_o[80], new_file_info[80];
508 FILE *fp;
509 int num_cols_offd = 0;
510#ifdef HYPRE_NO_GLOBAL_PARTITION
511 HYPRE_BigInt row_s, row_e, col_s, col_e;
512#endif
513 if (!matrix)
514 {
515 hypre_error_in_arg(1);
516 return hypre_error_flag;
517 }
518 comm = hypre_ParCSRMatrixComm(matrix);
519 global_num_rows = hypre_ParCSRMatrixGlobalNumRows(matrix);
520 global_num_cols = hypre_ParCSRMatrixGlobalNumCols(matrix);
521 col_map_offd = hypre_ParCSRMatrixColMapOffd(matrix);
522#ifndef HYPRE_NO_GLOBAL_PARTITION
523 row_starts = hypre_ParCSRMatrixRowStarts(matrix);
524 col_starts = hypre_ParCSRMatrixColStarts(matrix);
525#endif
526 if (hypre_ParCSRMatrixOffd(matrix))
527 num_cols_offd = hypre_CSRMatrixNumCols(hypre_ParCSRMatrixOffd(matrix));
528
529 MPI_Comm_rank(comm, &my_id);
530 MPI_Comm_size(comm, &num_procs);
531
532 sprintf(new_file_d,"%s.D.%d",file_name,my_id);
533 sprintf(new_file_o,"%s.O.%d",file_name,my_id);
534 sprintf(new_file_info,"%s.INFO.%d",file_name,my_id);
535 hypre_CSRMatrixPrint(hypre_ParCSRMatrixDiag(matrix),new_file_d);
536 if (num_cols_offd != 0)
537 hypre_CSRMatrixPrint(hypre_ParCSRMatrixOffd(matrix),new_file_o);
538
539 fp = fopen(new_file_info, "w");
540#ifdef HYPRE_LONG_LONG
541 fprintf(fp, "%lld\n", global_num_rows);
542 fprintf(fp, "%lld\n", global_num_cols);
543#else
544 fprintf(fp, "%d\n", global_num_rows);
545 fprintf(fp, "%d\n", global_num_cols);
546#endif
547 fprintf(fp, "%d\n", num_cols_offd);
548#ifdef HYPRE_NO_GLOBAL_PARTITION
549 row_s = hypre_ParCSRMatrixFirstRowIndex(matrix);
550 row_e = hypre_ParCSRMatrixLastRowIndex(matrix);
551 col_s = hypre_ParCSRMatrixFirstColDiag(matrix);
552 col_e = hypre_ParCSRMatrixLastColDiag(matrix);
553 /* add 1 to the ends because this is a starts partition */
554#ifdef HYPRE_LONG_LONG
555 fprintf(fp, "%lld %lld %lld %lld\n", row_s, row_e + 1, col_s, col_e + 1);
556#else
557 fprintf(fp, "%d %d %d %d\n", row_s, row_e + 1, col_s, col_e + 1);
558#endif
559#else
560 for (i=0; i < num_procs; i++)
561#ifdef HYPRE_LONG_LONG
562 fprintf(fp, "%lld %lld\n", row_starts[i], col_starts[i]);
563#else
564 fprintf(fp, "%d %d\n", row_starts[i], col_starts[i]);
565#endif
566#endif
567 for (i=0; i < num_cols_offd; i++)
568#ifdef HYPRE_LONG_LONG
569 fprintf(fp, "%lld\n", col_map_offd[i]);
570#else
571 fprintf(fp, "%d\n", col_map_offd[i]);
572#endif
573 fclose(fp);
574
575 return hypre_error_flag;
576}
577
578/*--------------------------------------------------------------------------
579 * hypre_ParCSRMatrixPrintIJ
580 *--------------------------------------------------------------------------*/
581
582int
583hypre_ParCSRMatrixPrintIJ( const hypre_ParCSRMatrix *matrix,
584 const int base_i,
585 const int base_j,
586 const char *filename )
587{
588 MPI_Comm comm;
589 HYPRE_BigInt global_num_rows;
590 HYPRE_BigInt global_num_cols;
591 HYPRE_BigInt first_row_index;
592 HYPRE_BigInt first_col_diag;
593 hypre_CSRMatrix *diag;
594 hypre_CSRMatrix *offd;
595 HYPRE_BigInt *col_map_offd;
596 int num_rows;
597 HYPRE_BigInt *row_starts;
598 HYPRE_BigInt *col_starts;
599 double *diag_data;
600 int *diag_i;
601 int *diag_j;
602 double *offd_data;
603 int *offd_i;
604 int *offd_j;
605 int myid, num_procs, i, j;
606 HYPRE_BigInt I, J;
607 char new_filename[255];
608 FILE *file;
609 int num_cols_offd, num_nonzeros_diag, num_nonzeros_offd;
610 int num_cols;
611 HYPRE_BigInt row, col;
612
613 if (!matrix)
614 {
615 hypre_error_in_arg(1);
616 return hypre_error_flag;
617 }
618 comm = hypre_ParCSRMatrixComm(matrix);
619 global_num_rows = hypre_ParCSRMatrixGlobalNumRows(matrix);
620 global_num_cols = hypre_ParCSRMatrixGlobalNumCols(matrix);
621 first_row_index = hypre_ParCSRMatrixFirstRowIndex(matrix);
622 first_col_diag = hypre_ParCSRMatrixFirstColDiag(matrix);
623 diag = hypre_ParCSRMatrixDiag(matrix);
624 offd = hypre_ParCSRMatrixOffd(matrix);
625 col_map_offd = hypre_ParCSRMatrixColMapOffd(matrix);
626 num_rows = hypre_ParCSRMatrixNumRows(matrix);
627 row_starts = hypre_ParCSRMatrixRowStarts(matrix);
628 col_starts = hypre_ParCSRMatrixColStarts(matrix);
629 MPI_Comm_rank(comm, &myid);
630 MPI_Comm_size(comm, &num_procs);
631
632 sprintf(new_filename,"%s.%05d", filename, myid);
633
634 if ((file = fopen(new_filename, "w")) == NULL)
635 {
636 printf("Error: can't open output file %s\n", new_filename);
637 hypre_error(HYPRE_ERROR_GENERIC);
638 return hypre_error_flag;
639 }
640
641 num_cols = hypre_CSRMatrixNumCols(diag);
642 num_cols_offd = hypre_CSRMatrixNumCols(offd);
643 num_nonzeros_offd = hypre_CSRMatrixNumNonzeros(offd);
644 num_nonzeros_diag = hypre_CSRMatrixNumNonzeros(diag);
645
646 diag_data = hypre_CSRMatrixData(diag);
647 diag_i = hypre_CSRMatrixI(diag);
648 diag_j = hypre_CSRMatrixJ(diag);
649 offd_i = hypre_CSRMatrixI(offd);
650 if (num_nonzeros_offd)
651 {
652 offd_data = hypre_CSRMatrixData(offd);
653 offd_j = hypre_CSRMatrixJ(offd);
654 }
655
656#ifdef HYPRE_LONG_LONG
657 fprintf(file, "%lld %lld\n", global_num_rows, global_num_cols);
658#else
659 fprintf(file, "%d %d\n", global_num_rows, global_num_cols);
660#endif
661 fprintf(file, "%d %d %d\n", num_rows, num_cols, num_cols_offd);
662 fprintf(file, "%d %d\n", num_nonzeros_diag, num_nonzeros_offd);
663
664#ifdef HYPRE_NO_GLOBAL_PARTITION
665 for (i=0; i <= 1; i++)
666 {
667 row = row_starts[i]+base_i;
668 col = col_starts[i]+base_j;
669#ifdef HYPRE_LONG_LONG
670 fprintf(file, "%lld %lld\n", row, col);
671#else
672 fprintf(file, "%d %d\n", row, col);
673#endif
674 }
675#else
676 for (i=0; i <= num_procs; i++)
677 {
678 row = row_starts[i]+base_i;
679 col = col_starts[i]+base_j;
680#ifdef HYPRE_LONG_LONG
681 fprintf(file, "%lld %lld\n", row, col);
682#else
683 fprintf(file, "%d %d\n", row, col);
684#endif
685 }
686#endif
687 for (i = 0; i < num_rows; i++)
688 {
689 I = first_row_index + i + base_i;
690
691 /* print diag columns */
692 for (j = diag_i[i]; j < diag_i[i+1]; j++)
693 {
694 J = first_col_diag + diag_j[j] + base_j;
695#ifdef HYPRE_LONG_LONG
696 if ( diag_data )
697 fprintf(file, "%lld %lld %.14e\n", I, J, diag_data[j]);
698 else
699 fprintf(file, "%lld %lld\n", I, J);
700#else
701 if ( diag_data )
702 fprintf(file, "%d %d %.14e\n", I, J, diag_data[j]);
703 else
704 fprintf(file, "%d %d\n", I, J);
705#endif
706 }
707
708 /* print offd columns */
709 if ( num_nonzeros_offd )
710 {
711 for (j = offd_i[i]; j < offd_i[i+1]; j++)
712 {
713 J = col_map_offd[offd_j[j]] + base_j;
714#ifdef HYPRE_LONG_LONG
715 if ( offd_data )
716 fprintf(file, "%lld %lld %.14e\n", I, J, offd_data[j]);
717 else
718 fprintf(file, "%lld %lld\n", I, J );
719#else
720 if ( offd_data )
721 fprintf(file, "%d %d %.14e\n", I, J, offd_data[j]);
722 else
723 fprintf(file, "%d %d\n", I, J );
724#endif
725 }
726 }
727 }
728
729 fclose(file);
730
731 return hypre_error_flag;
732}
733
734/*--------------------------------------------------------------------------
735 * hypre_ParCSRMatrixReadIJ
736 *--------------------------------------------------------------------------*/
737
738int
739hypre_ParCSRMatrixReadIJ( MPI_Comm comm,
740 const char *filename,
741 int *base_i_ptr,
742 int *base_j_ptr,
743 hypre_ParCSRMatrix **matrix_ptr)
744{
745 HYPRE_BigInt global_num_rows;
746 HYPRE_BigInt global_num_cols;
747 HYPRE_BigInt first_row_index;
748 HYPRE_BigInt first_col_diag;
749 HYPRE_BigInt last_col_diag;
750 hypre_ParCSRMatrix *matrix;
751 hypre_CSRMatrix *diag;
752 hypre_CSRMatrix *offd;
753 HYPRE_BigInt *col_map_offd;
754 HYPRE_BigInt *row_starts;
755 HYPRE_BigInt *col_starts;
756 int num_rows;
757 int base_i, base_j;
758 double *diag_data;
759 int *diag_i;
760 int *diag_j;
761 double *offd_data;
762 int *offd_i;
763 int *offd_j;
764 HYPRE_BigInt *aux_offd_j;
765 int myid, num_procs, i, j;
766 HYPRE_BigInt I, J;
767 char new_filename[255];
768 FILE *file;
769 int num_cols_offd, num_nonzeros_diag, num_nonzeros_offd;
770 int equal, i_col, num_cols;
771 int diag_cnt, offd_cnt, row_cnt;
772 double data;
773
774 MPI_Comm_size(comm, &num_procs);
775 MPI_Comm_rank(comm, &myid);
776
777 sprintf(new_filename,"%s.%05d", filename, myid);
778
779 if ((file = fopen(new_filename, "r")) == NULL)
780 {
781 printf("Error: can't open output file %s\n", new_filename);
782 hypre_error(HYPRE_ERROR_GENERIC);
783 return hypre_error_flag;
784 }
785
786#ifdef HYPRE_LONG_LONG
787 fscanf(file, "%lld %lld", &global_num_rows, &global_num_cols);
788#else
789 fscanf(file, "%d %d", &global_num_rows, &global_num_cols);
790#endif
791 fscanf(file, "%d %d %d", &num_rows, &num_cols, &num_cols_offd);
792 fscanf(file, "%d %d", &num_nonzeros_diag, &num_nonzeros_offd);
793
794 row_starts = hypre_CTAlloc(HYPRE_BigInt,num_procs+1);
795 col_starts = hypre_CTAlloc(HYPRE_BigInt,num_procs+1);
796
797 for (i = 0; i <= num_procs; i++)
798#ifdef HYPRE_LONG_LONG
799 fscanf(file, "%lld %lld", &row_starts[i], &col_starts[i]);
800#else
801 fscanf(file, "%d %d", &row_starts[i], &col_starts[i]);
802#endif
803
804 base_i = row_starts[0];
805 base_j = col_starts[0];
806
807 equal = 1;
808 for (i = 0; i <= num_procs; i++)
809 {
810 row_starts[i] -= base_i;
811 col_starts[i] -= base_j;
812 if (row_starts[i] != col_starts[i]) equal = 0;
813 }
814
815 if (equal)
816 {
817 hypre_TFree(col_starts);
818 col_starts = row_starts;
819 }
820 matrix = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
821 row_starts, col_starts,
822 num_cols_offd,
823 num_nonzeros_diag, num_nonzeros_offd);
824 hypre_ParCSRMatrixInitialize(matrix);
825
826 diag = hypre_ParCSRMatrixDiag(matrix);
827 offd = hypre_ParCSRMatrixOffd(matrix);
828
829 diag_data = hypre_CSRMatrixData(diag);
830 diag_i = hypre_CSRMatrixI(diag);
831 diag_j = hypre_CSRMatrixJ(diag);
832
833 offd_i = hypre_CSRMatrixI(offd);
834 if (num_nonzeros_offd)
835 {
836 offd_data = hypre_CSRMatrixData(offd);
837 offd_j = hypre_CSRMatrixJ(offd);
838 aux_offd_j = hypre_CTAlloc(HYPRE_BigInt, num_nonzeros_offd);
839 }
840
841 first_row_index = hypre_ParCSRMatrixFirstRowIndex(matrix);
842 first_col_diag = hypre_ParCSRMatrixFirstColDiag(matrix);
843 last_col_diag = first_col_diag+num_cols-1;
844
845 diag_cnt = 0;
846 offd_cnt = 0;
847 row_cnt = 0;
848 for (i = 0; i < num_nonzeros_diag+num_nonzeros_offd; i++)
849 {
850 /* read values */
851#ifdef HYPRE_LONG_LONG
852 fscanf(file, "%lld %lld %le", &I, &J, &data);
853#else
854 fscanf(file, "%d %d %le", &I, &J, &data);
855#endif
856 I = I-base_i-first_row_index;
857 J -= base_j;
858 if (I > row_cnt)
859 {
860 diag_i[I] = diag_cnt;
861 offd_i[I] = offd_cnt;
862 row_cnt++;
863 }
864 if (J < first_col_diag || J > last_col_diag)
865 {
866 aux_offd_j[offd_cnt] = J;
867 offd_data[offd_cnt++] = data;
868 }
869 else
870 {
871 diag_j[diag_cnt] = (int)(J - first_col_diag);
872 diag_data[diag_cnt++] = data;
873 }
874 }
875 diag_i[num_rows] = diag_cnt;
876 offd_i[num_rows] = offd_cnt;
877
878 fclose(file);
879
880 /* generate col_map_offd */
881 if (num_nonzeros_offd)
882 {
883 hypre_BigQsort0(aux_offd_j,0,num_nonzeros_offd-1);
884 col_map_offd = hypre_ParCSRMatrixColMapOffd(matrix);
885 col_map_offd[0] = aux_offd_j[0];
886 offd_cnt = 0;
887 for (i=1; i < num_nonzeros_offd; i++)
888 {
889 if (aux_offd_j[i] > col_map_offd[offd_cnt])
890 col_map_offd[++offd_cnt] = aux_offd_j[i];
891 }
892 for (i=0; i < num_nonzeros_offd; i++)
893 {
894 offd_j[i] = hypre_BigBinarySearch(col_map_offd, aux_offd_j[i], num_cols_offd);
895 }
896 hypre_TFree(aux_offd_j);
897 }
898
899 /* move diagonal element in first position in each row */
900 for (i=0; i < num_rows; i++)
901 {
902 i_col = diag_i[i];
903 for (j=i_col; j < diag_i[i+1]; j++)
904 {
905 if (diag_j[j] == i)
906 {
907 diag_j[j] = diag_j[i_col];
908 data = diag_data[j];
909 diag_data[j] = diag_data[i_col];
910 diag_data[i_col] = data;
911 diag_j[i_col] = i;
912 break;
913 }
914 }
915 }
916
917 *base_i_ptr = base_i;
918 *base_j_ptr = base_j;
919 *matrix_ptr = matrix;
920
921 return hypre_error_flag;
922}
923
924/*--------------------------------------------------------------------------
925 * hypre_ParCSRMatrixGetLocalRange
926 * returns the row numbers of the rows stored on this processor.
927 * "End" is actually the row number of the last row on this processor.
928 *--------------------------------------------------------------------------*/
929
930int
931hypre_ParCSRMatrixGetLocalRange( hypre_ParCSRMatrix *matrix,
932 HYPRE_BigInt *row_start,
933 HYPRE_BigInt *row_end,
934 HYPRE_BigInt *col_start,
935 HYPRE_BigInt *col_end )
936{
937
938#ifdef HYPRE_NO_GLOBAL_PARTITION
939 if (!matrix)
940 {
941 hypre_error_in_arg(1);
942 return hypre_error_flag;
943 }
944 *row_start = hypre_ParCSRMatrixFirstRowIndex(matrix);
945 *row_end = hypre_ParCSRMatrixLastRowIndex(matrix);
946 *col_start = hypre_ParCSRMatrixFirstColDiag(matrix);
947 *col_end = hypre_ParCSRMatrixLastColDiag(matrix);
948#else
949 int my_id;
950
951 if (!matrix)
952 {
953 hypre_error_in_arg(1);
954 return hypre_error_flag;
955 }
956
957 MPI_Comm_rank( hypre_ParCSRMatrixComm(matrix), &my_id );
958
959 *row_start = hypre_ParCSRMatrixRowStarts(matrix)[ my_id ];
960 *row_end = hypre_ParCSRMatrixRowStarts(matrix)[ my_id + 1 ]-1;
961 *col_start = hypre_ParCSRMatrixColStarts(matrix)[ my_id ];
962 *col_end = hypre_ParCSRMatrixColStarts(matrix)[ my_id + 1 ]-1;
963#endif
964
965 return hypre_error_flag;
966}
967
968/*--------------------------------------------------------------------------
969 * hypre_ParCSRMatrixGetRow
970 * Returns global column indices and/or values for a given row in the global
971 * matrix. Global row number is used, but the row must be stored locally or
972 * an error is returned. This implementation copies from the two matrices that
973 * store the local data, storing them in the hypre_ParCSRMatrix structure.
974 * Only a single row can be accessed via this function at any one time; the
975 * corresponding RestoreRow function must be called, to avoid bleeding memory,
976 * and to be able to look at another row.
977 * Either one of col_ind and values can be left null, and those values will
978 * not be returned.
979 * All indices are returned in 0-based indexing, no matter what is used under
980 * the hood. EXCEPTION: currently this only works if the local CSR matrices
981 * use 0-based indexing.
982 * This code, semantics, implementation, etc., are all based on PETSc's MPI_AIJ
983 * matrix code, adjusted for our data and software structures.
984 * AJC 4/99.
985 *--------------------------------------------------------------------------*/
986
987int
988hypre_ParCSRMatrixGetRow( hypre_ParCSRMatrix *mat,
989 HYPRE_BigInt row,
990 int *size,
991 HYPRE_BigInt **col_ind,
992 double **values )
993{
994 int my_id;
995 HYPRE_BigInt row_start, row_end;
996 hypre_CSRMatrix *Aa;
997 hypre_CSRMatrix *Ba;
998
999 if (!mat)
1000 {
1001 hypre_error_in_arg(1);
1002 return hypre_error_flag;
1003 }
1004 Aa = (hypre_CSRMatrix *) hypre_ParCSRMatrixDiag(mat);
1005 Ba = (hypre_CSRMatrix *) hypre_ParCSRMatrixOffd(mat);
1006
1007 if (hypre_ParCSRMatrixGetrowactive(mat)) return(-1);
1008
1009 MPI_Comm_rank( hypre_ParCSRMatrixComm(mat), &my_id );
1010
1011 hypre_ParCSRMatrixGetrowactive(mat) = 1;
1012#ifdef HYPRE_NO_GLOBAL_PARTITION
1013 row_start = hypre_ParCSRMatrixFirstRowIndex(mat);
1014 row_end = hypre_ParCSRMatrixLastRowIndex(mat) + 1;
1015#else
1016 row_end = hypre_ParCSRMatrixRowStarts(mat)[ my_id + 1 ];
1017 row_start = hypre_ParCSRMatrixRowStarts(mat)[ my_id ];
1018#endif
1019 if (row < row_start || row >= row_end) return(-1);
1020
1021 /* if buffer is not allocated and some information is requested,
1022 allocate buffer */
1023 if (!hypre_ParCSRMatrixRowvalues(mat) && ( col_ind || values ))
1024 {
1025 /*
1026 allocate enough space to hold information from the longest row.
1027 */
1028 int max = 1,tmp;
1029 int i;
1030 int m = row_end-row_start;
1031
1032 for ( i=0; i<m; i++ ) {
1033 tmp = hypre_CSRMatrixI(Aa)[i+1] - hypre_CSRMatrixI(Aa)[i] +
1034 hypre_CSRMatrixI(Ba)[i+1] - hypre_CSRMatrixI(Ba)[i];
1035 if (max < tmp) { max = tmp; }
1036 }
1037
1038 hypre_ParCSRMatrixRowvalues(mat) = (double *) hypre_CTAlloc
1039 ( double, max );
1040 hypre_ParCSRMatrixRowindices(mat) = (HYPRE_BigInt *) hypre_CTAlloc
1041 ( HYPRE_BigInt, max );
1042 }
1043
1044 /* Copy from dual sequential matrices into buffer */
1045 {
1046 double *vworkA, *vworkB, *v_p;
1047 int i, *cworkA, *cworkB;
1048 HYPRE_BigInt cstart = hypre_ParCSRMatrixFirstColDiag(mat);
1049 int nztot, nzA, nzB, lrow=(int)(row-row_start);
1050 HYPRE_BigInt *cmap, *idx_p;
1051
1052 nzA = hypre_CSRMatrixI(Aa)[lrow+1]-hypre_CSRMatrixI(Aa)[lrow];
1053 cworkA = &( hypre_CSRMatrixJ(Aa)[ hypre_CSRMatrixI(Aa)[lrow] ] );
1054 vworkA = &( hypre_CSRMatrixData(Aa)[ hypre_CSRMatrixI(Aa)[lrow] ] );
1055
1056 nzB = hypre_CSRMatrixI(Ba)[lrow+1]-hypre_CSRMatrixI(Ba)[lrow];
1057 cworkB = &( hypre_CSRMatrixJ(Ba)[ hypre_CSRMatrixI(Ba)[lrow] ] );
1058 vworkB = &( hypre_CSRMatrixData(Ba)[ hypre_CSRMatrixI(Ba)[lrow] ] );
1059
1060 nztot = nzA + nzB;
1061
1062 cmap = hypre_ParCSRMatrixColMapOffd(mat);
1063
1064 if (values || col_ind) {
1065 if (nztot) {
1066 /* Sort by increasing column numbers, assuming A and B already sorted */
1067 int imark = -1;
1068 if (values) {
1069 *values = v_p = hypre_ParCSRMatrixRowvalues(mat);
1070 for ( i=0; i<nzB; i++ ) {
1071 if (cmap[cworkB[i]] < cstart) v_p[i] = vworkB[i];
1072 else break;
1073 }
1074 imark = i;
1075 for ( i=0; i<nzA; i++ ) v_p[imark+i] = vworkA[i];
1076 for ( i=imark; i<nzB; i++ ) v_p[nzA+i] = vworkB[i];
1077 }
1078 if (col_ind) {
1079 *col_ind = idx_p = hypre_ParCSRMatrixRowindices(mat);
1080 if (imark > -1) {
1081 for ( i=0; i<imark; i++ ) {
1082 idx_p[i] = cmap[cworkB[i]];
1083 }
1084 } else {
1085 for ( i=0; i<nzB; i++ ) {
1086 if (cmap[cworkB[i]] < cstart) idx_p[i] = cmap[cworkB[i]];
1087 else break;
1088 }
1089 imark = i;
1090 }
1091 for ( i=0; i<nzA; i++ ) idx_p[imark+i] = cstart + cworkA[i];
1092 for ( i=imark; i<nzB; i++ ) idx_p[nzA+i] = cmap[cworkB[i]];
1093 }
1094 }
1095 else {
1096 if (col_ind) *col_ind = 0;
1097 if (values) *values = 0;
1098 }
1099 }
1100 *size = nztot;
1101
1102 } /* End of copy */
1103
1104
1105 return hypre_error_flag;
1106}
1107
1108/*--------------------------------------------------------------------------
1109 * hypre_ParCSRMatrixRestoreRow
1110 *--------------------------------------------------------------------------*/
1111
1112int
1113hypre_ParCSRMatrixRestoreRow( hypre_ParCSRMatrix *matrix,
1114 HYPRE_BigInt row,
1115 int *size,
1116 HYPRE_BigInt **col_ind,
1117 double **values )
1118{
1119
1120 if (!hypre_ParCSRMatrixGetrowactive(matrix)) {
1121 hypre_error(HYPRE_ERROR_GENERIC);
1122 return hypre_error_flag;
1123 }
1124
1125 hypre_ParCSRMatrixGetrowactive(matrix)=0;
1126
1127 return hypre_error_flag;
1128}
1129
1130/*--------------------------------------------------------------------------
1131 * hypre_CSRMatrixToParCSRMatrix:
1132 * generates a ParCSRMatrix distributed across the processors in comm
1133 * from a CSRMatrix on proc 0 .
1134 *
1135 * This shouldn't be used with the HYPRE_NO_GLOBAL_PARTITON option
1136 *
1137 *--------------------------------------------------------------------------*/
1138
1139hypre_ParCSRMatrix *
1140hypre_CSRMatrixToParCSRMatrix( MPI_Comm comm, hypre_CSRMatrix *A,
1141 HYPRE_BigInt *row_starts,
1142 HYPRE_BigInt *col_starts )
1143{
1144 HYPRE_BigInt *global_data;
1145 HYPRE_BigInt global_size;
1146 HYPRE_BigInt global_num_rows;
1147 HYPRE_BigInt global_num_cols;
1148 int *local_num_rows;
1149
1150 int num_procs, my_id;
1151 int *local_num_nonzeros;
1152 int num_nonzeros;
1153
1154 double *a_data;
1155 int *a_i;
1156 int *a_j;
1157
1158 hypre_CSRMatrix *local_A;
1159
1160 MPI_Request *requests;
1161 MPI_Status *status, status0;
1162 MPI_Datatype *csr_matrix_datatypes;
1163
1164 hypre_ParCSRMatrix *par_matrix;
1165
1166 HYPRE_BigInt first_col_diag;
1167 HYPRE_BigInt last_col_diag;
1168
1169 int i, j, ind;
1170
1171 MPI_Comm_rank(comm, &my_id);
1172 MPI_Comm_size(comm, &num_procs);
1173
1174 global_data = hypre_CTAlloc(HYPRE_BigInt, 2*num_procs+6);
1175 if (my_id == 0)
1176 {
1177 global_size = 3;
1178 if (row_starts)
1179 {
1180 if (col_starts)
1181 {
1182 if (col_starts != row_starts)
1183 {
1184 /* contains code for what to expect,
1185 if 0: row_starts = col_starts, only row_starts given
1186 if 1: only row_starts given, col_starts = NULL
1187 if 2: both row_starts and col_starts given
1188 if 3: only col_starts given, row_starts = NULL */
1189 global_data[3] = 2;
1190 global_size = 2*num_procs+6;
1191 for (i=0; i < num_procs+1; i++)
1192 global_data[i+4] = row_starts[i];
1193 for (i=0; i < num_procs+1; i++)
1194 global_data[i+num_procs+5] = col_starts[i];
1195 }
1196 else
1197 {
1198 global_data[3] = 0;
1199 global_size = num_procs+5;
1200 for (i=0; i < num_procs+1; i++)
1201 global_data[i+4] = row_starts[i];
1202 }
1203 }
1204 else
1205 {
1206 global_data[3] = 1;
1207 global_size = num_procs+5;
1208 for (i=0; i < num_procs+1; i++)
1209 global_data[i+4] = row_starts[i];
1210 }
1211 }
1212 else
1213 {
1214 if (col_starts)
1215 {
1216 global_data[3] = 3;
1217 global_size = num_procs+5;
1218 for (i=0; i < num_procs+1; i++)
1219 global_data[i+4] = col_starts[i];
1220 }
1221 }
1222 global_data[0] = hypre_CSRMatrixNumRows(A);
1223 global_data[1] = hypre_CSRMatrixNumCols(A);
1224 global_data[2] = global_size;
1225 a_data = hypre_CSRMatrixData(A);
1226 a_i = hypre_CSRMatrixI(A);
1227 a_j = hypre_CSRMatrixJ(A);
1228 }
1229 MPI_Bcast(global_data,3,MPI_HYPRE_BIG_INT,0,comm);
1230 global_num_rows = global_data[0];
1231 global_num_cols = global_data[1];
1232
1233 global_size = global_data[2];
1234 if (global_size > 3)
1235 {
1236 MPI_Bcast(&global_data[3],global_size-3,MPI_INT,0,comm);
1237 if (my_id > 0)
1238 {
1239 if (global_data[3] < 3)
1240 {
1241 row_starts = hypre_CTAlloc(HYPRE_BigInt, num_procs+1);
1242 for (i=0; i< num_procs+1; i++)
1243 {
1244 row_starts[i] = global_data[i+4];
1245 }
1246 if (global_data[3] == 0)
1247 col_starts = row_starts;
1248 if (global_data[3] == 2)
1249 {
1250 col_starts = hypre_CTAlloc(HYPRE_BigInt, num_procs+1);
1251 for (i=0; i < num_procs+1; i++)
1252 {
1253 col_starts[i] = global_data[i+num_procs+5];
1254 }
1255 }
1256 }
1257 else
1258 {
1259 col_starts = hypre_CTAlloc(HYPRE_BigInt, num_procs+1);
1260 for (i=0; i< num_procs+1; i++)
1261 {
1262 col_starts[i] = global_data[i+4];
1263 }
1264 }
1265 }
1266 }
1267 hypre_TFree(global_data);
1268
1269 local_num_rows = hypre_CTAlloc(int, num_procs);
1270 csr_matrix_datatypes = hypre_CTAlloc(MPI_Datatype, num_procs);
1271
1272 par_matrix = hypre_ParCSRMatrixCreate (comm, global_num_rows,
1273 global_num_cols,row_starts,col_starts,0,0,0);
1274
1275 row_starts = hypre_ParCSRMatrixRowStarts(par_matrix);
1276 col_starts = hypre_ParCSRMatrixColStarts(par_matrix);
1277
1278 for (i=0; i < num_procs; i++)
1279 local_num_rows[i] = (int) (row_starts[i+1] - row_starts[i]);
1280
1281 if (my_id == 0)
1282 {
1283 local_num_nonzeros = hypre_CTAlloc(int, num_procs);
1284 for (i=0; i < num_procs-1; i++)
1285 local_num_nonzeros[i] = a_i[(int)row_starts[i+1]]
1286 - a_i[(int)row_starts[i]];
1287 local_num_nonzeros[num_procs-1] = a_i[(int)global_num_rows]
1288 - a_i[(int)row_starts[num_procs-1]];
1289 }
1290 MPI_Scatter(local_num_nonzeros,1,MPI_INT,&num_nonzeros,1,MPI_INT,0,comm);
1291
1292 if (my_id == 0) num_nonzeros = local_num_nonzeros[0];
1293
1294 local_A = hypre_CSRMatrixCreate(local_num_rows[my_id], global_num_cols,
1295 num_nonzeros);
1296 if (my_id == 0)
1297 {
1298 requests = hypre_CTAlloc (MPI_Request, num_procs-1);
1299 status = hypre_CTAlloc(MPI_Status, num_procs-1);
1300 j=0;
1301 for (i=1; i < num_procs; i++)
1302 {
1303 ind = a_i[(int)row_starts[i]];
1304 hypre_BuildCSRMatrixMPIDataType(local_num_nonzeros[i],
1305 local_num_rows[i],
1306 &a_data[ind],
1307 &a_i[(int)row_starts[i]],
1308 &a_j[ind],
1309 &csr_matrix_datatypes[i]);
1310 MPI_Isend(MPI_BOTTOM, 1, csr_matrix_datatypes[i], i, 0, comm,
1311 &requests[j++]);
1312 MPI_Type_free(&csr_matrix_datatypes[i]);
1313 }
1314 hypre_CSRMatrixData(local_A) = a_data;
1315 hypre_CSRMatrixI(local_A) = a_i;
1316 hypre_CSRMatrixJ(local_A) = a_j;
1317 hypre_CSRMatrixOwnsData(local_A) = 0;
1318 MPI_Waitall(num_procs-1,requests,status);
1319 hypre_TFree(requests);
1320 hypre_TFree(status);
1321 hypre_TFree(local_num_nonzeros);
1322 }
1323 else
1324 {
1325 hypre_CSRMatrixInitialize(local_A);
1326 hypre_BuildCSRMatrixMPIDataType(num_nonzeros,
1327 local_num_rows[my_id],
1328 hypre_CSRMatrixData(local_A),
1329 hypre_CSRMatrixI(local_A),
1330 hypre_CSRMatrixJ(local_A),
1331 csr_matrix_datatypes);
1332 MPI_Recv(MPI_BOTTOM,1,csr_matrix_datatypes[0],0,0,comm,&status0);
1333 MPI_Type_free(csr_matrix_datatypes);
1334 }
1335
1336 first_col_diag = col_starts[my_id];
1337 last_col_diag = col_starts[my_id+1]-1;
1338
1339 GenerateDiagAndOffd(local_A, par_matrix, first_col_diag, last_col_diag);
1340
1341 /* set pointers back to NULL before destroying */
1342 if (my_id == 0)
1343 {
1344 hypre_CSRMatrixData(local_A) = NULL;
1345 hypre_CSRMatrixI(local_A) = NULL;
1346 hypre_CSRMatrixJ(local_A) = NULL;
1347 }
1348 hypre_CSRMatrixDestroy(local_A);
1349 hypre_TFree(local_num_rows);
1350 hypre_TFree(csr_matrix_datatypes);
1351
1352 return par_matrix;
1353}
1354
1355
1356
1357int
1358GenerateDiagAndOffd(hypre_CSRMatrix *A,
1359 hypre_ParCSRMatrix *matrix,
1360 HYPRE_BigInt first_col_diag,
1361 HYPRE_BigInt last_col_diag)
1362{
1363 int i, j;
1364 int jo, jd;
1365 int num_rows = hypre_CSRMatrixNumRows(A);
1366 int num_cols = hypre_CSRMatrixNumCols(A);
1367 double *a_data = hypre_CSRMatrixData(A);
1368 int *a_i = hypre_CSRMatrixI(A);
1369 int *a_j = hypre_CSRMatrixJ(A);
1370
1371 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(matrix);
1372 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(matrix);
1373
1374 HYPRE_BigInt *col_map_offd;
1375
1376 double *diag_data, *offd_data;
1377 int *diag_i, *offd_i;
1378 int *diag_j, *offd_j;
1379 int *marker;
1380 int num_cols_diag, num_cols_offd;
1381 int first_elmt = a_i[0];
1382 int num_nonzeros = a_i[num_rows]-first_elmt;
1383 int counter;
1384
1385 num_cols_diag = (int) (last_col_diag - first_col_diag) +1;
1386 num_cols_offd = 0;
1387
1388 if (num_cols - num_cols_diag)
1389 {
1390 hypre_CSRMatrixInitialize(diag);
1391 diag_i = hypre_CSRMatrixI(diag);
1392
1393 hypre_CSRMatrixInitialize(offd);
1394 offd_i = hypre_CSRMatrixI(offd);
1395 marker = hypre_CTAlloc(int,num_cols);
1396
1397 for (i=0; i < num_cols; i++)
1398 marker[i] = 0;
1399
1400 jo = 0;
1401 jd = 0;
1402 for (i=0; i < num_rows; i++)
1403 {
1404 offd_i[i] = jo;
1405 diag_i[i] = jd;
1406
1407 for (j=a_i[i]-first_elmt; j < a_i[i+1]-first_elmt; j++)
1408 if (a_j[j] < (int)first_col_diag || a_j[j] > (int)last_col_diag)
1409 {
1410 if (!marker[a_j[j]])
1411 {
1412 marker[a_j[j]] = 1;
1413 num_cols_offd++;
1414 }
1415 jo++;
1416 }
1417 else
1418 {
1419 jd++;
1420 }
1421 }
1422 offd_i[num_rows] = jo;
1423 diag_i[num_rows] = jd;
1424
1425 hypre_ParCSRMatrixColMapOffd(matrix) = hypre_CTAlloc(HYPRE_BigInt,num_cols_offd);
1426 col_map_offd = hypre_ParCSRMatrixColMapOffd(matrix);
1427
1428 counter = 0;
1429 for (i=0; i < num_cols; i++)
1430 if (marker[i])
1431 {
1432 col_map_offd[counter] = (HYPRE_BigInt) i;
1433 marker[i] = counter;
1434 counter++;
1435 }
1436
1437 hypre_CSRMatrixNumNonzeros(diag) = jd;
1438 hypre_CSRMatrixInitialize(diag);
1439 diag_data = hypre_CSRMatrixData(diag);
1440 diag_j = hypre_CSRMatrixJ(diag);
1441
1442 hypre_CSRMatrixNumNonzeros(offd) = jo;
1443 hypre_CSRMatrixNumCols(offd) = num_cols_offd;
1444 hypre_CSRMatrixInitialize(offd);
1445 offd_data = hypre_CSRMatrixData(offd);
1446 offd_j = hypre_CSRMatrixJ(offd);
1447
1448 jo = 0;
1449 jd = 0;
1450 for (i=0; i < num_rows; i++)
1451 {
1452 for (j=a_i[i]-first_elmt; j < a_i[i+1]-first_elmt; j++)
1453 if (a_j[j] < (int)first_col_diag || a_j[j] > (int)last_col_diag)
1454 {
1455 offd_data[jo] = a_data[j];
1456 offd_j[jo++] = marker[a_j[j]];
1457 }
1458 else
1459 {
1460 diag_data[jd] = a_data[j];
1461 diag_j[jd++] = a_j[j]-(int)first_col_diag;
1462 }
1463 }
1464 hypre_TFree(marker);
1465 }
1466 else
1467 {
1468 hypre_CSRMatrixNumNonzeros(diag) = num_nonzeros;
1469 hypre_CSRMatrixInitialize(diag);
1470 diag_data = hypre_CSRMatrixData(diag);
1471 diag_i = hypre_CSRMatrixI(diag);
1472 diag_j = hypre_CSRMatrixJ(diag);
1473
1474 for (i=0; i < num_nonzeros; i++)
1475 {
1476 diag_data[i] = a_data[i];
1477 diag_j[i] = a_j[i];
1478 }
1479 offd_i = hypre_CTAlloc(int, num_rows+1);
1480
1481 for (i=0; i < num_rows+1; i++)
1482 {
1483 diag_i[i] = a_i[i];
1484 offd_i[i] = 0;
1485 }
1486
1487 hypre_CSRMatrixNumCols(offd) = 0;
1488 hypre_CSRMatrixI(offd) = offd_i;
1489 }
1490
1491 return hypre_error_flag;
1492}
1493
1494hypre_CSRMatrix *
1495hypre_MergeDiagAndOffd(hypre_ParCSRMatrix *par_matrix)
1496{
1497 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(par_matrix);
1498 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(par_matrix);
1499 hypre_CSRMatrix *matrix;
1500
1501 int num_cols = hypre_ParCSRMatrixGlobalNumCols(par_matrix);
1502 HYPRE_BigInt first_col_diag = hypre_ParCSRMatrixFirstColDiag(par_matrix);
1503 HYPRE_BigInt *col_map_offd = hypre_ParCSRMatrixColMapOffd(par_matrix);
1504 int num_rows = hypre_CSRMatrixNumRows(diag);
1505
1506 int *diag_i = hypre_CSRMatrixI(diag);
1507 int *diag_j = hypre_CSRMatrixJ(diag);
1508 double *diag_data = hypre_CSRMatrixData(diag);
1509 int *offd_i = hypre_CSRMatrixI(offd);
1510 int *offd_j = hypre_CSRMatrixJ(offd);
1511 double *offd_data = hypre_CSRMatrixData(offd);
1512
1513 int *matrix_i;
1514 int *matrix_j;
1515 double *matrix_data;
1516
1517 int num_nonzeros, i, j;
1518 int count;
1519
1520/*
1521 if (!num_cols_offd)
1522 {
1523 matrix = hypre_CSRMatrixCreate(num_rows,num_cols,diag_i[num_rows]);
1524 hypre_CSRMatrixOwnsData(matrix) = 0;
1525 hypre_CSRMatrixI(matrix) = diag_i;
1526 hypre_CSRMatrixJ(matrix) = diag_j;
1527 hypre_CSRMatrixData(matrix) = diag_data;
1528 return matrix;
1529 }
1530*/
1531
1532 num_nonzeros = diag_i[num_rows] + offd_i[num_rows];
1533
1534 matrix = hypre_CSRMatrixCreate(num_rows,num_cols,num_nonzeros);
1535 hypre_CSRMatrixInitialize(matrix);
1536
1537 matrix_i = hypre_CSRMatrixI(matrix);
1538 matrix_j = hypre_CSRMatrixJ(matrix);
1539 matrix_data = hypre_CSRMatrixData(matrix);
1540
1541 count = 0;
1542 matrix_i[0] = 0;
1543 for (i=0; i < num_rows; i++)
1544 {
1545 for (j=diag_i[i]; j < diag_i[i+1]; j++)
1546 {
1547 matrix_data[count] = diag_data[j];
1548 matrix_j[count++] = diag_j[j]+(int)first_col_diag;
1549 }
1550 for (j=offd_i[i]; j < offd_i[i+1]; j++)
1551 {
1552 matrix_data[count] = offd_data[j];
1553 matrix_j[count++] = (int)col_map_offd[offd_j[j]];
1554 }
1555 matrix_i[i+1] = count;
1556 }
1557
1558 return matrix;
1559}
1560
1561/*--------------------------------------------------------------------------
1562 * hypre_ParCSRMatrixToCSRMatrixAll:
1563 * generates a CSRMatrix from a ParCSRMatrix on all processors that have
1564 * parts of the ParCSRMatrix
1565 *--------------------------------------------------------------------------*/
1566
1567hypre_CSRMatrix *
1568hypre_ParCSRMatrixToCSRMatrixAll(hypre_ParCSRMatrix *par_matrix)
1569{
1570 MPI_Comm comm = hypre_ParCSRMatrixComm(par_matrix);
1571 hypre_CSRMatrix *matrix;
1572 hypre_CSRMatrix *local_matrix;
1573 HYPRE_BigInt num_rows = hypre_ParCSRMatrixGlobalNumRows(par_matrix);
1574 HYPRE_BigInt num_cols = hypre_ParCSRMatrixGlobalNumCols(par_matrix);
1575#ifndef HYPRE_NO_GLOBAL_PARTITION
1576 HYPRE_BigInt *row_starts = hypre_ParCSRMatrixRowStarts(par_matrix);
1577#endif
1578 int *matrix_i;
1579 int *matrix_j;
1580 double *matrix_data;
1581
1582 int *local_matrix_i;
1583 int *local_matrix_j;
1584 double *local_matrix_data;
1585
1586 int i, j;
1587 int local_num_rows;
1588 int local_num_nonzeros;
1589 int num_nonzeros;
1590 int num_data;
1591 int num_requests;
1592 int vec_len, offset;
1593 int start_index;
1594 int proc_id;
1595 int num_procs, my_id;
1596 int num_types;
1597 int *used_procs;
1598
1599 MPI_Request *requests;
1600 MPI_Status *status;
1601 /* MPI_Datatype *data_type; */
1602
1603
1604#ifdef HYPRE_NO_GLOBAL_PARTITION
1605
1606 HYPRE_BigInt *new_vec_starts;
1607 int lower_j, upper_j;
1608
1609 int num_contacts;
1610 int contact_proc_list[1];
1611 HYPRE_BigInt contact_send_buf[1];
1612 int contact_send_buf_starts[2];
1613 int max_response_size;
1614 int *response_recv_buf=NULL;
1615 int *response_recv_buf_starts = NULL;
1616 hypre_DataExchangeResponse response_obj;
1617 hypre_ProcListElements send_proc_obj;
1618
1619 HYPRE_BigInt *send_info = NULL;
1620 MPI_Status status1;
1621 int count, tag1 = 11112, tag2 = 22223, tag3 = 33334;
1622 int start;
1623
1624#endif
1625
1626
1627
1628 MPI_Comm_size(comm, &num_procs);
1629 MPI_Comm_rank(comm, &my_id);
1630
1631#ifdef HYPRE_NO_GLOBAL_PARTITION
1632
1633 local_num_rows = (int)(hypre_ParCSRMatrixLastRowIndex(par_matrix) -
1634 hypre_ParCSRMatrixFirstRowIndex(par_matrix)) + 1;
1635
1636
1637 local_matrix = hypre_MergeDiagAndOffd(par_matrix); /* creates matrix */
1638 local_matrix_i = hypre_CSRMatrixI(local_matrix);
1639 local_matrix_j = hypre_CSRMatrixJ(local_matrix);
1640 local_matrix_data = hypre_CSRMatrixData(local_matrix);
1641
1642
1643/* determine procs that have vector data and store their ids in used_procs */
1644/* we need to do an exchange data for this. If I own row then I will contact
1645 processor 0 with the endpoint of my local range */
1646
1647 if (local_num_rows > 0)
1648 {
1649 num_contacts = 1;
1650 contact_proc_list[0] = 0;
1651 contact_send_buf[0] = hypre_ParCSRMatrixLastRowIndex(par_matrix);
1652 contact_send_buf_starts[0] = 0;
1653 contact_send_buf_starts[1] = 1;
1654
1655 }
1656 else
1657 {
1658 num_contacts = 0;
1659 contact_send_buf_starts[0] = 0;
1660 contact_send_buf_starts[1] = 0;
1661 }
1662 /*build the response object*/
1663 /*send_proc_obj will be for saving info from contacts */
1664 send_proc_obj.length = 0;
1665 send_proc_obj.storage_length = 10;
1666 send_proc_obj.id = hypre_CTAlloc(int, send_proc_obj.storage_length);
1667 send_proc_obj.vec_starts = hypre_CTAlloc(int, send_proc_obj.storage_length + 1);
1668 send_proc_obj.vec_starts[0] = 0;
1669 send_proc_obj.element_storage_length = 10;
1670 send_proc_obj.elements = hypre_CTAlloc(HYPRE_BigInt, send_proc_obj.element_storage_length);
1671
1672 max_response_size = 0; /* each response is null */
1673 response_obj.fill_response = hypre_FillResponseParToCSRMatrix;
1674 response_obj.data1 = NULL;
1675 response_obj.data2 = &send_proc_obj; /*this is where we keep info from contacts*/
1676
1677
1678 hypre_DataExchangeList(num_contacts,
1679 contact_proc_list, contact_send_buf,
1680 contact_send_buf_starts, sizeof(HYPRE_BigInt),
1681 0, &response_obj,
1682 max_response_size, 1,
1683 comm, (void**) &response_recv_buf,
1684 &response_recv_buf_starts);
1685
1686
1687
1688
1689 /* now processor 0 should have a list of ranges for processors that have rows -
1690 these are in send_proc_obj - it needs to create the new list of processors
1691 and also an array of vec starts - and send to those who own row*/
1692 if (my_id)
1693 {
1694 if (local_num_rows)
1695 {
1696 /* look for a message from processor 0 */
1697 MPI_Probe(0, tag1, comm, &status1);
1698 MPI_Get_count(&status1, MPI_HYPRE_BIG_INT, &count);
1699
1700 send_info = hypre_CTAlloc(HYPRE_BigInt, count);
1701 MPI_Recv(send_info, count, MPI_HYPRE_BIG_INT, 0, tag1, comm, &status1);
1702
1703 /* now unpack */
1704 num_types = (int) send_info[0];
1705 used_procs = hypre_CTAlloc(int, num_types);
1706 new_vec_starts = hypre_CTAlloc(HYPRE_BigInt, num_types+1);
1707
1708 for (i=1; i<= num_types; i++)
1709 {
1710 used_procs[i-1] = (int) send_info[i];
1711 }
1712 for (i=num_types+1; i< count; i++)
1713 {
1714 new_vec_starts[i-num_types-1] = send_info[i] ;
1715 }
1716 }
1717 else /* clean up and exit */
1718 {
1719 hypre_TFree(send_proc_obj.vec_starts);
1720 hypre_TFree(send_proc_obj.id);
1721 hypre_TFree(send_proc_obj.elements);
1722 if(response_recv_buf) hypre_TFree(response_recv_buf);
1723 if(response_recv_buf_starts) hypre_TFree(response_recv_buf_starts);
1724
1725
1726 if (hypre_CSRMatrixOwnsData(local_matrix))
1727 hypre_CSRMatrixDestroy(local_matrix);
1728 else
1729 hypre_TFree(local_matrix);
1730
1731
1732 return NULL;
1733 }
1734 }
1735 else /* my_id ==0 */
1736 {
1737 num_types = send_proc_obj.length;
1738 used_procs = hypre_CTAlloc(int, num_types);
1739 new_vec_starts = hypre_CTAlloc(HYPRE_BigInt, num_types+1);
1740
1741 new_vec_starts[0] = 0;
1742 for (i=0; i< num_types; i++)
1743 {
1744 used_procs[i] = send_proc_obj.id[i];
1745 new_vec_starts[i+1] = send_proc_obj.elements[i]+1;
1746 }
1747 qsort0(used_procs, 0, num_types-1);
1748
1749 hypre_BigQsort0(new_vec_starts, 0, num_types);
1750
1751 /*now we need to put into an array to send */
1752 count = 2*num_types+2;
1753 send_info = hypre_CTAlloc(HYPRE_BigInt, count);
1754 send_info[0] = (HYPRE_BigInt) num_types;
1755 for (i=1; i<= num_types; i++)
1756 {
1757 send_info[i] = (HYPRE_BigInt) used_procs[i-1];
1758 }
1759 for (i=num_types+1; i< count; i++)
1760 {
1761 send_info[i] = new_vec_starts[i-num_types-1];
1762 }
1763 requests = hypre_CTAlloc(MPI_Request, num_types);
1764 status = hypre_CTAlloc(MPI_Status, num_types);
1765
1766 /* don't send to myself - these are sorted so my id would be first*/
1767 start = 0;
1768 if (used_procs[0] == 0)
1769 {
1770 start = 1;
1771 }
1772
1773
1774 for (i=start; i < num_types; i++)
1775 {
1776 MPI_Isend(send_info, count, MPI_HYPRE_BIG_INT, used_procs[i], tag1, comm, &requests[i-start]);
1777 }
1778 MPI_Waitall(num_types-start, requests, status);
1779
1780 hypre_TFree(status);
1781 hypre_TFree(requests);
1782
1783
1784 }
1785 /* clean up */
1786 hypre_TFree(send_proc_obj.vec_starts);
1787 hypre_TFree(send_proc_obj.id);
1788 hypre_TFree(send_proc_obj.elements);
1789 hypre_TFree(send_info);
1790 if(response_recv_buf) hypre_TFree(response_recv_buf);
1791 if(response_recv_buf_starts) hypre_TFree(response_recv_buf_starts);
1792
1793 /* now proc 0 can exit if it has no rows */
1794 if (!local_num_rows)
1795 {
1796 if (hypre_CSRMatrixOwnsData(local_matrix))
1797 hypre_CSRMatrixDestroy(local_matrix);
1798 else
1799 hypre_TFree(local_matrix);
1800
1801 hypre_TFree(new_vec_starts);
1802 hypre_TFree(used_procs);
1803
1804 return NULL;
1805 }
1806
1807
1808 /* everyone left has rows and knows: new_vec_starts, num_types, and used_procs */
1809
1810 /* this matrix should be rather small */
1811 matrix_i = hypre_CTAlloc(int, num_rows+1);
1812
1813 num_requests = 4*num_types;
1814 requests = hypre_CTAlloc(MPI_Request, num_requests);
1815 status = hypre_CTAlloc(MPI_Status, num_requests);
1816
1817
1818 /* exchange contents of local_matrix_i - here we are sending to ourself also*/
1819
1820 j = 0;
1821 for (i = 0; i < num_types; i++)
1822 {
1823 proc_id = used_procs[i];
1824 vec_len = (int) (new_vec_starts[i+1] - new_vec_starts[i]);
1825 MPI_Irecv(&matrix_i[( (int) new_vec_starts[i]) + 1], vec_len, MPI_INT,
1826 proc_id, tag2, comm, &requests[j++]);
1827 }
1828 for (i = 0; i < num_types; i++)
1829 {
1830 proc_id = used_procs[i];
1831 MPI_Isend(&local_matrix_i[1], local_num_rows, MPI_INT,
1832 proc_id, tag2, comm, &requests[j++]);
1833 }
1834
1835 MPI_Waitall(j, requests, status);
1836
1837
1838
1839/* generate matrix_i from received data */
1840/* global numbering?*/
1841 offset = matrix_i[( (int) new_vec_starts[1])];
1842 for (i=1; i < num_types; i++)
1843 {
1844 lower_j = (int) new_vec_starts[i];
1845 upper_j = (int) new_vec_starts[i+1];
1846 for (j = lower_j; j < upper_j; j++)
1847 matrix_i[j+1] += offset;
1848 offset = matrix_i[upper_j];
1849 }
1850
1851 num_nonzeros = matrix_i[num_rows];
1852
1853 matrix = hypre_CSRMatrixCreate((int)num_rows, (int)num_cols, num_nonzeros);
1854 hypre_CSRMatrixI(matrix) = matrix_i;
1855 hypre_CSRMatrixInitialize(matrix);
1856 matrix_j = hypre_CSRMatrixJ(matrix);
1857 matrix_data = hypre_CSRMatrixData(matrix);
1858
1859/* generate datatypes for further data exchange and exchange remaining
1860 data, i.e. column info and actual data */
1861
1862 j = 0;
1863 for (i = 0; i < num_types; i++)
1864 {
1865 proc_id = used_procs[i];
1866 start_index = matrix_i[(int) new_vec_starts[i]];
1867 num_data = matrix_i[(int) new_vec_starts[i+1]] - start_index;
1868 MPI_Irecv(&matrix_data[start_index], num_data, MPI_DOUBLE,
1869 used_procs[i], tag1, comm, &requests[j++]);
1870 MPI_Irecv(&matrix_j[start_index], num_data, MPI_INT,
1871 used_procs[i], tag3, comm, &requests[j++]);
1872 }
1873 local_num_nonzeros = local_matrix_i[local_num_rows];
1874 for (i=0; i < num_types; i++)
1875 {
1876 MPI_Isend(local_matrix_data, local_num_nonzeros, MPI_DOUBLE,
1877 used_procs[i], tag1, comm, &requests[j++]);
1878 MPI_Isend(local_matrix_j, local_num_nonzeros, MPI_INT,
1879 used_procs[i], tag3, comm, &requests[j++]);
1880 }
1881
1882
1883 MPI_Waitall(num_requests, requests, status);
1884
1885 hypre_TFree(new_vec_starts);
1886
1887#else
1888 local_num_rows = (int)(row_starts[my_id+1] - row_starts[my_id]);
1889
1890/* if my_id contains no data, return NULL */
1891
1892 if (!local_num_rows)
1893 return NULL;
1894
1895 local_matrix = hypre_MergeDiagAndOffd(par_matrix);
1896 local_matrix_i = hypre_CSRMatrixI(local_matrix);
1897 local_matrix_j = hypre_CSRMatrixJ(local_matrix);
1898 local_matrix_data = hypre_CSRMatrixData(local_matrix);
1899
1900
1901 matrix_i = hypre_CTAlloc(int, num_rows+1);
1902
1903
1904/* determine procs that have vector data and store their ids in used_procs */
1905
1906 num_types = 0;
1907 for (i=0; i < num_procs; i++)
1908 if (row_starts[i+1]-row_starts[i] && i-my_id)
1909 num_types++;
1910 num_requests = 4*num_types;
1911
1912 used_procs = hypre_CTAlloc(int, num_types);
1913 j = 0;
1914 for (i=0; i < num_procs; i++)
1915 if (row_starts[i+1]-row_starts[i] && i-my_id)
1916 used_procs[j++] = i;
1917
1918 requests = hypre_CTAlloc(MPI_Request, num_requests);
1919 status = hypre_CTAlloc(MPI_Status, num_requests);
1920 /* data_type = hypre_CTAlloc(MPI_Datatype, num_types+1); */
1921
1922/* exchange contents of local_matrix_i */
1923
1924 j = 0;
1925 for (i = 0; i < num_types; i++)
1926 {
1927 proc_id = used_procs[i];
1928 vec_len = (int)(row_starts[proc_id+1] - row_starts[proc_id]);
1929 MPI_Irecv(&matrix_i[(int)row_starts[proc_id]+1], vec_len, MPI_INT,
1930 proc_id, 0, comm, &requests[j++]);
1931 }
1932 for (i = 0; i < num_types; i++)
1933 {
1934 proc_id = used_procs[i];
1935 MPI_Isend(&local_matrix_i[1], local_num_rows, MPI_INT,
1936 proc_id, 0, comm, &requests[j++]);
1937 }
1938
1939 vec_len = (int)(row_starts[my_id+1] - row_starts[my_id]);
1940 for (i=1; i <= vec_len; i++)
1941 matrix_i[(int)row_starts[my_id]+i] = local_matrix_i[i];
1942
1943 MPI_Waitall(j, requests, status);
1944
1945/* generate matrix_i from received data */
1946
1947 offset = matrix_i[(int)row_starts[1]];
1948 for (i=1; i < num_procs; i++)
1949 {
1950 for (j = (int)row_starts[i]; j < (int)row_starts[i+1]; j++)
1951 matrix_i[j+1] += offset;
1952 offset = matrix_i[(int)row_starts[i+1]];
1953 }
1954
1955 num_nonzeros = matrix_i[num_rows];
1956
1957 matrix = hypre_CSRMatrixCreate(num_rows, num_cols, num_nonzeros);
1958 hypre_CSRMatrixI(matrix) = matrix_i;
1959 hypre_CSRMatrixInitialize(matrix);
1960 matrix_j = hypre_CSRMatrixJ(matrix);
1961 matrix_data = hypre_CSRMatrixData(matrix);
1962
1963/* generate datatypes for further data exchange and exchange remaining
1964 data, i.e. column info and actual data */
1965
1966 j = 0;
1967 for (i = 0; i < num_types; i++)
1968 {
1969 proc_id = used_procs[i];
1970 start_index = matrix_i[(int)row_starts[proc_id]];
1971 num_data = matrix_i[(int)row_starts[proc_id+1]] - start_index;
1972 MPI_Irecv(&matrix_data[start_index], num_data, MPI_DOUBLE,
1973 used_procs[i], 0, comm, &requests[j++]);
1974 MPI_Irecv(&matrix_j[start_index], num_data, MPI_INT,
1975 used_procs[i], 0, comm, &requests[j++]);
1976 }
1977 local_num_nonzeros = local_matrix_i[local_num_rows];
1978 for (i=0; i < num_types; i++)
1979 {
1980 MPI_Isend(local_matrix_data, local_num_nonzeros, MPI_DOUBLE,
1981 used_procs[i], 0, comm, &requests[j++]);
1982 MPI_Isend(local_matrix_j, local_num_nonzeros, MPI_INT,
1983 used_procs[i], 0, comm, &requests[j++]);
1984 }
1985
1986 start_index = matrix_i[row_starts[my_id]];
1987 for (i=0; i < local_num_nonzeros; i++)
1988 {
1989 matrix_j[start_index+i] = local_matrix_j[i];
1990 matrix_data[start_index+i] = local_matrix_data[i];
1991 }
1992 MPI_Waitall(num_requests, requests, status);
1993/* for (i = 0; i < num_types; i++)
1994 {
1995 proc_id = used_procs[i];
1996 start_index = matrix_i[row_starts[proc_id]];
1997 num_data = matrix_i[row_starts[proc_id+1]] - start_index;
1998 hypre_BuildCSRJDataType(num_data,
1999 &matrix_data[start_index],
2000 &matrix_j[start_index],
2001 &data_type[i]);
2002 }
2003 local_num_nonzeros = local_matrix_i[local_num_rows];
2004 hypre_BuildCSRJDataType(local_num_nonzeros,
2005 local_matrix_data,
2006 local_matrix_j,
2007 &data_type[num_types]);
2008 j = 0;
2009 for (i=0; i < num_types; i++)
2010 MPI_Irecv(MPI_BOTTOM, 1, data_type[i],
2011 used_procs[i], 0, comm, &requests[j++]);
2012 for (i=0; i < num_types; i++)
2013 MPI_Isend(MPI_BOTTOM, 1, data_type[num_types],
2014 used_procs[i], 0, comm, &requests[j++]);
2015*/
2016 start_index = matrix_i[(int)row_starts[my_id]];
2017 for (i=0; i < local_num_nonzeros; i++)
2018 {
2019 matrix_j[start_index+i] = local_matrix_j[i];
2020 matrix_data[start_index+i] = local_matrix_data[i];
2021 }
2022 MPI_Waitall(num_requests, requests, status);
2023/*
2024 for (i=0; i <= num_types; i++)
2025 MPI_Type_free(&data_type[i]);
2026*/
2027
2028#endif
2029
2030 if (hypre_CSRMatrixOwnsData(local_matrix))
2031 hypre_CSRMatrixDestroy(local_matrix);
2032 else
2033 hypre_TFree(local_matrix);
2034
2035 if (num_requests)
2036 {
2037 hypre_TFree(requests);
2038 hypre_TFree(status);
2039 hypre_TFree(used_procs);
2040 }
2041 /* hypre_TFree(data_type); */
2042
2043 return matrix;
2044}
2045
2046/*--------------------------------------------------------------------------
2047 * hypre_ParCSRMatrixCopy,
2048 * copies B to A,
2049 * if copy_data = 0, only the structure of A is copied to B
2050 * the routine does not check whether the dimensions of A and B are compatible
2051 *--------------------------------------------------------------------------*/
2052
2053int
2054hypre_ParCSRMatrixCopy( hypre_ParCSRMatrix *A, hypre_ParCSRMatrix *B,
2055 int copy_data )
2056{
2057 hypre_CSRMatrix *A_diag;
2058 hypre_CSRMatrix *A_offd;
2059 HYPRE_BigInt *col_map_offd_A;
2060 hypre_CSRMatrix *B_diag;
2061 hypre_CSRMatrix *B_offd;
2062 HYPRE_BigInt *col_map_offd_B;
2063 int num_cols_offd;
2064 int i;
2065
2066 if (!A)
2067 {
2068 hypre_error_in_arg(1);
2069 return hypre_error_flag;
2070 }
2071 if (!B)
2072 {
2073 hypre_error_in_arg(1);
2074 return hypre_error_flag;
2075 }
2076 A_diag = hypre_ParCSRMatrixDiag(A);
2077 A_offd = hypre_ParCSRMatrixOffd(A);
2078 col_map_offd_A = hypre_ParCSRMatrixColMapOffd(A);
2079 B_diag = hypre_ParCSRMatrixDiag(B);
2080 B_offd = hypre_ParCSRMatrixOffd(B);
2081 col_map_offd_B = hypre_ParCSRMatrixColMapOffd(B);
2082 num_cols_offd = hypre_CSRMatrixNumCols(A_offd);
2083
2084 hypre_CSRMatrixCopy(A_diag, B_diag, copy_data);
2085 hypre_CSRMatrixCopy(A_offd, B_offd, copy_data);
2086 if (num_cols_offd && col_map_offd_B == NULL)
2087 {
2088 col_map_offd_B = hypre_CTAlloc(HYPRE_BigInt,num_cols_offd);
2089 hypre_ParCSRMatrixColMapOffd(B) = col_map_offd_B;
2090 }
2091 for (i = 0; i < num_cols_offd; i++)
2092 col_map_offd_B[i] = col_map_offd_A[i];
2093
2094 return hypre_error_flag;
2095}
2096/*--------------------------------------------------------------------
2097 * hypre_FillResponseParToCSRMatrix
2098 * Fill response function for determining the send processors
2099 * data exchange
2100 *--------------------------------------------------------------------*/
2101
2102int
2103hypre_FillResponseParToCSRMatrix(void *p_recv_contact_buf,
2104 int contact_size, int contact_proc, void *ro,
2105 MPI_Comm comm, void **p_send_response_buf,
2106 int *response_message_size )
2107{
2108 int myid;
2109 int i, index, count, elength;
2110
2111 HYPRE_BigInt *recv_contact_buf = (HYPRE_BigInt * ) p_recv_contact_buf;
2112
2113 hypre_DataExchangeResponse *response_obj = ro;
2114
2115 hypre_ProcListElements *send_proc_obj = response_obj->data2;
2116
2117
2118 MPI_Comm_rank(comm, &myid );
2119
2120
2121 /*check to see if we need to allocate more space in send_proc_obj for ids*/
2122 if (send_proc_obj->length == send_proc_obj->storage_length)
2123 {
2124 send_proc_obj->storage_length +=10; /*add space for 10 more processors*/
2125 send_proc_obj->id = hypre_TReAlloc(send_proc_obj->id,int,
2126 send_proc_obj->storage_length);
2127 send_proc_obj->vec_starts = hypre_TReAlloc(send_proc_obj->vec_starts,int,
2128 send_proc_obj->storage_length + 1);
2129 }
2130
2131 /*initialize*/
2132 count = send_proc_obj->length;
2133 index = send_proc_obj->vec_starts[count]; /*this is the number of elements*/
2134
2135 /*send proc*/
2136 send_proc_obj->id[count] = contact_proc;
2137
2138 /*do we need more storage for the elements?*/
2139 if (send_proc_obj->element_storage_length < index + contact_size)
2140 {
2141 elength = hypre_max(contact_size, 10);
2142 elength += index;
2143 send_proc_obj->elements = hypre_TReAlloc(send_proc_obj->elements,
2144 HYPRE_BigInt, elength);
2145 send_proc_obj->element_storage_length = elength;
2146 }
2147 /*populate send_proc_obj*/
2148 for (i=0; i< contact_size; i++)
2149 {
2150 send_proc_obj->elements[index++] = recv_contact_buf[i];
2151 }
2152 send_proc_obj->vec_starts[count+1] = index;
2153 send_proc_obj->length++;
2154
2155
2156 /*output - no message to return (confirmation) */
2157 *response_message_size = 0;
2158
2159
2160 return hypre_error_flag;
2161
2162}
2163
2164/*--------------------------------------------------------------------------
2165 * hypre_ParCSRMatrixCompleteClone
2166 * Creates and returns a new copy of the argument, A.
2167 * Data is not copied, only structural information is reproduced.
2168 * The following variables are not copied because they will be constructed
2169 * later if needed: CommPkg, CommPkgT, rowindices, rowvalues
2170 *--------------------------------------------------------------------------*/
2171/* This differs from Hypre_ParCSRMatrixClone in parcsr_ls/par_gsmg.c, because
2172 that Clone function makes a matrix with different global parameters. */
2173
2174hypre_ParCSRMatrix * hypre_ParCSRMatrixCompleteClone( hypre_ParCSRMatrix * A )
2175{
2176 hypre_ParCSRMatrix * B = hypre_CTAlloc(hypre_ParCSRMatrix, 1);
2177 int i, ncols_offd;
2178
2179 hypre_ParCSRMatrixComm( B ) = hypre_ParCSRMatrixComm( A );
2180 hypre_ParCSRMatrixGlobalNumRows( B ) = hypre_ParCSRMatrixGlobalNumRows( A );
2181 hypre_ParCSRMatrixGlobalNumCols( B ) = hypre_ParCSRMatrixGlobalNumCols( A );
2182 hypre_ParCSRMatrixFirstRowIndex( B ) = hypre_ParCSRMatrixFirstRowIndex( A );
2183 hypre_ParCSRMatrixFirstColDiag( B ) = hypre_ParCSRMatrixFirstColDiag( A );
2184 hypre_ParCSRMatrixLastRowIndex( B ) = hypre_ParCSRMatrixLastRowIndex( A );
2185 hypre_ParCSRMatrixLastColDiag( B ) = hypre_ParCSRMatrixLastColDiag( A );
2186 hypre_ParCSRMatrixDiag( B ) = hypre_CSRMatrixClone( hypre_ParCSRMatrixDiag( A ) );
2187 hypre_ParCSRMatrixOffd( B ) = hypre_CSRMatrixClone( hypre_ParCSRMatrixOffd( A ) );
2188 hypre_ParCSRMatrixRowStarts( B ) = hypre_ParCSRMatrixRowStarts( A );
2189 hypre_ParCSRMatrixColStarts( B ) = hypre_ParCSRMatrixColStarts( A );
2190 /* note that B doesn't own row_starts and col_starts, this isn't a fully deep copy */
2191 hypre_ParCSRMatrixCommPkg( B ) = NULL;
2192 hypre_ParCSRMatrixCommPkgT( B ) = NULL;
2193 hypre_ParCSRMatrixOwnsData( B ) = 1;
2194 hypre_ParCSRMatrixOwnsRowStarts( B ) = 0;
2195 hypre_ParCSRMatrixOwnsColStarts( B ) = 0;
2196 hypre_ParCSRMatrixNumNonzeros( B ) = hypre_ParCSRMatrixNumNonzeros( A );
2197 hypre_ParCSRMatrixDNumNonzeros( B ) = hypre_ParCSRMatrixNumNonzeros( A );
2198 hypre_ParCSRMatrixRowindices( B ) = NULL;
2199 hypre_ParCSRMatrixRowvalues( B ) = NULL;
2200 hypre_ParCSRMatrixGetrowactive( B ) = 0;
2201 ncols_offd = hypre_CSRMatrixNumCols( hypre_ParCSRMatrixOffd( B ) );
2202
2203 hypre_ParCSRMatrixColMapOffd( B ) = hypre_CTAlloc( HYPRE_BigInt, ncols_offd );
2204 for ( i=0; i<ncols_offd; ++i )
2205 hypre_ParCSRMatrixColMapOffd( B )[i] = hypre_ParCSRMatrixColMapOffd( A )[i];
2206
2207 return B;
2208}
2209
2210/*--------------------------------------------------------------------------
2211 * hypre_ParCSRMatrixUnion
2212 * Creates and returns a new matrix whose elements are the union of A and B.
2213 * Data is not copied, only structural information is created.
2214 * A and B must have the same communicator, numbers and distributions of rows
2215 * and columns (they can differ in which row-column pairs are nonzero, thus
2216 * in which columns are in a offd block)
2217 *--------------------------------------------------------------------------*/
2218
2219hypre_ParCSRMatrix * hypre_ParCSRMatrixUnion( hypre_ParCSRMatrix * A,
2220 hypre_ParCSRMatrix * B )
2221{
2222 hypre_ParCSRMatrix * C;
2223 HYPRE_BigInt * col_map_offd_C = NULL;
2224 int num_procs, my_id, p;
2225 MPI_Comm comm = hypre_ParCSRMatrixComm( A );
2226
2227 MPI_Comm_rank(comm,&my_id);
2228 MPI_Comm_size(comm,&num_procs);
2229
2230 C = hypre_CTAlloc( hypre_ParCSRMatrix, 1 );
2231 hypre_ParCSRMatrixComm( C ) = hypre_ParCSRMatrixComm( A );
2232 hypre_ParCSRMatrixGlobalNumRows( C ) = hypre_ParCSRMatrixGlobalNumRows( A );
2233 hypre_ParCSRMatrixGlobalNumCols( C ) = hypre_ParCSRMatrixGlobalNumCols( A );
2234 hypre_ParCSRMatrixFirstRowIndex( C ) = hypre_ParCSRMatrixFirstRowIndex( A );
2235 hypre_assert( hypre_ParCSRMatrixFirstRowIndex( B )
2236 == hypre_ParCSRMatrixFirstRowIndex( A ) );
2237 hypre_ParCSRMatrixRowStarts( C ) = hypre_ParCSRMatrixRowStarts( A );
2238 hypre_ParCSRMatrixOwnsRowStarts( C ) = 0;
2239 /* hypre_ParCSRMatrixColStarts( C ) = hypre_CTAlloc( int, num_procs+1 ); */
2240 hypre_ParCSRMatrixColStarts( C ) = hypre_ParCSRMatrixColStarts( A );
2241 hypre_ParCSRMatrixOwnsColStarts( C ) = 0;
2242 for ( p=0; p<=num_procs; ++p )
2243 hypre_assert( hypre_ParCSRMatrixColStarts(A)
2244 == hypre_ParCSRMatrixColStarts(B) );
2245 hypre_ParCSRMatrixFirstColDiag( C ) = hypre_ParCSRMatrixFirstColDiag( A );
2246 hypre_ParCSRMatrixLastRowIndex( C ) = hypre_ParCSRMatrixLastRowIndex( A );
2247 hypre_ParCSRMatrixLastColDiag( C ) = hypre_ParCSRMatrixLastColDiag( A );
2248
2249 hypre_ParCSRMatrixDiag( C ) =
2250 hypre_CSRMatrixUnion( hypre_ParCSRMatrixDiag(A), hypre_ParCSRMatrixDiag(B),
2251 0, 0, 0 );
2252 hypre_ParCSRMatrixOffd( C ) =
2253 hypre_CSRMatrixUnion( hypre_ParCSRMatrixOffd(A), hypre_ParCSRMatrixOffd(B),
2254 hypre_ParCSRMatrixColMapOffd(A),
2255 hypre_ParCSRMatrixColMapOffd(B), &col_map_offd_C );
2256 hypre_ParCSRMatrixColMapOffd( C ) = col_map_offd_C;
2257 hypre_ParCSRMatrixCommPkg( C ) = NULL;
2258 hypre_ParCSRMatrixCommPkgT( C ) = NULL;
2259 hypre_ParCSRMatrixOwnsData( C ) = 1;
2260 /* SetNumNonzeros, SetDNumNonzeros are global, need MPI_Allreduce.
2261 I suspect, but don't know, that other parts of hypre do not assume that
2262 the correct values have been set.
2263 hypre_ParCSRMatrixSetNumNonzeros( C );
2264 hypre_ParCSRMatrixSetDNumNonzeros( C );*/
2265 hypre_ParCSRMatrixNumNonzeros( C ) = 0;
2266 hypre_ParCSRMatrixDNumNonzeros( C ) = 0.0;
2267 hypre_ParCSRMatrixRowindices( C ) = NULL;
2268 hypre_ParCSRMatrixRowvalues( C ) = NULL;
2269 hypre_ParCSRMatrixGetrowactive( C ) = 0;
2270
2271 return C;
2272}
Note: See TracBrowser for help on using the repository browser.