Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
48void 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}
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)
double smsrres(int m, int n, double *val, int *indx, double *xlocal, double *x, double *b)
Definition: smsrres.c:47