IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
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"
51 void
52 Vec_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;
61 END_FUNC_DH}
62 
63 #undef __FUNC__
64 #define __FUNC__ "Vec_dhDestroy"
65 void
66 Vec_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;
73 END_FUNC_DH}
74 
75 
76 #undef __FUNC__
77 #define __FUNC__ "Vec_dhInit"
78 void
79 Vec_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;
84 END_FUNC_DH}
85 
86 #undef __FUNC__
87 #define __FUNC__ "Vec_dhCopy"
88 void
89 Vec_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));
98 END_FUNC_DH}
99 
100 
101 #undef __FUNC__
102 #define __FUNC__ "Vec_dhDuplicate"
103 void
104 Vec_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;
116 END_FUNC_DH}
117 
118 #undef __FUNC__
119 #define __FUNC__ "Vec_dhSet"
120 void
121 Vec_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;
129 END_FUNC_DH}
130 
131 #undef __FUNC__
132 #define __FUNC__ "Vec_dhSetRand"
133 void
134 Vec_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;
158 END_FUNC_DH}
159 
160 
161 #undef __FUNC__
162 #define __FUNC__ "Vec_dhPrint"
163 void
164 Vec_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  }
264 END_FUNC_DH}
265 
266 
267 #undef __FUNC__
268 #define __FUNC__ "Vec_dhPrintBIN"
269 void
270 Vec_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;
283 END_FUNC_DH}
284 
285 #define MAX_JUNK 200
286 
287 #undef __FUNC__
288 #define __FUNC__ "Vec_dhRead"
289 void
290 Vec_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;
364 END_FUNC_DH}
365 
366 #undef __FUNC__
367 #define __FUNC__ "Vec_dhReadBIN"
368 extern void
369 Vec_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;
378 END_FUNC_DH}
Definition: Vec_dh.h:52