source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_jacobi_interp.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: 21.8 KB
RevLine 
[2aa6644]1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13#include "headers.h"
14#include "par_amg.h"
15
16/* #define HYPRE_JACINT_PRINT_ROW_SUMS*/
17/* #define HYPRE_JACINT_PRINT_SOME_ROWS */
18/* #define HYPRE_JACINT_PRINT_MATRICES*/
19#define HYPRE_MAX_PRINTABLE_MATRIX 125
20/*#define HYPRE_JACINT_PRINT_DIAGNOSTICS*/
21
22void hypre_BoomerAMGJacobiInterp( hypre_ParCSRMatrix * A,
23 hypre_ParCSRMatrix ** P,
24 hypre_ParCSRMatrix * S,
25 int num_functions, int * dof_func,
26 int * CF_marker, int level,
27 double truncation_threshold,
28 double truncation_threshold_minus )
29/* nji steps of Jacobi interpolation, with nji presently just set in the code.*/
30{
31 double weight_AF = 1.0; /* weight multiplied by A's fine row elements */
32 int * dof_func_offd = NULL;
33 int nji = 1;
34 int iji;
35
36 hypre_ParCSRMatrix_dof_func_offd( A,
37 num_functions,
38 dof_func,
39 &dof_func_offd );
40
41 for ( iji=0; iji<nji; ++iji )
42 {
43 hypre_BoomerAMGJacobiInterp_1( A, P, S, CF_marker, level,
44 truncation_threshold, truncation_threshold_minus,
45 dof_func, dof_func_offd,
46 weight_AF );
47 }
48
49 if ( dof_func_offd != NULL )
50 hypre_TFree( dof_func_offd );
51}
52
53void hypre_BoomerAMGJacobiInterp_1( hypre_ParCSRMatrix * A,
54 hypre_ParCSRMatrix ** P,
55 hypre_ParCSRMatrix * S,
56 int * CF_marker, int level,
57 double truncation_threshold,
58 double truncation_threshold_minus,
59 int * dof_func, int * dof_func_offd,
60 double weight_AF
61 )
62 /*A is the linear system.
63 P is an interpolation matrix, input and output
64 CF_marker identifies coarse and fine points
65 If we imagine P and A as split into coarse and fine submatrices,
66
67 [ AFF AFC ] [ AF ] [ IFC ]
68 A = [ ] = [ ] , P = [ ]
69 [ ACF ACC ] [ AC ] [ ICC ]
70 (note that ICC is an identity matrix, applied to coarse points only)
71 then this function computes
72
73 IFCnew = IFCold - DFF(-1) * ( AFF*IFCold + AFC )
74 = IFCold - DFF(-1) * AF * Pold)
75 where DFF is the diagonal of AFF, (-1) represents the inverse, and
76 where "old" denotes a value on entry to this function, "new" a returned value.
77
78 */
79{
80 hypre_ParCSRMatrix * Pnew;
81 hypre_ParCSRMatrix * C;
82 hypre_CSRMatrix *P_diag = hypre_ParCSRMatrixDiag(*P);
83 hypre_CSRMatrix *P_offd = hypre_ParCSRMatrixOffd(*P);
84 double *P_diag_data = hypre_CSRMatrixData(P_diag);
85 int *P_diag_i = hypre_CSRMatrixI(P_diag);
86 int *P_diag_j = hypre_CSRMatrixJ(P_diag);
87 double *P_offd_data = hypre_CSRMatrixData(P_offd);
88 int *P_offd_i = hypre_CSRMatrixI(P_offd);
89 hypre_CSRMatrix *C_diag;
90 hypre_CSRMatrix *C_offd;
91 hypre_CSRMatrix *Pnew_diag;
92 hypre_CSRMatrix *Pnew_offd;
93 int num_rows_diag_P = hypre_CSRMatrixNumRows(P_diag);
94 int i;
95 int Jnochanges=0, Jchanges, Pnew_num_nonzeros;
96 int CF_coarse=0;
97 int * J_marker = hypre_CTAlloc( int, num_rows_diag_P );
98 int nc, ncmax, ncmin, nc1;
99 int num_procs, my_id;
100 MPI_Comm comm = hypre_ParCSRMatrixComm( A );
101#ifdef HYPRE_JACINT_PRINT_ROW_SUMS
102 int m, nmav, npav;
103 double PIi, PIimax, PIimin, PIimav, PIipav, randthresh;
104 double eps = 1.0e-17;
105#endif
106#ifdef HYPRE_JACINT_PRINT_MATRICES
107 char filename[80];
108 int i_dummy, j_dummy;
109 int *base_i_ptr = &i_dummy;
110 int *base_j_ptr = &j_dummy;
111#endif
112#ifdef HYPRE_JACINT_PRINT_SOME_ROWS
113 int sample_rows[50], n_sample_rows=0, isamp;
114#endif
115
116 MPI_Comm_size(comm, &num_procs);
117 MPI_Comm_rank(comm,&my_id);
118
119
120 for ( i=0; i<num_rows_diag_P; ++i )
121 {
122 J_marker[i] = CF_marker[i];
123 if (CF_marker[i]>=0) ++CF_coarse;
124 }
125#ifdef HYPRE_JACINT_PRINT_DIAGNOSTICS
126 printf("%i %i Jacobi_Interp_1, P has %i+%i=%i nonzeros, local sum %e\n", my_id, level,
127 hypre_CSRMatrixNumNonzeros(P_diag), hypre_CSRMatrixNumNonzeros(P_offd),
128 hypre_CSRMatrixNumNonzeros(P_diag)+hypre_CSRMatrixNumNonzeros(P_offd),
129 hypre_ParCSRMatrixLocalSumElts(*P) );
130#endif
131
132 /* row sum computations, for output */
133#ifdef HYPRE_JACINT_PRINT_ROW_SUMS
134 PIimax=-1.0e12, PIimin=1.0e12, PIimav=0, PIipav=0;
135 nmav=0, npav=0;
136 for ( i=0; i<num_rows_diag_P; ++i )
137 {
138 PIi = 0; /* i-th value of P*1, i.e. sum of row i of P */
139 for ( m=P_diag_i[i]; m<P_diag_i[i+1]; ++m )
140 PIi += P_diag_data[m];
141 for ( m=P_offd_i[i]; m<P_offd_i[i+1]; ++m )
142 PIi += P_offd_data[m];
143 if (CF_marker[i]<0)
144 {
145 PIimax = hypre_max( PIimax, PIi );
146 PIimin = hypre_min( PIimin, PIi );
147 if (PIi<=1-eps) { PIimav+=PIi; ++nmav; };
148 if (PIi>=1+eps) { PIipav+=PIi; ++npav; };
149 }
150 }
151 if ( nmav>0 ) PIimav = PIimav/nmav;
152 if ( npav>0 ) PIipav = PIipav/npav;
153 printf("%i %i P in max,min row sums %e %e\n", my_id, level, PIimax, PIimin );
154#endif
155
156 ncmax=0; ncmin=num_rows_diag_P; nc1=0;
157 for ( i=0; i<num_rows_diag_P; ++i )
158 if (CF_marker[i]<0)
159 {
160 nc = P_diag_i[i+1] - P_diag_i[i];
161 if (nc<=1)
162 {
163 ++nc1;
164 }
165 ncmax = hypre_max( nc, ncmax );
166 ncmin = hypre_min( nc, ncmin );
167 }
168#if 0
169 /* a very agressive reduction in how much the Jacobi step does: */
170 for ( i=0; i<num_rows_diag_P; ++i )
171 if (CF_marker[i]<0)
172 {
173 nc = P_diag_i[i+1] - P_diag_i[i];
174 if (nc>ncmin+1)
175 /*if ( nc > ncmin + 0.5*(ncmax-ncmin) )*/
176 {
177 J_marker[i] = 1;
178 ++Jnochanges;
179 }
180 }
181#endif
182 Jchanges = num_rows_diag_P - Jnochanges - CF_coarse;
183
184#ifdef HYPRE_JACINT_PRINT_SOME_ROWS
185 printf("some rows to be changed: ");
186 randthresh = 15/(double)Jchanges;
187 for ( i=0; i<num_rows_diag_P; ++i )
188 {
189 if ( J_marker[i]<0 )
190 {
191 if ( ((double)rand())/RAND_MAX < randthresh )
192 {
193 printf( "%i: ", i );
194 for ( m=P_diag_i[i]; m<P_diag_i[i+1]; ++m )
195 printf( " %i %f, ", P_diag_j[m], P_diag_data[m] );
196 printf("; ");
197 sample_rows[n_sample_rows] = i;
198 ++n_sample_rows;
199 }
200 }
201 }
202 printf("\n");
203#endif
204#ifdef HYPRE_JACINT_PRINT_DIAGNOSTICS
205 printf("%i %i P has %i rows, %i changeable, %i don't change-good, %i coarse\n",
206 my_id, level, num_rows_diag_P, Jchanges, Jnochanges, CF_coarse );
207 printf("%i %i min,max diag cols per row: %i, %i; no.rows w.<=1 col: %i\n", my_id, level, ncmin, ncmax, nc1 );
208#endif
209#ifdef HYPRE_JACINT_PRINT_MATRICES
210 if ( num_rows_diag_P <= HYPRE_MAX_PRINTABLE_MATRIX )
211 {
212 sprintf( filename, "Ain%i", level );
213 hypre_ParCSRMatrixPrintIJ( A,0,0,filename);
214 sprintf( filename, "Sin%i", level );
215 hypre_ParCSRMatrixPrintIJ( S,0,0,filename);
216 sprintf( filename, "Pin%i", level );
217 hypre_ParCSRMatrixPrintIJ( *P,0,0,filename);
218 }
219#endif
220
221 C = hypre_ParMatmul_FC( A, *P, J_marker, dof_func, dof_func_offd );
222 /* hypre_parMatmul_FC creates and returns C, a variation of the
223 matrix product A*P in which only the "Fine"-designated rows have
224 been computed. (all columns are Coarse because all columns of P
225 are). "Fine" is defined solely by the marker array, and for
226 example could be a proper subset of the fine points of a
227 multigrid hierarchy.
228 As a matrix, C is the size of A*P. But only the marked rows have
229 been computed.
230 */
231#ifdef HYPRE_JACINT_PRINT_MATRICES
232 sprintf( filename, "C%i", level );
233 if ( num_rows_diag_P <= HYPRE_MAX_PRINTABLE_MATRIX ) hypre_ParCSRMatrixPrintIJ( C,0,0,filename);
234#endif
235 C_diag = hypre_ParCSRMatrixDiag(C);
236 C_offd = hypre_ParCSRMatrixOffd(C);
237#ifdef HYPRE_JACINT_PRINT_DIAGNOSTICS
238 printf("%i %i Jacobi_Interp_1 after matmul, C has %i+%i=%i nonzeros, local sum %e\n",
239 my_id, level, hypre_CSRMatrixNumNonzeros(C_diag),
240 hypre_CSRMatrixNumNonzeros(C_offd),
241 hypre_CSRMatrixNumNonzeros(C_diag)+hypre_CSRMatrixNumNonzeros(C_offd),
242 hypre_ParCSRMatrixLocalSumElts(C) );
243#endif
244
245 hypre_ParMatScaleDiagInv_F( C, A, weight_AF, J_marker );
246 /* hypre_ParMatScaleDiagInv scales of its first argument by premultiplying with
247 a submatrix of the inverse of the diagonal of its second argument.
248 The marker array determines which diagonal elements are used. The marker
249 array should select exactly the right number of diagonal elements (the number
250 of rows of AP_FC).
251 */
252#ifdef HYPRE_JACINT_PRINT_MATRICES
253 sprintf( filename, "Cout%i", level );
254 if ( num_rows_diag_P <= HYPRE_MAX_PRINTABLE_MATRIX ) hypre_ParCSRMatrixPrintIJ( C,0,0,filename);
255#endif
256
257 Pnew = hypre_ParMatMinus_F( *P, C, J_marker );
258 /* hypre_ParMatMinus_F subtracts rows of its second argument from selected rows
259 of its first argument. The marker array determines which rows of the first
260 argument are affected, and they should exactly correspond to all the rows
261 of the second argument.
262 */
263 Pnew_diag = hypre_ParCSRMatrixDiag(Pnew);
264 Pnew_offd = hypre_ParCSRMatrixOffd(Pnew);
265 Pnew_num_nonzeros = hypre_CSRMatrixNumNonzeros(Pnew_diag)+hypre_CSRMatrixNumNonzeros(Pnew_offd);
266#ifdef HYPRE_JACINT_PRINT_DIAGNOSTICS
267 printf("%i %i Jacobi_Interp_1 after MatMinus, Pnew has %i+%i=%i nonzeros, local sum %e\n",
268 my_id, level, hypre_CSRMatrixNumNonzeros(Pnew_diag),
269 hypre_CSRMatrixNumNonzeros(Pnew_offd), Pnew_num_nonzeros,
270 hypre_ParCSRMatrixLocalSumElts(Pnew) );
271#endif
272
273 /* Transfer ownership of col_starts from P to Pnew ... */
274 if ( hypre_ParCSRMatrixColStarts(*P) &&
275 hypre_ParCSRMatrixColStarts(*P)==hypre_ParCSRMatrixColStarts(Pnew) )
276 {
277 if ( hypre_ParCSRMatrixOwnsColStarts(*P) && !hypre_ParCSRMatrixOwnsColStarts(Pnew) )
278 {
279 hypre_ParCSRMatrixSetColStartsOwner(*P,0);
280 hypre_ParCSRMatrixSetColStartsOwner(Pnew,1);
281 }
282 }
283
284 hypre_ParCSRMatrixDestroy( C );
285 hypre_ParCSRMatrixDestroy( *P );
286
287 /* Note that I'm truncating all the fine rows, not just the J-marked ones. */
288#if 0
289 if ( Pnew_num_nonzeros < 10000 ) /* a fixed number like this makes it no.procs.-depdendent */
290 { /* ad-hoc attempt to reduce zero-matrix problems seen in testing..*/
291 truncation_threshold = 1.0e-6 * truncation_threshold;
292 truncation_threshold_minus = 1.0e-6 * truncation_threshold_minus;
293 }
294#endif
295 hypre_BoomerAMGTruncateInterp( Pnew, truncation_threshold,
296 truncation_threshold_minus, CF_marker );
297
298 hypre_MatvecCommPkgCreate ( Pnew );
299
300
301 *P = Pnew;
302
303 P_diag = hypre_ParCSRMatrixDiag(*P);
304 P_offd = hypre_ParCSRMatrixOffd(*P);
305 P_diag_data = hypre_CSRMatrixData(P_diag);
306 P_diag_i = hypre_CSRMatrixI(P_diag);
307 P_diag_j = hypre_CSRMatrixJ(P_diag);
308 P_offd_data = hypre_CSRMatrixData(P_offd);
309 P_offd_i = hypre_CSRMatrixI(P_offd);
310
311 /* row sum computations, for output */
312#ifdef HYPRE_JACINT_PRINT_ROW_SUMS
313 PIimax=-1.0e12, PIimin=1.0e12, PIimav=0, PIipav=0;
314 nmav=0, npav=0;
315 for ( i=0; i<num_rows_diag_P; ++i )
316 {
317 PIi = 0; /* i-th value of P*1, i.e. sum of row i of P */
318 for ( m=P_diag_i[i]; m<P_diag_i[i+1]; ++m )
319 PIi += P_diag_data[m];
320 for ( m=P_offd_i[i]; m<P_offd_i[i+1]; ++m )
321 PIi += P_offd_data[m];
322 if (CF_marker[i]<0)
323 {
324 PIimax = hypre_max( PIimax, PIi );
325 PIimin = hypre_min( PIimin, PIi );
326 if (PIi<=1-eps) { PIimav+=PIi; ++nmav; };
327 if (PIi>=1+eps) { PIipav+=PIi; ++npav; };
328 }
329 }
330 if ( nmav>0 ) PIimav = PIimav/nmav;
331 if ( npav>0 ) PIipav = PIipav/npav;
332 printf("%i %i P out max,min row sums %e %e\n", my_id, level, PIimax, PIimin );
333#endif
334
335#ifdef HYPRE_JACINT_PRINT_SOME_ROWS
336 printf("some changed rows: ");
337 for ( isamp=0; isamp<n_sample_rows; ++isamp )
338 {
339 i = sample_rows[isamp];
340 printf( "%i: ", i );
341 for ( m=P_diag_i[i]; m<P_diag_i[i+1]; ++m )
342 printf( " %i %f, ", P_diag_j[m], P_diag_data[m] );
343 printf("; ");
344 }
345 printf("\n");
346#endif
347 ncmax=0; ncmin=num_rows_diag_P; nc1=0;
348 for ( i=0; i<num_rows_diag_P; ++i )
349 if (CF_marker[i]<0)
350 {
351 nc = P_diag_i[i+1] - P_diag_i[i];
352 if (nc<=1) ++nc1;
353 ncmax = hypre_max( nc, ncmax );
354 ncmin = hypre_min( nc, ncmin );
355 }
356#ifdef HYPRE_JACINT_PRINT_DIAGNOSTICS
357 printf("%i %i P has %i rows, %i changeable, %i too good, %i coarse\n",
358 my_id, level, num_rows_diag_P, num_rows_diag_P-Jnochanges-CF_coarse, Jnochanges, CF_coarse );
359 printf("%i %i min,max diag cols per row: %i, %i; no.rows w.<=1 col: %i\n", my_id, level, ncmin, ncmax, nc1 );
360
361 printf("%i %i Jacobi_Interp_1 after truncation (%e), Pnew has %i+%i=%i nonzeros, local sum %e\n",
362 my_id, level, truncation_threshold,
363 hypre_CSRMatrixNumNonzeros(Pnew_diag), hypre_CSRMatrixNumNonzeros(Pnew_offd),
364 hypre_CSRMatrixNumNonzeros(Pnew_diag)+hypre_CSRMatrixNumNonzeros(Pnew_offd),
365 hypre_ParCSRMatrixLocalSumElts(Pnew) );
366#endif
367
368 /* Programming Notes:
369 1. Judging by around line 299 of par_interp.c, they typical use of CF_marker
370 is that CF_marker>=0 means Coarse, CF_marker<0 means Fine.
371 */
372#ifdef HYPRE_JACINT_PRINT_MATRICES
373 sprintf( filename, "Pout%i", level );
374 if ( num_rows_diag_P <= HYPRE_MAX_PRINTABLE_MATRIX ) hypre_ParCSRMatrixPrintIJ( *P,0,0,filename);
375#endif
376
377 hypre_TFree( J_marker );
378
379}
380
381void hypre_BoomerAMGTruncateInterp( hypre_ParCSRMatrix *P,
382 double eps, double dlt,
383 int * CF_marker )
384/* Truncate the interpolation matrix P, but only in rows for which the
385 marker is <0. Truncation means that an element P(i,j) is set to 0 if
386 P(i,j)>0 and P(i,j)<eps*max( P(i,j) ) or if
387 P(i,j)>0 and P(i,j)<dlt*max( -P(i,j) ) or if
388 P(i,j)<0 and P(i,j)>dlt*min( -P(i,j) ) or if
389 P(i,j)<0 and P(i,j)>eps*min( P(i,j) )
390 ( 0<eps,dlt<1, typically 0.1=dlt<eps=0.2, )
391 The min and max are only computed locally, as I'm guessing that there isn't
392 usually much to be gained (in the way of improved performance) by getting
393 them perfectly right.
394*/
395
396/* The function hypre_BoomerAMGInterpTruncation in par_interp.c is
397 very similar. It looks at fabs(value) rather than separately
398 dealing with value<0 and value>0 as recommended by Klaus Stuben,
399 thus as this function does. In this function, only "marked" rows
400 are affected. Lastly, in hypre_BoomerAMGInterpTruncation, if any
401 element gets discarded, it reallocates arrays to the new size.
402*/
403{
404 hypre_CSRMatrix *P_diag = hypre_ParCSRMatrixDiag(P);
405 hypre_CSRMatrix *P_offd = hypre_ParCSRMatrixOffd(P);
406 double *P_diag_data = hypre_CSRMatrixData(P_diag);
407 int *P_diag_i = hypre_CSRMatrixI(P_diag);
408 int *P_diag_j = hypre_CSRMatrixJ(P_diag);
409 double *P_offd_data = hypre_CSRMatrixData(P_offd);
410 int *P_offd_i = hypre_CSRMatrixI(P_offd);
411 int *P_offd_j = hypre_CSRMatrixJ(P_offd);
412 int *new_P_diag_i;
413 int *new_P_offd_i;
414 int num_rows_diag_P = hypre_CSRMatrixNumRows(P_diag);
415 int num_rows_offd_P = hypre_CSRMatrixNumRows(P_offd);
416 int num_nonzeros_diag = hypre_CSRMatrixNumNonzeros(P_diag);
417 int num_nonzeros_offd = hypre_CSRMatrixNumNonzeros(P_offd);
418#if 0
419 MPI_Comm comm = hypre_ParCSRMatrixComm( P );
420 double vmax1, vmin1;
421#endif
422 double vmax = 0.0;
423 double vmin = 0.0;
424 double v, old_sum, new_sum, scale, wmax, wmin;
425 int i1, m, m1d, m1o;
426
427 /* compute vmax = eps*max(P(i,j)), vmin = eps*min(P(i,j)) */
428 for ( i1 = 0; i1 < num_rows_diag_P; i1++ )
429 {
430 for ( m=P_diag_i[i1]; m<P_diag_i[i1+1]; ++m )
431 {
432 v = P_diag_data[m];
433 vmax = hypre_max( v, vmax );
434 vmin = hypre_min( v, vmin );
435 }
436 for ( m=P_offd_i[i1]; m<P_offd_i[i1+1]; ++m )
437 {
438 v = P_offd_data[m];
439 vmax = hypre_max( v, vmax );
440 vmin = hypre_min( v, vmin );
441 }
442 }
443#if 0
444 /* This can make max,min global so results don't depend on no. processors
445 We don't want this except for testing, or maybe this could be put
446 someplace better. I don't like adding communication here, for a minor reason.
447 */
448 vmax1 = vmax; vmin1 = vmin;
449 MPI_Allreduce( &vmax1, &vmax, 1, MPI_DOUBLE, MPI_MAX, comm );
450 MPI_Allreduce( &vmin1, &vmin, 1, MPI_DOUBLE, MPI_MIN, comm );
451#endif
452 if ( vmax <= 0.0 ) vmax = 1.0; /* make sure no v is v>vmax if no v is v>0 */
453 if ( vmin >= 0.0 ) vmin = -1.0; /* make sure no v is v<vmin if no v is v<0 */
454 wmax = - dlt * vmin;
455 wmin = - dlt * vmax;
456 vmax *= eps;
457 vmin *= eps;
458
459 /* Repack the i,j,and data arrays so as to discard the small elements of P.
460 Elements of Coarse rows (CF_marker>=0) are always kept.
461 The arrays are not re-allocated, so there will generally be unused space
462 at the ends of the arrays. */
463 new_P_diag_i = hypre_CTAlloc( int, num_rows_diag_P+1 );
464 new_P_offd_i = hypre_CTAlloc( int, num_rows_offd_P+1 );
465 m1d = P_diag_i[0];
466 m1o = P_offd_i[0];
467 for ( i1 = 0; i1 < num_rows_diag_P; i1++ )
468 {
469 old_sum = 0;
470 new_sum = 0;
471 for ( m=P_diag_i[i1]; m<P_diag_i[i1+1]; ++m )
472 {
473 v = P_diag_data[m];
474 old_sum += v;
475 if ( CF_marker[i1]>=0 || ( v>=vmax && v>=wmax ) || ( v<=vmin && v<=wmin ) )
476 { /* keep v */
477 new_sum += v;
478 P_diag_j[m1d] = P_diag_j[m];
479 P_diag_data[m1d] = P_diag_data[m];
480 ++m1d;
481 }
482 else
483 { /* discard v */
484 --num_nonzeros_diag;
485 }
486 }
487 for ( m=P_offd_i[i1]; m<P_offd_i[i1+1]; ++m )
488 {
489 v = P_offd_data[m];
490 old_sum += v;
491 if ( CF_marker[i1]>=0 || ( v>=vmax && v>=wmax ) || ( v<=vmin && v<=wmin ) )
492 { /* keep v */
493 new_sum += v;
494 P_offd_j[m1o] = P_offd_j[m];
495 P_offd_data[m1o] = P_offd_data[m];
496 ++m1o;
497 }
498 else
499 { /* discard v */
500 --num_nonzeros_offd;
501 }
502 }
503
504 new_P_diag_i[i1+1] = m1d;
505 if ( i1<num_rows_offd_P ) new_P_offd_i[i1+1] = m1o;
506
507 /* rescale to keep row sum the same */
508 if (new_sum!=0) scale = old_sum/new_sum; else scale = 1.0;
509 for ( m=new_P_diag_i[i1]; m<new_P_diag_i[i1+1]; ++m )
510 P_diag_data[m] *= scale;
511 if ( i1<num_rows_offd_P ) /* this test fails when there is no offd block */
512 for ( m=new_P_offd_i[i1]; m<new_P_offd_i[i1+1]; ++m )
513 P_offd_data[m] *= scale;
514
515 }
516
517 for ( i1 = 1; i1 <= num_rows_diag_P; i1++ )
518 {
519 P_diag_i[i1] = new_P_diag_i[i1];
520 if ( i1<=num_rows_offd_P && num_nonzeros_offd>0 ) P_offd_i[i1] = new_P_offd_i[i1];
521 }
522 hypre_TFree( new_P_diag_i );
523 if ( num_rows_offd_P>0 ) hypre_TFree( new_P_offd_i );
524
525 hypre_CSRMatrixNumNonzeros(P_diag) = num_nonzeros_diag;
526 hypre_CSRMatrixNumNonzeros(P_offd) = num_nonzeros_offd;
527 hypre_ParCSRMatrixSetDNumNonzeros( P );
528 hypre_ParCSRMatrixSetNumNonzeros( P );
529
530}
531
532
533
534/*
535 hypre_ParCSRMatrix_dof_func_offd allocates, computes and returns dof_func_offd.
536 The caller is responsible for freeing dof_func_offd.
537 This function has code copied from hypre_BoomerAMGCreateS and hypre_BoomerAMGCreateSabs
538 They should be retrofitted to call this function. Or, better, call this function separately
539 and pass the result into them through an argument (less communication, less computation).
540*/
541
542int
543hypre_ParCSRMatrix_dof_func_offd(
544 hypre_ParCSRMatrix *A,
545 int num_functions,
546 int *dof_func,
547 int **dof_func_offd )
548{
549 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
550 hypre_ParCSRCommHandle *comm_handle;
551 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
552
553 int num_cols_offd = 0;
554 int Solve_err_flag = 0;
555 int num_sends;
556 int *int_buf_data;
557 int index, start, i, j;
558
559 num_cols_offd = hypre_CSRMatrixNumCols(A_offd);
560 *dof_func_offd = NULL;
561 if (num_cols_offd)
562 {
563 if (num_functions > 1)
564 *dof_func_offd = hypre_CTAlloc(int, num_cols_offd);
565 }
566
567
568 /*-------------------------------------------------------------------
569 * Get the dof_func data for the off-processor columns
570 *-------------------------------------------------------------------*/
571
572 if (!comm_pkg)
573 {
574#ifdef HYPRE_NO_GLOBAL_PARTITION
575 hypre_NewCommPkgCreate(A);
576#else
577 hypre_MatvecCommPkgCreate(A);
578#endif
579 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
580 }
581
582 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
583 if (num_functions > 1)
584 {
585 int_buf_data = hypre_CTAlloc(int,hypre_ParCSRCommPkgSendMapStart(comm_pkg,
586 num_sends));
587 index = 0;
588 for (i = 0; i < num_sends; i++)
589 {
590 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
591 for (j=start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
592 int_buf_data[index++]
593 = dof_func[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
594 }
595
596 comm_handle = hypre_ParCSRCommHandleCreate( 11, comm_pkg, int_buf_data,
597 *dof_func_offd);
598
599 hypre_ParCSRCommHandleDestroy(comm_handle);
600 hypre_TFree(int_buf_data);
601 }
602
603 return(Solve_err_flag);
604}
Note: See TracBrowser for help on using the repository browser.