source: CIVL/examples/mpi-omp/AMG2013/IJ_mv/IJMatrix_parcsr.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: 88.0 KB
Line 
1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14
15/******************************************************************************
16 *
17 * IJMatrix_ParCSR interface
18 *
19 *****************************************************************************/
20
21#include "headers.h"
22
23#include "../HYPRE.h"
24
25/******************************************************************************
26 *
27 * hypre_IJMatrixCreateParCSR
28 *
29 *****************************************************************************/
30
31int
32hypre_IJMatrixCreateParCSR(hypre_IJMatrix *matrix)
33{
34 MPI_Comm comm = hypre_IJMatrixComm(matrix);
35 HYPRE_BigInt *row_partitioning = hypre_IJMatrixRowPartitioning(matrix);
36 HYPRE_BigInt *col_partitioning = hypre_IJMatrixColPartitioning(matrix);
37 hypre_ParCSRMatrix *par_matrix;
38 HYPRE_BigInt *row_starts;
39 HYPRE_BigInt *col_starts;
40 int num_procs;
41 int i;
42
43 MPI_Comm_size(comm,&num_procs);
44
45#ifdef HYPRE_NO_GLOBAL_PARTITION
46 row_starts = hypre_CTAlloc(HYPRE_BigInt,2);
47 if (hypre_IJMatrixGlobalFirstRow(matrix))
48 for (i=0; i < 2; i++)
49 row_starts[i] = row_partitioning[i]- hypre_IJMatrixGlobalFirstRow(matrix);
50 else
51 for (i=0; i < 2; i++)
52 row_starts[i] = row_partitioning[i];
53 if (row_partitioning != col_partitioning)
54 {
55 col_starts = hypre_CTAlloc(HYPRE_BigInt,2);
56 if (hypre_IJMatrixGlobalFirstCol(matrix))
57 for (i=0; i < 2; i++)
58 col_starts[i] = col_partitioning[i]-hypre_IJMatrixGlobalFirstCol(matrix);
59 else
60 for (i=0; i < 2; i++)
61 col_starts[i] = col_partitioning[i];
62 }
63 else
64 col_starts = row_starts;
65
66 par_matrix = hypre_ParCSRMatrixCreate(comm, hypre_IJMatrixGlobalNumRows(matrix),
67 hypre_IJMatrixGlobalNumCols(matrix),row_starts, col_starts, 0, 0, 0);
68
69#else
70 row_starts = hypre_CTAlloc(HYPRE_BigInt,num_procs+1);
71 if (row_partitioning[0])
72 for (i=0; i < num_procs+1; i++)
73 row_starts[i] = row_partitioning[i]-row_partitioning[0];
74 else
75 for (i=0; i < num_procs+1; i++)
76 row_starts[i] = row_partitioning[i];
77 if (row_partitioning != col_partitioning)
78 {
79 col_starts = hypre_CTAlloc(HYPRE_BigInt,num_procs+1);
80 if (col_partitioning[0])
81 for (i=0; i < num_procs+1; i++)
82 col_starts[i] = col_partitioning[i]-col_partitioning[0];
83 else
84 for (i=0; i < num_procs+1; i++)
85 col_starts[i] = col_partitioning[i];
86 }
87 else
88 col_starts = row_starts;
89 par_matrix = hypre_ParCSRMatrixCreate(comm,row_starts[num_procs],
90 col_starts[num_procs],row_starts, col_starts, 0, 0, 0);
91#endif
92
93 hypre_IJMatrixObject(matrix) = par_matrix;
94
95
96 return hypre_error_flag;
97
98}
99
100
101/******************************************************************************
102 *
103 * hypre_IJMatrixSetRowSizesParCSR
104 *
105 *****************************************************************************/
106int
107hypre_IJMatrixSetRowSizesParCSR(hypre_IJMatrix *matrix,
108 const int *sizes)
109{
110 int local_num_rows, local_num_cols;
111 int i, my_id;
112 int *row_space;
113 HYPRE_BigInt *row_partitioning = hypre_IJMatrixRowPartitioning(matrix);
114 HYPRE_BigInt *col_partitioning = hypre_IJMatrixColPartitioning(matrix);
115 hypre_AuxParCSRMatrix *aux_matrix;
116 MPI_Comm comm = hypre_IJMatrixComm(matrix);
117
118 MPI_Comm_rank(comm,&my_id);
119#ifdef HYPRE_NO_GLOBAL_PARTITION
120 local_num_rows = (int) (row_partitioning[1]-row_partitioning[0]);
121 local_num_cols = (int)(col_partitioning[1]-col_partitioning[0]);
122#else
123 local_num_rows = (int)(row_partitioning[my_id+1]-row_partitioning[my_id]);
124 local_num_cols = (int)(col_partitioning[my_id+1]-col_partitioning[my_id]);
125#endif
126 aux_matrix = hypre_IJMatrixTranslator(matrix);
127 row_space = NULL;
128 if (aux_matrix)
129 row_space = hypre_AuxParCSRMatrixRowSpace(aux_matrix);
130 if (!row_space)
131 row_space = hypre_CTAlloc(int, local_num_rows);
132 for (i = 0; i < local_num_rows; i++)
133 row_space[i] = sizes[i];
134 if (!aux_matrix)
135 {
136 hypre_AuxParCSRMatrixCreate(&aux_matrix, local_num_rows,
137 local_num_cols, row_space);
138 hypre_IJMatrixTranslator(matrix) = aux_matrix;
139 }
140 hypre_AuxParCSRMatrixRowSpace(aux_matrix) = row_space;
141
142 return hypre_error_flag;
143
144}
145
146/******************************************************************************
147 *
148 * hypre_IJMatrixSetDiagOffdSizesParCSR
149 * sets diag_i inside the diag part of the ParCSRMatrix
150 * and offd_i inside the offd part,
151 * requires exact row sizes for diag and offd
152 *
153 *****************************************************************************/
154int
155hypre_IJMatrixSetDiagOffdSizesParCSR(hypre_IJMatrix *matrix,
156 const int *diag_sizes,
157 const int *offdiag_sizes)
158{
159 int local_num_rows;
160 int i;
161 hypre_ParCSRMatrix *par_matrix = hypre_IJMatrixObject(matrix);
162 hypre_AuxParCSRMatrix *aux_matrix = hypre_IJMatrixTranslator(matrix);
163 hypre_CSRMatrix *diag;
164 hypre_CSRMatrix *offd;
165 int *diag_i;
166 int *offd_i;
167
168 if (!par_matrix)
169 {
170 hypre_IJMatrixCreateParCSR(matrix);
171 par_matrix = hypre_IJMatrixObject(matrix);
172 }
173
174 diag = hypre_ParCSRMatrixDiag(par_matrix);
175 diag_i = hypre_CSRMatrixI(diag);
176 local_num_rows = hypre_CSRMatrixNumRows(diag);
177 if (!diag_i)
178 diag_i = hypre_CTAlloc(int, local_num_rows+1);
179 for (i = 0; i < local_num_rows; i++)
180 diag_i[i+1] = diag_i[i] + diag_sizes[i];
181 hypre_CSRMatrixI(diag) = diag_i;
182 hypre_CSRMatrixNumNonzeros(diag) = diag_i[local_num_rows];
183 offd = hypre_ParCSRMatrixOffd(par_matrix);
184 offd_i = hypre_CSRMatrixI(offd);
185 if (!offd_i)
186 offd_i = hypre_CTAlloc(int, local_num_rows+1);
187 for (i = 0; i < local_num_rows; i++)
188 offd_i[i+1] = offd_i[i] + offdiag_sizes[i];
189 hypre_CSRMatrixI(offd) = offd_i;
190 hypre_CSRMatrixNumNonzeros(offd) = offd_i[local_num_rows];
191 if (!aux_matrix)
192 {
193 hypre_AuxParCSRMatrixCreate(&aux_matrix, local_num_rows,
194 hypre_CSRMatrixNumCols(diag), NULL);
195 hypre_IJMatrixTranslator(matrix) = aux_matrix;
196 }
197 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
198 if (offd_i[local_num_rows])
199 hypre_AuxParCSRMatrixAuxOffdJ(aux_matrix) =
200 hypre_CTAlloc(HYPRE_BigInt, offd_i[local_num_rows]);
201
202 return hypre_error_flag;
203
204}
205
206/******************************************************************************
207 *
208 * hypre_IJMatrixSetMaxOffProcElmtsParCSR
209 *
210 *****************************************************************************/
211
212int
213hypre_IJMatrixSetMaxOffProcElmtsParCSR(hypre_IJMatrix *matrix,
214 int max_off_proc_elmts)
215{
216 hypre_AuxParCSRMatrix *aux_matrix;
217 int local_num_rows, local_num_cols, my_id;
218 HYPRE_BigInt *row_partitioning = hypre_IJMatrixRowPartitioning(matrix);
219 HYPRE_BigInt *col_partitioning = hypre_IJMatrixColPartitioning(matrix);
220 MPI_Comm comm = hypre_IJMatrixComm(matrix);
221
222 MPI_Comm_rank(comm,&my_id);
223 aux_matrix = hypre_IJMatrixTranslator(matrix);
224 if (!aux_matrix)
225 {
226#ifdef HYPRE_NO_GLOBAL_PARTITION
227 local_num_rows = (int)( row_partitioning[1]-row_partitioning[0]);
228 local_num_cols = (int)( col_partitioning[1]-col_partitioning[0]);
229#else
230 local_num_rows = (int)( row_partitioning[my_id+1]-row_partitioning[my_id]);
231 local_num_cols = (int)( col_partitioning[my_id+1]-col_partitioning[my_id]);
232#endif
233 hypre_AuxParCSRMatrixCreate(&aux_matrix, local_num_rows,
234 local_num_cols, NULL);
235 hypre_IJMatrixTranslator(matrix) = aux_matrix;
236 }
237 hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix) = max_off_proc_elmts;
238
239 return hypre_error_flag;
240}
241
242/******************************************************************************
243 *
244 * hypre_IJMatrixInitializeParCSR
245 *
246 * initializes AuxParCSRMatrix and ParCSRMatrix as necessary
247 *
248 *****************************************************************************/
249
250int
251hypre_IJMatrixInitializeParCSR(hypre_IJMatrix *matrix)
252{
253 hypre_ParCSRMatrix *par_matrix = hypre_IJMatrixObject(matrix);
254 hypre_AuxParCSRMatrix *aux_matrix = hypre_IJMatrixTranslator(matrix);
255 int local_num_rows;
256
257 if (hypre_IJMatrixAssembleFlag(matrix) == 0)
258 {
259 if (!par_matrix)
260 {
261 hypre_IJMatrixCreateParCSR(matrix);
262 par_matrix = hypre_IJMatrixObject(matrix);
263 }
264 local_num_rows =
265 hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(par_matrix));
266 if (!aux_matrix)
267 {
268 hypre_AuxParCSRMatrixCreate(&aux_matrix, local_num_rows,
269 hypre_CSRMatrixNumCols(hypre_ParCSRMatrixDiag(par_matrix)),
270 NULL);
271 hypre_IJMatrixTranslator(matrix) = aux_matrix;
272 }
273
274 hypre_ParCSRMatrixInitialize(par_matrix);
275 hypre_AuxParCSRMatrixInitialize(aux_matrix);
276 if (! hypre_AuxParCSRMatrixNeedAux(aux_matrix))
277 {
278 int i, *indx_diag, *indx_offd, *diag_i, *offd_i;
279 diag_i = hypre_CSRMatrixI(hypre_ParCSRMatrixDiag(par_matrix));
280 offd_i = hypre_CSRMatrixI(hypre_ParCSRMatrixOffd(par_matrix));
281 indx_diag = hypre_AuxParCSRMatrixIndxDiag(aux_matrix);
282
283 indx_offd = hypre_AuxParCSRMatrixIndxOffd(aux_matrix);
284 for (i=0; i < local_num_rows; i++)
285 {
286 indx_diag[i] = diag_i[i];
287 indx_offd[i] = offd_i[i];
288 }
289 }
290 }
291 else /* AB 4/06 - the assemble routine destroys the aux matrix - so we need
292 to recreate if initialize is called again*/
293 {
294 if (!aux_matrix)
295 {
296 local_num_rows =
297 hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(par_matrix));
298 hypre_AuxParCSRMatrixCreate(&aux_matrix, local_num_rows,
299 hypre_CSRMatrixNumCols(hypre_ParCSRMatrixDiag(par_matrix)),
300 NULL);
301 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
302 hypre_IJMatrixTranslator(matrix) = aux_matrix;
303 }
304
305 }
306 return hypre_error_flag;
307}
308
309/******************************************************************************
310 *
311 * hypre_IJMatrixGetRowCountsParCSR
312 *
313 * gets the number of columns for rows specified by the user
314 *
315 *****************************************************************************/
316
317int hypre_IJMatrixGetRowCountsParCSR( hypre_IJMatrix *matrix,
318 int nrows,
319 HYPRE_BigInt *rows,
320 int *ncols)
321{
322 HYPRE_BigInt row_index;
323 MPI_Comm comm = hypre_IJMatrixComm(matrix);
324 hypre_ParCSRMatrix *par_matrix = hypre_IJMatrixObject(matrix);
325
326 HYPRE_BigInt *row_partitioning = hypre_IJMatrixRowPartitioning(matrix);
327
328 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(par_matrix);
329 int *diag_i = hypre_CSRMatrixI(diag);
330
331 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(par_matrix);
332 int *offd_i = hypre_CSRMatrixI(offd);
333
334 int i, my_id;
335
336 MPI_Comm_rank(comm,&my_id);
337
338 for (i=0; i < nrows; i++)
339 {
340 row_index = rows[i];
341#ifdef HYPRE_NO_GLOBAL_PARTITION
342 if (row_index >= row_partitioning[0] &&
343 row_index < row_partitioning[1])
344 {
345 /* compute local row number */
346 row_index -= row_partitioning[0];
347#else
348 if (row_index >= row_partitioning[my_id] &&
349 row_index < row_partitioning[my_id+1])
350 {
351 /* compute local row number */
352 row_index -= row_partitioning[my_id];
353#endif
354 ncols[i] = diag_i[row_index+1]-diag_i[row_index]+offd_i[row_index+1]
355 -offd_i[row_index];
356 }
357 else
358 {
359 ncols[i] = 0;
360#ifdef HYPRE_LONG_LONG
361 printf ("Warning! Row %lld is not on Proc. %d!\n", row_index, my_id);
362#else
363 printf ("Warning! Row %d is not on Proc. %d!\n", row_index, my_id);
364#endif
365 }
366 }
367
368 return hypre_error_flag;
369
370}
371/******************************************************************************
372 *
373 * hypre_IJMatrixGetValuesParCSR
374 *
375 * gets values of an IJMatrix
376 *
377 *****************************************************************************/
378int
379hypre_IJMatrixGetValuesParCSR( hypre_IJMatrix *matrix,
380 int nrows,
381 int *ncols,
382 HYPRE_BigInt *rows,
383 HYPRE_BigInt *cols,
384 double *values)
385{
386 MPI_Comm comm = hypre_IJMatrixComm(matrix);
387 hypre_ParCSRMatrix *par_matrix = hypre_IJMatrixObject(matrix);
388 int assemble_flag = hypre_IJMatrixAssembleFlag(matrix);
389
390 hypre_CSRMatrix *diag;
391 int *diag_i;
392 int *diag_j;
393 double *diag_data;
394
395 hypre_CSRMatrix *offd;
396 int *offd_i;
397 int *offd_j;
398 double *offd_data;
399
400 HYPRE_BigInt *col_map_offd;
401 HYPRE_BigInt *col_starts = hypre_ParCSRMatrixColStarts(par_matrix);
402
403 HYPRE_BigInt *row_partitioning = hypre_IJMatrixRowPartitioning(matrix);
404
405#ifndef HYPRE_NO_GLOBAL_PARTITION
406 HYPRE_BigInt *col_partitioning = hypre_IJMatrixColPartitioning(matrix);
407#endif
408
409 int i, j, n, ii, indx;
410 HYPRE_BigInt col_indx;
411 int num_procs, my_id;
412 HYPRE_BigInt col_0, col_n, row;
413 int row_local, row_size;
414 int warning = 0;
415 int *counter;
416
417 MPI_Comm_size(comm,&num_procs);
418 MPI_Comm_rank(comm,&my_id);
419
420 if (assemble_flag == 0)
421 {
422 hypre_error_in_arg(1);
423 printf("Error! Matrix not assembled yet! HYPRE_IJMatrixGetValues\n");
424 }
425
426#ifdef HYPRE_NO_GLOBAL_PARTITION
427 col_0 = col_starts[0];
428 col_n = col_starts[1]-1;
429#else
430 col_0 = col_starts[my_id];
431 col_n = col_starts[my_id+1]-1;
432#endif
433
434 diag = hypre_ParCSRMatrixDiag(par_matrix);
435 diag_i = hypre_CSRMatrixI(diag);
436 diag_j = hypre_CSRMatrixJ(diag);
437 diag_data = hypre_CSRMatrixData(diag);
438
439 offd = hypre_ParCSRMatrixOffd(par_matrix);
440 offd_i = hypre_CSRMatrixI(offd);
441 if (num_procs > 1)
442 {
443 offd_j = hypre_CSRMatrixJ(offd);
444 offd_data = hypre_CSRMatrixData(offd);
445 col_map_offd = hypre_ParCSRMatrixColMapOffd(par_matrix);
446 }
447
448 if (nrows < 0)
449 {
450 nrows = -nrows;
451
452 counter = hypre_CTAlloc(int,nrows+1);
453 counter[0] = 0;
454 for (i=0; i < nrows; i++)
455 counter[i+1] = counter[i]+ncols[i];
456
457 indx = 0;
458 for (i=0; i < nrows; i++)
459 {
460 row = rows[i];
461#ifdef HYPRE_NO_GLOBAL_PARTITION
462 if (row >= row_partitioning[0] && row < row_partitioning[1])
463 {
464 row_local = (int)(row - row_partitioning[0]);
465#else
466 if (row >= row_partitioning[my_id] && row < row_partitioning[my_id+1])
467 {
468 row_local = (int)(row - row_partitioning[my_id]);
469#endif
470 row_size = diag_i[row_local+1]-diag_i[row_local]+
471 offd_i[row_local+1]-offd_i[row_local];
472 if (counter[i]+row_size > counter[nrows])
473 {
474 hypre_error_in_arg(1);
475 printf ("Error! Not enough memory! HYPRE_IJMatrixGetValues\n");
476 }
477 if (ncols[i] < row_size)
478 warning = 1;
479 for (j = diag_i[row_local]; j < diag_i[row_local+1]; j++)
480 {
481 cols[indx] = (HYPRE_BigInt) diag_j[j] + col_0;
482 values[indx++] = diag_data[j];
483 }
484 for (j = offd_i[row_local]; j < offd_i[row_local+1]; j++)
485 {
486 cols[indx] = col_map_offd[offd_j[j]];
487 values[indx++] = offd_data[j];
488 }
489 counter[i+1] = indx;
490 }
491 else
492#ifdef HYPRE_LONG_LONG
493 printf ("Warning! Row %lld is not on Proc. %d!\n", row, my_id);
494#else
495 printf ("Warning! Row %d is not on Proc. %d!\n", row, my_id);
496#endif
497 }
498 if (warning)
499 {
500 for (i=0; i < nrows; i++)
501 ncols[i] = counter[i+1] - counter[i];
502 printf ("Warning! ncols has been changed!\n");
503 }
504 hypre_TFree(counter);
505 }
506 else
507 {
508 indx = 0;
509 for (ii=0; ii < nrows; ii++)
510 {
511 row = rows[ii];
512 n = ncols[ii];
513#ifdef HYPRE_NO_GLOBAL_PARTITION
514 if (row >= row_partitioning[0] && row < row_partitioning[1])
515 {
516 row_local = (int)(row - row_partitioning[0]);
517 /* compute local row number */
518 for (i=0; i < n; i++)
519 {
520 col_indx = cols[indx] - hypre_IJMatrixGlobalFirstCol(matrix);
521
522#else
523 if (row >= row_partitioning[my_id] && row < row_partitioning[my_id+1])
524 {
525 row_local = (int)(row - row_partitioning[my_id]);
526 /* compute local row number */
527 for (i=0; i < n; i++)
528 {
529 col_indx = cols[indx] - col_partitioning[0];
530#endif
531
532
533
534 values[indx] = 0.0;
535 if (col_indx < col_0 || col_indx > col_n)
536 /* search in offd */
537 {
538 for (j=offd_i[row_local]; j < offd_i[row_local+1]; j++)
539 {
540 if (col_map_offd[offd_j[j]] == col_indx)
541 {
542 values[indx] = offd_data[j];
543 break;
544 }
545 }
546 }
547 else /* search in diag */
548 {
549 col_indx = col_indx - col_0;
550 for (j=diag_i[row_local]; j < diag_i[row_local+1]; j++)
551 {
552 if (diag_j[j] == (int)col_indx)
553 {
554 values[indx] = diag_data[j];
555 break;
556 }
557 }
558 }
559 indx++;
560 }
561 }
562 else
563#ifdef HYPRE_LONG_LONG
564 printf ("Warning! Row %lld is not on Proc. %d!\n", row, my_id);
565#else
566 printf ("Warning! Row %d is not on Proc. %d!\n", row, my_id);
567#endif
568 }
569 }
570
571 return hypre_error_flag;
572
573}
574/******************************************************************************
575 *
576 * hypre_IJMatrixSetValuesParCSR
577 *
578 * sets or adds row values to an IJMatrix before assembly,
579 *
580 *****************************************************************************/
581int
582hypre_IJMatrixSetValuesParCSR( hypre_IJMatrix *matrix,
583 int nrows,
584 int *ncols,
585 const HYPRE_BigInt *rows,
586 const HYPRE_BigInt *cols,
587 const double *values)
588{
589 hypre_ParCSRMatrix *par_matrix;
590 hypre_CSRMatrix *diag, *offd;
591 hypre_AuxParCSRMatrix *aux_matrix;
592 HYPRE_BigInt *row_partitioning;
593 HYPRE_BigInt *col_partitioning;
594 MPI_Comm comm = hypre_IJMatrixComm(matrix);
595 int num_procs, my_id;
596 int row_local;
597 HYPRE_BigInt col_0, col_n, row;
598 int i, ii, j, n, not_found;
599 HYPRE_BigInt **aux_j;
600 HYPRE_BigInt *local_j;
601 HYPRE_BigInt *tmp_j;
602 double **aux_data;
603 double *local_data;
604 double *tmp_data;
605 int diag_space, offd_space;
606 int *row_length, *row_space;
607 int need_aux;
608 int tmp_indx, indx;
609 int space, size, old_size;
610 int cnt, cnt_diag, cnt_offd;
611 int pos_diag, pos_offd;
612 int len_diag, len_offd;
613 int offd_indx, diag_indx;
614 int *diag_i;
615 int *diag_j;
616 double *diag_data;
617 int *offd_i;
618 int *offd_j;
619 HYPRE_BigInt *aux_offd_j;
620 double *offd_data;
621 HYPRE_BigInt first;
622 int current_num_elmts;
623 int max_off_proc_elmts;
624 int off_proc_i_indx;
625 HYPRE_BigInt *off_proc_i;
626 HYPRE_BigInt *off_proc_j;
627 double *off_proc_data;
628 MPI_Comm_size(comm, &num_procs);
629 MPI_Comm_rank(comm, &my_id);
630 par_matrix = hypre_IJMatrixObject( matrix );
631 row_partitioning = hypre_IJMatrixRowPartitioning(matrix);
632 col_partitioning = hypre_IJMatrixColPartitioning(matrix);
633
634#ifdef HYPRE_NO_GLOBAL_PARTITION
635 col_0 = col_partitioning[0];
636 col_n = col_partitioning[1]-1;
637 first = hypre_IJMatrixGlobalFirstCol(matrix);
638#else
639 col_0 = col_partitioning[my_id];
640 col_n = col_partitioning[my_id+1]-1;
641 first = col_partitioning[0];
642#endif
643 if (nrows < 0)
644 {
645 hypre_error_in_arg(2);
646 printf("Error! nrows negative! HYPRE_IJMatrixSetValues\n");
647 }
648
649 if (hypre_IJMatrixAssembleFlag(matrix))
650 {
651 HYPRE_BigInt *col_map_offd;
652 int num_cols_offd;
653 int j_offd;
654 indx = 0;
655 for (ii=0; ii < nrows; ii++)
656 {
657 row = rows[ii];
658 n = ncols[ii];
659
660 /* processor owns the row */
661
662#ifdef HYPRE_NO_GLOBAL_PARTITION
663 if (row >= row_partitioning[0] && row < row_partitioning[1])
664 {
665 row_local = (int)(row - row_partitioning[0]);
666#else
667 if (row >= row_partitioning[my_id] && row < row_partitioning[my_id+1])
668 {
669 row_local = (int)(row - row_partitioning[my_id]);
670#endif
671
672 /* compute local row number */
673 diag = hypre_ParCSRMatrixDiag(par_matrix);
674 diag_i = hypre_CSRMatrixI(diag);
675 diag_j = hypre_CSRMatrixJ(diag);
676 diag_data = hypre_CSRMatrixData(diag);
677 offd = hypre_ParCSRMatrixOffd(par_matrix);
678 offd_i = hypre_CSRMatrixI(offd);
679 num_cols_offd = hypre_CSRMatrixNumCols(offd);
680 if (num_cols_offd)
681 {
682 col_map_offd = hypre_ParCSRMatrixColMapOffd(par_matrix);
683 offd_j = hypre_CSRMatrixJ(offd);
684 offd_data = hypre_CSRMatrixData(offd);
685 }
686 size = diag_i[row_local+1] - diag_i[row_local]
687 + offd_i[row_local+1] - offd_i[row_local];
688
689 if (n > size)
690 {
691 hypre_error(HYPRE_ERROR_GENERIC);
692#ifdef HYPRE_LONG_LONG
693 printf (" row %lld too long! \n", row);
694#else
695 printf (" row %d too long! \n", row);
696#endif
697 return hypre_error_flag;
698 }
699
700 pos_diag = diag_i[row_local];
701 pos_offd = offd_i[row_local];
702 len_diag = diag_i[row_local+1];
703 len_offd = offd_i[row_local+1];
704 not_found = 1;
705
706 for (i=0; i < n; i++)
707 {
708 if (cols[indx] < col_0 || cols[indx] > col_n)
709 /* insert into offd */
710 {
711 j_offd = hypre_BigBinarySearch(col_map_offd,cols[indx]-first,
712 num_cols_offd);
713 if (j_offd == -1)
714 {
715 hypre_error(HYPRE_ERROR_GENERIC);
716#ifdef HYPRE_LONG_LONG
717 printf (" Error, element %lld %lld does not exist\n",
718 row, cols[indx]);
719#else
720 printf (" Error, element %d %d does not exist\n",
721 row, cols[indx]);
722#endif
723 return hypre_error_flag;
724 }
725 for (j=pos_offd; j < len_offd; j++)
726 {
727 if (offd_j[j] == j_offd)
728 {
729 offd_data[j] = values[indx];
730 not_found = 0;
731 break;
732 }
733 }
734 if (not_found)
735 {
736 hypre_error(HYPRE_ERROR_GENERIC);
737#ifdef HYPRE_LONG_LONG
738 printf (" Error, element %lld %lld does not exist\n",
739 row, cols[indx]);
740#else
741 printf (" Error, element %d %d does not exist\n",
742 row, cols[indx]);
743#endif
744 return hypre_error_flag;
745 }
746 not_found = 1;
747 }
748 /* diagonal element */
749 else if (cols[indx] == row)
750 {
751 if (diag_j[pos_diag] != row_local)
752 {
753 hypre_error(HYPRE_ERROR_GENERIC);
754#ifdef HYPRE_LONG_LONG
755 printf (" Error, element %lld %lld does not exist\n",
756 row, cols[indx]);
757#else
758 printf (" Error, element %d %d does not exist\n",
759 row, cols[indx]);
760#endif
761 /* return -1;*/
762 return hypre_error_flag;
763 }
764 diag_data[pos_diag] = values[indx];
765 }
766 else /* insert into diag */
767 {
768 for (j=pos_diag; j < len_diag; j++)
769 {
770 if (diag_j[j] == (int)(cols[indx]-col_0))
771 {
772 diag_data[j] = values[indx];
773 not_found = 0;
774 break;
775 }
776 }
777 if (not_found)
778 {
779 hypre_error(HYPRE_ERROR_GENERIC);
780#ifdef HYPRE_LONG_LONG
781 printf (" Error, element %lld %lld does not exist\n",
782 row, cols[indx]);
783#else
784 printf (" Error, element %d %d does not exist\n",
785 row, cols[indx]);
786#endif
787 /* return -1; */
788 return hypre_error_flag;
789 }
790 }
791 indx++;
792 }
793 }
794
795 /* processor does not own the row */
796
797 else
798 {
799 aux_matrix = hypre_IJMatrixTranslator(matrix);
800 if (!aux_matrix)
801 {
802#ifdef HYPRE_NO_GLOBAL_PARTITION
803 size = (int)(row_partitioning[1]-row_partitioning[0]);
804#else
805 size = (int)(row_partitioning[my_id+1]-row_partitioning[my_id]);
806#endif
807 hypre_AuxParCSRMatrixCreate(&aux_matrix, size, size, NULL);
808 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
809 hypre_IJMatrixTranslator(matrix) = aux_matrix;
810 }
811 current_num_elmts
812 = hypre_AuxParCSRMatrixCurrentNumElmts(aux_matrix);
813 max_off_proc_elmts
814 = hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix);
815 off_proc_i_indx = hypre_AuxParCSRMatrixOffProcIIndx(aux_matrix);
816 off_proc_i = hypre_AuxParCSRMatrixOffProcI(aux_matrix);
817 off_proc_j = hypre_AuxParCSRMatrixOffProcJ(aux_matrix);
818 off_proc_data = hypre_AuxParCSRMatrixOffProcData(aux_matrix);
819
820 if (!max_off_proc_elmts)
821 {
822 max_off_proc_elmts = hypre_max(n,1000);
823 hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix) =
824 max_off_proc_elmts;
825 hypre_AuxParCSRMatrixOffProcI(aux_matrix)
826 = hypre_CTAlloc(HYPRE_BigInt,2*max_off_proc_elmts);
827 hypre_AuxParCSRMatrixOffProcJ(aux_matrix)
828 = hypre_CTAlloc(HYPRE_BigInt,max_off_proc_elmts);
829 hypre_AuxParCSRMatrixOffProcData(aux_matrix)
830 = hypre_CTAlloc(double,max_off_proc_elmts);
831 off_proc_i = hypre_AuxParCSRMatrixOffProcI(aux_matrix);
832 off_proc_j = hypre_AuxParCSRMatrixOffProcJ(aux_matrix);
833 off_proc_data = hypre_AuxParCSRMatrixOffProcData(aux_matrix);
834 }
835 else if (current_num_elmts + n > max_off_proc_elmts)
836 {
837 max_off_proc_elmts += 3*n;
838 off_proc_i = hypre_TReAlloc(off_proc_i,HYPRE_BigInt,2*max_off_proc_elmts);
839 off_proc_j = hypre_TReAlloc(off_proc_j,HYPRE_BigInt,max_off_proc_elmts);
840 off_proc_data = hypre_TReAlloc(off_proc_data,double,
841 max_off_proc_elmts);
842 hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix)
843 = max_off_proc_elmts;
844 hypre_AuxParCSRMatrixOffProcI(aux_matrix) = off_proc_i;
845 hypre_AuxParCSRMatrixOffProcJ(aux_matrix) = off_proc_j;
846 hypre_AuxParCSRMatrixOffProcData(aux_matrix) = off_proc_data;
847 }
848 off_proc_i[off_proc_i_indx++] = row;
849 off_proc_i[off_proc_i_indx++] = n;
850 for (i=0; i < n; i++)
851 {
852 off_proc_j[current_num_elmts] = cols[indx];
853 off_proc_data[current_num_elmts++] = values[indx++];
854 }
855 hypre_AuxParCSRMatrixOffProcIIndx(aux_matrix) = off_proc_i_indx;
856 hypre_AuxParCSRMatrixCurrentNumElmts(aux_matrix)
857 = current_num_elmts;
858 }
859 } /* end of for loop */
860 }
861 else
862 {
863 aux_matrix = hypre_IJMatrixTranslator(matrix);
864 row_space = hypre_AuxParCSRMatrixRowSpace(aux_matrix);
865 row_length = hypre_AuxParCSRMatrixRowLength(aux_matrix);
866 need_aux = hypre_AuxParCSRMatrixNeedAux(aux_matrix);
867 indx = 0;
868 for (ii=0; ii < nrows; ii++)
869 {
870 row = rows[ii];
871 n = ncols[ii];
872 /* processor owns the row */
873#ifdef HYPRE_NO_GLOBAL_PARTITION
874 if (row >= row_partitioning[0] && row < row_partitioning[1])
875 {
876 row_local = (int)(row - row_partitioning[0]);
877#else
878 if (row >= row_partitioning[my_id] && row < row_partitioning[my_id+1])
879 {
880 row_local = (int)(row - row_partitioning[my_id]);
881
882#endif
883 /* compute local row number */
884 if (need_aux)
885 {
886 aux_j = hypre_AuxParCSRMatrixAuxJ(aux_matrix);
887 aux_data = hypre_AuxParCSRMatrixAuxData(aux_matrix);
888 local_j = aux_j[row_local];
889 local_data = aux_data[row_local];
890 space = row_space[row_local];
891 old_size = row_length[row_local];
892 size = space - old_size;
893 if (size < n)
894 {
895 size = n - size;
896 tmp_j = hypre_CTAlloc(HYPRE_BigInt,size);
897 tmp_data = hypre_CTAlloc(double,size);
898 }
899 else
900 {
901 tmp_j = NULL;
902 }
903 tmp_indx = 0;
904 not_found = 1;
905 size = old_size;
906 for (i=0; i < n; i++)
907 {
908 for (j=0; j < old_size; j++)
909 {
910 if (local_j[j] == cols[indx])
911 {
912 local_data[j] = values[indx];
913 not_found = 0;
914 break;
915 }
916 }
917 if (not_found)
918 {
919 if (size < space)
920 {
921 local_j[size] = cols[indx];
922 local_data[size++] = values[indx];
923 }
924 else
925 {
926 tmp_j[tmp_indx] = cols[indx];
927 tmp_data[tmp_indx++] = values[indx];
928 }
929 }
930 not_found = 1;
931 indx++;
932 }
933
934 row_length[row_local] = size+tmp_indx;
935
936 if (tmp_indx)
937 {
938 aux_j[row_local] = hypre_TReAlloc(aux_j[row_local],HYPRE_BigInt,
939 size+tmp_indx);
940 aux_data[row_local] = hypre_TReAlloc(aux_data[row_local],
941 double,size+tmp_indx);
942 row_space[row_local] = size+tmp_indx;
943 local_j = aux_j[row_local];
944 local_data = aux_data[row_local];
945 }
946
947 cnt = size;
948
949 for (i=0; i < tmp_indx; i++)
950 {
951 local_j[cnt] = tmp_j[i];
952 local_data[cnt++] = tmp_data[i];
953 }
954
955 if (tmp_j)
956 {
957 hypre_TFree(tmp_j);
958 hypre_TFree(tmp_data);
959 }
960 }
961 else /* insert immediately into data in ParCSRMatrix structure */
962 {
963 int col_j;
964 offd_indx =hypre_AuxParCSRMatrixIndxOffd(aux_matrix)[row_local];
965 diag_indx =hypre_AuxParCSRMatrixIndxDiag(aux_matrix)[row_local];
966 diag = hypre_ParCSRMatrixDiag(par_matrix);
967 diag_i = hypre_CSRMatrixI(diag);
968 diag_j = hypre_CSRMatrixJ(diag);
969 diag_data = hypre_CSRMatrixData(diag);
970 offd = hypre_ParCSRMatrixOffd(par_matrix);
971 offd_i = hypre_CSRMatrixI(offd);
972 if (num_procs > 1)
973 {
974 offd_j = hypre_CSRMatrixJ(offd);
975 offd_data = hypre_CSRMatrixData(offd);
976 aux_offd_j = hypre_AuxParCSRMatrixAuxOffdJ(aux_matrix);
977 }
978 cnt_diag = diag_indx;
979 cnt_offd = offd_indx;
980 diag_space = diag_i[row_local+1];
981 offd_space = offd_i[row_local+1];
982 not_found = 1;
983 for (i=0; i < n; i++)
984 {
985 if (cols[indx] < col_0 || cols[indx] > col_n)
986 /* insert into offd */
987 {
988 for (j=offd_i[row_local]; j < offd_indx; j++)
989 {
990 if (aux_offd_j[j] == cols[indx])
991 {
992 offd_data[j] = values[indx];
993 not_found = 0;
994 break;
995 }
996 }
997 if (not_found)
998 {
999 if (cnt_offd < offd_space)
1000 {
1001 aux_offd_j[cnt_offd] = cols[indx];
1002 offd_data[cnt_offd++] = values[indx];
1003 }
1004 else
1005 {
1006 hypre_error(HYPRE_ERROR_GENERIC);
1007#ifdef HYPRE_LONG_LONG
1008 printf("Error in row %lld ! Too many elements!\n",
1009 row);
1010#else
1011 printf("Error in row %d ! Too many elements!\n",
1012 row);
1013#endif
1014 /* return 1; */
1015 return hypre_error_flag;
1016 }
1017 }
1018 not_found = 1;
1019 }
1020 else /* insert into diag */
1021 {
1022 for (j=diag_i[row_local]; j < diag_indx; j++)
1023 {
1024 col_j = (int)(cols[indx]-col_0);
1025 if (diag_j[j] == col_j)
1026 {
1027 diag_data[j] = values[indx];
1028 not_found = 0;
1029 break;
1030 }
1031 }
1032 if (not_found)
1033 {
1034 if (cnt_diag < diag_space)
1035 {
1036 diag_j[cnt_diag] = col_j;
1037 diag_data[cnt_diag++] = values[indx];
1038 }
1039 else
1040 {
1041 hypre_error(HYPRE_ERROR_GENERIC);
1042#ifdef HYPRE_LONG_LONG
1043 printf("Error in row %lld ! Too many elements !\n",
1044 row);
1045#else
1046 printf("Error in row %d ! Too many elements !\n",
1047 row);
1048#endif
1049 /* return 1; */
1050 return hypre_error_flag;
1051 }
1052 }
1053 not_found = 1;
1054 }
1055 indx++;
1056 }
1057
1058 hypre_AuxParCSRMatrixIndxDiag(aux_matrix)[row_local] = cnt_diag;
1059 hypre_AuxParCSRMatrixIndxOffd(aux_matrix)[row_local] = cnt_offd;
1060
1061 }
1062 }
1063
1064 /* processor does not own the row */
1065 else
1066 {
1067
1068 current_num_elmts
1069 = hypre_AuxParCSRMatrixCurrentNumElmts(aux_matrix);
1070 max_off_proc_elmts
1071 = hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix);
1072 off_proc_i_indx = hypre_AuxParCSRMatrixOffProcIIndx(aux_matrix);
1073 off_proc_i = hypre_AuxParCSRMatrixOffProcI(aux_matrix);
1074 off_proc_j = hypre_AuxParCSRMatrixOffProcJ(aux_matrix);
1075 off_proc_data = hypre_AuxParCSRMatrixOffProcData(aux_matrix);
1076
1077 if (!max_off_proc_elmts)
1078 {
1079 max_off_proc_elmts = hypre_max(n,1000);
1080 hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix) =
1081 max_off_proc_elmts;
1082 hypre_AuxParCSRMatrixOffProcI(aux_matrix)
1083 = hypre_CTAlloc(HYPRE_BigInt,2*max_off_proc_elmts);
1084 hypre_AuxParCSRMatrixOffProcJ(aux_matrix)
1085 = hypre_CTAlloc(HYPRE_BigInt,max_off_proc_elmts);
1086 hypre_AuxParCSRMatrixOffProcData(aux_matrix)
1087 = hypre_CTAlloc(double,max_off_proc_elmts);
1088 off_proc_i = hypre_AuxParCSRMatrixOffProcI(aux_matrix);
1089 off_proc_j = hypre_AuxParCSRMatrixOffProcJ(aux_matrix);
1090 off_proc_data = hypre_AuxParCSRMatrixOffProcData(aux_matrix);
1091 }
1092 else if (current_num_elmts + n > max_off_proc_elmts)
1093 {
1094 max_off_proc_elmts += 3*n;
1095 off_proc_i = hypre_TReAlloc(off_proc_i,HYPRE_BigInt,2*max_off_proc_elmts);
1096 off_proc_j = hypre_TReAlloc(off_proc_j,HYPRE_BigInt,max_off_proc_elmts);
1097 off_proc_data = hypre_TReAlloc(off_proc_data,double,
1098 max_off_proc_elmts);
1099 hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix)
1100 = max_off_proc_elmts;
1101 hypre_AuxParCSRMatrixOffProcI(aux_matrix) = off_proc_i;
1102 hypre_AuxParCSRMatrixOffProcJ(aux_matrix) = off_proc_j;
1103 hypre_AuxParCSRMatrixOffProcData(aux_matrix) = off_proc_data;
1104 }
1105 off_proc_i[off_proc_i_indx++] = row;
1106 off_proc_i[off_proc_i_indx++] = n;
1107 for (i=0; i < n; i++)
1108 {
1109 off_proc_j[current_num_elmts] = cols[indx];
1110 off_proc_data[current_num_elmts++] = values[indx++];
1111 }
1112 hypre_AuxParCSRMatrixOffProcIIndx(aux_matrix) = off_proc_i_indx;
1113 hypre_AuxParCSRMatrixCurrentNumElmts(aux_matrix)
1114 = current_num_elmts;
1115 }
1116 }
1117 }
1118 return hypre_error_flag;
1119
1120}
1121
1122/******************************************************************************
1123 *
1124 * hypre_IJMatrixAddToValuesParCSR
1125 *
1126 * adds row values to an IJMatrix
1127 *
1128 *****************************************************************************/
1129int
1130hypre_IJMatrixAddToValuesParCSR( hypre_IJMatrix *matrix,
1131 int nrows,
1132 int *ncols,
1133 const HYPRE_BigInt *rows,
1134 const HYPRE_BigInt *cols,
1135 const double *values)
1136{
1137 hypre_ParCSRMatrix *par_matrix;
1138 hypre_CSRMatrix *diag, *offd;
1139 hypre_AuxParCSRMatrix *aux_matrix;
1140 HYPRE_BigInt *row_partitioning;
1141 HYPRE_BigInt *col_partitioning;
1142 MPI_Comm comm = hypre_IJMatrixComm(matrix);
1143 int num_procs, my_id;
1144 int row_local;
1145 HYPRE_BigInt row;
1146 HYPRE_BigInt col_0, col_n;
1147 int i, ii, j, n, not_found, col_j;
1148 HYPRE_BigInt **aux_j;
1149 HYPRE_BigInt *local_j;
1150 HYPRE_BigInt *tmp_j;
1151 double **aux_data;
1152 double *local_data;
1153 double *tmp_data;
1154 int diag_space, offd_space;
1155 int *row_length, *row_space;
1156 int need_aux;
1157 int tmp_indx, indx;
1158 int space, size, old_size;
1159 int cnt, cnt_diag, cnt_offd;
1160 int pos_diag, pos_offd;
1161 int len_diag, len_offd;
1162 int offd_indx, diag_indx;
1163 int first;
1164 int *diag_i;
1165 int *diag_j;
1166 double *diag_data;
1167 int *offd_i;
1168 int *offd_j;
1169 HYPRE_BigInt *aux_offd_j;
1170 double *offd_data;
1171 int current_num_elmts;
1172 int max_off_proc_elmts;
1173 int off_proc_i_indx;
1174 HYPRE_BigInt *off_proc_i;
1175 HYPRE_BigInt *off_proc_j;
1176 double *off_proc_data;
1177
1178 MPI_Comm_size(comm, &num_procs);
1179 MPI_Comm_rank(comm, &my_id);
1180 par_matrix = hypre_IJMatrixObject( matrix );
1181 row_partitioning = hypre_IJMatrixRowPartitioning(matrix);
1182 col_partitioning = hypre_IJMatrixColPartitioning(matrix);
1183#ifdef HYPRE_NO_GLOBAL_PARTITION
1184 col_0 = col_partitioning[0];
1185 col_n = col_partitioning[1]-1;
1186 first = hypre_IJMatrixGlobalFirstCol(matrix);
1187#else
1188 col_0 = col_partitioning[my_id];
1189 col_n = col_partitioning[my_id+1]-1;
1190 first = col_partitioning[0];
1191#endif
1192 if (hypre_IJMatrixAssembleFlag(matrix))
1193 {
1194
1195
1196 int num_cols_offd;
1197 HYPRE_BigInt *col_map_offd;
1198 int j_offd;
1199 indx = 0;
1200
1201
1202 /* AB - 4/06 - need to get this object*/
1203 aux_matrix = hypre_IJMatrixTranslator(matrix);
1204
1205 for (ii=0; ii < nrows; ii++)
1206 {
1207 row = rows[ii];
1208 n = ncols[ii];
1209#ifdef HYPRE_NO_GLOBAL_PARTITION
1210 if (row >= row_partitioning[0] && row < row_partitioning[1])
1211 {
1212 row_local = (int)(row - row_partitioning[0]);
1213#else
1214 if (row >= row_partitioning[my_id] && row < row_partitioning[my_id+1])
1215 {
1216 row_local = (int)(row - row_partitioning[my_id]);
1217#endif
1218 /* compute local row number */
1219 diag = hypre_ParCSRMatrixDiag(par_matrix);
1220 diag_i = hypre_CSRMatrixI(diag);
1221 diag_j = hypre_CSRMatrixJ(diag);
1222 diag_data = hypre_CSRMatrixData(diag);
1223 offd = hypre_ParCSRMatrixOffd(par_matrix);
1224 offd_i = hypre_CSRMatrixI(offd);
1225 num_cols_offd = hypre_CSRMatrixNumCols(offd);
1226 if (num_cols_offd)
1227 {
1228 col_map_offd = hypre_ParCSRMatrixColMapOffd(par_matrix);
1229 offd_j = hypre_CSRMatrixJ(offd);
1230 offd_data = hypre_CSRMatrixData(offd);
1231 }
1232 size = diag_i[row_local+1] - diag_i[row_local]
1233 + offd_i[row_local+1] - offd_i[row_local];
1234
1235 if (n > size)
1236 {
1237 hypre_error(HYPRE_ERROR_GENERIC);
1238#ifdef HYPRE_LONG_LONG
1239 printf (" row %lld too long! \n", row);
1240#else
1241 printf (" row %d too long! \n", row);
1242#endif
1243 /* return -1; */
1244 return hypre_error_flag;
1245 }
1246
1247 pos_diag = diag_i[row_local];
1248 pos_offd = offd_i[row_local];
1249 len_diag = diag_i[row_local+1];
1250 len_offd = offd_i[row_local+1];
1251 not_found = 1;
1252
1253 for (i=0; i < n; i++)
1254 {
1255 if (cols[indx] < col_0 || cols[indx] > col_n)
1256 /* insert into offd */
1257 {
1258 j_offd = hypre_BigBinarySearch(col_map_offd,cols[indx]-first,
1259 num_cols_offd);
1260 if (j_offd == -1)
1261 {
1262 hypre_error(HYPRE_ERROR_GENERIC);
1263#ifdef HYPRE_LONG_LONG
1264 printf (" Error, element %lld %lld does not exist\n",
1265 row, cols[indx]);
1266#else
1267 printf (" Error, element %d %d does not exist\n",
1268 row, cols[indx]);
1269#endif
1270 return hypre_error_flag;
1271 /* return -1; */
1272 }
1273 for (j=pos_offd; j < len_offd; j++)
1274 {
1275 if (offd_j[j] == j_offd)
1276 {
1277 offd_data[j] += values[indx];
1278 not_found = 0;
1279 break;
1280 }
1281 }
1282 if (not_found)
1283 {
1284 hypre_error(HYPRE_ERROR_GENERIC);
1285#ifdef HYPRE_LONG_LONG
1286 printf (" Error, element %lld %lld does not exist\n",
1287 row, cols[indx]);
1288#else
1289 printf (" Error, element %d %d does not exist\n",
1290 row, cols[indx]);
1291#endif
1292 /* return -1;*/
1293 return hypre_error_flag;
1294 }
1295 not_found = 1;
1296 }
1297 /* diagonal element */
1298 else if (cols[indx] == row)
1299 {
1300 if (diag_j[pos_diag] != row_local)
1301 {
1302 hypre_error(HYPRE_ERROR_GENERIC);
1303#ifdef HYPRE_LONG_LONG
1304 printf (" Error, element %lld %lld does not exist\n",
1305 row, cols[indx]);
1306#else
1307 printf (" Error, element %d %d does not exist\n",
1308 row, cols[indx]);
1309#endif
1310 /* return -1; */
1311 return hypre_error_flag;
1312 }
1313 diag_data[pos_diag] += values[indx];
1314 }
1315 else /* insert into diag */
1316 {
1317 for (j=pos_diag; j < len_diag; j++)
1318 {
1319 if (diag_j[j] == (int)(cols[indx]-col_0))
1320 {
1321 diag_data[j] += values[indx];
1322 not_found = 0;
1323 break;
1324 }
1325 }
1326 if (not_found)
1327 {
1328 hypre_error(HYPRE_ERROR_GENERIC);
1329#ifdef HYPRE_LONG_LONG
1330 printf (" Error, element %lld %lld does not exist\n",
1331 row, cols[indx]);
1332#else
1333 printf (" Error, element %d %d does not exist\n",
1334 row, cols[indx]);
1335#endif
1336 /* return -1;*/
1337 return hypre_error_flag;
1338 }
1339 }
1340 indx++;
1341 }
1342 }
1343 /* not my row */
1344 else
1345 {
1346 if (!aux_matrix)
1347 {
1348#ifdef HYPRE_NO_GLOBAL_PARTITION
1349 size = (int)( row_partitioning[1]-row_partitioning[0]);
1350#else
1351 size = (int)( row_partitioning[my_id+1]-row_partitioning[my_id]);
1352#endif
1353 hypre_AuxParCSRMatrixCreate(&aux_matrix, size, size, NULL);
1354 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
1355 hypre_IJMatrixTranslator(matrix) = aux_matrix;
1356 }
1357 current_num_elmts
1358 = hypre_AuxParCSRMatrixCurrentNumElmts(aux_matrix);
1359 max_off_proc_elmts
1360 = hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix);
1361 off_proc_i_indx = hypre_AuxParCSRMatrixOffProcIIndx(aux_matrix);
1362 off_proc_i = hypre_AuxParCSRMatrixOffProcI(aux_matrix);
1363 off_proc_j = hypre_AuxParCSRMatrixOffProcJ(aux_matrix);
1364 off_proc_data = hypre_AuxParCSRMatrixOffProcData(aux_matrix);
1365
1366 if (!max_off_proc_elmts)
1367 {
1368 max_off_proc_elmts = hypre_max(n,1000);
1369 hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix) =
1370 max_off_proc_elmts;
1371 hypre_AuxParCSRMatrixOffProcI(aux_matrix)
1372 = hypre_CTAlloc(HYPRE_BigInt,2*max_off_proc_elmts);
1373 hypre_AuxParCSRMatrixOffProcJ(aux_matrix)
1374 = hypre_CTAlloc(HYPRE_BigInt,max_off_proc_elmts);
1375 hypre_AuxParCSRMatrixOffProcData(aux_matrix)
1376 = hypre_CTAlloc(double,max_off_proc_elmts);
1377 off_proc_i = hypre_AuxParCSRMatrixOffProcI(aux_matrix);
1378 off_proc_j = hypre_AuxParCSRMatrixOffProcJ(aux_matrix);
1379 off_proc_data = hypre_AuxParCSRMatrixOffProcData(aux_matrix);
1380 }
1381 else if (current_num_elmts + n > max_off_proc_elmts)
1382 {
1383 max_off_proc_elmts += 3*n;
1384 off_proc_i = hypre_TReAlloc(off_proc_i,HYPRE_BigInt,2*max_off_proc_elmts);
1385 off_proc_j = hypre_TReAlloc(off_proc_j,HYPRE_BigInt,max_off_proc_elmts);
1386 off_proc_data = hypre_TReAlloc(off_proc_data,double,
1387 max_off_proc_elmts);
1388 hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix)
1389 = max_off_proc_elmts;
1390 hypre_AuxParCSRMatrixOffProcI(aux_matrix) = off_proc_i;
1391 hypre_AuxParCSRMatrixOffProcJ(aux_matrix) = off_proc_j;
1392 hypre_AuxParCSRMatrixOffProcData(aux_matrix) = off_proc_data;
1393 }
1394
1395 /* AB - 4/6 - the row should be negative to indicate an add */
1396 /* off_proc_i[off_proc_i_indx++] = row; */
1397 off_proc_i[off_proc_i_indx++] = -row-1;
1398
1399 off_proc_i[off_proc_i_indx++] = n;
1400 for (i=0; i < n; i++)
1401 {
1402 off_proc_j[current_num_elmts] = cols[indx];
1403 off_proc_data[current_num_elmts++] = values[indx++];
1404 }
1405 hypre_AuxParCSRMatrixOffProcIIndx(aux_matrix) = off_proc_i_indx;
1406 hypre_AuxParCSRMatrixCurrentNumElmts(aux_matrix)
1407 = current_num_elmts;
1408 }
1409 }
1410 }
1411
1412/* not assembled */
1413 else
1414 {
1415 aux_matrix = hypre_IJMatrixTranslator(matrix);
1416 row_space = hypre_AuxParCSRMatrixRowSpace(aux_matrix);
1417 row_length = hypre_AuxParCSRMatrixRowLength(aux_matrix);
1418 need_aux = hypre_AuxParCSRMatrixNeedAux(aux_matrix);
1419 indx = 0;
1420 for (ii=0; ii < nrows; ii++)
1421 {
1422 row = rows[ii];
1423 n = ncols[ii];
1424#ifdef HYPRE_NO_GLOBAL_PARTITION
1425 if (row >= row_partitioning[0] && row < row_partitioning[1])
1426 {
1427 row_local = (int)(row - row_partitioning[0]);
1428#else
1429 if (row >= row_partitioning[my_id] && row < row_partitioning[my_id+1])
1430 {
1431 row_local = (int)(row - row_partitioning[my_id]);
1432#endif
1433 /* compute local row number */
1434 if (need_aux)
1435 {
1436 aux_j = hypre_AuxParCSRMatrixAuxJ(aux_matrix);
1437 aux_data = hypre_AuxParCSRMatrixAuxData(aux_matrix);
1438 local_j = aux_j[row_local];
1439 local_data = aux_data[row_local];
1440 space = row_space[row_local];
1441 old_size = row_length[row_local];
1442 size = space - old_size;
1443 if (size < n)
1444 {
1445 size = n - size;
1446 tmp_j = hypre_CTAlloc(HYPRE_BigInt,size);
1447 tmp_data = hypre_CTAlloc(double,size);
1448 }
1449 else
1450 {
1451 tmp_j = NULL;
1452 }
1453 tmp_indx = 0;
1454 not_found = 1;
1455 size = old_size;
1456 for (i=0; i < n; i++)
1457 {
1458 for (j=0; j < old_size; j++)
1459 {
1460 if (local_j[j] == cols[indx])
1461 {
1462 local_data[j] += values[indx];
1463 not_found = 0;
1464 break;
1465 }
1466 }
1467 if (not_found)
1468 {
1469 if (size < space)
1470 {
1471 local_j[size] = cols[indx];
1472 local_data[size++] = values[indx];
1473 }
1474 else
1475 {
1476 tmp_j[tmp_indx] = cols[indx];
1477 tmp_data[tmp_indx++] = values[indx];
1478 }
1479 }
1480 not_found = 1;
1481 indx++;
1482 }
1483
1484 row_length[row_local] = size+tmp_indx;
1485
1486 if (tmp_indx)
1487 {
1488 aux_j[row_local] = hypre_TReAlloc(aux_j[row_local],HYPRE_BigInt,
1489 size+tmp_indx);
1490 aux_data[row_local] = hypre_TReAlloc(aux_data[row_local],
1491 double,size+tmp_indx);
1492 row_space[row_local] = size+tmp_indx;
1493 local_j = aux_j[row_local];
1494 local_data = aux_data[row_local];
1495 }
1496
1497 cnt = size;
1498
1499 for (i=0; i < tmp_indx; i++)
1500 {
1501 local_j[cnt] = tmp_j[i];
1502 local_data[cnt++] = tmp_data[i];
1503 }
1504
1505 if (tmp_j)
1506 {
1507 hypre_TFree(tmp_j);
1508 hypre_TFree(tmp_data);
1509 }
1510 }
1511 else /* insert immediately into data in ParCSRMatrix structure */
1512 {
1513 offd_indx = hypre_AuxParCSRMatrixIndxOffd(aux_matrix)[row_local];
1514 diag_indx = hypre_AuxParCSRMatrixIndxDiag(aux_matrix)[row_local];
1515 diag = hypre_ParCSRMatrixDiag(par_matrix);
1516 diag_i = hypre_CSRMatrixI(diag);
1517 diag_j = hypre_CSRMatrixJ(diag);
1518 diag_data = hypre_CSRMatrixData(diag);
1519 offd = hypre_ParCSRMatrixOffd(par_matrix);
1520 offd_i = hypre_CSRMatrixI(offd);
1521 if (num_procs > 1)
1522 {
1523 offd_j = hypre_CSRMatrixJ(offd);
1524 offd_data = hypre_CSRMatrixData(offd);
1525 aux_offd_j = hypre_AuxParCSRMatrixAuxOffdJ(aux_matrix);
1526 }
1527 cnt_diag = diag_indx;
1528 cnt_offd = offd_indx;
1529 diag_space = diag_i[row_local+1];
1530 offd_space = offd_i[row_local+1];
1531 not_found = 1;
1532 for (i=0; i < n; i++)
1533 {
1534 if (cols[indx] < col_0 || cols[indx] > col_n)
1535 /* insert into offd */
1536 {
1537 for (j=offd_i[row_local]; j < offd_indx; j++)
1538 {
1539 if (aux_offd_j[j] == cols[indx])
1540 {
1541 offd_data[j] += values[indx];
1542 not_found = 0;
1543 break;
1544 }
1545 }
1546 if (not_found)
1547 {
1548 if (cnt_offd < offd_space)
1549 {
1550 aux_offd_j[cnt_offd] = cols[indx];
1551 offd_data[cnt_offd++] = values[indx];
1552 }
1553 else
1554 {
1555 hypre_error(HYPRE_ERROR_GENERIC);
1556#ifdef HYPRE_LONG_LONG
1557 printf("Error in row %lld ! Too many elements!\n",
1558 row);
1559#else
1560 printf("Error in row %d ! Too many elements!\n",
1561 row);
1562#endif
1563 /* return 1;*/
1564 return hypre_error_flag;
1565 }
1566 }
1567 not_found = 1;
1568 }
1569 else /* insert into diag */
1570 {
1571 for (j=diag_i[row_local]; j < diag_indx; j++)
1572 {
1573 col_j = (int)(cols[indx]-col_0);
1574 if (diag_j[j] == col_j)
1575 {
1576 diag_data[j] += values[indx];
1577 not_found = 0;
1578 break;
1579 }
1580 }
1581 if (not_found)
1582 {
1583 if (cnt_diag < diag_space)
1584 {
1585 diag_j[cnt_diag] = col_j;
1586 diag_data[cnt_diag++] = values[indx];
1587 }
1588 else
1589 {
1590 hypre_error(HYPRE_ERROR_GENERIC);
1591#ifdef HYPRE_LONG_LONG
1592 printf("Error in row %lld ! Too many elements !\n",
1593 row);
1594#else
1595 printf("Error in row %d ! Too many elements !\n",
1596 row);
1597#endif
1598 /* return 1; */
1599 return hypre_error_flag;
1600 }
1601 }
1602 not_found = 1;
1603 }
1604 indx++;
1605 }
1606
1607 hypre_AuxParCSRMatrixIndxDiag(aux_matrix)[row_local] = cnt_diag;
1608 hypre_AuxParCSRMatrixIndxOffd(aux_matrix)[row_local] = cnt_offd;
1609
1610 }
1611 }
1612 /* not my row */
1613 else
1614 {
1615 current_num_elmts
1616 = hypre_AuxParCSRMatrixCurrentNumElmts(aux_matrix);
1617 max_off_proc_elmts
1618 = hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix);
1619 off_proc_i_indx = hypre_AuxParCSRMatrixOffProcIIndx(aux_matrix);
1620 off_proc_i = hypre_AuxParCSRMatrixOffProcI(aux_matrix);
1621 off_proc_j = hypre_AuxParCSRMatrixOffProcJ(aux_matrix);
1622 off_proc_data = hypre_AuxParCSRMatrixOffProcData(aux_matrix);
1623
1624 if (!max_off_proc_elmts)
1625 {
1626 max_off_proc_elmts = hypre_max(n,1000);
1627 hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix) =
1628 max_off_proc_elmts;
1629 hypre_AuxParCSRMatrixOffProcI(aux_matrix)
1630 = hypre_CTAlloc(HYPRE_BigInt,2*max_off_proc_elmts);
1631 hypre_AuxParCSRMatrixOffProcJ(aux_matrix)
1632 = hypre_CTAlloc(HYPRE_BigInt,max_off_proc_elmts);
1633 hypre_AuxParCSRMatrixOffProcData(aux_matrix)
1634 = hypre_CTAlloc(double,max_off_proc_elmts);
1635 off_proc_i = hypre_AuxParCSRMatrixOffProcI(aux_matrix);
1636 off_proc_j = hypre_AuxParCSRMatrixOffProcJ(aux_matrix);
1637 off_proc_data = hypre_AuxParCSRMatrixOffProcData(aux_matrix);
1638 }
1639 else if (current_num_elmts + n > max_off_proc_elmts)
1640 {
1641 max_off_proc_elmts += 3*n;
1642 off_proc_i = hypre_TReAlloc(off_proc_i,HYPRE_BigInt,2*max_off_proc_elmts);
1643 off_proc_j = hypre_TReAlloc(off_proc_j,HYPRE_BigInt,max_off_proc_elmts);
1644 off_proc_data = hypre_TReAlloc(off_proc_data,double,
1645 max_off_proc_elmts);
1646 hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix)
1647 = max_off_proc_elmts;
1648 hypre_AuxParCSRMatrixOffProcI(aux_matrix) = off_proc_i;
1649 hypre_AuxParCSRMatrixOffProcJ(aux_matrix) = off_proc_j;
1650 hypre_AuxParCSRMatrixOffProcData(aux_matrix) = off_proc_data;
1651 }
1652 off_proc_i[off_proc_i_indx++] = -row-1;
1653 off_proc_i[off_proc_i_indx++] = n;
1654 for (i=0; i < n; i++)
1655 {
1656 off_proc_j[current_num_elmts] = cols[indx];
1657 off_proc_data[current_num_elmts++] = values[indx++];
1658 }
1659 hypre_AuxParCSRMatrixOffProcIIndx(aux_matrix) = off_proc_i_indx;
1660 hypre_AuxParCSRMatrixCurrentNumElmts(aux_matrix)
1661 = current_num_elmts;
1662 }
1663 }
1664 }
1665 return hypre_error_flag;
1666
1667}
1668
1669/******************************************************************************
1670 *
1671 * hypre_IJMatrixAssembleParCSR
1672 *
1673 * assembles IJMatrix from AuxParCSRMatrix auxiliary structure
1674 *****************************************************************************/
1675int
1676hypre_IJMatrixAssembleParCSR(hypre_IJMatrix *matrix)
1677{
1678 MPI_Comm comm = hypre_IJMatrixComm(matrix);
1679 hypre_ParCSRMatrix *par_matrix = hypre_IJMatrixObject(matrix);
1680 hypre_AuxParCSRMatrix *aux_matrix = hypre_IJMatrixTranslator(matrix);
1681 HYPRE_BigInt *row_partitioning = hypre_IJMatrixRowPartitioning(matrix);
1682 HYPRE_BigInt *col_partitioning = hypre_IJMatrixColPartitioning(matrix);
1683
1684 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(par_matrix);
1685 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(par_matrix);
1686 int *diag_i = hypre_CSRMatrixI(diag);
1687 int *offd_i = hypre_CSRMatrixI(offd);
1688 int *diag_j;
1689 int *offd_j = NULL;
1690 double *diag_data;
1691 double *offd_data = NULL;
1692 int i, j, j0;
1693 int num_cols_offd;
1694 int *diag_pos;
1695 HYPRE_BigInt *col_map_offd;
1696 int *row_length;
1697 HYPRE_BigInt **aux_j;
1698 double **aux_data;
1699 int my_id, num_procs;
1700 int num_rows;
1701 int i_diag, i_offd;
1702 HYPRE_BigInt *local_j;
1703 double *local_data;
1704 HYPRE_BigInt col_0, col_n;
1705 int nnz_offd;
1706 HYPRE_BigInt *aux_offd_j;
1707 HYPRE_BigInt *tmp_j;
1708 double temp;
1709#ifdef HYPRE_NO_GLOBAL_PARTITION
1710 HYPRE_BigInt base = hypre_IJMatrixGlobalFirstCol(matrix);
1711#else
1712 HYPRE_BigInt base = col_partitioning[0];
1713#endif
1714 int off_proc_i_indx;
1715 int max_off_proc_elmts;
1716 int current_num_elmts;
1717 HYPRE_BigInt *off_proc_i;
1718 HYPRE_BigInt *off_proc_j;
1719 double *off_proc_data;
1720 int offd_proc_elmts;
1721
1722 if (aux_matrix)
1723 {
1724 off_proc_i_indx = hypre_AuxParCSRMatrixOffProcIIndx(aux_matrix);
1725 MPI_Allreduce(&off_proc_i_indx,&offd_proc_elmts,1,MPI_INT, MPI_SUM,comm);
1726 if (offd_proc_elmts)
1727 {
1728 max_off_proc_elmts=hypre_AuxParCSRMatrixMaxOffProcElmts(aux_matrix);
1729 current_num_elmts=hypre_AuxParCSRMatrixCurrentNumElmts(aux_matrix);
1730 off_proc_i=hypre_AuxParCSRMatrixOffProcI(aux_matrix);
1731 off_proc_j=hypre_AuxParCSRMatrixOffProcJ(aux_matrix);
1732 off_proc_data=hypre_AuxParCSRMatrixOffProcData(aux_matrix);
1733 hypre_IJMatrixAssembleOffProcValsParCSR(matrix,off_proc_i_indx,
1734 max_off_proc_elmts, current_num_elmts, off_proc_i,
1735 off_proc_j, off_proc_data);
1736 }
1737 }
1738
1739 if (hypre_IJMatrixAssembleFlag(matrix) == 0)
1740 {
1741 MPI_Comm_size(comm, &num_procs);
1742 MPI_Comm_rank(comm, &my_id);
1743#ifdef HYPRE_NO_GLOBAL_PARTITION
1744 num_rows = (int)(row_partitioning[1] - row_partitioning[0]);
1745 col_0 = col_partitioning[0];
1746 col_n = col_partitioning[1]-1;
1747#else
1748 num_rows = (int)(row_partitioning[my_id+1] - row_partitioning[my_id]);
1749 col_0 = col_partitioning[my_id];
1750 col_n = col_partitioning[my_id+1]-1;
1751#endif
1752 /* Got till here !! Warning!!! Check this carefully!! Problems!! ***/
1753
1754 /* move data into ParCSRMatrix if not there already */
1755 if (hypre_AuxParCSRMatrixNeedAux(aux_matrix))
1756 {
1757 aux_j = hypre_AuxParCSRMatrixAuxJ(aux_matrix);
1758 aux_data = hypre_AuxParCSRMatrixAuxData(aux_matrix);
1759 row_length = hypre_AuxParCSRMatrixRowLength(aux_matrix);
1760 diag_pos = hypre_CTAlloc(int, num_rows);
1761 i_diag = 0;
1762 i_offd = 0;
1763 for (i=0; i < num_rows; i++)
1764 {
1765 local_j = aux_j[i];
1766 local_data = aux_data[i];
1767 diag_pos[i] = -1;
1768 for (j=0; j < row_length[i]; j++)
1769 {
1770 if (local_j[j] < col_0 || local_j[j] > col_n)
1771 i_offd++;
1772 else
1773 {
1774 i_diag++;
1775 if ((int)(local_j[j]-col_0) == i) diag_pos[i] = j;
1776 }
1777 }
1778 diag_i[i+1] = i_diag;
1779 offd_i[i+1] = i_offd;
1780 }
1781
1782 if (hypre_CSRMatrixJ(diag))
1783 hypre_TFree(hypre_CSRMatrixJ(diag));
1784 if (hypre_CSRMatrixData(diag))
1785 hypre_TFree(hypre_CSRMatrixData(diag));
1786 if (hypre_CSRMatrixJ(offd))
1787 hypre_TFree(hypre_CSRMatrixJ(offd));
1788 if (hypre_CSRMatrixData(offd))
1789 hypre_TFree(hypre_CSRMatrixData(offd));
1790 diag_j = hypre_CTAlloc(int,i_diag);
1791 diag_data = hypre_CTAlloc(double,i_diag);
1792 if (i_offd > 0)
1793 {
1794 offd_j = hypre_CTAlloc(int,i_offd);
1795 offd_data = hypre_CTAlloc(double,i_offd);
1796 aux_offd_j = hypre_CTAlloc(HYPRE_BigInt, i_offd);
1797 tmp_j = hypre_CTAlloc(HYPRE_BigInt, i_offd);
1798 }
1799
1800 i_diag = 0;
1801 i_offd = 0;
1802 for (i=0; i < num_rows; i++)
1803 {
1804 local_j = aux_j[i];
1805 local_data = aux_data[i];
1806 if (diag_pos[i] > -1)
1807 {
1808 diag_j[i_diag] = (int)(local_j[diag_pos[i]] - col_0);
1809 diag_data[i_diag++] = local_data[diag_pos[i]];
1810 }
1811 for (j=0; j < row_length[i]; j++)
1812 {
1813 if (local_j[j] < col_0 || local_j[j] > col_n)
1814 {
1815 aux_offd_j[i_offd] = local_j[j];
1816 tmp_j[i_offd] = aux_offd_j[i_offd];
1817 offd_data[i_offd++] = local_data[j];
1818 }
1819 else if (j != diag_pos[i])
1820 {
1821 diag_j[i_diag] = (int)(local_j[j] - col_0);
1822 diag_data[i_diag++] = local_data[j];
1823 }
1824 }
1825 }
1826 hypre_CSRMatrixJ(diag) = diag_j;
1827 hypre_CSRMatrixData(diag) = diag_data;
1828 hypre_CSRMatrixNumNonzeros(diag) = diag_i[num_rows];
1829 if (i_offd > 0)
1830 {
1831 hypre_CSRMatrixJ(offd) = offd_j;
1832 hypre_CSRMatrixData(offd) = offd_data;
1833 }
1834 hypre_CSRMatrixNumNonzeros(offd) = offd_i[num_rows];
1835 hypre_TFree(diag_pos);
1836 }
1837 else
1838 {
1839 /* move diagonal element into first space */
1840
1841 diag_j = hypre_CSRMatrixJ(diag);
1842 diag_data = hypre_CSRMatrixData(diag);
1843 for (i=0; i < num_rows; i++)
1844 {
1845 j0 = diag_i[i];
1846 for (j=j0; j < diag_i[i+1]; j++)
1847 {
1848 if (diag_j[j] == i)
1849 {
1850 temp = diag_data[j0];
1851 diag_data[j0] = diag_data[j];
1852 diag_data[j] = temp;
1853 diag_j[j] = diag_j[j0];
1854 diag_j[j0] = i;
1855 }
1856 }
1857 }
1858
1859 offd_j = hypre_CSRMatrixJ(offd);
1860 }
1861
1862 /* generate the nonzero rows inside offd and diag by calling */
1863
1864 hypre_CSRMatrixSetRownnz(diag);
1865 hypre_CSRMatrixSetRownnz(offd);
1866 nnz_offd = offd_i[num_rows];
1867
1868 /* generate col_map_offd */
1869 if (nnz_offd)
1870 {
1871 hypre_BigQsort0(tmp_j,0,nnz_offd-1);
1872 num_cols_offd = 1;
1873 for (i=0; i < nnz_offd-1; i++)
1874 {
1875 if (tmp_j[i+1] > tmp_j[i])
1876 {
1877 tmp_j[num_cols_offd] = tmp_j[i+1];
1878 num_cols_offd++;
1879 }
1880 }
1881 col_map_offd = hypre_CTAlloc(HYPRE_BigInt,num_cols_offd);
1882 for (i=0; i < num_cols_offd; i++)
1883 col_map_offd[i] = tmp_j[i];
1884
1885 for (i=0; i < nnz_offd; i++)
1886 {
1887 offd_j[i]=hypre_BigBinarySearch(col_map_offd,aux_offd_j[i],num_cols_offd);
1888 }
1889 if (base)
1890 {
1891 for (i=0; i < num_cols_offd; i++)
1892 col_map_offd[i] -= base;
1893 }
1894 hypre_ParCSRMatrixColMapOffd(par_matrix) = col_map_offd;
1895 hypre_CSRMatrixNumCols(offd) = num_cols_offd;
1896 hypre_TFree(aux_offd_j);
1897 hypre_TFree(tmp_j);
1898 }
1899 hypre_AuxParCSRMatrixDestroy(aux_matrix);
1900 hypre_IJMatrixTranslator(matrix) = NULL;
1901 hypre_IJMatrixAssembleFlag(matrix) = 1;
1902 }
1903 return hypre_error_flag;
1904}
1905
1906/******************************************************************************
1907 *
1908 * hypre_IJMatrixDestroyParCSR
1909 *
1910 * frees an IJMatrix
1911 *
1912 *****************************************************************************/
1913int
1914hypre_IJMatrixDestroyParCSR(hypre_IJMatrix *matrix)
1915{
1916 hypre_ParCSRMatrixDestroy(hypre_IJMatrixObject(matrix));
1917 hypre_AuxParCSRMatrixDestroy(hypre_IJMatrixTranslator(matrix));
1918
1919 return hypre_error_flag;
1920}
1921
1922
1923
1924
1925
1926
1927
1928/******************************************************************************
1929 *
1930 * hypre_IJMatrixAssembleOffProcValsParCSR
1931 *
1932 * This is for handling set and get values calls to off-proc. entries -
1933 * it is called from matrix assemble. There is an alternate version for
1934 * when the assumed partition is being used.
1935 *
1936 *****************************************************************************/
1937
1938#ifndef HYPRE_NO_GLOBAL_PARTITION
1939
1940int
1941hypre_IJMatrixAssembleOffProcValsParCSR( hypre_IJMatrix *matrix,
1942 int off_proc_i_indx,
1943 int max_off_proc_elmts,
1944 int current_num_elmts,
1945 HYPRE_BigInt *off_proc_i,
1946 HYPRE_BigInt *off_proc_j,
1947 double *off_proc_data)
1948{
1949 MPI_Comm comm = hypre_IJMatrixComm(matrix);
1950 MPI_Request *requests;
1951 MPI_Status *status;
1952 int i, ii, j, j2, jj, n;
1953 HYPRE_BigInt row;
1954 int iii, iid, indx, ip;
1955 int proc_id, num_procs, my_id;
1956 int num_sends, num_sends3;
1957 int num_recvs;
1958 int num_requests;
1959 int vec_start, vec_len;
1960 int *send_procs;
1961 int *chunks;
1962 HYPRE_BigInt *send_i;
1963 int *send_map_starts;
1964 int *dbl_send_map_starts;
1965 int *recv_procs;
1966 int *recv_chunks;
1967 HYPRE_BigInt *recv_i;
1968 int *recv_vec_starts;
1969 int *dbl_recv_vec_starts;
1970 int *info;
1971 int *int_buffer;
1972 int *proc_id_mem;
1973 HYPRE_BigInt *partitioning;
1974 int *displs;
1975 int *recv_buf;
1976 double *send_data;
1977 double *recv_data;
1978
1979 MPI_Comm_size(comm,&num_procs);
1980 MPI_Comm_rank(comm, &my_id);
1981 partitioning = hypre_IJMatrixRowPartitioning(matrix);
1982
1983 info = hypre_CTAlloc(int,num_procs);
1984 chunks = hypre_CTAlloc(int,num_procs);
1985 proc_id_mem = hypre_CTAlloc(int,off_proc_i_indx/2);
1986 j=0;
1987 for (i=0; i < off_proc_i_indx; i++)
1988 {
1989 row = off_proc_i[i++];
1990 if (row < 0) row = -row-1;
1991 n = (int) off_proc_i[i];
1992 proc_id = hypre_FindProc(partitioning,row,num_procs);
1993 proc_id_mem[j++] = proc_id;
1994 info[proc_id] += n;
1995 chunks[proc_id]++;
1996 }
1997
1998 /* determine send_procs and amount of data to be sent */
1999 num_sends = 0;
2000 for (i=0; i < num_procs; i++)
2001 {
2002 if (info[i])
2003 {
2004 num_sends++;
2005 }
2006 }
2007 send_procs = hypre_CTAlloc(int,num_sends);
2008 send_map_starts = hypre_CTAlloc(int,num_sends+1);
2009 dbl_send_map_starts = hypre_CTAlloc(int,num_sends+1);
2010 num_sends3 = 3*num_sends;
2011 int_buffer = hypre_CTAlloc(int,3*num_sends);
2012 j = 0;
2013 j2 = 0;
2014 send_map_starts[0] = 0;
2015 dbl_send_map_starts[0] = 0;
2016 for (i=0; i < num_procs; i++)
2017 {
2018 if (info[i])
2019 {
2020 send_procs[j++] = i;
2021 send_map_starts[j] = send_map_starts[j-1]+2*chunks[i]+info[i];
2022 dbl_send_map_starts[j] = dbl_send_map_starts[j-1]+info[i];
2023 int_buffer[j2++] = i;
2024 int_buffer[j2++] = chunks[i];
2025 int_buffer[j2++] = info[i];
2026 }
2027 }
2028
2029 hypre_TFree(chunks);
2030
2031 MPI_Allgather(&num_sends3,1,MPI_INT,info,1,MPI_INT,comm);
2032
2033 displs = hypre_CTAlloc(int, num_procs+1);
2034 displs[0] = 0;
2035 for (i=1; i < num_procs+1; i++)
2036 displs[i] = displs[i-1]+info[i-1];
2037 recv_buf = hypre_CTAlloc(int, displs[num_procs]);
2038
2039 MPI_Allgatherv(int_buffer,num_sends3,MPI_INT,recv_buf,info,displs,
2040 MPI_INT,comm);
2041
2042 hypre_TFree(int_buffer);
2043 hypre_TFree(info);
2044
2045 /* determine recv procs and amount of data to be received */
2046 num_recvs = 0;
2047 for (j=0; j < displs[num_procs]; j+=3)
2048 {
2049 if (recv_buf[j] == my_id)
2050 num_recvs++;
2051 }
2052
2053 recv_procs = hypre_CTAlloc(int,num_recvs);
2054 recv_chunks = hypre_CTAlloc(int,num_recvs);
2055 recv_vec_starts = hypre_CTAlloc(int,num_recvs+1);
2056 dbl_recv_vec_starts = hypre_CTAlloc(int,num_recvs+1);
2057
2058 j2 = 0;
2059 recv_vec_starts[0] = 0;
2060 dbl_recv_vec_starts[0] = 0;
2061 for (i=0; i < num_procs; i++)
2062 {
2063 for (j=displs[i]; j < displs[i+1]; j+=3)
2064 {
2065 if (recv_buf[j] == my_id)
2066 {
2067 recv_procs[j2] = i;
2068 recv_chunks[j2++] = recv_buf[j+1];
2069 recv_vec_starts[j2] = recv_vec_starts[j2-1]+2*recv_buf[j+1]
2070 +recv_buf[j+2];
2071 dbl_recv_vec_starts[j2] = dbl_recv_vec_starts[j2-1]+recv_buf[j+2];
2072 }
2073 if (j2 == num_recvs) break;
2074 }
2075 }
2076 hypre_TFree(recv_buf);
2077 hypre_TFree(displs);
2078
2079 /* set up data to be sent to send procs */
2080 /* send_i contains for each send proc : row no., no. of elmts and column
2081 indices, send_data contains corresponding values */
2082
2083 send_i = hypre_CTAlloc(HYPRE_BigInt,send_map_starts[num_sends]);
2084 send_data = hypre_CTAlloc(double,dbl_send_map_starts[num_sends]);
2085 recv_i = hypre_CTAlloc(HYPRE_BigInt,recv_vec_starts[num_recvs]);
2086 recv_data = hypre_CTAlloc(double,dbl_recv_vec_starts[num_recvs]);
2087
2088 j=0;
2089 jj=0;
2090 for (i=0; i < off_proc_i_indx; i++)
2091 {
2092 row = off_proc_i[i++];
2093 n = (int) off_proc_i[i];
2094 proc_id = proc_id_mem[i/2];
2095 indx = hypre_BinarySearch(send_procs,proc_id,num_sends);
2096 iii = send_map_starts[indx];
2097 iid = dbl_send_map_starts[indx];
2098 send_i[iii++] = row;
2099 send_i[iii++] = (HYPRE_BigInt) n;
2100 for (ii = 0; ii < n; ii++)
2101 {
2102 send_i[iii++] = off_proc_j[jj];
2103 send_data[iid++] = off_proc_data[jj++];
2104 }
2105 send_map_starts[indx] = iii;
2106 dbl_send_map_starts[indx] = iid;
2107 }
2108
2109 hypre_TFree(proc_id_mem);
2110
2111 for (i=num_sends; i > 0; i--)
2112 {
2113 send_map_starts[i] = send_map_starts[i-1];
2114 dbl_send_map_starts[i] = dbl_send_map_starts[i-1];
2115 }
2116 send_map_starts[0] = 0;
2117 dbl_send_map_starts[0] = 0;
2118
2119 num_requests = num_recvs+num_sends;
2120
2121 requests = hypre_CTAlloc(MPI_Request, num_requests);
2122 status = hypre_CTAlloc(MPI_Status, num_requests);
2123
2124 j=0;
2125 for (i=0; i < num_recvs; i++)
2126 {
2127 vec_start = recv_vec_starts[i];
2128 vec_len = recv_vec_starts[i+1] - vec_start;
2129 ip = recv_procs[i];
2130 MPI_Irecv(&recv_i[vec_start], vec_len, MPI_HYPRE_BIG_INT, ip, 0, comm,
2131 &requests[j++]);
2132 }
2133
2134 for (i=0; i < num_sends; i++)
2135 {
2136 vec_start = send_map_starts[i];
2137 vec_len = send_map_starts[i+1] - vec_start;
2138 ip = send_procs[i];
2139 MPI_Isend(&send_i[vec_start], vec_len, MPI_HYPRE_BIG_INT, ip, 0, comm,
2140 &requests[j++]);
2141 }
2142
2143 if (num_requests)
2144 {
2145 MPI_Waitall(num_requests, requests, status);
2146 }
2147
2148 j=0;
2149 for (i=0; i < num_recvs; i++)
2150 {
2151 vec_start = dbl_recv_vec_starts[i];
2152 vec_len = dbl_recv_vec_starts[i+1] - vec_start;
2153 ip = recv_procs[i];
2154 MPI_Irecv(&recv_data[vec_start], vec_len, MPI_DOUBLE, ip, 0, comm,
2155 &requests[j++]);
2156 }
2157
2158 for (i=0; i < num_sends; i++)
2159 {
2160 vec_start = dbl_send_map_starts[i];
2161 vec_len = dbl_send_map_starts[i+1] - vec_start;
2162 ip = send_procs[i];
2163 MPI_Isend(&send_data[vec_start], vec_len, MPI_DOUBLE, ip, 0, comm,
2164 &requests[j++]);
2165 }
2166
2167 if (num_requests)
2168 {
2169 MPI_Waitall(num_requests, requests, status);
2170 hypre_TFree(requests);
2171 hypre_TFree(status);
2172 }
2173
2174 hypre_TFree(send_i);
2175 hypre_TFree(send_data);
2176 hypre_TFree(send_procs);
2177 hypre_TFree(send_map_starts);
2178 hypre_TFree(dbl_send_map_starts);
2179 hypre_TFree(recv_procs);
2180 hypre_TFree(recv_vec_starts);
2181 hypre_TFree(dbl_recv_vec_starts);
2182
2183 j = 0;
2184 j2 = 0;
2185 for (i=0; i < num_recvs; i++)
2186 {
2187 for (ii=0; ii < recv_chunks[i]; ii++)
2188 {
2189 row = recv_i[j];
2190 if (row < 0)
2191 {
2192 row = -row-1;
2193 hypre_IJMatrixAddToValuesParCSR(matrix,1,(int*)&recv_i[j+1],&row,
2194 &recv_i[j+2],&recv_data[j2]);
2195 j2 += recv_i[j+1];
2196 j += recv_i[j+1]+2;
2197 }
2198 else
2199 {
2200 hypre_IJMatrixSetValuesParCSR(matrix,1,(int*)&recv_i[j+1],&row,
2201 &recv_i[j+2],&recv_data[j2]);
2202 j2 += recv_i[j+1];
2203 j += recv_i[j+1]+2;
2204 }
2205 }
2206 }
2207 hypre_TFree(recv_chunks);
2208 hypre_TFree(recv_i);
2209 hypre_TFree(recv_data);
2210
2211 return hypre_error_flag;
2212
2213}
2214
2215#else
2216
2217
2218/* assumed partition version */
2219
2220int
2221hypre_IJMatrixAssembleOffProcValsParCSR( hypre_IJMatrix *matrix,
2222 int off_proc_i_indx,
2223 int max_off_proc_elmts,
2224 int current_num_elmts,
2225 HYPRE_BigInt *off_proc_i,
2226 HYPRE_BigInt *off_proc_j,
2227 double *off_proc_data)
2228{
2229
2230 MPI_Comm comm = hypre_IJMatrixComm(matrix);
2231
2232 int i, j, k, in_i;
2233 int myid;
2234
2235 int proc_id, last_proc, prev_id, tmp_id;
2236 int max_response_size;
2237 HYPRE_BigInt global_num_cols;
2238 int ex_num_contacts = 0, num_rows = 0;
2239 HYPRE_BigInt range_start, range_end;
2240 int num_elements;
2241 int storage;
2242 int indx;
2243 HYPRE_BigInt row;
2244 int num_ranges;
2245 int num_recvs;
2246 int counter;
2247 HYPRE_BigInt upper_bound;
2248 int num_real_procs;
2249
2250
2251 HYPRE_BigInt *row_list=NULL;
2252 int *row_list_num_elements=NULL;
2253 int *a_proc_id=NULL, *orig_order=NULL;
2254 int *real_proc_id = NULL, *us_real_proc_id = NULL;
2255 int *ex_contact_procs = NULL, *ex_contact_vec_starts = NULL;
2256 HYPRE_BigInt *ex_contact_buf = NULL;
2257 int *recv_starts=NULL;
2258 HYPRE_BigInt *response_buf = NULL;
2259 int *response_buf_starts=NULL;
2260 int *num_rows_per_proc = NULL, *num_elements_total = NULL;
2261
2262 int obj_size_bytes, int_size, double_size, big_int_size;
2263 int tmp_int;
2264 HYPRE_BigInt tmp_big_int;
2265 HYPRE_BigInt *col_ptr;
2266 HYPRE_BigInt *big_int_data = NULL;
2267 int big_int_data_size = 0, double_data_size = 0;
2268
2269 void *void_contact_buf = NULL;
2270 void *index_ptr;
2271 void *recv_data_ptr;
2272
2273 double tmp_double;
2274 double *col_data_ptr;
2275 double *double_data = NULL;
2276
2277 hypre_DataExchangeResponse response_obj1, response_obj2;
2278 hypre_ProcListElements send_proc_obj;
2279
2280 hypre_IJAssumedPart *apart;
2281
2282 hypre_ParCSRMatrix *par_matrix = hypre_IJMatrixObject(matrix);
2283
2284 MPI_Comm_rank(comm, &myid);
2285 global_num_cols = hypre_IJMatrixGlobalNumCols(matrix);
2286
2287 /* verify that we have created the assumed partition */
2288 if (hypre_ParCSRMatrixAssumedPartition(par_matrix) == NULL)
2289 {
2290 hypre_ParCSRMatrixCreateAssumedPartition(par_matrix);
2291 }
2292
2293 apart = hypre_ParCSRMatrixAssumedPartition(par_matrix);
2294
2295 num_rows = off_proc_i_indx/2;
2296
2297 row_list = hypre_CTAlloc(HYPRE_BigInt, num_rows);
2298 row_list_num_elements = hypre_CTAlloc(int, num_rows);
2299 a_proc_id = hypre_CTAlloc(int, num_rows);
2300 orig_order = hypre_CTAlloc(int, num_rows);
2301 real_proc_id = hypre_CTAlloc(int, num_rows);
2302
2303 /* get the assumed processor id for each row */
2304 if (num_rows > 0 )
2305 {
2306 for (i=0; i < num_rows; i++)
2307 {
2308 row = off_proc_i[i*2];
2309 if (row < 0) row = -row-1;
2310 row_list[i] = row;
2311 row_list_num_elements[i] = (int) off_proc_i[i*2+1];
2312
2313 hypre_GetAssumedPartitionProcFromRow (comm, row, global_num_cols, &proc_id);
2314 a_proc_id[i] = proc_id;
2315 orig_order[i] = i;
2316 }
2317
2318 /* now we need to find the actual order of each row - sort on row -
2319 this will result in proc ids sorted also...*/
2320
2321 hypre_BigQsortb2i(row_list, a_proc_id, orig_order, 0, num_rows -1);
2322
2323
2324 /* calculate the number of contacts */
2325 ex_num_contacts = 1;
2326 last_proc = a_proc_id[0];
2327 for (i=1; i < num_rows; i++)
2328 {
2329 if (a_proc_id[i] > last_proc)
2330 {
2331 ex_num_contacts++;
2332 last_proc = a_proc_id[i];
2333 }
2334 }
2335
2336 }
2337
2338
2339 /* now we will go through a create a contact list - need to contact
2340 assumed processors and find out who the actual row owner is - we
2341 will contact with a range (2 numbers) */
2342
2343 ex_contact_procs = hypre_CTAlloc(int, ex_num_contacts);
2344 ex_contact_vec_starts = hypre_CTAlloc(int, ex_num_contacts+1);
2345 ex_contact_buf = hypre_CTAlloc(HYPRE_BigInt, ex_num_contacts*2);
2346
2347 counter = 0;
2348 range_end = -1;
2349 for (i=0; i< num_rows; i++)
2350 {
2351 if (row_list[i] > range_end)
2352 {
2353 /* assumed proc */
2354 proc_id = a_proc_id[i];
2355
2356 /* end of prev. range */
2357 if (counter > 0) ex_contact_buf[counter*2 - 1] = row_list[i-1];
2358
2359 /*start new range*/
2360 ex_contact_procs[counter] = proc_id;
2361 ex_contact_vec_starts[counter] = counter*2;
2362 ex_contact_buf[counter*2] = row_list[i];
2363 counter++;
2364
2365 hypre_GetAssumedPartitionRowRange(comm, proc_id, global_num_cols,
2366 &range_start, &range_end);
2367
2368
2369 }
2370 }
2371 /*finish the starts*/
2372 ex_contact_vec_starts[counter] = counter*2;
2373 /*finish the last range*/
2374 if (counter > 0)
2375 ex_contact_buf[counter*2 - 1] = row_list[num_rows - 1];
2376
2377
2378 /*don't allocate space for responses */
2379
2380
2381 /*create response object - can use same fill response as used in the commpkg
2382 routine */
2383 response_obj1.fill_response = hypre_RangeFillResponseIJDetermineRecvProcs;
2384 response_obj1.data1 = apart; /* this is necessary so we can fill responses*/
2385 response_obj1.data2 = NULL;
2386
2387
2388 max_response_size = 6; /* 6 means we can fit 3 ranges*/
2389
2390
2391 hypre_DataExchangeList(ex_num_contacts, ex_contact_procs,
2392 ex_contact_buf, ex_contact_vec_starts, sizeof(HYPRE_BigInt),
2393 sizeof(HYPRE_BigInt), &response_obj1, max_response_size, 1,
2394 comm, (void**) &response_buf, &response_buf_starts);
2395
2396
2397 /* now response_buf contains
2398 a proc_id followed by an upper bound for the range. */
2399
2400 hypre_TFree(ex_contact_procs);
2401 hypre_TFree(ex_contact_buf);
2402 hypre_TFree(ex_contact_vec_starts);
2403
2404 hypre_TFree(a_proc_id);
2405
2406
2407 /*how many ranges were returned?*/
2408 num_ranges = response_buf_starts[ex_num_contacts];
2409 num_ranges = num_ranges/2;
2410
2411
2412 prev_id = -1;
2413 j = 0;
2414 counter = 0;
2415 num_real_procs = 0;
2416
2417
2418 /* loop through ranges - create a list of actual processor ids*/
2419 for (i=0; i<num_ranges; i++)
2420 {
2421 upper_bound = response_buf[i*2+1];
2422 counter = 0;
2423 tmp_id = response_buf[i*2];
2424
2425 /* loop through row_list entries - counting how many are in the range */
2426 while (row_list[j] <= upper_bound && j < num_rows)
2427 {
2428 real_proc_id[j] = tmp_id;
2429 j++;
2430 counter++;
2431 }
2432 if (counter > 0 && tmp_id != prev_id)
2433 {
2434 num_real_procs++;
2435 }
2436
2437 prev_id = tmp_id;
2438
2439 }
2440
2441
2442
2443 /* now we have the list of real procesors ids (real_proc_id) - and
2444 the number of distinct ones - so now we can set up data to be
2445 sent - we have int data and double data. that we will need to
2446 pack together */
2447
2448
2449 /*first find out how many rows and
2450 elements we need to send per proc - so we can do storage */
2451
2452 ex_contact_procs = hypre_CTAlloc(int, num_real_procs);
2453 num_rows_per_proc = hypre_CTAlloc(int, num_real_procs);
2454 num_elements_total = hypre_CTAlloc(int, num_real_procs);
2455
2456 counter = 0;
2457
2458 if (num_real_procs > 0 )
2459 {
2460 ex_contact_procs[0] = real_proc_id[0];
2461 num_rows_per_proc[0] = 1;
2462 num_elements_total[0] = row_list_num_elements[orig_order[0]];
2463
2464 for (i=1; i < num_rows; i++) /* loop through real procs - these are sorted
2465 (row_list is sorted also)*/
2466 {
2467 if (real_proc_id[i] == ex_contact_procs[counter]) /* same processor */
2468 {
2469 num_rows_per_proc[counter] += 1; /*another row */
2470 num_elements_total[counter] += row_list_num_elements[orig_order[i]];
2471 }
2472 else /* new processor */
2473 {
2474 counter++;
2475 ex_contact_procs[counter] = real_proc_id[i];
2476 num_rows_per_proc[counter] = 1;
2477 num_elements_total[counter] = row_list_num_elements[orig_order[i]];
2478
2479 }
2480 }
2481 }
2482
2483 /* to pack together, we need to use the largest obj. size of
2484 (int) and (double) - if these are much different, then we are
2485 wasting some storage, but I do not think that it will be a
2486 large amount since this function should not be used on really
2487 large amounts of data anyway*/
2488 int_size = sizeof(int);
2489 double_size = sizeof(double);
2490 big_int_size = sizeof(HYPRE_BigInt);
2491
2492
2493 obj_size_bytes = hypre_max(big_int_size, double_size);
2494
2495 /* set up data to be sent to send procs */
2496 /* for each proc, ex_contact_buf contains #rows, row #,
2497 no. elements, col indicies, col data, row #, no. elements, col
2498 indicies, col data, etc. */
2499
2500 /* first calculate total storage and make vec_starts arrays */
2501 storage = 0;
2502 ex_contact_vec_starts = hypre_CTAlloc(int, num_real_procs + 1);
2503 ex_contact_vec_starts[0] = -1;
2504
2505 for (i=0; i < num_real_procs; i++)
2506 {
2507 storage += 1 + 2 * num_rows_per_proc[i] + 2* num_elements_total[i];
2508 ex_contact_vec_starts[i+1] = -storage-1; /* need negative for next loop */
2509 }
2510
2511 hypre_TFree(num_elements_total);
2512
2513 void_contact_buf = hypre_MAlloc(storage*obj_size_bytes);
2514 index_ptr = void_contact_buf; /* step through with this index */
2515
2516
2517 /* for each proc: #rows, row #, no. elements,
2518 col indicies, col data, row #, no. elements, col indicies, col data, etc. */
2519
2520 /* un-sort real_proc_id - we want to access data arrays in order, so
2521 cheaper to do this*/
2522 us_real_proc_id = hypre_CTAlloc(int, num_rows);
2523 for (i=0; i < num_rows; i++)
2524 {
2525 us_real_proc_id[orig_order[i]] = real_proc_id[i];
2526 }
2527 hypre_TFree(real_proc_id);
2528
2529 counter = 0; /* index into data arrays */
2530 prev_id = -1;
2531 for (i=0; i < num_rows; i++)
2532 {
2533 proc_id = us_real_proc_id[i];
2534 row = off_proc_i[i*2]; /*can't use row list[i] - you loose the negative signs that
2535 differentiate add/set values */
2536 num_elements = row_list_num_elements[i];
2537 /* find position of this processor */
2538 indx = hypre_BinarySearch(ex_contact_procs, proc_id, num_real_procs);
2539 in_i = ex_contact_vec_starts[indx];
2540
2541 index_ptr = (void *) ((char *) void_contact_buf + in_i*obj_size_bytes);
2542
2543 if (in_i < 0) /* first time for this processor - add the number of rows to the buffer */
2544 {
2545 in_i = -in_i - 1;
2546 /* re-calc. index_ptr since in_i was negative */
2547 index_ptr = (void *) ((char *) void_contact_buf + in_i*obj_size_bytes);
2548
2549 tmp_int = num_rows_per_proc[indx];
2550 memcpy( index_ptr, &tmp_int, int_size);
2551 index_ptr = (void *) ((char *) index_ptr + obj_size_bytes);
2552
2553 in_i++;
2554 }
2555 /* add row # */
2556 memcpy( index_ptr, &row, big_int_size);
2557 index_ptr = (void *) ((char *) index_ptr + obj_size_bytes);
2558 in_i++;
2559
2560 /* add number of elements */
2561 memcpy( index_ptr, &num_elements, int_size);
2562 index_ptr = (void *) ((char *) index_ptr + obj_size_bytes);
2563 in_i++;
2564
2565 /* now add col indices */
2566 for (j=0; j< num_elements; j++)
2567 {
2568 tmp_big_int = off_proc_j[counter+j]; /* col number */
2569
2570 memcpy( index_ptr, &tmp_big_int, big_int_size);
2571 index_ptr = (void *) ((char *) index_ptr + obj_size_bytes);
2572 in_i ++;
2573
2574 }
2575
2576 /* now add data */
2577 for (j=0; j< num_elements; j++)
2578 {
2579 tmp_double = off_proc_data[counter++]; /* value */
2580
2581 memcpy( index_ptr, &tmp_double, double_size);
2582 index_ptr = (void *) ((char *) index_ptr + obj_size_bytes);
2583 in_i++;
2584
2585 }
2586
2587
2588 /* increment the indexes to keep track of where we are - we
2589 * adjust below to be actual starts*/
2590 ex_contact_vec_starts[indx] = in_i;
2591
2592 }
2593
2594 /* some clean up */
2595
2596 hypre_TFree(response_buf);
2597 hypre_TFree(response_buf_starts);
2598
2599 hypre_TFree(us_real_proc_id);
2600 hypre_TFree(orig_order);
2601 hypre_TFree(row_list);
2602 hypre_TFree(row_list_num_elements);
2603 hypre_TFree(num_rows_per_proc);
2604
2605
2606 for (i=num_real_procs; i > 0; i--)
2607 {
2608 ex_contact_vec_starts[i] = ex_contact_vec_starts[i-1];
2609 }
2610
2611 ex_contact_vec_starts[0] = 0;
2612
2613
2614 /* now send the data */
2615
2616 /***********************************/
2617
2618 /* first get the interger info in send_proc_obj */
2619
2620 /* the response we expect is just a confirmation*/
2621 response_buf = NULL;
2622 response_buf_starts = NULL;
2623
2624
2625 /*build the response object*/
2626
2627 /* use the send_proc_obj for the info kept from contacts */
2628 /*estimate inital storage allocation */
2629 send_proc_obj.length = 0;
2630 send_proc_obj.storage_length = num_real_procs + 5;
2631 send_proc_obj.id = NULL; /* don't care who send it to us */
2632 send_proc_obj.vec_starts = hypre_CTAlloc(int, send_proc_obj.storage_length + 1);
2633 send_proc_obj.vec_starts[0] = 0;
2634 send_proc_obj.element_storage_length = storage + 20;
2635 send_proc_obj.v_elements = hypre_MAlloc(obj_size_bytes*send_proc_obj.element_storage_length);
2636
2637 response_obj2.fill_response = hypre_FillResponseIJOffProcVals;
2638 response_obj2.data1 = NULL;
2639 response_obj2.data2 = &send_proc_obj;
2640
2641 max_response_size = 0;
2642
2643 hypre_DataExchangeList(num_real_procs, ex_contact_procs,
2644 void_contact_buf, ex_contact_vec_starts, obj_size_bytes,
2645 0, &response_obj2, max_response_size, 2,
2646 comm, (void **) &response_buf, &response_buf_starts);
2647
2648
2649
2650
2651
2652 hypre_TFree(response_buf);
2653 hypre_TFree(response_buf_starts);
2654
2655 hypre_TFree(ex_contact_procs);
2656 hypre_TFree(void_contact_buf);
2657 hypre_TFree(ex_contact_vec_starts);
2658
2659
2660 /* Now we can unpack the send_proc_objects and call set
2661 and add to values functions */
2662
2663 num_recvs = send_proc_obj.length;
2664
2665 /* alias */
2666 recv_data_ptr = send_proc_obj.v_elements;
2667 recv_starts = send_proc_obj.vec_starts;
2668
2669 for (i=0; i < num_recvs; i++)
2670 {
2671
2672 indx = recv_starts[i];
2673
2674 /* get the number of rows for this recv */
2675 memcpy( &num_rows, recv_data_ptr, int_size);
2676 recv_data_ptr = (void *) ((char *)recv_data_ptr + obj_size_bytes);
2677 indx++;
2678
2679
2680 for (j=0; j < num_rows; j++) /* for each row: unpack info */
2681 {
2682 /* row # */
2683 memcpy( &row, recv_data_ptr, big_int_size);
2684 recv_data_ptr = (void *) ((char *)recv_data_ptr + obj_size_bytes);
2685 indx++;
2686
2687 /* num elements for this row */
2688 memcpy( &num_elements, recv_data_ptr, int_size);
2689 recv_data_ptr = (void *) ((char *)recv_data_ptr + obj_size_bytes);
2690 indx++;
2691
2692 /* col indices */
2693 if (big_int_size == obj_size_bytes)
2694 {
2695 col_ptr = (HYPRE_BigInt *) recv_data_ptr;
2696 recv_data_ptr = (void *) ((char *)recv_data_ptr + num_elements*obj_size_bytes);
2697 }
2698 else /* copy data */
2699 {
2700 if (big_int_data_size < num_elements)
2701 {
2702 big_int_data = hypre_TReAlloc(big_int_data, HYPRE_BigInt, num_elements + 10);
2703 }
2704 for (k=0; k< num_elements; k++)
2705 {
2706 memcpy( &big_int_data[k], recv_data_ptr, big_int_size);
2707 recv_data_ptr = (void *) ((char *)recv_data_ptr + obj_size_bytes);
2708 }
2709 col_ptr = big_int_data;
2710 }
2711
2712 /* col data */
2713 if (double_size == obj_size_bytes)
2714 {
2715 col_data_ptr = (double *) recv_data_ptr;
2716 recv_data_ptr = (void *) ((char *)recv_data_ptr + num_elements*obj_size_bytes);
2717 }
2718 else /* copy data */
2719 {
2720 if (double_data_size < num_elements)
2721 {
2722 double_data = hypre_TReAlloc(double_data, double, num_elements + 10);
2723 }
2724 for (k=0; k< num_elements; k++)
2725 {
2726 memcpy( &double_data[k], recv_data_ptr, double_size);
2727 recv_data_ptr = (void *) ((char *)recv_data_ptr + obj_size_bytes);
2728 }
2729 col_data_ptr = double_data;
2730
2731 }
2732
2733 if (row < 0) /* Add */
2734 {
2735 row = -row-1;
2736 hypre_IJMatrixAddToValuesParCSR(matrix,1,&num_elements,&row,
2737 col_ptr,col_data_ptr);
2738 }
2739 else /* Set */
2740 {
2741 hypre_IJMatrixSetValuesParCSR(matrix,1,&num_elements,&row,
2742 col_ptr,col_data_ptr);
2743 }
2744 indx += (num_elements*2);
2745
2746 }
2747 }
2748 hypre_TFree(send_proc_obj.v_elements);
2749 hypre_TFree(send_proc_obj.vec_starts);
2750
2751 if (big_int_data) hypre_TFree(big_int_data);
2752 if (double_data) hypre_TFree(double_data);
2753
2754
2755
2756 return hypre_error_flag;
2757}
2758
2759#endif
2760
2761
2762/*--------------------------------------------------------------------
2763 * hypre_FillResponseIJOffProcVals
2764 * Fill response function for the previous function (2nd data exchange)
2765 *--------------------------------------------------------------------*/
2766
2767int
2768hypre_FillResponseIJOffProcVals(void *p_recv_contact_buf,
2769 int contact_size, int contact_proc, void *ro,
2770 MPI_Comm comm, void **p_send_response_buf,
2771 int *response_message_size )
2772
2773
2774{
2775 int myid;
2776 int index, count, elength;
2777
2778 int object_size;
2779 void *index_ptr;
2780
2781 hypre_DataExchangeResponse *response_obj = ro;
2782
2783 hypre_ProcListElements *send_proc_obj = response_obj->data2;
2784
2785 object_size = hypre_max(sizeof(HYPRE_BigInt), sizeof(double));
2786
2787 MPI_Comm_rank(comm, &myid );
2788
2789
2790 /*check to see if we need to allocate more space in send_proc_obj for vec starts*/
2791 if (send_proc_obj->length == send_proc_obj->storage_length)
2792 {
2793 send_proc_obj->storage_length +=20; /*add space for 20 more contact*/
2794 send_proc_obj->vec_starts = hypre_TReAlloc(send_proc_obj->vec_starts,int,
2795 send_proc_obj->storage_length + 1);
2796 }
2797
2798 /*initialize*/
2799 count = send_proc_obj->length;
2800 index = send_proc_obj->vec_starts[count]; /*this is the current number of elements*/
2801
2802 /*do we need more storage for the elements?*/
2803 if (send_proc_obj->element_storage_length < index + contact_size)
2804 {
2805 elength = hypre_max(contact_size, 100);
2806 elength += index;
2807 send_proc_obj->v_elements = hypre_ReAlloc(send_proc_obj->v_elements,
2808 elength*object_size);
2809 send_proc_obj->element_storage_length = elength;
2810 }
2811 /*populate send_proc_obj*/
2812 index_ptr = (void *) ((char *) send_proc_obj->v_elements + index*object_size);
2813
2814 memcpy(index_ptr, p_recv_contact_buf , object_size*contact_size);
2815
2816
2817 send_proc_obj->vec_starts[count+1] = index + contact_size;
2818 send_proc_obj->length++;
2819
2820
2821 /*output - no message to return (confirmation) */
2822 *response_message_size = 0;
2823
2824
2825 return hypre_error_flag;
2826
2827}
2828
2829
2830
2831
2832/*--------------------------------------------------------------------*/
2833
2834int hypre_FindProc(HYPRE_BigInt *list, HYPRE_BigInt value, int list_length)
2835{
2836 int low, high, m;
2837
2838 low = 0;
2839 high = list_length;
2840 if (value >= list[high] || value < list[low])
2841 return -1;
2842 else
2843 {
2844 while (low+1 < high)
2845 {
2846 m = (low + high) / 2;
2847 if (value < list[m])
2848 {
2849 high = m;
2850 }
2851 else if (value >= list[m])
2852 {
2853 low = m;
2854 }
2855 }
2856 return low;
2857 }
2858}
2859
Note: See TracBrowser for help on using the repository browser.