194 struct ITLIN_IO *itlin_ioctl = NULL;
196 int zibnum_fwalloc(
int size,
double **ptr,
char vname[])
199 fprintf(stdout,
"\n allocation of %s follows, size=%i\n",vname,size);
201 *ptr = (
double *) malloc((
int) size*
sizeof(double)) ;
203 fprintf(stdout,
"\n allocation of %s done\n",vname);
206 {
if (ERRORLEVEL>0 && ERROR)
207 fprintf(ERROR,
"\n allocation of %s failed!\n",vname);
210 {
for(i=0;i<size;i++) (*ptr)[i]=0.0;
return 0;};
213 int zibnum_iwalloc(
int size,
int **ptr,
char vname[])
216 fprintf(stdout,
"\n allocation of %s follows, size=%i\n",vname,size);
218 *ptr = (
int *) malloc((
int) size*
sizeof(int)) ;
220 fprintf(stdout,
"\n allocation of %s done\n",vname);
223 {
if (ERRORLEVEL>0 && ERROR)
224 fprintf(ERROR,
"\n allocation of %s failed!\n",vname);
227 {
for(i=0;i<size;i++) (*ptr)[i]=0;
return 0;};
230 int zibnum_pfwalloc(
int size,
double ***ptr,
char vname[])
233 fprintf(stdout,
"\n allocation of %s follows, size=%i\n",vname,size);
235 *ptr = (
double **) malloc((
int) size*
sizeof(
double*)) ;
237 fprintf(stdout,
"\n allocation of %s done\n",vname);
240 {
if (ERRORLEVEL>0 && ERROR)
241 fprintf(ERROR,
"\n allocation of %s failed!\n",vname);
244 {
for(i=0;i<size;i++) (*ptr)[i]=NULL;
return 0;};
247 double zibnum_scaled_norm2(
int n,
double *v,
double *scale)
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 );
254 double zibnum_scaled_sprod(
int n,
double *v1,
double *v2,
double *scale)
256 double t1, t2, rval = 0.0;
258 {t1=v1[i]/scale[i]; t2=v2[i]/scale[i]; rval += t1*t2;};
259 return rval / (double)n ;
262 double zibnum_norm2(
int n,
double *v)
265 for (i=0;i<n;i++) rval += v[i]*v[i];
266 return sqrt( rval / (
double)n );
269 void zibnum_scale(
int n,
double *v1,
double *v2,
double *scale)
271 for (i=0;i<n;i++) v2[i]=v1[i]/scale[i];
275 void zibnum_descale(
int n,
double *v1,
double *v2,
double *scale)
277 for (i=0;i<n;i++) v2[i]=v1[i]*scale[i];
281 void itlin_noprecon(
int n,
double *x,
double *z)
283 for (j=0;j<n;j++) z[j]=x[j];
return;
286 void itlin_dataout(
int k,
int n,
double *x,
struct ITLIN_DATA *data)
288 double *res = data->res;
290 { fprintf(FITER,
"%5i",k);
291 for (i=0;i<n;i++) fprintf(FITER,
" % 14.10e",x[i]);
295 { fprintf(FRES,
"%5i",k);
296 for (i=0;i<n;i++) fprintf(FRES,
" % 14.10e",res[i]);
300 fprintf(FMISC,
"%5i %1i % 14.10e % 14.10e % 14.10e % 14.10e\n",
301 k,data->codeid, data->residuum, data->normdx, data->tau,
303 if ( DATALEVEL==0 )
return;
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");
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);
318 fprintf(DATA,
" % 14.10e % 14.10e % 14.10e\n",
321 { fprintf(DATA,
" % 14.10e",x[i]); i++;
322 if (i<n) fprintf(DATA,
" % 14.10e",x[i]);
325 fprintf(DATA,
" % 14.10e % 14.10e\n",data->residuum,data->normdx);
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
334 double large = 1.0/SMALL, default_scale;
336 {
if ( ERRORLEVEL>0 )
337 fprintf(ERROR,
"\n Error - Routine matvec must be supplied\n");
341 {
if ( ERRORLEVEL>0 )
342 fprintf(ERROR,
"\n Error - Number of Equations/unknowns must be >0\n");
346 {
if ( ERRORLEVEL>0 )
347 fprintf(ERROR,
"\n Error - opt->tol must be positive\n");
351 {
if ( opt->tol > TOLMAX )
355 "\n User prescribed RTOL decreased to reasonable largest value RTOL=%e\n",
358 else if ( opt->tol < TOLMIN )
362 "\n User prescribed RTOL increased to reasonable smallest value RTOL=%e\n",
366 default_scale = opt->tol;
367 if ( opt->scale && itlin_code == 1 )
369 {
if ( (opt->scale)[i] < 0.0 )
370 {
if ( ERRORLEVEL>0 )
372 "\n Error - negative value (opt->scale)[%i] supplied\n",i);
375 else if ( (opt->scale)[i] == 0.0 ) (opt->scale)[i] = default_scale;
376 else if ( (opt->scale)[i] < SMALL )
377 {
if ( ERRORLEVEL>1 )
379 "\n Warning: (opt->scale)[%i] too small - increased to %e\n",
381 (opt->scale)[i] = SMALL;
383 else if ( (opt->scale)[i] > large )
384 {
if ( ERRORLEVEL>1 )
386 "\n Warning: (opt->scale)[%i] too large - decreased to %e\n",
388 (opt->scale)[i] = large;
392 if ( itlin_code == 0 )
394 if ( opt->i_max <= 0 )
395 { opt->i_max = MIN(10,n);
396 if ( ERRORLEVEL > 1 )
398 " Warning: opt->i_max not set; is reset to i_max=%i\n",opt->i_max);
400 else if ( opt->i_max > n )
401 { opt->i_max = MIN(10,n);
402 if ( ERRORLEVEL > 1 )
404 " Warning: opt->i_max is greater than n; is reset to i_max=%i\n",opt->i_max);
407 else if ( itlin_code == 1 )
409 if ( opt->i_max <= 0 )
411 if ( ERRORLEVEL > 1 )
413 " Warning: opt->i_max not set; is reset to i_max=%i\n",opt->i_max);
415 if ( opt->rho <= 0.0 )
417 if ( ERRORLEVEL > 1 )
419 " Warning: opt->rho not set; is reset to rho=%e\n",opt->rho);
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",
430 if ( itlin_code == GMRES || itlin_code == GBIT )
431 fprintf(MONITOR,
" The maximum number of iterations before restart is: %i\n",
433 if ( itlin_code == 1 )
435 fprintf(MONITOR,
" The safety factor is rho = %6e\n",opt->rho);
437 else if ( itlin_code == 2 )
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");
void matvec(int n, double *x, double *b)
Used for gmres.