source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_rap.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: 66.3 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#include "headers.h"
17
18/*--------------------------------------------------------------------------
19 * OLD NOTES:
20 * Sketch of John's code to build RAP
21 *
22 * Uses two integer arrays icg and ifg as marker arrays
23 *
24 * icg needs to be of size n_fine; size of ia.
25 * A negative value of icg(i) indicates i is a f-point, otherwise
26 * icg(i) is the converts from fine to coarse grid orderings.
27 * Note that I belive the code assumes that if i<j and both are
28 * c-points, then icg(i) < icg(j).
29 * ifg needs to be of size n_coarse; size of irap
30 * I don't think it has meaning as either input or output.
31 *
32 * In the code, both the interpolation and restriction operator
33 * are stored row-wise in the array b. If i is a f-point,
34 * ib(i) points the row of the interpolation operator for point
35 * i. If i is a c-point, ib(i) points the row of the restriction
36 * operator for point i.
37 *
38 * In the CSR storage for rap, its guaranteed that the rows will
39 * be ordered ( i.e. ic<jc -> irap(ic) < irap(jc)) but I don't
40 * think there is a guarantee that the entries within a row will
41 * be ordered in any way except that the diagonal entry comes first.
42 *
43 * As structured now, the code requires that the size of rap be
44 * predicted up front. To avoid this, one could execute the code
45 * twice, the first time would only keep track of icg ,ifg and ka.
46 * Then you would know how much memory to allocate for rap and jrap.
47 * The second time would fill in these arrays. Actually you might
48 * be able to include the filling in of jrap into the first pass;
49 * just overestimate its size (its an integer array) and cut it
50 * back before the second time through. This would avoid some if tests
51 * in the second pass.
52 *
53 * Questions
54 * 1) parallel (PetSc) version?
55 * 2) what if we don't store R row-wise and don't
56 * even want to store a copy of it in this form
57 * temporarily?
58 *--------------------------------------------------------------------------*/
59
60hypre_BigCSRMatrix *
61hypre_ExchangeRAPData( hypre_BigCSRMatrix *RAP_int,
62 hypre_ParCSRCommPkg *comm_pkg_RT)
63{
64 int *RAP_int_i;
65 HYPRE_BigInt *RAP_int_j = NULL;
66 double *RAP_int_data = NULL;
67 int num_cols = 0;
68
69 MPI_Comm comm = hypre_ParCSRCommPkgComm(comm_pkg_RT);
70 int num_recvs = hypre_ParCSRCommPkgNumRecvs(comm_pkg_RT);
71 int *recv_procs = hypre_ParCSRCommPkgRecvProcs(comm_pkg_RT);
72 int *recv_vec_starts = hypre_ParCSRCommPkgRecvVecStarts(comm_pkg_RT);
73 int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg_RT);
74 int *send_procs = hypre_ParCSRCommPkgSendProcs(comm_pkg_RT);
75 int *send_map_starts = hypre_ParCSRCommPkgSendMapStarts(comm_pkg_RT);
76
77 hypre_BigCSRMatrix *RAP_ext;
78
79 int *RAP_ext_i;
80 HYPRE_BigInt *RAP_ext_j = NULL;
81 double *RAP_ext_data = NULL;
82
83 hypre_ParCSRCommHandle *comm_handle;
84 hypre_ParCSRCommPkg *tmp_comm_pkg;
85
86 int *jdata_recv_vec_starts;
87 int *jdata_send_map_starts;
88
89 int num_rows;
90 int num_nonzeros;
91 int i, j;
92 int num_procs;
93
94 MPI_Comm_size(comm,&num_procs);
95
96
97 RAP_ext_i = hypre_CTAlloc(int, send_map_starts[num_sends]+1);
98
99 jdata_recv_vec_starts = hypre_CTAlloc(int, num_recvs+1);
100 jdata_send_map_starts = hypre_CTAlloc(int, num_sends+1);
101
102/*--------------------------------------------------------------------------
103 * recompute RAP_int_i so that RAP_int_i[j+1] contains the number of
104 * elements of row j (to be determined through send_map_elmnts on the
105 * receiving end)
106 *--------------------------------------------------------------------------*/
107
108 if (num_recvs)
109 {
110 RAP_int_i = hypre_BigCSRMatrixI(RAP_int);
111 RAP_int_j = hypre_BigCSRMatrixJ(RAP_int);
112 RAP_int_data = hypre_BigCSRMatrixData(RAP_int);
113 num_cols = hypre_BigCSRMatrixNumCols(RAP_int);
114 }
115
116 jdata_recv_vec_starts[0] = 0;
117 for (i=0; i < num_recvs; i++)
118 {
119 jdata_recv_vec_starts[i+1] = RAP_int_i[recv_vec_starts[i+1]];
120 }
121
122 for (i=num_recvs; i > 0; i--)
123 for (j = recv_vec_starts[i]; j > recv_vec_starts[i-1]; j--)
124 RAP_int_i[j] -= RAP_int_i[j-1];
125
126/*--------------------------------------------------------------------------
127 * initialize communication
128 *--------------------------------------------------------------------------*/
129
130 if (num_recvs && num_sends)
131 comm_handle = hypre_ParCSRCommHandleCreate(12,comm_pkg_RT,
132 &RAP_int_i[1], &RAP_ext_i[1]);
133 else if (num_recvs)
134 comm_handle = hypre_ParCSRCommHandleCreate(12,comm_pkg_RT,
135 &RAP_int_i[1], NULL);
136 else if (num_sends)
137 comm_handle = hypre_ParCSRCommHandleCreate(12,comm_pkg_RT,
138 NULL, &RAP_ext_i[1]);
139
140 tmp_comm_pkg = hypre_CTAlloc(hypre_ParCSRCommPkg, 1);
141 hypre_ParCSRCommPkgComm(tmp_comm_pkg) = comm;
142 hypre_ParCSRCommPkgNumSends(tmp_comm_pkg) = num_recvs;
143 hypre_ParCSRCommPkgNumRecvs(tmp_comm_pkg) = num_sends;
144 hypre_ParCSRCommPkgSendProcs(tmp_comm_pkg) = recv_procs;
145 hypre_ParCSRCommPkgRecvProcs(tmp_comm_pkg) = send_procs;
146
147 hypre_ParCSRCommHandleDestroy(comm_handle);
148 comm_handle = NULL;
149
150/*--------------------------------------------------------------------------
151 * compute num_nonzeros for RAP_ext
152 *--------------------------------------------------------------------------*/
153
154 for (i=0; i < num_sends; i++)
155 for (j = send_map_starts[i]; j < send_map_starts[i+1]; j++)
156 RAP_ext_i[j+1] += RAP_ext_i[j];
157
158 num_rows = send_map_starts[num_sends];
159 num_nonzeros = RAP_ext_i[num_rows];
160 if (num_nonzeros)
161 {
162
163 RAP_ext_j = hypre_CTAlloc(HYPRE_BigInt, num_nonzeros);
164 RAP_ext_data = hypre_CTAlloc(double, num_nonzeros);
165 }
166
167 for (i=0; i < num_sends+1; i++)
168 {
169 jdata_send_map_starts[i] = RAP_ext_i[send_map_starts[i]];
170 }
171
172 hypre_ParCSRCommPkgRecvVecStarts(tmp_comm_pkg) = jdata_send_map_starts;
173 hypre_ParCSRCommPkgSendMapStarts(tmp_comm_pkg) = jdata_recv_vec_starts;
174
175 comm_handle = hypre_ParCSRCommHandleCreate(1,tmp_comm_pkg,RAP_int_data,
176 RAP_ext_data);
177 hypre_ParCSRCommHandleDestroy(comm_handle);
178 comm_handle = NULL;
179
180 comm_handle = hypre_ParCSRCommHandleCreate(21,tmp_comm_pkg,RAP_int_j,
181 RAP_ext_j);
182 RAP_ext = hypre_BigCSRMatrixCreate(num_rows,num_cols,num_nonzeros);
183
184 hypre_BigCSRMatrixI(RAP_ext) = RAP_ext_i;
185 if (num_nonzeros)
186 {
187 hypre_BigCSRMatrixJ(RAP_ext) = RAP_ext_j;
188 hypre_BigCSRMatrixData(RAP_ext) = RAP_ext_data;
189 }
190
191 hypre_ParCSRCommHandleDestroy(comm_handle);
192 comm_handle = NULL;
193
194 hypre_TFree(jdata_recv_vec_starts);
195 hypre_TFree(jdata_send_map_starts);
196 hypre_TFree(tmp_comm_pkg);
197
198 return RAP_ext;
199}
200
201/*--------------------------------------------------------------------------
202 * hypre_BoomerAMGBuildCoarseOperator
203 *--------------------------------------------------------------------------*/
204
205int
206hypre_BoomerAMGBuildCoarseOperator( hypre_ParCSRMatrix *RT,
207 hypre_ParCSRMatrix *A,
208 hypre_ParCSRMatrix *P,
209 hypre_ParCSRMatrix **RAP_ptr )
210
211{
212 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
213
214 hypre_CSRMatrix *RT_diag = hypre_ParCSRMatrixDiag(RT);
215 hypre_CSRMatrix *RT_offd = hypre_ParCSRMatrixOffd(RT);
216 int num_cols_diag_RT = hypre_CSRMatrixNumCols(RT_diag);
217 int num_cols_offd_RT = hypre_CSRMatrixNumCols(RT_offd);
218 int num_rows_offd_RT = hypre_CSRMatrixNumRows(RT_offd);
219 hypre_ParCSRCommPkg *comm_pkg_RT = hypre_ParCSRMatrixCommPkg(RT);
220 int num_recvs_RT = 0;
221 int num_sends_RT = 0;
222 int *send_map_starts_RT;
223 int *send_map_elmts_RT;
224
225 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
226
227 double *A_diag_data = hypre_CSRMatrixData(A_diag);
228 int *A_diag_i = hypre_CSRMatrixI(A_diag);
229 int *A_diag_j = hypre_CSRMatrixJ(A_diag);
230
231 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
232
233 double *A_offd_data = hypre_CSRMatrixData(A_offd);
234 int *A_offd_i = hypre_CSRMatrixI(A_offd);
235 int *A_offd_j = hypre_CSRMatrixJ(A_offd);
236
237 int num_cols_diag_A = hypre_CSRMatrixNumCols(A_diag);
238 int num_cols_offd_A = hypre_CSRMatrixNumCols(A_offd);
239
240 hypre_CSRMatrix *P_diag = hypre_ParCSRMatrixDiag(P);
241
242 double *P_diag_data = hypre_CSRMatrixData(P_diag);
243 int *P_diag_i = hypre_CSRMatrixI(P_diag);
244 int *P_diag_j = hypre_CSRMatrixJ(P_diag);
245
246 hypre_CSRMatrix *P_offd = hypre_ParCSRMatrixOffd(P);
247 HYPRE_BigInt *col_map_offd_P = hypre_ParCSRMatrixColMapOffd(P);
248
249 double *P_offd_data = hypre_CSRMatrixData(P_offd);
250 int *P_offd_i = hypre_CSRMatrixI(P_offd);
251 int *P_offd_j = hypre_CSRMatrixJ(P_offd);
252
253 HYPRE_BigInt first_col_diag_P = hypre_ParCSRMatrixFirstColDiag(P);
254 HYPRE_BigInt last_col_diag_P;
255 int num_cols_diag_P = hypre_CSRMatrixNumCols(P_diag);
256 int num_cols_offd_P = hypre_CSRMatrixNumCols(P_offd);
257 HYPRE_BigInt *coarse_partitioning = hypre_ParCSRMatrixColStarts(P);
258 HYPRE_BigInt *RT_partitioning = hypre_ParCSRMatrixColStarts(RT);
259
260 hypre_ParCSRMatrix *RAP;
261 HYPRE_BigInt *col_map_offd_RAP;
262 HYPRE_BigInt *new_col_map_offd_RAP;
263
264 hypre_BigCSRMatrix *RAP_int = NULL;
265
266 double *RAP_int_data;
267 int *RAP_int_i;
268 HYPRE_BigInt *RAP_int_j;
269
270 hypre_BigCSRMatrix *RAP_ext;
271
272 double *RAP_ext_data;
273 int *RAP_ext_i;
274 HYPRE_BigInt *RAP_ext_j;
275 int *int_RAP_ext_j = NULL;
276
277 hypre_CSRMatrix *RAP_diag;
278
279 double *RAP_diag_data;
280 int *RAP_diag_i;
281 int *RAP_diag_j;
282
283 hypre_CSRMatrix *RAP_offd;
284
285 double *RAP_offd_data;
286 int *RAP_offd_i;
287 int *RAP_offd_j;
288
289 int RAP_size;
290 int RAP_ext_size;
291 int RAP_diag_size;
292 int RAP_offd_size;
293 int P_ext_diag_size;
294 int P_ext_offd_size;
295 HYPRE_BigInt first_col_diag_RAP;
296 HYPRE_BigInt last_col_diag_RAP;
297 int num_cols_offd_RAP = 0;
298
299 hypre_CSRMatrix *R_diag;
300
301 double *R_diag_data;
302 int *R_diag_i;
303 int *R_diag_j;
304
305 hypre_CSRMatrix *R_offd;
306
307 double *R_offd_data;
308 int *R_offd_i;
309 int *R_offd_j;
310
311 hypre_BigCSRMatrix *Ps_ext;
312
313 double *Ps_ext_data;
314 int *Ps_ext_i;
315 HYPRE_BigInt *Ps_ext_j;
316
317 double *P_ext_diag_data;
318 int *P_ext_diag_i;
319 int *P_ext_diag_j;
320
321 double *P_ext_offd_data;
322 int *P_ext_offd_i;
323 int *P_ext_offd_j;
324
325 HYPRE_BigInt *col_map_offd_Pext;
326 int *map_P_to_Pext;
327 int *map_P_to_RAP;
328 int *map_Pext_to_RAP;
329
330 int *P_marker;
331 int **P_mark_array;
332 int **A_mark_array;
333 int *A_marker;
334 HYPRE_BigInt *temp;
335
336 HYPRE_BigInt n_coarse, n_coarse_RT;
337 HYPRE_BigInt value;
338 int square = 1;
339 int num_cols_offd_Pext = 0;
340
341 int ic, i, j, k;
342 int i1, i2, i3, ii, ns, ne, size, rest;
343 int cnt, cnt_offd, cnt_diag;
344 int jj1, jj2, jj3, jcol;
345
346 int *jj_count, *jj_cnt_diag, *jj_cnt_offd;
347 int jj_counter, jj_count_diag, jj_count_offd;
348 int jj_row_begining, jj_row_begin_diag, jj_row_begin_offd;
349 int start_indexing = 0; /* start indexing for RAP_data at 0 */
350 int num_nz_cols_A;
351 int num_procs;
352 int num_threads;
353
354 double r_entry;
355 double r_a_product;
356 double r_a_p_product;
357
358 double zero = 0.0;
359
360 int tmp_ii, tmp_ii_1, tmp_ii_2, tmp_ii_10, tmp_ii_11;
361
362 /*-----------------------------------------------------------------------
363 * Copy ParCSRMatrix RT into CSRMatrix R so that we have row-wise access
364 * to restriction .
365 *-----------------------------------------------------------------------*/
366
367 MPI_Comm_size(comm,&num_procs);
368 num_threads = hypre_NumThreads();
369
370 if (comm_pkg_RT)
371 {
372 num_recvs_RT = hypre_ParCSRCommPkgNumRecvs(comm_pkg_RT);
373 num_sends_RT = hypre_ParCSRCommPkgNumSends(comm_pkg_RT);
374 send_map_starts_RT =hypre_ParCSRCommPkgSendMapStarts(comm_pkg_RT);
375 send_map_elmts_RT = hypre_ParCSRCommPkgSendMapElmts(comm_pkg_RT);
376 }
377
378 hypre_CSRMatrixTranspose(RT_diag,&R_diag,1);
379 if (num_cols_offd_RT)
380 {
381 hypre_CSRMatrixTranspose(RT_offd,&R_offd,1);
382 R_offd_data = hypre_CSRMatrixData(R_offd);
383 R_offd_i = hypre_CSRMatrixI(R_offd);
384 R_offd_j = hypre_CSRMatrixJ(R_offd);
385 }
386
387 /*-----------------------------------------------------------------------
388 * Access the CSR vectors for R. Also get sizes of fine and
389 * coarse grids.
390 *-----------------------------------------------------------------------*/
391
392 R_diag_data = hypre_CSRMatrixData(R_diag);
393 R_diag_i = hypre_CSRMatrixI(R_diag);
394 R_diag_j = hypre_CSRMatrixJ(R_diag);
395
396 n_coarse = hypre_ParCSRMatrixGlobalNumCols(P);
397 num_nz_cols_A = num_cols_diag_A + num_cols_offd_A;
398
399 n_coarse_RT = hypre_ParCSRMatrixGlobalNumCols(RT);
400 if (n_coarse != n_coarse_RT)
401 square = 0;
402
403 /*-----------------------------------------------------------------------
404 * Generate Ps_ext, i.e. portion of P that is stored on neighbor procs
405 * and needed locally for triple matrix product
406 *-----------------------------------------------------------------------*/
407
408 if (num_procs > 1)
409 {
410 Ps_ext = hypre_ParCSRMatrixExtractBigExt(P,A,1);
411 Ps_ext_data = hypre_BigCSRMatrixData(Ps_ext);
412 Ps_ext_i = hypre_BigCSRMatrixI(Ps_ext);
413 Ps_ext_j = hypre_BigCSRMatrixJ(Ps_ext);
414 }
415 P_ext_diag_i = hypre_CTAlloc(int,num_cols_offd_A+1);
416 P_ext_offd_i = hypre_CTAlloc(int,num_cols_offd_A+1);
417
418
419 P_ext_diag_size = 0;
420 P_ext_offd_size = 0;
421 last_col_diag_P = first_col_diag_P + (HYPRE_BigInt)num_cols_diag_P - 1;
422
423 for (i=0; i < num_cols_offd_A; i++)
424 {
425 for (j=Ps_ext_i[i]; j < Ps_ext_i[i+1]; j++)
426 if (Ps_ext_j[j] < first_col_diag_P || Ps_ext_j[j] > last_col_diag_P)
427 P_ext_offd_size++;
428 else
429 P_ext_diag_size++;
430 P_ext_diag_i[i+1] = P_ext_diag_size;
431 P_ext_offd_i[i+1] = P_ext_offd_size;
432 }
433
434 if (P_ext_diag_size)
435 {
436 P_ext_diag_j = hypre_CTAlloc(int, P_ext_diag_size);
437 P_ext_diag_data = hypre_CTAlloc(double, P_ext_diag_size);
438 }
439 if (P_ext_offd_size)
440 {
441 P_ext_offd_j = hypre_CTAlloc(int, P_ext_offd_size);
442 P_ext_offd_data = hypre_CTAlloc(double, P_ext_offd_size);
443 }
444
445 if (P_ext_offd_size || num_cols_offd_P)
446 temp = hypre_CTAlloc(HYPRE_BigInt, P_ext_offd_size+num_cols_offd_P);
447
448 cnt_offd = 0;
449 cnt_diag = 0;
450 cnt = 0;
451 for (i=0; i < num_cols_offd_A; i++)
452 {
453 for (j=Ps_ext_i[i]; j < Ps_ext_i[i+1]; j++)
454 {
455 value = Ps_ext_j[j];
456 if (value < first_col_diag_P || value > last_col_diag_P)
457 {
458 Ps_ext_j[cnt_offd] = value;
459 temp[cnt_offd] = value;
460 P_ext_offd_data[cnt_offd++] = Ps_ext_data[j];
461 }
462 else
463 {
464 P_ext_diag_j[cnt_diag] = (int)(value - first_col_diag_P);
465 P_ext_diag_data[cnt_diag++] = Ps_ext_data[j];
466 }
467 }
468 }
469 if (P_ext_offd_size || num_cols_offd_P)
470 {
471 cnt = P_ext_offd_size;
472 for (i=0; i < num_cols_offd_P; i++)
473 temp[cnt++] = col_map_offd_P[i];
474 }
475 if (cnt)
476 {
477 hypre_BigQsort0(temp, 0, cnt-1);
478
479 num_cols_offd_Pext = 1;
480 value = temp[0];
481 for (i=1; i < cnt; i++)
482 {
483 if (temp[i] > value)
484 {
485 value = temp[i];
486 temp[num_cols_offd_Pext++] = value;
487 }
488 }
489 }
490
491 if (num_cols_offd_Pext)
492 col_map_offd_Pext = hypre_CTAlloc(HYPRE_BigInt,num_cols_offd_Pext);
493
494 for (i=0; i < num_cols_offd_Pext; i++)
495 col_map_offd_Pext[i] = temp[i];
496
497 for (i=0 ; i < P_ext_offd_size; i++)
498 P_ext_offd_j[i] = hypre_BigBinarySearch(col_map_offd_Pext,
499 Ps_ext_j[i],
500 num_cols_offd_Pext);
501
502 if (P_ext_offd_size || num_cols_offd_P)
503 hypre_TFree(temp);
504 if (num_procs > 1)
505 {
506 hypre_BigCSRMatrixDestroy(Ps_ext);
507 Ps_ext = NULL;
508 }
509
510
511 if (num_cols_offd_P)
512 {
513 map_P_to_Pext = hypre_CTAlloc(int,num_cols_offd_P);
514
515 cnt = 0;
516 for (i=0; i < num_cols_offd_Pext; i++)
517 if (col_map_offd_Pext[i] == col_map_offd_P[cnt])
518 {
519 map_P_to_Pext[cnt++] = i;
520 if (cnt == num_cols_offd_P) break;
521 }
522 }
523
524 /*-----------------------------------------------------------------------
525 * First Pass: Determine size of RAP_int and set up RAP_int_i if there
526 * are more than one processor and nonzero elements in R_offd
527 *-----------------------------------------------------------------------*/
528
529 P_mark_array = hypre_CTAlloc(int *, num_threads);
530 A_mark_array = hypre_CTAlloc(int *, num_threads);
531
532 if (num_cols_offd_RT)
533 {
534 jj_count = hypre_CTAlloc(int, num_threads);
535
536#define HYPRE_SMP_PRIVATE i,ii,ic,i1,i2,i3,jj1,jj2,jj3,ns,ne,size,rest,jj_counter,jj_row_begining,A_marker,P_marker
537#include "../utilities/hypre_smp_forloop.h"
538 for (ii = 0; ii < num_threads; ii++)
539 {
540 size = num_cols_offd_RT/num_threads;
541 rest = num_cols_offd_RT - size*num_threads;
542 if (ii < rest)
543 {
544 ns = ii*size+ii;
545 ne = (ii+1)*size+ii+1;
546 }
547 else
548 {
549 ns = ii*size+rest;
550 ne = (ii+1)*size+rest;
551 }
552
553 /*-----------------------------------------------------------------------
554 * Allocate marker arrays.
555 *-----------------------------------------------------------------------*/
556
557 if (num_cols_offd_Pext || num_cols_diag_P)
558 {
559 P_mark_array[ii] = hypre_CTAlloc(int, num_cols_diag_P+num_cols_offd_Pext);
560 P_marker = P_mark_array[ii];
561 }
562 A_mark_array[ii] = hypre_CTAlloc(int, num_nz_cols_A);
563 A_marker = A_mark_array[ii];
564 /*-----------------------------------------------------------------------
565 * Initialize some stuff.
566 *-----------------------------------------------------------------------*/
567
568 jj_counter = start_indexing;
569 for (ic = 0; ic < num_cols_diag_P+num_cols_offd_Pext; ic++)
570 {
571 P_marker[ic] = -1;
572 }
573 for (i = 0; i < num_nz_cols_A; i++)
574 {
575 A_marker[i] = -1;
576 }
577
578 /*-----------------------------------------------------------------------
579 * Loop over exterior c-points
580 *-----------------------------------------------------------------------*/
581
582 for (ic = ns; ic < ne; ic++)
583 {
584
585 jj_row_begining = jj_counter;
586
587 /*--------------------------------------------------------------------
588 * Loop over entries in row ic of R_offd.
589 *--------------------------------------------------------------------*/
590
591 for (jj1 = R_offd_i[ic]; jj1 < R_offd_i[ic+1]; jj1++)
592 {
593 i1 = R_offd_j[jj1];
594
595 /*-----------------------------------------------------------------
596 * Loop over entries in row i1 of A_offd.
597 *-----------------------------------------------------------------*/
598
599 for (jj2 = A_offd_i[i1]; jj2 < A_offd_i[i1+1]; jj2++)
600 {
601 i2 = A_offd_j[jj2];
602
603 /*--------------------------------------------------------------
604 * Check A_marker to see if point i2 has been previously
605 * visited. New entries in RAP only occur from unmarked points.
606 *--------------------------------------------------------------*/
607
608 if (A_marker[i2] != ic)
609 {
610
611 /*-----------------------------------------------------------
612 * Mark i2 as visited.
613 *-----------------------------------------------------------*/
614
615 A_marker[i2] = ic;
616
617 /*-----------------------------------------------------------
618 * Loop over entries in row i2 of P_ext.
619 *-----------------------------------------------------------*/
620
621 for (jj3 = P_ext_diag_i[i2]; jj3 < P_ext_diag_i[i2+1]; jj3++)
622 {
623 i3 = P_ext_diag_j[jj3];
624
625 /*--------------------------------------------------------
626 * Check P_marker to see that RAP_{ic,i3} has not already
627 * been accounted for. If it has not, mark it and increment
628 * counter.
629 *--------------------------------------------------------*/
630
631 if (P_marker[i3] < jj_row_begining)
632 {
633 P_marker[i3] = jj_counter;
634 jj_counter++;
635 }
636 }
637 for (jj3 = P_ext_offd_i[i2]; jj3 < P_ext_offd_i[i2+1]; jj3++)
638 {
639 i3 = P_ext_offd_j[jj3] + num_cols_diag_P;
640
641 /*--------------------------------------------------------
642 * Check P_marker to see that RAP_{ic,i3} has not already
643 * been accounted for. If it has not, mark it and increment
644 * counter.
645 *--------------------------------------------------------*/
646
647 if (P_marker[i3] < jj_row_begining)
648 {
649 P_marker[i3] = jj_counter;
650 jj_counter++;
651 }
652 }
653 }
654 }
655 /*-----------------------------------------------------------------
656 * Loop over entries in row i1 of A_diag.
657 *-----------------------------------------------------------------*/
658
659 for (jj2 = A_diag_i[i1]; jj2 < A_diag_i[i1+1]; jj2++)
660 {
661 i2 = A_diag_j[jj2];
662
663 /*--------------------------------------------------------------
664 * Check A_marker to see if point i2 has been previously
665 * visited. New entries in RAP only occur from unmarked points.
666 *--------------------------------------------------------------*/
667
668 if (A_marker[i2+num_cols_offd_A] != ic)
669 {
670
671 /*-----------------------------------------------------------
672 * Mark i2 as visited.
673 *-----------------------------------------------------------*/
674
675 A_marker[i2+num_cols_offd_A] = ic;
676
677 /*-----------------------------------------------------------
678 * Loop over entries in row i2 of P_diag.
679 *-----------------------------------------------------------*/
680
681 for (jj3 = P_diag_i[i2]; jj3 < P_diag_i[i2+1]; jj3++)
682 {
683 i3 = P_diag_j[jj3];
684
685 /*--------------------------------------------------------
686 * Check P_marker to see that RAP_{ic,i3} has not already
687 * been accounted for. If it has not, mark it and increment
688 * counter.
689 *--------------------------------------------------------*/
690
691 if (P_marker[i3] < jj_row_begining)
692 {
693 P_marker[i3] = jj_counter;
694 jj_counter++;
695 }
696 }
697 /*-----------------------------------------------------------
698 * Loop over entries in row i2 of P_offd.
699 *-----------------------------------------------------------*/
700
701 for (jj3 = P_offd_i[i2]; jj3 < P_offd_i[i2+1]; jj3++)
702 {
703 i3 = map_P_to_Pext[P_offd_j[jj3]] + num_cols_diag_P;
704
705 /*--------------------------------------------------------
706 * Check P_marker to see that RAP_{ic,i3} has not already
707 * been accounted for. If it has not, mark it and increment
708 * counter.
709 *--------------------------------------------------------*/
710
711 if (P_marker[i3] < jj_row_begining)
712 {
713 P_marker[i3] = jj_counter;
714 jj_counter++;
715 }
716 }
717 }
718 }
719 }
720 }
721
722 jj_count[ii] = jj_counter;
723
724 }
725
726 /*-----------------------------------------------------------------------
727 * Allocate RAP_int_data and RAP_int_j arrays.
728 *-----------------------------------------------------------------------*/
729 for (i = 0; i < num_threads-1; i++)
730 jj_count[i+1] += jj_count[i];
731
732 RAP_size = jj_count[num_threads-1];
733 RAP_int_i = hypre_CTAlloc(int, num_cols_offd_RT+1);
734 RAP_int_data = hypre_CTAlloc(double, RAP_size);
735 RAP_int_j = hypre_CTAlloc(HYPRE_BigInt, RAP_size);
736
737 RAP_int_i[num_cols_offd_RT] = RAP_size;
738
739 /*-----------------------------------------------------------------------
740 * Second Pass: Fill in RAP_int_data and RAP_int_j.
741 *-----------------------------------------------------------------------*/
742
743#define HYPRE_SMP_PRIVATE i,ii,ic,i1,i2,i3,jj1,jj2,jj3,ns,ne,size,rest,jj_counter,jj_row_begining,A_marker,P_marker,r_entry,r_a_product,r_a_p_product
744#include "../utilities/hypre_smp_forloop.h"
745 for (ii = 0; ii < num_threads; ii++)
746 {
747 size = num_cols_offd_RT/num_threads;
748 rest = num_cols_offd_RT - size*num_threads;
749 if (ii < rest)
750 {
751 ns = ii*size+ii;
752 ne = (ii+1)*size+ii+1;
753 }
754 else
755 {
756 ns = ii*size+rest;
757 ne = (ii+1)*size+rest;
758 }
759
760 /*-----------------------------------------------------------------------
761 * Initialize some stuff.
762 *-----------------------------------------------------------------------*/
763 if (num_cols_offd_Pext || num_cols_diag_P)
764 P_marker = P_mark_array[ii];
765 A_marker = A_mark_array[ii];
766
767 jj_counter = start_indexing;
768 if (ii > 0) jj_counter = jj_count[ii-1];
769
770 for (ic = 0; ic < num_cols_diag_P+num_cols_offd_Pext; ic++)
771 {
772 P_marker[ic] = -1;
773 }
774 for (i = 0; i < num_nz_cols_A; i++)
775 {
776 A_marker[i] = -1;
777 }
778
779 /*-----------------------------------------------------------------------
780 * Loop over exterior c-points.
781 *-----------------------------------------------------------------------*/
782
783 for (ic = ns; ic < ne; ic++)
784 {
785
786 jj_row_begining = jj_counter;
787 RAP_int_i[ic] = jj_counter;
788
789 /*--------------------------------------------------------------------
790 * Loop over entries in row ic of R_offd.
791 *--------------------------------------------------------------------*/
792
793 for (jj1 = R_offd_i[ic]; jj1 < R_offd_i[ic+1]; jj1++)
794 {
795 i1 = R_offd_j[jj1];
796 r_entry = R_offd_data[jj1];
797
798 /*-----------------------------------------------------------------
799 * Loop over entries in row i1 of A_offd.
800 *-----------------------------------------------------------------*/
801
802 for (jj2 = A_offd_i[i1]; jj2 < A_offd_i[i1+1]; jj2++)
803 {
804 i2 = A_offd_j[jj2];
805 r_a_product = r_entry * A_offd_data[jj2];
806
807 /*--------------------------------------------------------------
808 * Check A_marker to see if point i2 has been previously
809 * visited. New entries in RAP only occur from unmarked points.
810 *--------------------------------------------------------------*/
811
812 if (A_marker[i2] != ic)
813 {
814
815 /*-----------------------------------------------------------
816 * Mark i2 as visited.
817 *-----------------------------------------------------------*/
818
819 A_marker[i2] = ic;
820
821 /*-----------------------------------------------------------
822 * Loop over entries in row i2 of P_ext.
823 *-----------------------------------------------------------*/
824
825 for (jj3 = P_ext_diag_i[i2]; jj3 < P_ext_diag_i[i2+1]; jj3++)
826 {
827 i3 = P_ext_diag_j[jj3];
828 r_a_p_product = r_a_product * P_ext_diag_data[jj3];
829
830 /*--------------------------------------------------------
831 * Check P_marker to see that RAP_{ic,i3} has not already
832 * been accounted for. If it has not, create a new entry.
833 * If it has, add new contribution.
834 *--------------------------------------------------------*/
835
836 if (P_marker[i3] < jj_row_begining)
837 {
838 P_marker[i3] = jj_counter;
839 RAP_int_data[jj_counter] = r_a_p_product;
840 RAP_int_j[jj_counter] = (HYPRE_BigInt)i3 + first_col_diag_P;
841 jj_counter++;
842 }
843 else
844 {
845 RAP_int_data[P_marker[i3]] += r_a_p_product;
846 }
847 }
848
849 for (jj3 = P_ext_offd_i[i2]; jj3 < P_ext_offd_i[i2+1]; jj3++)
850 {
851 i3 = P_ext_offd_j[jj3] + num_cols_diag_P;
852 r_a_p_product = r_a_product * P_ext_offd_data[jj3];
853
854 /*--------------------------------------------------------
855 * Check P_marker to see that RAP_{ic,i3} has not already
856 * been accounted for. If it has not, create a new entry.
857 * If it has, add new contribution.
858 *--------------------------------------------------------*/
859
860 if (P_marker[i3] < jj_row_begining)
861 {
862 P_marker[i3] = jj_counter;
863 RAP_int_data[jj_counter] = r_a_p_product;
864 RAP_int_j[jj_counter]
865 = col_map_offd_Pext[i3-num_cols_diag_P];
866 jj_counter++;
867 }
868 else
869 {
870 RAP_int_data[P_marker[i3]] += r_a_p_product;
871 }
872 }
873 }
874
875 /*--------------------------------------------------------------
876 * If i2 is previously visited ( A_marker[12]=ic ) it yields
877 * no new entries in RAP and can just add new contributions.
878 *--------------------------------------------------------------*/
879
880 else
881 {
882 for (jj3 = P_ext_diag_i[i2]; jj3 < P_ext_diag_i[i2+1]; jj3++)
883 {
884 i3 = P_ext_diag_j[jj3];
885 r_a_p_product = r_a_product * P_ext_diag_data[jj3];
886 RAP_int_data[P_marker[i3]] += r_a_p_product;
887 }
888 for (jj3 = P_ext_offd_i[i2]; jj3 < P_ext_offd_i[i2+1]; jj3++)
889 {
890 i3 = P_ext_offd_j[jj3] + num_cols_diag_P;
891 r_a_p_product = r_a_product * P_ext_offd_data[jj3];
892 RAP_int_data[P_marker[i3]] += r_a_p_product;
893 }
894 }
895 }
896
897 /*-----------------------------------------------------------------
898 * Loop over entries in row i1 of A_diag.
899 *-----------------------------------------------------------------*/
900
901 for (jj2 = A_diag_i[i1]; jj2 < A_diag_i[i1+1]; jj2++)
902 {
903 i2 = A_diag_j[jj2];
904 r_a_product = r_entry * A_diag_data[jj2];
905
906 /*--------------------------------------------------------------
907 * Check A_marker to see if point i2 has been previously
908 * visited. New entries in RAP only occur from unmarked points.
909 *--------------------------------------------------------------*/
910
911 if (A_marker[i2+num_cols_offd_A] != ic)
912 {
913
914 /*-----------------------------------------------------------
915 * Mark i2 as visited.
916 *-----------------------------------------------------------*/
917
918 A_marker[i2+num_cols_offd_A] = ic;
919
920 /*-----------------------------------------------------------
921 * Loop over entries in row i2 of P_diag.
922 *-----------------------------------------------------------*/
923
924 for (jj3 = P_diag_i[i2]; jj3 < P_diag_i[i2+1]; jj3++)
925 {
926 i3 = P_diag_j[jj3];
927 r_a_p_product = r_a_product * P_diag_data[jj3];
928
929 /*--------------------------------------------------------
930 * Check P_marker to see that RAP_{ic,i3} has not already
931 * been accounted for. If it has not, create a new entry.
932 * If it has, add new contribution.
933 *--------------------------------------------------------*/
934
935 if (P_marker[i3] < jj_row_begining)
936 {
937 P_marker[i3] = jj_counter;
938 RAP_int_data[jj_counter] = r_a_p_product;
939 RAP_int_j[jj_counter] = (HYPRE_BigInt)i3 + first_col_diag_P;
940 jj_counter++;
941 }
942 else
943 {
944 RAP_int_data[P_marker[i3]] += r_a_p_product;
945 }
946 }
947 for (jj3 = P_offd_i[i2]; jj3 < P_offd_i[i2+1]; jj3++)
948 {
949 i3 = map_P_to_Pext[P_offd_j[jj3]] + num_cols_diag_P;
950 r_a_p_product = r_a_product * P_offd_data[jj3];
951
952 /*--------------------------------------------------------
953 * Check P_marker to see that RAP_{ic,i3} has not already
954 * been accounted for. If it has not, create a new entry.
955 * If it has, add new contribution.
956 *--------------------------------------------------------*/
957
958 if (P_marker[i3] < jj_row_begining)
959 {
960 P_marker[i3] = jj_counter;
961 RAP_int_data[jj_counter] = r_a_p_product;
962 RAP_int_j[jj_counter] =
963 col_map_offd_Pext[i3-num_cols_diag_P];
964 jj_counter++;
965 }
966 else
967 {
968 RAP_int_data[P_marker[i3]] += r_a_p_product;
969 }
970 }
971 }
972
973 /*--------------------------------------------------------------
974 * If i2 is previously visited ( A_marker[12]=ic ) it yields
975 * no new entries in RAP and can just add new contributions.
976 *--------------------------------------------------------------*/
977
978 else
979 {
980 for (jj3 = P_diag_i[i2]; jj3 < P_diag_i[i2+1]; jj3++)
981 {
982 i3 = P_diag_j[jj3];
983 r_a_p_product = r_a_product * P_diag_data[jj3];
984 RAP_int_data[P_marker[i3]] += r_a_p_product;
985 }
986 for (jj3 = P_offd_i[i2]; jj3 < P_offd_i[i2+1]; jj3++)
987 {
988 i3 = map_P_to_Pext[P_offd_j[jj3]] + num_cols_diag_P;
989 r_a_p_product = r_a_product * P_offd_data[jj3];
990 RAP_int_data[P_marker[i3]] += r_a_p_product;
991 }
992 }
993 }
994 }
995 }
996 if (num_cols_offd_Pext || num_cols_diag_P)
997 hypre_TFree(P_mark_array[ii]);
998 hypre_TFree(A_mark_array[ii]);
999 }
1000
1001 RAP_int = hypre_BigCSRMatrixCreate(num_cols_offd_RT,num_rows_offd_RT,RAP_size);
1002 hypre_BigCSRMatrixI(RAP_int) = RAP_int_i;
1003 hypre_BigCSRMatrixJ(RAP_int) = RAP_int_j;
1004 hypre_BigCSRMatrixData(RAP_int) = RAP_int_data;
1005 hypre_TFree(jj_count);
1006 }
1007
1008 RAP_ext_size = 0;
1009 if (num_sends_RT || num_recvs_RT)
1010 {
1011 RAP_ext = hypre_ExchangeRAPData(RAP_int,comm_pkg_RT);
1012 RAP_ext_i = hypre_BigCSRMatrixI(RAP_ext);
1013 RAP_ext_j = hypre_BigCSRMatrixJ(RAP_ext);
1014 RAP_ext_data = hypre_BigCSRMatrixData(RAP_ext);
1015 RAP_ext_size = RAP_ext_i[hypre_BigCSRMatrixNumRows(RAP_ext)];
1016 }
1017 if (num_cols_offd_RT)
1018 {
1019 hypre_BigCSRMatrixDestroy(RAP_int);
1020 RAP_int = NULL;
1021 }
1022 RAP_diag_i = hypre_CTAlloc(int, num_cols_diag_RT+1);
1023 RAP_offd_i = hypre_CTAlloc(int, num_cols_diag_RT+1);
1024
1025 first_col_diag_RAP = first_col_diag_P;
1026 last_col_diag_RAP = first_col_diag_P + (HYPRE_BigInt)num_cols_diag_P - 1;
1027
1028 /*-----------------------------------------------------------------------
1029 * check for new nonzero columns in RAP_offd generated through RAP_ext
1030 *-----------------------------------------------------------------------*/
1031
1032 if (RAP_ext_size || num_cols_offd_Pext)
1033 {
1034 temp = hypre_CTAlloc(HYPRE_BigInt,RAP_ext_size+num_cols_offd_Pext);
1035 cnt = 0;
1036 for (i=0; i < RAP_ext_size; i++)
1037 if (RAP_ext_j[i] < first_col_diag_RAP
1038 || RAP_ext_j[i] > last_col_diag_RAP)
1039 temp[cnt++] = RAP_ext_j[i];
1040 for (i=0; i < num_cols_offd_Pext; i++)
1041 temp[cnt++] = col_map_offd_Pext[i];
1042
1043
1044 if (cnt)
1045 {
1046 hypre_BigQsort0(temp,0,cnt-1);
1047 value = temp[0];
1048 num_cols_offd_RAP = 1;
1049 for (i=1; i < cnt; i++)
1050 {
1051 if (temp[i] > value)
1052 {
1053 value = temp[i];
1054 temp[num_cols_offd_RAP++] = value;
1055 }
1056 }
1057 }
1058
1059 /* now evaluate col_map_offd_RAP */
1060 if (num_cols_offd_RAP)
1061 col_map_offd_RAP = hypre_CTAlloc(HYPRE_BigInt, num_cols_offd_RAP);
1062
1063 for (i=0 ; i < num_cols_offd_RAP; i++)
1064 col_map_offd_RAP[i] = temp[i];
1065
1066 hypre_TFree(temp);
1067 }
1068
1069 if (num_cols_offd_P)
1070 {
1071 map_P_to_RAP = hypre_CTAlloc(int,num_cols_offd_P);
1072
1073 cnt = 0;
1074 for (i=0; i < num_cols_offd_RAP; i++)
1075 if (col_map_offd_RAP[i] == col_map_offd_P[cnt])
1076 {
1077 map_P_to_RAP[cnt++] = i;
1078 if (cnt == num_cols_offd_P) break;
1079 }
1080 }
1081
1082 if (num_cols_offd_Pext)
1083 {
1084 map_Pext_to_RAP = hypre_CTAlloc(int,num_cols_offd_Pext);
1085
1086 cnt = 0;
1087 for (i=0; i < num_cols_offd_RAP; i++)
1088 if (col_map_offd_RAP[i] == col_map_offd_Pext[cnt])
1089 {
1090 map_Pext_to_RAP[cnt++] = i;
1091 if (cnt == num_cols_offd_Pext) break;
1092 }
1093 }
1094
1095 /*-----------------------------------------------------------------------
1096 * Convert RAP_ext column indices
1097 *-----------------------------------------------------------------------*/
1098
1099 int_RAP_ext_j = hypre_CTAlloc(int, RAP_ext_size);
1100
1101 for (i=0; i < RAP_ext_size; i++)
1102 if (RAP_ext_j[i] < first_col_diag_RAP
1103 || RAP_ext_j[i] > last_col_diag_RAP)
1104 int_RAP_ext_j[i] = num_cols_diag_P
1105 + hypre_BigBinarySearch(col_map_offd_RAP,
1106 RAP_ext_j[i],num_cols_offd_RAP);
1107 else
1108 int_RAP_ext_j[i] = (int)(RAP_ext_j[i]-first_col_diag_RAP);
1109
1110/* need to allocate new P_marker etc. and make further changes */
1111 /*-----------------------------------------------------------------------
1112 * Initialize some stuff.
1113 *-----------------------------------------------------------------------*/
1114 jj_cnt_diag = hypre_CTAlloc(int, num_threads);
1115 jj_cnt_offd = hypre_CTAlloc(int, num_threads);
1116 size = num_cols_diag_RT/num_threads;
1117 rest = num_cols_diag_RT - size*num_threads;
1118
1119#define HYPRE_SMP_PRIVATE i,j,k,jcol,ii,ic,i1,i2,i3,jj1,jj2,jj3,ns,ne,jj_count_diag,jj_count_offd,jj_row_begin_diag,jj_row_begin_offd,A_marker,P_marker, tmp_ii_10, tmp_ii, tmp_ii_1
1120#include "../utilities/hypre_smp_forloop.h"
1121 for (ii = 0; ii < num_threads; ii++)
1122 {
1123 /*size = num_cols_diag_RT/num_threads;
1124 rest = num_cols_diag_RT - size*num_threads;*/
1125 if (ii < rest)
1126 {
1127 ns = ii*size+ii;
1128 ne = (ii+1)*size+ii+1;
1129 }
1130 else
1131 {
1132 ns = ii*size+rest;
1133 ne = (ii+1)*size+rest;
1134 }
1135
1136 P_mark_array[ii] = hypre_CTAlloc(int, num_cols_diag_P+num_cols_offd_RAP);
1137 A_mark_array[ii] = hypre_CTAlloc(int, num_nz_cols_A);
1138 P_marker = P_mark_array[ii];
1139 A_marker = A_mark_array[ii];
1140 jj_count_diag = start_indexing;
1141 jj_count_offd = start_indexing;
1142
1143 for (ic = 0; ic < num_cols_diag_P+num_cols_offd_RAP; ic++)
1144 {
1145 P_marker[ic] = -1;
1146 }
1147 for (i = 0; i < num_nz_cols_A; i++)
1148 {
1149 A_marker[i] = -1;
1150 }
1151
1152 /*-----------------------------------------------------------------------
1153 * Loop over interior c-points.
1154 *-----------------------------------------------------------------------*/
1155
1156 for (ic = ns; ic < ne; ic++)
1157 {
1158
1159 /*--------------------------------------------------------------------
1160 * Set marker for diagonal entry, RAP_{ic,ic}. and for all points
1161 * being added to row ic of RAP_diag and RAP_offd through RAP_ext
1162 *--------------------------------------------------------------------*/
1163
1164 jj_row_begin_diag = jj_count_diag;
1165 jj_row_begin_offd = jj_count_offd;
1166
1167 if (square)
1168 P_marker[ic] = jj_count_diag++;
1169
1170 for (i=0; i < num_sends_RT; i++)
1171 {
1172 tmp_ii_10 = send_map_starts_RT[i+1];
1173 for (j = send_map_starts_RT[i]; j < tmp_ii_10; j++)
1174 if (send_map_elmts_RT[j] == ic)
1175 {
1176 for (k=RAP_ext_i[j]; k < RAP_ext_i[j+1]; k++)
1177 {
1178 jcol = int_RAP_ext_j[k];
1179 if (jcol < num_cols_diag_P)
1180 {
1181 if (P_marker[jcol] < jj_row_begin_diag)
1182 {
1183 P_marker[jcol] = jj_count_diag;
1184 jj_count_diag++;
1185 }
1186 }
1187 else
1188 {
1189 if (P_marker[jcol] < jj_row_begin_offd)
1190 {
1191 P_marker[jcol] = jj_count_offd;
1192 jj_count_offd++;
1193 }
1194 }
1195 }
1196 break;
1197 }
1198 }
1199 /*--------------------------------------------------------------------
1200 * Loop over entries in row ic of R_diag.
1201 *--------------------------------------------------------------------*/
1202
1203 for (jj1 = R_diag_i[ic]; jj1 < R_diag_i[ic+1]; jj1++)
1204 {
1205 i1 = R_diag_j[jj1];
1206
1207 /*-----------------------------------------------------------------
1208 * Loop over entries in row i1 of A_offd.
1209 *-----------------------------------------------------------------*/
1210
1211 if (num_cols_offd_A)
1212 {
1213 for (jj2 = A_offd_i[i1]; jj2 < A_offd_i[i1+1]; jj2++)
1214 {
1215 i2 = A_offd_j[jj2];
1216
1217 /*--------------------------------------------------------------
1218 * Check A_marker to see if point i2 has been previously
1219 * visited. New entries in RAP only occur from unmarked points.
1220 *--------------------------------------------------------------*/
1221
1222 if (A_marker[i2] != ic)
1223 {
1224
1225 /*-----------------------------------------------------------
1226 * Mark i2 as visited.
1227 *-----------------------------------------------------------*/
1228
1229 A_marker[i2] = ic;
1230
1231 /*-----------------------------------------------------------
1232 * Loop over entries in row i2 of P_ext.
1233 *-----------------------------------------------------------*/
1234
1235 for (jj3 = P_ext_diag_i[i2]; jj3 < P_ext_diag_i[i2+1]; jj3++)
1236 {
1237 i3 = P_ext_diag_j[jj3];
1238
1239 /*--------------------------------------------------------
1240 * Check P_marker to see that RAP_{ic,i3} has not already
1241 * been accounted for. If it has not, mark it and increment
1242 * counter.
1243 *--------------------------------------------------------*/
1244
1245 if (P_marker[i3] < jj_row_begin_diag)
1246 {
1247 P_marker[i3] = jj_count_diag;
1248 jj_count_diag++;
1249 }
1250 }
1251 for (jj3 = P_ext_offd_i[i2]; jj3 < P_ext_offd_i[i2+1]; jj3++)
1252 {
1253 i3 = map_Pext_to_RAP[P_ext_offd_j[jj3]]+num_cols_diag_P;
1254
1255 /*--------------------------------------------------------
1256 * Check P_marker to see that RAP_{ic,i3} has not already
1257 * been accounted for. If it has not, mark it and increment
1258 * counter.
1259 *--------------------------------------------------------*/
1260
1261 if (P_marker[i3] < jj_row_begin_offd)
1262 {
1263 P_marker[i3] = jj_count_offd;
1264 jj_count_offd++;
1265 }
1266 }
1267 }
1268 }
1269 }
1270 /*-----------------------------------------------------------------
1271 * Loop over entries in row i1 of A_diag.
1272 *-----------------------------------------------------------------*/
1273
1274 tmp_ii = A_diag_i[i1+1];
1275 for (jj2 = A_diag_i[i1]; jj2 < tmp_ii; jj2++)
1276 {
1277 i2 = A_diag_j[jj2];
1278
1279 /*--------------------------------------------------------------
1280 * Check A_marker to see if point i2 has been previously
1281 * visited. New entries in RAP only occur from unmarked points.
1282 *--------------------------------------------------------------*/
1283
1284 if (A_marker[i2+num_cols_offd_A] != ic)
1285 {
1286 /*-----------------------------------------------------------
1287 * Mark i2 as visited.
1288 *-----------------------------------------------------------*/
1289
1290 A_marker[i2+num_cols_offd_A] = ic;
1291
1292 /*-----------------------------------------------------------
1293 * Loop over entries in row i2 of P_diag.
1294 *-----------------------------------------------------------*/
1295
1296 tmp_ii_1 = P_diag_i[i2+1];
1297 for (jj3 = P_diag_i[i2]; jj3 < tmp_ii_1; jj3++)
1298 {
1299 i3 = P_diag_j[jj3];
1300
1301 /*--------------------------------------------------------
1302 * Check P_marker to see that RAP_{ic,i3} has not already
1303 * been accounted for. If it has not, mark it and increment
1304 * counter.
1305 *--------------------------------------------------------*/
1306
1307 if (P_marker[i3] < jj_row_begin_diag)
1308 {
1309 P_marker[i3] = jj_count_diag;
1310 jj_count_diag++;
1311 }
1312 }
1313 /*-----------------------------------------------------------
1314 * Loop over entries in row i2 of P_offd.
1315 *-----------------------------------------------------------*/
1316
1317 if (num_cols_offd_P)
1318 {
1319 for (jj3 = P_offd_i[i2]; jj3 < P_offd_i[i2+1]; jj3++)
1320 {
1321 i3 = map_P_to_RAP[P_offd_j[jj3]] + num_cols_diag_P;
1322
1323 /*--------------------------------------------------------
1324 * Check P_marker to see that RAP_{ic,i3} has not already
1325 * been accounted for. If it has not, mark it and increment
1326 * counter.
1327 *--------------------------------------------------------*/
1328
1329 if (P_marker[i3] < jj_row_begin_offd)
1330 {
1331 P_marker[i3] = jj_count_offd;
1332 jj_count_offd++;
1333 }
1334 }
1335 }
1336 }
1337 }
1338 }
1339
1340 /*--------------------------------------------------------------------
1341 * Set RAP_diag_i and RAP_offd_i for this row.
1342 *--------------------------------------------------------------------*/
1343/*
1344 RAP_diag_i[ic] = jj_row_begin_diag;
1345 RAP_offd_i[ic] = jj_row_begin_offd;
1346*/
1347 }
1348 jj_cnt_diag[ii] = jj_count_diag;
1349 jj_cnt_offd[ii] = jj_count_offd;
1350 }
1351
1352 for (i=0; i < num_threads-1; i++)
1353 {
1354 jj_cnt_diag[i+1] += jj_cnt_diag[i];
1355 jj_cnt_offd[i+1] += jj_cnt_offd[i];
1356 }
1357
1358 jj_count_diag = jj_cnt_diag[num_threads-1];
1359 jj_count_offd = jj_cnt_offd[num_threads-1];
1360
1361 RAP_diag_i[num_cols_diag_RT] = jj_count_diag;
1362 RAP_offd_i[num_cols_diag_RT] = jj_count_offd;
1363
1364 /*-----------------------------------------------------------------------
1365 * Allocate RAP_diag_data and RAP_diag_j arrays.
1366 * Allocate RAP_offd_data and RAP_offd_j arrays.
1367 *-----------------------------------------------------------------------*/
1368
1369 RAP_diag_size = jj_count_diag;
1370 /*printf ("RAP_diag_size %d\n", RAP_diag_size);
1371 fflush(NULL);*/
1372 if (RAP_diag_size)
1373 {
1374 RAP_diag_data = hypre_CTAlloc(double, RAP_diag_size);
1375 RAP_diag_j = hypre_CTAlloc(int, RAP_diag_size);
1376 }
1377
1378 RAP_offd_size = jj_count_offd;
1379 if (RAP_offd_size)
1380 {
1381 RAP_offd_data = hypre_CTAlloc(double, RAP_offd_size);
1382 RAP_offd_j = hypre_CTAlloc(int, RAP_offd_size);
1383 }
1384
1385 if (RAP_offd_size == 0 && num_cols_offd_RAP != 0)
1386 {
1387 num_cols_offd_RAP = 0;
1388 hypre_TFree(col_map_offd_RAP);
1389 }
1390 /*-----------------------------------------------------------------------
1391 * Second Pass: Fill in RAP_diag_data and RAP_diag_j.
1392 * Second Pass: Fill in RAP_offd_data and RAP_offd_j.
1393 *-----------------------------------------------------------------------*/
1394
1395#define HYPRE_SMP_PRIVATE i,j,k,jcol,ii,ic,i1,i2,i3,jj1,jj2,jj3,ns,ne,size,rest,jj_count_diag,jj_count_offd,jj_row_begin_diag,jj_row_begin_offd,A_marker,P_marker,r_entry,r_a_product,r_a_p_product, tmp_ii_11, tmp_ii, tmp_ii_1, tmp_ii_2
1396#include "../utilities/hypre_smp_forloop.h"
1397 for (ii = 0; ii < num_threads; ii++)
1398 {
1399 size = num_cols_diag_RT/num_threads;
1400 rest = num_cols_diag_RT - size*num_threads;
1401 if (ii < rest)
1402 {
1403 ns = ii*size+ii;
1404 ne = (ii+1)*size+ii+1;
1405 }
1406 else
1407 {
1408 ns = ii*size+rest;
1409 ne = (ii+1)*size+rest;
1410 }
1411
1412 /*-----------------------------------------------------------------------
1413 * Initialize some stuff.
1414 *-----------------------------------------------------------------------*/
1415
1416 P_marker = P_mark_array[ii];
1417 A_marker = A_mark_array[ii];
1418 for (ic = 0; ic < num_cols_diag_P+num_cols_offd_RAP; ic++)
1419 {
1420 P_marker[ic] = -1;
1421 }
1422 for (i = 0; i < num_nz_cols_A ; i++)
1423 {
1424 A_marker[i] = -1;
1425 }
1426
1427 jj_count_diag = start_indexing;
1428 jj_count_offd = start_indexing;
1429 if (ii > 0)
1430 {
1431 jj_count_diag = jj_cnt_diag[ii-1];
1432 jj_count_offd = jj_cnt_offd[ii-1];
1433 }
1434
1435 /*-----------------------------------------------------------------------
1436 * Loop over interior c-points.
1437 *-----------------------------------------------------------------------*/
1438
1439 for (ic = ns; ic < ne; ic++)
1440 {
1441
1442 /*--------------------------------------------------------------------
1443 * Create diagonal entry, RAP_{ic,ic} and add entries of RAP_ext
1444 *--------------------------------------------------------------------*/
1445
1446 jj_row_begin_diag = jj_count_diag;
1447 jj_row_begin_offd = jj_count_offd;
1448 RAP_diag_i[ic] = jj_row_begin_diag;
1449 RAP_offd_i[ic] = jj_row_begin_offd;
1450
1451 if (square)
1452 {
1453 P_marker[ic] = jj_count_diag;
1454 /*if (jj_count_diag > RAP_diag_size-1)
1455 { printf(" warning jj_count_diag %d\n", jj_count_diag);
1456 fflush(NULL);}*/
1457 RAP_diag_data[jj_count_diag] = zero;
1458 RAP_diag_j[jj_count_diag] = ic;
1459 jj_count_diag++;
1460 }
1461
1462 for (i=0; i < num_sends_RT; i++)
1463 {
1464 tmp_ii_11 = send_map_starts_RT[i+1];
1465 for (j = send_map_starts_RT[i]; j < tmp_ii_11; j++)
1466 if (send_map_elmts_RT[j] == ic)
1467 {
1468 for (k=RAP_ext_i[j]; k < RAP_ext_i[j+1]; k++)
1469 {
1470 jcol = int_RAP_ext_j[k];
1471 if (jcol < num_cols_diag_P)
1472 {
1473 if (P_marker[jcol] < jj_row_begin_diag)
1474 {
1475 P_marker[jcol] = jj_count_diag;
1476 RAP_diag_data[jj_count_diag]
1477 = RAP_ext_data[k];
1478 RAP_diag_j[jj_count_diag] = jcol;
1479 jj_count_diag++;
1480 }
1481 else
1482 RAP_diag_data[P_marker[jcol]]
1483 += RAP_ext_data[k];
1484 }
1485 else
1486 {
1487 if (P_marker[jcol] < jj_row_begin_offd)
1488 {
1489 P_marker[jcol] = jj_count_offd;
1490 RAP_offd_data[jj_count_offd]
1491 = RAP_ext_data[k];
1492 RAP_offd_j[jj_count_offd]
1493 = jcol-num_cols_diag_P;
1494 jj_count_offd++;
1495 }
1496 else
1497 RAP_offd_data[P_marker[jcol]]
1498 += RAP_ext_data[k];
1499 }
1500 }
1501 break;
1502 }
1503 }
1504 /*--------------------------------------------------------------------
1505 * Loop over entries in row ic of R_diag.
1506 *--------------------------------------------------------------------*/
1507
1508 for (jj1 = R_diag_i[ic]; jj1 < R_diag_i[ic+1]; jj1++)
1509 {
1510 i1 = R_diag_j[jj1];
1511 r_entry = R_diag_data[jj1];
1512
1513 /*-----------------------------------------------------------------
1514 * Loop over entries in row i1 of A_offd.
1515 *-----------------------------------------------------------------*/
1516
1517 if (num_cols_offd_A)
1518 {
1519 for (jj2 = A_offd_i[i1]; jj2 < A_offd_i[i1+1]; jj2++)
1520 {
1521 i2 = A_offd_j[jj2];
1522 r_a_product = r_entry * A_offd_data[jj2];
1523
1524 /*--------------------------------------------------------------
1525 * Check A_marker to see if point i2 has been previously
1526 * visited. New entries in RAP only occur from unmarked points.
1527 *--------------------------------------------------------------*/
1528
1529 if (A_marker[i2] != ic)
1530 {
1531
1532 /*-----------------------------------------------------------
1533 * Mark i2 as visited.
1534 *-----------------------------------------------------------*/
1535
1536 A_marker[i2] = ic;
1537
1538 /*-----------------------------------------------------------
1539 * Loop over entries in row i2 of P_ext.
1540 *-----------------------------------------------------------*/
1541
1542 for (jj3 = P_ext_diag_i[i2]; jj3 < P_ext_diag_i[i2+1]; jj3++)
1543 {
1544 i3 = P_ext_diag_j[jj3];
1545 r_a_p_product = r_a_product * P_ext_diag_data[jj3];
1546
1547 /*--------------------------------------------------------
1548 * Check P_marker to see that RAP_{ic,i3} has not already
1549 * been accounted for. If it has not, create a new entry.
1550 * If it has, add new contribution.
1551 *--------------------------------------------------------*/
1552 if (P_marker[i3] < jj_row_begin_diag)
1553 {
1554 P_marker[i3] = jj_count_diag;
1555 RAP_diag_data[jj_count_diag] = r_a_p_product;
1556 RAP_diag_j[jj_count_diag] = i3;
1557 jj_count_diag++;
1558 }
1559 else
1560 RAP_diag_data[P_marker[i3]] += r_a_p_product;
1561 }
1562 for (jj3 = P_ext_offd_i[i2]; jj3 < P_ext_offd_i[i2+1]; jj3++)
1563 {
1564 i3 = map_Pext_to_RAP[P_ext_offd_j[jj3]] + num_cols_diag_P;
1565 r_a_p_product = r_a_product * P_ext_offd_data[jj3];
1566
1567 /*--------------------------------------------------------
1568 * Check P_marker to see that RAP_{ic,i3} has not already
1569 * been accounted for. If it has not, create a new entry.
1570 * If it has, add new contribution.
1571 *--------------------------------------------------------*/
1572 if (P_marker[i3] < jj_row_begin_offd)
1573 {
1574 P_marker[i3] = jj_count_offd;
1575 RAP_offd_data[jj_count_offd] = r_a_p_product;
1576 RAP_offd_j[jj_count_offd] = i3 - num_cols_diag_P;
1577 jj_count_offd++;
1578 }
1579 else
1580 RAP_offd_data[P_marker[i3]] += r_a_p_product;
1581 }
1582 }
1583
1584 /*--------------------------------------------------------------
1585 * If i2 is previously visited ( A_marker[12]=ic ) it yields
1586 * no new entries in RAP and can just add new contributions.
1587 *--------------------------------------------------------------*/
1588 else
1589 {
1590 for (jj3 = P_ext_diag_i[i2]; jj3 < P_ext_diag_i[i2+1]; jj3++)
1591 {
1592 i3 = P_ext_diag_j[jj3];
1593 r_a_p_product = r_a_product * P_ext_diag_data[jj3];
1594 RAP_diag_data[P_marker[i3]] += r_a_p_product;
1595 }
1596 for (jj3 = P_ext_offd_i[i2]; jj3 < P_ext_offd_i[i2+1]; jj3++)
1597 {
1598 i3 = map_Pext_to_RAP[P_ext_offd_j[jj3]] + num_cols_diag_P;
1599 r_a_p_product = r_a_product * P_ext_offd_data[jj3];
1600 RAP_offd_data[P_marker[i3]] += r_a_p_product;
1601 }
1602 }
1603 }
1604 }
1605
1606 /*-----------------------------------------------------------------
1607 * Loop over entries in row i1 of A_diag.
1608 *-----------------------------------------------------------------*/
1609
1610 /* for (jj2 = A_diag_i[i1]; jj2 < A_diag_i[i1+1]; jj2++) */
1611
1612 tmp_ii = A_diag_i[i1+1];
1613 for (jj2 = A_diag_i[i1]; jj2 < tmp_ii; jj2++)
1614 {
1615 i2 = A_diag_j[jj2];
1616 r_a_product = r_entry * A_diag_data[jj2];
1617
1618 /*--------------------------------------------------------------
1619 * Check A_marker to see if point i2 has been previously
1620 * visited. New entries in RAP only occur from unmarked points.
1621 *--------------------------------------------------------------*/
1622
1623 if (A_marker[i2+num_cols_offd_A] != ic)
1624 {
1625 /*-----------------------------------------------------------
1626 * Mark i2 as visited.
1627 *-----------------------------------------------------------*/
1628
1629 A_marker[i2+num_cols_offd_A] = ic;
1630
1631 /*-----------------------------------------------------------
1632 * Loop over entries in row i2 of P_diag.
1633 *-----------------------------------------------------------*/
1634
1635 tmp_ii_1 = P_diag_i[i2+1];
1636 for (jj3 = P_diag_i[i2]; jj3 < tmp_ii_1; jj3++)
1637 {
1638 i3 = P_diag_j[jj3];
1639 r_a_p_product = r_a_product * P_diag_data[jj3];
1640
1641 /*--------------------------------------------------------
1642 * Check P_marker to see that RAP_{ic,i3} has not already
1643 * been accounted for. If it has not, create a new entry.
1644 * If it has, add new contribution.
1645 *--------------------------------------------------------*/
1646
1647 if (P_marker[i3] < jj_row_begin_diag)
1648 {
1649 P_marker[i3] = jj_count_diag;
1650 RAP_diag_data[jj_count_diag] = r_a_p_product;
1651 RAP_diag_j[jj_count_diag] = P_diag_j[jj3];
1652 jj_count_diag++;
1653 }
1654 else
1655 {
1656 RAP_diag_data[P_marker[i3]] += r_a_p_product;
1657 }
1658 }
1659 if (num_cols_offd_P)
1660 {
1661 for (jj3 = P_offd_i[i2]; jj3 < P_offd_i[i2+1]; jj3++)
1662 {
1663 i3 = map_P_to_RAP[P_offd_j[jj3]] + num_cols_diag_P;
1664 r_a_p_product = r_a_product * P_offd_data[jj3];
1665
1666 /*--------------------------------------------------------
1667 * Check P_marker to see that RAP_{ic,i3} has not already
1668 * been accounted for. If it has not, create a new entry.
1669 * If it has, add new contribution.
1670 *--------------------------------------------------------*/
1671
1672 if (P_marker[i3] < jj_row_begin_offd)
1673 {
1674 P_marker[i3] = jj_count_offd;
1675 RAP_offd_data[jj_count_offd] = r_a_p_product;
1676 RAP_offd_j[jj_count_offd] = i3 - num_cols_diag_P;
1677 jj_count_offd++;
1678 }
1679 else
1680 {
1681 RAP_offd_data[P_marker[i3]] += r_a_p_product;
1682 }
1683 }
1684 }
1685 }
1686
1687 /*--------------------------------------------------------------
1688 * If i2 is previously visited ( A_marker[12]=ic ) it yields
1689 * no new entries in RAP and can just add new contributions.
1690 *--------------------------------------------------------------*/
1691
1692 else
1693 {
1694 tmp_ii_2 = P_diag_i[i2+1];
1695
1696 for (jj3 = P_diag_i[i2]; jj3 < tmp_ii_2; jj3++)
1697 {
1698 i3 = P_diag_j[jj3];
1699 /* __prefetch_by_load((P_marker + i3)); */
1700 r_a_p_product = r_a_product * P_diag_data[jj3];
1701 RAP_diag_data[P_marker[i3]] += r_a_p_product;
1702 }
1703 if (num_cols_offd_P)
1704 {
1705 for (jj3 = P_offd_i[i2]; jj3 < P_offd_i[i2+1]; jj3++)
1706 {
1707 i3 = map_P_to_RAP[P_offd_j[jj3]] + num_cols_diag_P;
1708 r_a_p_product = r_a_product * P_offd_data[jj3];
1709 RAP_offd_data[P_marker[i3]] += r_a_p_product;
1710 }
1711 }
1712 }
1713 }
1714 }
1715 }
1716 hypre_TFree(P_mark_array[ii]);
1717 hypre_TFree(A_mark_array[ii]);
1718 }
1719
1720 /* check if really all off-diagonal entries occurring in col_map_offd_RAP
1721 are represented and eliminate if necessary */
1722
1723 P_marker = hypre_CTAlloc(int,num_cols_offd_RAP);
1724 for (i=0; i < num_cols_offd_RAP; i++)
1725 P_marker[i] = -1;
1726
1727 jj_count_offd = 0;
1728 for (i=0; i < RAP_offd_size; i++)
1729 {
1730 i3 = RAP_offd_j[i];
1731 if (P_marker[i3])
1732 {
1733 P_marker[i3] = 0;
1734 jj_count_offd++;
1735 }
1736 }
1737
1738 if (jj_count_offd < num_cols_offd_RAP)
1739 {
1740 new_col_map_offd_RAP = hypre_CTAlloc(HYPRE_BigInt,jj_count_offd);
1741 jj_counter = 0;
1742 for (i=0; i < num_cols_offd_RAP; i++)
1743 if (!P_marker[i])
1744 {
1745 P_marker[i] = jj_counter;
1746 new_col_map_offd_RAP[jj_counter++] = col_map_offd_RAP[i];
1747 }
1748
1749 for (i=0; i < RAP_offd_size; i++)
1750 {
1751 i3 = RAP_offd_j[i];
1752 RAP_offd_j[i] = P_marker[i3];
1753 }
1754
1755 num_cols_offd_RAP = jj_count_offd;
1756 hypre_TFree(col_map_offd_RAP);
1757 col_map_offd_RAP = new_col_map_offd_RAP;
1758 }
1759 hypre_TFree(P_marker);
1760
1761 RAP = hypre_ParCSRMatrixCreate(comm, n_coarse_RT, n_coarse,
1762 RT_partitioning, coarse_partitioning,
1763 num_cols_offd_RAP, RAP_diag_size,
1764 RAP_offd_size);
1765
1766/* Have RAP own coarse_partitioning instead of P */
1767 hypre_ParCSRMatrixSetColStartsOwner(P,0);
1768 hypre_ParCSRMatrixSetColStartsOwner(RT,0);
1769
1770 RAP_diag = hypre_ParCSRMatrixDiag(RAP);
1771 hypre_CSRMatrixI(RAP_diag) = RAP_diag_i;
1772 if (RAP_diag_size)
1773 {
1774 hypre_CSRMatrixData(RAP_diag) = RAP_diag_data;
1775 hypre_CSRMatrixJ(RAP_diag) = RAP_diag_j;
1776 }
1777
1778 RAP_offd = hypre_ParCSRMatrixOffd(RAP);
1779 hypre_CSRMatrixI(RAP_offd) = RAP_offd_i;
1780 if (num_cols_offd_RAP)
1781 {
1782 hypre_CSRMatrixData(RAP_offd) = RAP_offd_data;
1783 hypre_CSRMatrixJ(RAP_offd) = RAP_offd_j;
1784 hypre_ParCSRMatrixColMapOffd(RAP) = col_map_offd_RAP;
1785 }
1786 if (num_procs > 1)
1787 {
1788 /* hypre_GenerateRAPCommPkg(RAP, A); */
1789#ifdef HYPRE_NO_GLOBAL_PARTITION
1790 hypre_NewCommPkgCreate(RAP);
1791#else
1792 hypre_MatvecCommPkgCreate(RAP);
1793#endif
1794 }
1795
1796 *RAP_ptr = RAP;
1797
1798 /*-----------------------------------------------------------------------
1799 * Free R, P_ext and marker arrays.
1800 *-----------------------------------------------------------------------*/
1801
1802 hypre_CSRMatrixDestroy(R_diag);
1803 R_diag = NULL;
1804
1805 if (num_cols_offd_RT)
1806 {
1807 hypre_CSRMatrixDestroy(R_offd);
1808 R_offd = NULL;
1809 }
1810
1811 if (num_sends_RT || num_recvs_RT)
1812 {
1813 hypre_BigCSRMatrixDestroy(RAP_ext);
1814 RAP_ext = NULL;
1815 }
1816 hypre_TFree(P_mark_array);
1817 hypre_TFree(A_mark_array);
1818 hypre_TFree(P_ext_diag_i);
1819 hypre_TFree(P_ext_offd_i);
1820 hypre_TFree(jj_cnt_diag);
1821 hypre_TFree(jj_cnt_offd);
1822 hypre_TFree(int_RAP_ext_j);
1823 if (num_cols_offd_P)
1824 {
1825 hypre_TFree(map_P_to_Pext);
1826 hypre_TFree(map_P_to_RAP);
1827 }
1828 if (num_cols_offd_Pext)
1829 {
1830 hypre_TFree(col_map_offd_Pext);
1831 hypre_TFree(map_Pext_to_RAP);
1832 }
1833 if (P_ext_diag_size)
1834 {
1835 hypre_TFree(P_ext_diag_data);
1836 hypre_TFree(P_ext_diag_j);
1837 }
1838 if (P_ext_offd_size)
1839 {
1840 hypre_TFree(P_ext_offd_data);
1841 hypre_TFree(P_ext_offd_j);
1842 }
1843 return(0);
1844
1845}
1846
Note: See TracBrowser for help on using the repository browser.