Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
distrib_vbr_matrix.c
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#include <stdlib.h>
44#include <stdio.h>
45#include "paz_aztec.h"
46#include "prototypes.h"
47
48void distrib_vbr_matrix(int *proc_config,
49 int *N_global, int *N_blk_global,
50 int *n_nonzeros, int *n_blk_nonzeros,
51 int *N_update, int **update,
52 double **val, int **indx, int **rpntr, int **cpntr,
53 int **bpntr, int **bindx,
54 double **x, double **b, double **bt, double **xexact)
55#undef DEBUG
56
57{
58 int i, n_entries, N_columns, n_global_nonzeros, n_global_blk_nonzeros;
59 int N_local;
60 int ii, j, row, have_xexact = 0 ;
61 int kk = 0;
62 int max_ii = 0, max_jj = 0;
63 int ione = 1;
64 double value;
65 double *cnt;
66 int *rpntr1, *bindx1, *bpntr1, *indx1;
67 double *val1, *b1, *bt1, *x1, *xexact1;
68
69 printf("Processor %d of %d entering distrib_matrix.\n",
70 proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
71
72 /*************** Distribute global matrix to all processors ************/
73
74 if(proc_config[PAZ_node] == 0)
75 {
76 if ((*xexact) != NULL) have_xexact = 1;
77 printf("Broadcasting exact solution\n");
78 }
79
80 if(proc_config[PAZ_N_procs] > 1)
81 {
82
83 PAZ_broadcast((char *) N_global, sizeof(int), proc_config, PAZ_PACK);
84 PAZ_broadcast((char *) N_blk_global, sizeof(int), proc_config, PAZ_PACK);
85 PAZ_broadcast((char *) n_nonzeros, sizeof(int), proc_config, PAZ_PACK);
86 PAZ_broadcast((char *) n_blk_nonzeros, sizeof(int), proc_config, PAZ_PACK);
87 PAZ_broadcast((char *) &have_xexact, sizeof(int), proc_config, PAZ_PACK);
88 PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
89
90 printf("Processor %d of %d done with global parameter broadcast.\n",
91 proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
92
93 if(proc_config[PAZ_node] != 0)
94 {
95 *bpntr = (int *) calloc(*N_blk_global+1,sizeof(int)) ;
96 *rpntr = (int *) calloc(*N_blk_global+1,sizeof(int)) ;
97 *bindx = (int *) calloc(*n_blk_nonzeros+1,sizeof(int)) ;
98 *indx = (int *) calloc(*n_blk_nonzeros+1,sizeof(int)) ;
99 *val = (double *) calloc(*n_nonzeros+1,sizeof(double)) ;
100 printf("Processor %d of %d done with global calloc.\n",
101 proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
102}
103
104 PAZ_broadcast((char *) (*bpntr), sizeof(int) *(*N_blk_global+1),
105 proc_config, PAZ_PACK);
106 PAZ_broadcast((char *) (*rpntr), sizeof(int) *(*N_blk_global+1),
107 proc_config, PAZ_PACK);
108 PAZ_broadcast((char *) (*bindx), sizeof(int) *(*n_blk_nonzeros+1),
109 proc_config, PAZ_PACK);
110 PAZ_broadcast((char *) (*indx), sizeof(int) *(*n_blk_nonzeros+1),
111 proc_config, PAZ_PACK);
112 PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
113 PAZ_broadcast((char *) (*val), sizeof(double)*(*n_nonzeros+1),
114 proc_config, PAZ_PACK);
115 PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
116
117 printf("Processor %d of %d done with matrix broadcast.\n",
118 proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
119
120 /* Set rhs and initialize guess */
121 if(proc_config[PAZ_node] != 0)
122 {
123 (*b) = (double *) calloc(*N_global,sizeof(double)) ;
124 (*bt) = (double *) calloc(*N_global,sizeof(double)) ;
125 (*x) = (double *) calloc(*N_global,sizeof(double)) ;
126 if (have_xexact)
127 (*xexact) = (double *) calloc(*N_global,sizeof(double)) ;
128 }
129
130 PAZ_broadcast((char *) (*x), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
131 PAZ_broadcast((char *) (*b), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
132 PAZ_broadcast((char *) (*bt), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
133 if (have_xexact)
134 PAZ_broadcast((char *)
135 (*xexact), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
136 PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
137 printf("Processor %d of %d done with rhs/guess broadcast.\n",
138 proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
139
140 }
141
142 /********************** Generate update map *************************/
143
144 PAZ_read_update(N_update, update, proc_config, *N_blk_global,
145 1, PAZ_linear) ;
146
147 printf("Processor %d of %d has %d rows of %d total block rows.\n",
148 proc_config[PAZ_node],proc_config[PAZ_N_procs],*N_update,*N_blk_global) ;
149
150 /*************** Construct local matrix from global matrix ************/
151
152 /* The local matrix is a copy of the rows assigned to this processor.
153 It is stored in MSR format and still has global indices (PAZ_transform
154 will complete conversion to local indices.
155 */
156
157 if(proc_config[PAZ_N_procs] > 1)
158 {
159 n_global_nonzeros = *n_nonzeros;
160 n_global_blk_nonzeros = *n_blk_nonzeros;
161
162 *n_nonzeros = 0;
163 *n_blk_nonzeros = 0;
164 N_local = 0;
165
166 for (i=0; i<*N_update; i++)
167 {
168 row = (*update)[i];
169 *n_nonzeros += (*indx)[(*bpntr)[row+1]] - (*indx)[(*bpntr)[row]];
170 *n_blk_nonzeros += (*bpntr)[row+1] - (*bpntr)[row];
171 N_local += (*rpntr)[row+1] - (*rpntr)[row];
172
173 }
174
175 printf("Processor %d of %d has %d nonzeros of %d total nonzeros.\n",
176 proc_config[PAZ_node],proc_config[PAZ_N_procs],
177 *n_nonzeros,n_global_nonzeros) ;
178
179 printf("Processor %d of %d has %d block nonzeros of %d total block nonzeros.\n",
180 proc_config[PAZ_node],proc_config[PAZ_N_procs],
181 *n_blk_nonzeros,n_global_blk_nonzeros) ;
182
183 printf("Processor %d of %d has %d equations of %d total equations.\n",
184 proc_config[PAZ_node],proc_config[PAZ_N_procs],
185 N_local,*N_global) ;
186
187#ifdef DEBUG
188 { double sum1 = 0.0;
189 for (i=0;i<*N_global; i++) sum1 += (*b)[i];
190
191 printf("Processor %d of %d has sum of b = %12.4g.\n",
192 proc_config[PAZ_node],proc_config[PAZ_N_procs],sum1) ;
193 }
194#endif /* DEBUG */
195
196 /* Allocate memory for local matrix */
197
198 bpntr1 = (int *) calloc(*N_update+1,sizeof(int)) ;
199 rpntr1 = (int *) calloc(*N_update+1,sizeof(int)) ;
200 bindx1 = (int *) calloc(*n_blk_nonzeros+1,sizeof(int)) ;
201 indx1 = (int *) calloc(*n_blk_nonzeros+1,sizeof(int)) ;
202 val1 = (double *) calloc(*n_nonzeros+1,sizeof(double)) ;
203 b1 = (double *) calloc(N_local,sizeof(double)) ;
204 bt1 = (double *) calloc(N_local,sizeof(double)) ;
205 x1 = (double *) calloc(N_local,sizeof(double)) ;
206 if (have_xexact)
207 xexact1 = (double *) calloc(N_local,sizeof(double)) ;
208
209 {
210 int cur_blk_size, indx_offset, len_val, row_offset, row_offset1;
211 double *val_ptr, *val1_ptr;
212
213 bpntr1[0] = 0;
214 indx1[0] = 0;
215 rpntr1[0] = 0;
216 for (i=0; i<*N_update; i++)
217 {
218 row = (*update)[i];
219 cur_blk_size = (*rpntr)[row+1] - (*rpntr)[row];
220 rpntr1[i+1] = rpntr1[i] + cur_blk_size;
221 row_offset = (*rpntr)[row];
222 row_offset1 = rpntr1[i];
223 for (j = 0; j<cur_blk_size; j++)
224 {
225 b1[row_offset1+j] = (*b)[row_offset+j];
226 x1[row_offset1+j] = (*x)[row_offset+j];
227 if (have_xexact) xexact1[row_offset1+j] = (*xexact)[row_offset+j];
228 }
229 bpntr1[i+1] = bpntr1[i];
230
231#ifdef DEBUG
232 printf("Proc %d of %d: Global row = %d: Local row = %d:
233 b = %12.4g: x = %12.4g: bindx = %d: val = %12.4g \n",
234 proc_config[PAZ_node],proc_config[PAZ_N_procs],
235 row, i, b1[i], x1[i], bindx1[i], val1[i]) ;
236#endif
237 indx_offset = (*indx)[(*bpntr)[row]] - indx1[bpntr1[i]];
238 for (j = (*bpntr)[row]; j < (*bpntr)[row+1]; j++)
239 {
240 indx1[bpntr1 [i+1] + 1] = (*indx)[j+1] - indx_offset;
241 bindx1[bpntr1 [i+1] ] = (*bindx)[j];
242 bpntr1[i+1] ++;
243 }
244 len_val = indx1[bpntr1[i+1]] - indx1[bpntr1[i]];
245 val_ptr = (*val)+(*indx)[(*bpntr)[row]];
246 val1_ptr = val1+indx1[bpntr1[i]];
247 for (j = 0; j<len_val; j++)
248 {
249 *val1_ptr = *val_ptr;
250 val_ptr++; val1_ptr++;
251 }
252 }
253 }
254 printf("Processor %d of %d done with extracting local operators.\n",
255 proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
256
257 if (have_xexact)
258 {
259 printf(
260 "The residual using VBR format and exact solution on processor %d is %12.4g\n",
261 proc_config[PAZ_node],
262 svbrres (N_local, *N_global, *N_update, val1, indx1, bindx1,
263 rpntr1, (*rpntr), bpntr1, bpntr1+1,
264 (*xexact), b1));
265 }
266
267 /* Release memory for global matrix, rhs and solution */
268
269 free ((void *) (*val));
270 free ((void *) (*indx));
271 free ((void *) (*bindx));
272 free ((void *) (*bpntr));
273 free ((void *) (*rpntr));
274 free ((void *) (*b));
275 free ((void *) (*bt));
276 free ((void *) (*x));
277 if (have_xexact) free((void *) *xexact);
278
279 /* Return local matrix through same pointers. */
280
281 *val = val1;
282 *indx = indx1;
283 *bindx = bindx1;
284 *bpntr = bpntr1;
285 *rpntr = rpntr1;
286 *b = b1;
287 *bt = bt1;
288 *x = x1;
289 if (have_xexact) *xexact = xexact1;
290
291 }
292 if (have_xexact && proc_config[PAZ_N_procs] == 1)
293 {
294 printf(
295 "The residual using VBR format and exact solution on processor %d is %12.4g\n",
296 proc_config[PAZ_node],
297 svbrres (*N_global, *N_global, *N_update, (*val), (*indx), (*bindx),
298 (*rpntr), (*rpntr), (*bpntr), (*bpntr)+1,
299 (*xexact), (*b)));
300 }
301
302
303 printf("Processor %d of %d leaving distrib_matrix.\n",
304 proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
305
306 /* end distrib_matrix */
307}
void distrib_vbr_matrix(int *proc_config, int *N_global, int *N_blk_global, int *n_nonzeros, int *n_blk_nonzeros, int *N_update, int **update, double **val, int **indx, int **rpntr, int **cpntr, int **bpntr, int **bindx, double **x, double **b, double **bt, double **xexact)
double svbrres(int m, int n, int m_blk, double *val, int *indx, int *bindx, int *rpntr, int *cpntr, int *bpntrb, int *bpntre, double *x, double *b)
Definition: svbrres.c:49