Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
distrib_msr_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_msr_matrix(int *proc_config, int *N_global,
49  int *n_nonzeros, int *N_update, int **update,
50  double **val, int **bindx,
51  double **x, double **b, double **bt, double **xexact)
52 #undef DEBUG
53 
54 {
55  int i, n_entries, N_columns, n_global_nonzeros;
56  int ii, j, row, have_xexact = 0 ;
57  int kk = 0;
58  int max_ii = 0, max_jj = 0;
59  int ione = 1;
60  double value;
61  double *cnt;
62  int *pntr, *bindx1, *pntr1;
63  double *val1, *b1, *bt1, *x1, *xexact1;
64 
65  printf("Processor %d of %d entering distrib_matrix.\n",
66  proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
67 
68  /*************** Distribute global matrix to all processors ************/
69 
70  if(proc_config[PAZ_node] == 0)
71  {
72  if ((*xexact) != NULL) have_xexact = 1;
73  printf("Broadcasting exact solution\n");
74  }
75 
76  if(proc_config[PAZ_N_procs] > 1) {
77 
78  PAZ_broadcast((char *) N_global, sizeof(int), proc_config, PAZ_PACK);
79  PAZ_broadcast((char *) n_nonzeros, sizeof(int), proc_config, PAZ_PACK);
80  PAZ_broadcast((char *) &have_xexact, sizeof(int), proc_config, PAZ_PACK);
81  PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
82 
83  if(proc_config[PAZ_node] != 0)
84  {
85  (*bindx) = (int *) calloc(*n_nonzeros+1,sizeof(int)) ;
86  (*val) = (double *) calloc(*n_nonzeros+1,sizeof(double)) ;
87  }
88 
89  PAZ_broadcast((char *) (*bindx), sizeof(int) *(*n_nonzeros+1),
90  proc_config, PAZ_PACK);
91  PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
92  PAZ_broadcast((char *) (*val), sizeof(double)*(*n_nonzeros+1),
93  proc_config, PAZ_PACK);
94  PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
95 
96  printf("Processor %d of %d done with matrix broadcast.\n",
97  proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
98 
99  /* Set rhs and initialize guess */
100  if(proc_config[PAZ_node] != 0)
101  {
102  (*b) = (double *) calloc(*N_global,sizeof(double)) ;
103  (*bt) = (double *) calloc(*N_global,sizeof(double)) ;
104  (*x) = (double *) calloc(*N_global,sizeof(double)) ;
105  if (have_xexact)
106  (*xexact) = (double *) calloc(*N_global,sizeof(double)) ;
107  }
108 
109  PAZ_broadcast((char *) (*x), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
110  PAZ_broadcast((char *) (*b), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
111  PAZ_broadcast((char *) (*bt), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
112  if (have_xexact)
113  PAZ_broadcast((char *)
114  (*xexact), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
115  PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
116  printf("Processor %d of %d done with rhs/guess broadcast.\n",
117  proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
118 
119  }
120 
121  /********************** Generate update map *************************/
122 
123  PAZ_read_update(N_update, update, proc_config, (*N_global),
124  1, PAZ_linear) ;
125 
126  printf("Processor %d of %d has %d rows of %d total rows.\n",
127  proc_config[PAZ_node],proc_config[PAZ_N_procs],*N_update,(*N_global)) ;
128 
129  /*************** Construct local matrix from global matrix ************/
130 
131  /* The local matrix is a copy of the rows assigned to this processor.
132  It is stored in MSR format and still has global indices (PAZ_transform
133  will complete conversion to local indices.
134  */
135 
136  if(proc_config[PAZ_N_procs] > 1) {
137  n_global_nonzeros = *n_nonzeros;
138 
139  *n_nonzeros = *N_update;
140 
141  for (i=0; i<*N_update; i++)
142  *n_nonzeros += (*bindx)[(*update)[i]+1] - (*bindx)[(*update)[i]];
143 
144  printf("Processor %d of %d has %d nonzeros of %d total nonzeros.\n",
145  proc_config[PAZ_node],proc_config[PAZ_N_procs],
146  *n_nonzeros,n_global_nonzeros) ;
147 
148 #ifdef DEBUG
149  { double sum1 = 0.0;
150  for (i=0;i<(*N_global); i++) sum1 += (*b)[i];
151 
152  printf("Processor %d of %d has sum of b = %12.4g.\n",
153  proc_config[PAZ_node],proc_config[PAZ_N_procs],sum1) ;
154  }
155 #endif /* DEBUG */
156 
157  /* Allocate memory for local matrix */
158 
159  bindx1 = (int *) calloc(*n_nonzeros+1,sizeof(int)) ;
160  val1 = (double *) calloc(*n_nonzeros+1,sizeof(double)) ;
161  b1 = (double *) calloc(*N_update,sizeof(double)) ;
162  bt1 = (double *) calloc(*N_update,sizeof(double)) ;
163  x1 = (double *) calloc(*N_update,sizeof(double)) ;
164  if (have_xexact)
165  xexact1 = (double *) calloc(*N_update,sizeof(double)) ;
166 
167  bindx1[0] = *N_update+1;
168 
169  for (i=0; i<*N_update; i++)
170  {
171  row = (*update)[i];
172  b1[i] = (*b)[row];
173  bt1[i] = (*bt)[row];
174  x1[i] = (*x)[row];
175  if (have_xexact) xexact1[i] = (*xexact)[row];
176  val1[i] = (*val)[row];
177  bindx1[i+1] = bindx1[i];
178 
179 #ifdef DEBUG
180  printf("Proc %d of %d: Global row = %d: Local row = %d:
181  b = %12.4g: x = %12.4g: bindx = %d: val = %12.4g \n",
182  proc_config[PAZ_node],proc_config[PAZ_N_procs],
183  row, i, b1[i], x1[i], bindx1[i], val1[i]) ;
184 #endif
185 
186  for (j = (*bindx)[row]; j < (*bindx)[row+1]; j++)
187  {
188  val1[ bindx1 [i+1] ] = (*val)[j];
189  bindx1[bindx1 [i+1] ] = (*bindx)[j];
190  bindx1[i+1] ++;
191  }
192  }
193 
194  printf("Processor %d of %d done with extracting local operators.\n",
195  proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
196 
197  if (have_xexact)
198  {
199  printf(
200  "The residual using MSR format and exact solution on processor %d is %12.4g\n",
201  proc_config[PAZ_node],
202  smsrres (*N_update, (*N_global), val1, bindx1, xexact1, (*xexact), b1));
203  }
204 
205  /* Release memory for global matrix, rhs and solution */
206 
207  free ((void *) (*val));
208  free ((void *) (*bindx));
209  free ((void *) (*b));
210  free ((void *) (*bt));
211  free ((void *) (*x));
212  if (have_xexact) free((void *) *xexact);
213 
214  /* Return local matrix through same pointers. */
215 
216  *val = val1;
217  *bindx = bindx1;
218  *b = b1;
219  *bt = bt1;
220  *x = x1;
221  if (have_xexact) *xexact = xexact1;
222 
223  }
224  if (have_xexact && proc_config[PAZ_N_procs] == 1)
225  {
226  printf(
227  "The residual using MSR format and exact solution on processor %d is %12.4g\n",
228  proc_config[PAZ_node],
229  smsrres (*N_update, (*N_global), (*val), (*bindx),
230  (*xexact), (*xexact), (*b)));
231  }
232 
233 
234  printf("Processor %d of %d leaving distrib_matrix.\n",
235  proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
236 
237  /* end distrib_matrix */
238 }
double smsrres(int m, int n, double *val, int *indx, double *xlocal, double *x, double *b)
Definition: smsrres.c:47
void distrib_msr_matrix(int *proc_config, int *N_global, int *n_nonzeros, int *N_update, int **update, double **val, int **bindx, double **x, double **b, double **bt, double **xexact)