Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
48 void 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 }
svbrres
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
distrib_vbr_matrix
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)
Definition: distrib_vbr_matrix.c:48
Generated by   doxygen 1.8.5