VERB_code_2.3
utils.c
1 /*
2  Common utility routines for the NewtonLib package.
3 
4  * Written by L. Weimann
5  * Purpose Performing certain common tasks of NewtonLib codes,
6  part itlin.
7  * Category ???. - Utilities
8  * Keywords Memory allocation, scaled norm, scaled scalarproduct,
9  data output, parameter check.
10  * Version 1.0
11  * Revision May 2006
12  * Latest Change May 2006
13  * Library NewtonLib
14  * Code C, Double Precision
15  * Environment Standard C environment on PC's,
16  workstations and hosts.
17  * Copyright (c) Konrad-Zuse-Zentrum fuer
18  Informationstechnik Berlin (ZIB)
19  Takustrasse 7, D-14195 Berlin-Dahlem
20  phone : + 49/30/84185-0
21  fax : + 49/30/84185-125
22  * Contact Lutz Weimann
23  ZIB, Division Scientific Computing,
24  Department Numerical Analysis and Modelling
25  phone : + 49/30/84185-185
26  fax : + 49/30/84185-107
27  e-mail: weimann@zib.de
28 
29  ---------------------------------------------------------------
30 
31  * Licence
32  You may use or modify this code for your own non commercial
33  purposes for an unlimited time.
34  In any case you should not deliver this code without a special
35  permission of ZIB.
36  In case you intend to use the code commercially, we oblige you
37  to sign an according licence agreement with ZIB.
38 
39  * Warranty
40  This code has been tested up to a certain level. Defects and
41  weaknesses, which may be included in the code, do not establish
42  any warranties by ZIB. ZIB does not take over any liabilities
43  which may follow from acquisition or application of this code.
44 
45  * Software status
46  This code is under care of ZIB and belongs to ZIB software class 2.
47 
48  ------------------------------------------------------------
49 
50 
51  This file contains the following routines and functions:
52  --------------------------------------------------------
53 
54  zibnum_fwalloc - allocate memory for a double precision array
55  zibnum_iwalloc - allocate memory for an integer array
56  zibnum_pfwalloc - allocate memory for an array of pointers to
57  double precision arrays
58  zibnum_scaled_norm2 - compute the scaled Euclidian-norm of a vector
59  zibnum_scaled_sprod - compute a scaled scalar product
60  zibnum_norm2 - compute the unscaled Euclidian-norm of a vector
61  zibnum_scale - compute a scaled vector from an unscaled vector
62  zibnum_descale - compute the unscaled vector from a scaled vector
63  zibnum_dataout - write data output
64  zibnum_parcheck_and_print - check and print parameter and options settings
65 
66  The calls of the routines/functions are as follows:
67  ---------------------------------------------------
68 
69  int zibnum_fwalloc(int size, double **ptr, char vname[])
70 
71  This function allocates memory for a double array via the malloc function.
72  The parameters are:
73  int size (input) The number of elements of type double to be allocated.
74  double **ptr (output) The pointer to memory, of type (double *), which has
75  been returned by the malloc function.
76  char vname[] (input) An identifying character string of the memory portion
77  to be allocated, used for print within a possible error
78  message.
79  The int return code is 0, if the memory allocation was successful, or -999
80  if the allocation failed.
81  ---
82  int zibnum_iwalloc(int size, int **ptr, char vname[])
83 
84  This function allocates memory for a int array via the malloc function.
85  The parameters are:
86  int size (input) The number of elements of type int to be allocated.
87  int **ptr (output) The pointer to memory, of type (int *), which has
88  been returned by the malloc function.
89  char vname[] (input) An identifying character string of the memory portion
90  to be allocated, used for print within a possible error
91  message.
92  The int return code is 0, if the memory allocation was successful, or -998
93  if the allocation failed.
94  ---
95  int zibnum_pfwalloc(int size, double ***ptr, char vname[])
96 
97  This function allocates memory for a pointer array via the malloc function.
98  The parameters are:
99  int size (input) The number of elements of type pointer to be allocated.
100  int ***ptr (output) The pointer to memory, of type (double **) (i.e.
101  to pointers which are pointing to memory allocated
102  for double arrays), which has been returned by the
103  malloc function.
104  char vname[] (input) An identifying character string of the memory portion
105  to be allocated, used for print within a possible error
106  message.
107  The int return code is 0, if the memory allocation was successful, or -997
108  if the allocation failed.
109  ---
110  Note: In order to activate some debug output on dynamic memory allocation
111  set the C preprocessor flag DEBUG, i.e. set the option -DDEBUG if
112  you use the GNU C-compiler (gcc).
113  ---
114  double zibnum_scaled_norm2(int n, double *v, double *scale)
115 
116  This function computes the scaled norm of the (double) vector v of
117  size n, using the (double) vector scale for scaling, as below:
118  result := Sqrt ( ( Sum(i=0 to n-1) (v[i]/scale[i])^2 ) / n ) .
119  ---
120  double zibnum_scaled_sprod(int n, double *v1, double *v2, double *scale)
121 
122  This function computes the scaled scalar product of the (double) vectors
123  v1 and v2 of size n, using the (double) vector scale for scaling, as below:
124  result := ( Sum(i=0 to n-1) (v1[i]/scale[i])*(v2[i]/scale[i]) ) / n .
125  ---
126  double zibnum_norm2(int n, double *v)
127 
128  This function computes the (ordinary) norm of the (double) vector v of
129  size n, as below:
130  result := Sqrt ( ( Sum(i=0 to n-1) v[i]^2 ) / n ) .
131  ---
132  void zibnum_scale(int n, double *v1, double *v2, double *scale)
133 
134  This routine computes the scaled vector of the vector v1 of size n and
135  stores the result to the vector v2, as below:
136  v2[i] = v1[i]/scale[i] for i=0,...,n-1 .
137  ---
138  void zibnum_descale(int n, double *v1, double *v2, double *scale)
139 
140  This routine computes the descaled vector of the vector v1 of size n and
141  stores the result to the vector v2, as below:
142  v2[i] = v1[i]*scale[i] for i=0,...,n-1 .
143  ---
144  void itlin_dataout(int k, int n, double *x, struct ITLIN_DATA *data)
145 
146  This routine is designed to print out computed data within each iteration
147  step. It may be replaced or extended for special purposes.
148  The parameters are (all input arguments):
149 
150  int k The current iteration step number.
151  int n The size of any double arrays mentioned below.
152  double *x An array which holds the current iterate x^k .
153 
154  The fields of the struct ITLIN_DATA are:
155 
156  double* data->res An array which holds the residuum r=b-A*x^k.
157  double data->residuum
158  double data->normdx The norm of the latest solution correction.
159  Not used by GMRES.
160  double data->tau The GBIT parameter tau or the gamma_sum of PCG
161  which is used in the PCG termination criterium.
162  Not used by GMRES.
163  double data->t The GBIT parameter t or the quantity gamma_i^2 of PCG.
164  Not used by GMRES.
165  enum data->mode The mode of the current iterate:
166  Initial (=1): This is the first call of itlin_dataout.
167  Intermediate (=2): This is an intermediate call of
168  itlin_dataout.
169  Solution (=3): This is a final call of itlin_dataout,
170  and *x holds an approximate solution.
171  Final (=4) : This is a final call of itlin_dataout,
172  but *x does not hold a solution.
173  ---
174  int itlin_parcheck_and_print(int n, MATVEC *matvec,
175  struct ITLIN_OPT *opt, int itlin_code)
176 
177  This function checks and prints parameter and options settings.
178  The parameters are:
179 
180  int n The number of equations and unknowns.
181  struct ITLIN_OPT *opt A pointer to an options data structure.
182  int itlin_code The identification number of the calling code.
183  Valid identifications are:
184  0 : GMRES
185  1 : GBIT
186  2 : PCG
187 
188 */
189 #include <stdlib.h>
190 #include <stdio.h>
191 #include <math.h>
192 #include "itlin.h"
193 
194 struct ITLIN_IO *itlin_ioctl = NULL;
195 
196 int zibnum_fwalloc(int size, double **ptr, char vname[])
197 { int i;
198 #ifdef DEBUG
199  fprintf(stdout,"\n allocation of %s follows, size=%i\n",vname,size);
200 #endif
201  *ptr = (double *) malloc((int) size*sizeof(double)) ;
202 #ifdef DEBUG
203  fprintf(stdout,"\n allocation of %s done\n",vname);
204 #endif
205  if (*ptr == NULL)
206  { if (ERRORLEVEL>0 && ERROR)
207  fprintf(ERROR,"\n allocation of %s failed!\n",vname);
208  return -999;}
209  else
210  {for(i=0;i<size;i++) (*ptr)[i]=0.0; return 0;};
211 }
212 
213 int zibnum_iwalloc(int size, int **ptr, char vname[])
214 { int i;
215 #ifdef DEBUG
216  fprintf(stdout,"\n allocation of %s follows, size=%i\n",vname,size);
217 #endif
218  *ptr = (int *) malloc((int) size*sizeof(int)) ;
219 #ifdef DEBUG
220  fprintf(stdout,"\n allocation of %s done\n",vname);
221 #endif
222  if (*ptr == NULL)
223  { if (ERRORLEVEL>0 && ERROR)
224  fprintf(ERROR,"\n allocation of %s failed!\n",vname);
225  return -998;}
226  else
227  {for(i=0;i<size;i++) (*ptr)[i]=0; return 0;};
228 }
229 
230 int zibnum_pfwalloc(int size, double ***ptr, char vname[])
231 { int i;
232 #ifdef DEBUG
233  fprintf(stdout,"\n allocation of %s follows, size=%i\n",vname,size);
234 #endif
235  *ptr = (double **) malloc((int) size*sizeof(double*)) ;
236 #ifdef DEBUG
237  fprintf(stdout,"\n allocation of %s done\n",vname);
238 #endif
239  if (*ptr == NULL)
240  { if (ERRORLEVEL>0 && ERROR)
241  fprintf(ERROR,"\n allocation of %s failed!\n",vname);
242  return -997;}
243  else
244  {for(i=0;i<size;i++) (*ptr)[i]=NULL; return 0;};
245 }
246 
247 double zibnum_scaled_norm2(int n, double *v, double *scale)
248 { int i;
249  double t, rval = 0.0;
250  for (i=0;i<n;i++) {t=v[i]/scale[i]; rval += t*t;};
251  return sqrt( rval / (double)n );
252 }
253 
254 double zibnum_scaled_sprod(int n, double *v1, double *v2, double *scale)
255 { int i;
256  double t1, t2, rval = 0.0;
257  for (i=0;i<n;i++)
258  {t1=v1[i]/scale[i]; t2=v2[i]/scale[i]; rval += t1*t2;};
259  return rval / (double)n ;
260 }
261 
262 double zibnum_norm2(int n, double *v)
263 { int i;
264  double rval = 0.0;
265  for (i=0;i<n;i++) rval += v[i]*v[i];
266  return sqrt( rval / (double)n );
267 }
268 
269 void zibnum_scale(int n, double *v1, double *v2, double *scale)
270 { int i;
271  for (i=0;i<n;i++) v2[i]=v1[i]/scale[i];
272  return ;
273 }
274 
275 void zibnum_descale(int n, double *v1, double *v2, double *scale)
276 { int i;
277  for (i=0;i<n;i++) v2[i]=v1[i]*scale[i];
278  return ;
279 }
280 
281 void itlin_noprecon(int n, double *x, double *z)
282 { int j;
283  for (j=0;j<n;j++) z[j]=x[j]; return;
284 }
285 
286 void itlin_dataout(int k, int n, double *x, struct ITLIN_DATA *data)
287 { int i;
288  double *res = data->res;
289  if (FITER)
290  { fprintf(FITER,"%5i",k);
291  for (i=0;i<n;i++) fprintf(FITER," % 14.10e",x[i]);
292  fprintf(FITER,"\n");
293  };
294  if (FRES)
295  { fprintf(FRES,"%5i",k);
296  for (i=0;i<n;i++) fprintf(FRES," % 14.10e",res[i]);
297  fprintf(FRES,"\n");
298  };
299  if (FMISC)
300  fprintf(FMISC,"%5i %1i % 14.10e % 14.10e % 14.10e % 14.10e\n",
301  k,data->codeid, data->residuum, data->normdx, data->tau,
302  data->t);
303  if ( DATALEVEL==0 ) return;
304  if ( k==0 )
305  { fprintf(DATA," Start data:\n N = %i\n\n",n);
306  fprintf(DATA," Format: iteration-number, (x(i),i=1,...N), Residuum , Normdx\n\n");
307  fprintf(DATA," Initial data:\n\n");
308  }
309  else if ( data->mode==Solution )
310  fprintf(DATA,"\n Solution data:\n\n");
311  else if ( data->mode==Final )
312  fprintf(DATA,"\n Final data:\n\n");
313  else if ( k==1 && DATALEVEL>1 )
314  fprintf(DATA,"\n Intermediate data:\n\n");
315  if ( k==0 || data->mode==Solution || data->mode==Final || DATALEVEL > 1 )
316  { fprintf(DATA," %4i\n",k);
317  for (i=0;i<n-2;i+=3)
318  fprintf(DATA," % 14.10e % 14.10e % 14.10e\n",
319  x[i],x[i+1],x[i+2]);
320  if (i<n)
321  { fprintf(DATA," % 14.10e",x[i]); i++;
322  if (i<n) fprintf(DATA," % 14.10e",x[i]);
323  fprintf(DATA,"\n");
324  };
325  fprintf(DATA," % 14.10e % 14.10e\n",data->residuum,data->normdx);
326  };
327 }
328 
329 int itlin_parcheck_and_print(int n, MATVEC *matvec,
330  struct ITLIN_OPT *opt, int itlin_code)
331 #define TOLMIN 1.0e-15
332 #define TOLMAX 1.0e-1
333 { int i, fail=0;
334  double large = 1.0/SMALL, default_scale;
335  if (!matvec)
336  { if ( ERRORLEVEL>0 )
337  fprintf(ERROR,"\n Error - Routine matvec must be supplied\n");
338  return -99;
339  };
340  if ( n<=0 )
341  { if ( ERRORLEVEL>0 )
342  fprintf(ERROR,"\n Error - Number of Equations/unknowns must be >0\n");
343  return 20;
344  };
345  if ( opt->tol <= 0 )
346  { if ( ERRORLEVEL>0 )
347  fprintf(ERROR,"\n Error - opt->tol must be positive\n");
348  return 21;
349  }
350  else
351  { if ( opt->tol > TOLMAX )
352  { opt->tol = TOLMAX;
353  if ( ERRORLEVEL>1 )
354  fprintf(ERROR,
355  "\n User prescribed RTOL decreased to reasonable largest value RTOL=%e\n",
356  opt->tol);
357  }
358  else if ( opt->tol < TOLMIN )
359  { opt->tol = TOLMIN;
360  if ( ERRORLEVEL>1 )
361  fprintf(ERROR,
362  "\n User prescribed RTOL increased to reasonable smallest value RTOL=%e\n",
363  opt->tol);
364  };
365  };
366  default_scale = opt->tol;
367  if ( opt->scale && itlin_code == 1 )
368  { for (i=0;i<n;i++)
369  { if ( (opt->scale)[i] < 0.0 )
370  { if ( ERRORLEVEL>0 )
371  fprintf(ERROR,
372  "\n Error - negative value (opt->scale)[%i] supplied\n",i);
373  return 22;
374  }
375  else if ( (opt->scale)[i] == 0.0 ) (opt->scale)[i] = default_scale;
376  else if ( (opt->scale)[i] < SMALL )
377  { if ( ERRORLEVEL>1 )
378  fprintf(ERROR,
379  "\n Warning: (opt->scale)[%i] too small - increased to %e\n",
380  i,SMALL);
381  (opt->scale)[i] = SMALL;
382  }
383  else if ( (opt->scale)[i] > large )
384  { if ( ERRORLEVEL>1 )
385  fprintf(ERROR,
386  "\n Warning: (opt->scale)[%i] too large - decreased to %e\n",
387  i,large);
388  (opt->scale)[i] = large;
389  };
390  };
391  };
392  if ( itlin_code == 0 )
393  {
394  if ( opt->i_max <= 0 )
395  { opt->i_max = MIN(10,n);
396  if ( ERRORLEVEL > 1 )
397  fprintf(ERROR,
398  " Warning: opt->i_max not set; is reset to i_max=%i\n",opt->i_max);
399  }
400  else if ( opt->i_max > n )
401  { opt->i_max = MIN(10,n);
402  if ( ERRORLEVEL > 1 )
403  fprintf(ERROR,
404  " Warning: opt->i_max is greater than n; is reset to i_max=%i\n",opt->i_max);
405  };
406  }
407  else if ( itlin_code == 1 )
408  {
409  if ( opt->i_max <= 0 )
410  { opt->i_max = 10;
411  if ( ERRORLEVEL > 1 )
412  fprintf(ERROR,
413  " Warning: opt->i_max not set; is reset to i_max=%i\n",opt->i_max);
414  };
415  if ( opt->rho <= 0.0 )
416  { opt->rho = 1.0;
417  if ( ERRORLEVEL > 1 )
418  fprintf(ERROR,
419  " Warning: opt->rho not set; is reset to rho=%e\n",opt->rho);
420  };
421  };
422  if ( MONITORLEVEL==0 ) return 0;
423  fprintf(MONITOR,"\n Problem dimension: n = %i\n",n);
424  if ( itlin_code == 0 )
425  fprintf(MONITOR,"\n Prescribed residuum threshold: %e\n",opt->tol);
426  else if ( itlin_code == 1 )
427  fprintf(MONITOR,"\n Prescribed relative precision: %e\n",opt->tol);
428  fprintf(MONITOR," The maximum permitted number of iteration steps is: %i\n",
429  opt->maxiter);
430  if ( itlin_code == GMRES || itlin_code == GBIT )
431  fprintf(MONITOR," The maximum number of iterations before restart is: %i\n",
432  opt->i_max);
433  if ( itlin_code == 1 )
434  {
435  fprintf(MONITOR," The safety factor is rho = %6e\n",opt->rho);
436  }
437  else if ( itlin_code == 2 )
438  {
439  if ( opt->convcheck == Absolute )
440  fprintf(MONITOR," The absolute error will be checked for termination\n");
441  else if ( opt->convcheck == Relative )
442  fprintf(MONITOR," The relative error will be checked for termination\n");
443  };
444  return 0;
445 }
void matvec(int n, double *x, double *b)
Used for gmres.