IFPACK Development
Loading...
Searching...
No Matches
Vec_dh.c
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 "Vec_dh.h"
45#include "Mem_dh.h"
46#include "SubdomainGraph_dh.h"
47#include "io_dh.h"
48
49#undef __FUNC__
50#define __FUNC__ "Vec_dhCreate"
51void
52Vec_dhCreate (Vec_dh * v)
53{
54 START_FUNC_DH
55 struct _vec_dh *tmp =
56 (struct _vec_dh *) MALLOC_DH (sizeof (struct _vec_dh));
57 CHECK_V_ERROR;
58 *v = tmp;
59 tmp->n = 0;
60 tmp->vals = NULL;
61END_FUNC_DH}
62
63#undef __FUNC__
64#define __FUNC__ "Vec_dhDestroy"
65void
66Vec_dhDestroy (Vec_dh v)
67{
68 START_FUNC_DH if (v->vals != NULL)
69 FREE_DH (v->vals);
70 CHECK_V_ERROR;
71 FREE_DH (v);
72 CHECK_V_ERROR;
73END_FUNC_DH}
74
75
76#undef __FUNC__
77#define __FUNC__ "Vec_dhInit"
78void
79Vec_dhInit (Vec_dh v, int size)
80{
81 START_FUNC_DH v->n = size;
82 v->vals = (double *) MALLOC_DH (size * sizeof (double));
83 CHECK_V_ERROR;
84END_FUNC_DH}
85
86#undef __FUNC__
87#define __FUNC__ "Vec_dhCopy"
88void
89Vec_dhCopy (Vec_dh x, Vec_dh y)
90{
91 START_FUNC_DH if (x->vals == NULL)
92 SET_V_ERROR ("x->vals is NULL");
93 if (y->vals == NULL)
94 SET_V_ERROR ("y->vals is NULL");
95 if (x->n != y->n)
96 SET_V_ERROR ("x and y are different lengths");
97 memcpy (y->vals, x->vals, x->n * sizeof (double));
98END_FUNC_DH}
99
100
101#undef __FUNC__
102#define __FUNC__ "Vec_dhDuplicate"
103void
104Vec_dhDuplicate (Vec_dh v, Vec_dh * out)
105{
106 START_FUNC_DH Vec_dh tmp;
107 int size = v->n;
108 if (v->vals == NULL)
109 SET_V_ERROR ("v->vals is NULL");
110 Vec_dhCreate (out);
111 CHECK_V_ERROR;
112 tmp = *out;
113 tmp->n = size;
114 tmp->vals = (double *) MALLOC_DH (size * sizeof (double));
115 CHECK_V_ERROR;
116END_FUNC_DH}
117
118#undef __FUNC__
119#define __FUNC__ "Vec_dhSet"
120void
121Vec_dhSet (Vec_dh v, double value)
122{
123 START_FUNC_DH int i, m = v->n;
124 double *vals = v->vals;
125 if (v->vals == NULL)
126 SET_V_ERROR ("v->vals is NULL");
127 for (i = 0; i < m; ++i)
128 vals[i] = value;
129END_FUNC_DH}
130
131#undef __FUNC__
132#define __FUNC__ "Vec_dhSetRand"
133void
134Vec_dhSetRand (Vec_dh v)
135{
136 START_FUNC_DH int i, m = v->n;
137 double max = 0.0;
138 double *vals = v->vals;
139
140 if (v->vals == NULL)
141 SET_V_ERROR ("v->vals is NULL");
142
143#ifdef WIN32
144 for (i = 0; i < m; ++i)
145 vals[i] = rand ();
146#else
147 for (i = 0; i < m; ++i)
148 vals[i] = rand ();
149#endif
150
151 /* find largest value in vector, and scale vector,
152 * so all values are in [0.0,1.0]
153 */
154 for (i = 0; i < m; ++i)
155 max = MAX (max, vals[i]);
156 for (i = 0; i < m; ++i)
157 vals[i] = vals[i] / max;
158END_FUNC_DH}
159
160
161#undef __FUNC__
162#define __FUNC__ "Vec_dhPrint"
163void
164Vec_dhPrint (Vec_dh v, SubdomainGraph_dh sg, char *filename)
165{
166 START_FUNC_DH double *vals = v->vals;
167 int pe, i, m = v->n;
168 FILE *fp;
169
170 if (v->vals == NULL)
171 SET_V_ERROR ("v->vals is NULL");
172
173 /*--------------------------------------------------------
174 * case 1: no permutation information
175 *--------------------------------------------------------*/
176 if (sg == NULL)
177 {
178 for (pe = 0; pe < np_dh; ++pe)
179 {
180 MPI_Barrier (comm_dh);
181 if (pe == myid_dh)
182 {
183 if (pe == 0)
184 {
185 fp = openFile_dh (filename, "w");
186 CHECK_V_ERROR;
187 }
188 else
189 {
190 fp = openFile_dh (filename, "a");
191 CHECK_V_ERROR;
192 }
193
194 for (i = 0; i < m; ++i)
195 fprintf (fp, "%g\n", vals[i]);
196
197 closeFile_dh (fp);
198 CHECK_V_ERROR;
199 }
200 }
201 }
202
203 /*--------------------------------------------------------
204 * case 2: single mpi task, multiple subdomains
205 *--------------------------------------------------------*/
206 else if (np_dh == 1)
207 {
208 int i, j;
209
210 fp = openFile_dh (filename, "w");
211 CHECK_V_ERROR;
212
213 for (i = 0; i < sg->blocks; ++i)
214 {
215 int oldBlock = sg->n2o_sub[i];
216 int beg_row = sg->beg_rowP[oldBlock];
217 int end_row = beg_row + sg->row_count[oldBlock];
218
219 printf ("seq: block= %i beg= %i end= %i\n", oldBlock, beg_row,
220 end_row);
221
222
223 for (j = beg_row; j < end_row; ++j)
224 {
225 fprintf (fp, "%g\n", vals[j]);
226 }
227 }
228 }
229
230 /*--------------------------------------------------------
231 * case 3: multiple mpi tasks, one subdomain per task
232 *--------------------------------------------------------*/
233 else
234 {
235 int id = sg->o2n_sub[myid_dh];
236 for (pe = 0; pe < np_dh; ++pe)
237 {
238 MPI_Barrier (comm_dh);
239 if (id == pe)
240 {
241 if (pe == 0)
242 {
243 fp = openFile_dh (filename, "w");
244 CHECK_V_ERROR;
245 }
246 else
247 {
248 fp = openFile_dh (filename, "a");
249 CHECK_V_ERROR;
250 }
251
252 fprintf (stderr, "par: block= %i\n", id);
253
254 for (i = 0; i < m; ++i)
255 {
256 fprintf (fp, "%g\n", vals[i]);
257 }
258
259 closeFile_dh (fp);
260 CHECK_V_ERROR;
261 }
262 }
263 }
264END_FUNC_DH}
265
266
267#undef __FUNC__
268#define __FUNC__ "Vec_dhPrintBIN"
269void
270Vec_dhPrintBIN (Vec_dh v, SubdomainGraph_dh sg, char *filename)
271{
272 START_FUNC_DH if (np_dh > 1)
273 {
274 SET_V_ERROR ("only implemented for a single MPI task");
275 }
276 if (sg != NULL)
277 {
278 SET_V_ERROR ("not implemented for reordered vector; ensure sg=NULL");
279 }
280
281 io_dh_print_ebin_vec_private (v->n, 0, v->vals, NULL, NULL, NULL, filename);
282 CHECK_V_ERROR;
283END_FUNC_DH}
284
285#define MAX_JUNK 200
286
287#undef __FUNC__
288#define __FUNC__ "Vec_dhRead"
289void
290Vec_dhRead (Vec_dh * vout, int ignore, char *filename)
291{
292 START_FUNC_DH Vec_dh tmp;
293 FILE *fp;
294 int items, n, i;
295 double *v, w;
296 char junk[MAX_JUNK];
297
298 Vec_dhCreate (&tmp);
299 CHECK_V_ERROR;
300 *vout = tmp;
301
302 if (np_dh > 1)
303 {
304 SET_V_ERROR ("only implemented for a single MPI task");
305 }
306
307 fp = openFile_dh (filename, "w");
308 CHECK_V_ERROR;
309
310 /* skip over file lines */
311 if (ignore)
312 {
313 printf ("Vec_dhRead:: ignoring following header lines:\n");
314 printf
315 ("--------------------------------------------------------------\n");
316 for (i = 0; i < ignore; ++i)
317 {
318 fgets (junk, MAX_JUNK, fp);
319 printf ("%s", junk);
320 }
321 printf
322 ("--------------------------------------------------------------\n");
323 }
324
325 /* count floating point entries in file */
326 n = 0;
327 while (!feof (fp))
328 {
329 items = fscanf (fp, "%lg", &w);
330 if (items != 1)
331 {
332 break;
333 }
334 ++n;
335 }
336
337 printf ("Vec_dhRead:: n= %i\n", n);
338
339 /* allocate storage */
340 tmp->n = n;
341 v = tmp->vals = (double *) MALLOC_DH (n * sizeof (double));
342 CHECK_V_ERROR;
343
344 /* reset file, and skip over header again */
345 rewind (fp);
346 rewind (fp);
347 for (i = 0; i < ignore; ++i)
348 {
349 fgets (junk, MAX_JUNK, fp);
350 }
351
352 /* read values */
353 for (i = 0; i < n; ++i)
354 {
355 items = fscanf (fp, "%lg", v + i);
356 if (items != 1)
357 {
358 sprintf (msgBuf_dh, "failed to read value %i of %i", i + 1, n);
359 }
360 }
361
362 closeFile_dh (fp);
363 CHECK_V_ERROR;
364END_FUNC_DH}
365
366#undef __FUNC__
367#define __FUNC__ "Vec_dhReadBIN"
368extern void
369Vec_dhReadBIN (Vec_dh * vout, char *filename)
370{
371 START_FUNC_DH Vec_dh tmp;
372
373 Vec_dhCreate (&tmp);
374 CHECK_V_ERROR;
375 *vout = tmp;
376 io_dh_read_ebin_vec_private (&tmp->n, &tmp->vals, filename);
377 CHECK_V_ERROR;
378END_FUNC_DH}
Definition: Vec_dh.h:53