00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 #include <stdlib.h>
00190 #include <stdio.h>
00191 #include <math.h>
00192 #include "itlin.h"
00193
00194 struct ITLIN_IO *itlin_ioctl = NULL;
00195
00196 int zibnum_fwalloc(int size, double **ptr, char vname[])
00197 { int i;
00198 #ifdef DEBUG
00199 fprintf(stdout,"\n allocation of %s follows, size=%i\n",vname,size);
00200 #endif
00201 *ptr = (double *) malloc((int) size*sizeof(double)) ;
00202 #ifdef DEBUG
00203 fprintf(stdout,"\n allocation of %s done\n",vname);
00204 #endif
00205 if (*ptr == NULL)
00206 { if (ERRORLEVEL>0 && ERROR)
00207 fprintf(ERROR,"\n allocation of %s failed!\n",vname);
00208 return -999;}
00209 else
00210 {for(i=0;i<size;i++) (*ptr)[i]=0.0; return 0;};
00211 }
00212
00213 int zibnum_iwalloc(int size, int **ptr, char vname[])
00214 { int i;
00215 #ifdef DEBUG
00216 fprintf(stdout,"\n allocation of %s follows, size=%i\n",vname,size);
00217 #endif
00218 *ptr = (int *) malloc((int) size*sizeof(int)) ;
00219 #ifdef DEBUG
00220 fprintf(stdout,"\n allocation of %s done\n",vname);
00221 #endif
00222 if (*ptr == NULL)
00223 { if (ERRORLEVEL>0 && ERROR)
00224 fprintf(ERROR,"\n allocation of %s failed!\n",vname);
00225 return -998;}
00226 else
00227 {for(i=0;i<size;i++) (*ptr)[i]=0; return 0;};
00228 }
00229
00230 int zibnum_pfwalloc(int size, double ***ptr, char vname[])
00231 { int i;
00232 #ifdef DEBUG
00233 fprintf(stdout,"\n allocation of %s follows, size=%i\n",vname,size);
00234 #endif
00235 *ptr = (double **) malloc((int) size*sizeof(double*)) ;
00236 #ifdef DEBUG
00237 fprintf(stdout,"\n allocation of %s done\n",vname);
00238 #endif
00239 if (*ptr == NULL)
00240 { if (ERRORLEVEL>0 && ERROR)
00241 fprintf(ERROR,"\n allocation of %s failed!\n",vname);
00242 return -997;}
00243 else
00244 {for(i=0;i<size;i++) (*ptr)[i]=NULL; return 0;};
00245 }
00246
00247 double zibnum_scaled_norm2(int n, double *v, double *scale)
00248 { int i;
00249 double t, rval = 0.0;
00250 for (i=0;i<n;i++) {t=v[i]/scale[i]; rval += t*t;};
00251 return sqrt( rval / (double)n );
00252 }
00253
00254 double zibnum_scaled_sprod(int n, double *v1, double *v2, double *scale)
00255 { int i;
00256 double t1, t2, rval = 0.0;
00257 for (i=0;i<n;i++)
00258 {t1=v1[i]/scale[i]; t2=v2[i]/scale[i]; rval += t1*t2;};
00259 return rval / (double)n ;
00260 }
00261
00262 double zibnum_norm2(int n, double *v)
00263 { int i;
00264 double rval = 0.0;
00265 for (i=0;i<n;i++) rval += v[i]*v[i];
00266 return sqrt( rval / (double)n );
00267 }
00268
00269 void zibnum_scale(int n, double *v1, double *v2, double *scale)
00270 { int i;
00271 for (i=0;i<n;i++) v2[i]=v1[i]/scale[i];
00272 return ;
00273 }
00274
00275 void zibnum_descale(int n, double *v1, double *v2, double *scale)
00276 { int i;
00277 for (i=0;i<n;i++) v2[i]=v1[i]*scale[i];
00278 return ;
00279 }
00280
00281 void itlin_noprecon(int n, double *x, double *z)
00282 { int j;
00283 for (j=0;j<n;j++) z[j]=x[j]; return;
00284 }
00285
00286 void itlin_dataout(int k, int n, double *x, struct ITLIN_DATA *data)
00287 { int i;
00288 double *res = data->res;
00289 if (FITER)
00290 { fprintf(FITER,"%5i",k);
00291 for (i=0;i<n;i++) fprintf(FITER," % 14.10e",x[i]);
00292 fprintf(FITER,"\n");
00293 };
00294 if (FRES)
00295 { fprintf(FRES,"%5i",k);
00296 for (i=0;i<n;i++) fprintf(FRES," % 14.10e",res[i]);
00297 fprintf(FRES,"\n");
00298 };
00299 if (FMISC)
00300 fprintf(FMISC,"%5i %1i % 14.10e % 14.10e % 14.10e % 14.10e\n",
00301 k,data->codeid, data->residuum, data->normdx, data->tau,
00302 data->t);
00303 if ( DATALEVEL==0 ) return;
00304 if ( k==0 )
00305 { fprintf(DATA," Start data:\n N = %i\n\n",n);
00306 fprintf(DATA," Format: iteration-number, (x(i),i=1,...N), Residuum , Normdx\n\n");
00307 fprintf(DATA," Initial data:\n\n");
00308 }
00309 else if ( data->mode==Solution )
00310 fprintf(DATA,"\n Solution data:\n\n");
00311 else if ( data->mode==Final )
00312 fprintf(DATA,"\n Final data:\n\n");
00313 else if ( k==1 && DATALEVEL>1 )
00314 fprintf(DATA,"\n Intermediate data:\n\n");
00315 if ( k==0 || data->mode==Solution || data->mode==Final || DATALEVEL > 1 )
00316 { fprintf(DATA," %4i\n",k);
00317 for (i=0;i<n-2;i+=3)
00318 fprintf(DATA," % 14.10e % 14.10e % 14.10e\n",
00319 x[i],x[i+1],x[i+2]);
00320 if (i<n)
00321 { fprintf(DATA," % 14.10e",x[i]); i++;
00322 if (i<n) fprintf(DATA," % 14.10e",x[i]);
00323 fprintf(DATA,"\n");
00324 };
00325 fprintf(DATA," % 14.10e % 14.10e\n",data->residuum,data->normdx);
00326 };
00327 }
00328
00329 int itlin_parcheck_and_print(int n, MATVEC *matvec,
00330 struct ITLIN_OPT *opt, int itlin_code)
00331 #define TOLMIN 1.0e-15
00332 #define TOLMAX 1.0e-1
00333 { int i, fail=0;
00334 double large = 1.0/SMALL, default_scale;
00335 if (!matvec)
00336 { if ( ERRORLEVEL>0 )
00337 fprintf(ERROR,"\n Error - Routine matvec must be supplied\n");
00338 return -99;
00339 };
00340 if ( n<=0 )
00341 { if ( ERRORLEVEL>0 )
00342 fprintf(ERROR,"\n Error - Number of Equations/unknowns must be >0\n");
00343 return 20;
00344 };
00345 if ( opt->tol <= 0 )
00346 { if ( ERRORLEVEL>0 )
00347 fprintf(ERROR,"\n Error - opt->tol must be positive\n");
00348 return 21;
00349 }
00350 else
00351 { if ( opt->tol > TOLMAX )
00352 { opt->tol = TOLMAX;
00353 if ( ERRORLEVEL>1 )
00354 fprintf(ERROR,
00355 "\n User prescribed RTOL decreased to reasonable largest value RTOL=%e\n",
00356 opt->tol);
00357 }
00358 else if ( opt->tol < TOLMIN )
00359 { opt->tol = TOLMIN;
00360 if ( ERRORLEVEL>1 )
00361 fprintf(ERROR,
00362 "\n User prescribed RTOL increased to reasonable smallest value RTOL=%e\n",
00363 opt->tol);
00364 };
00365 };
00366 default_scale = opt->tol;
00367 if ( opt->scale && itlin_code == 1 )
00368 { for (i=0;i<n;i++)
00369 { if ( (opt->scale)[i] < 0.0 )
00370 { if ( ERRORLEVEL>0 )
00371 fprintf(ERROR,
00372 "\n Error - negative value (opt->scale)[%i] supplied\n",i);
00373 return 22;
00374 }
00375 else if ( (opt->scale)[i] == 0.0 ) (opt->scale)[i] = default_scale;
00376 else if ( (opt->scale)[i] < SMALL )
00377 { if ( ERRORLEVEL>1 )
00378 fprintf(ERROR,
00379 "\n Warning: (opt->scale)[%i] too small - increased to %e\n",
00380 i,SMALL);
00381 (opt->scale)[i] = SMALL;
00382 }
00383 else if ( (opt->scale)[i] > large )
00384 { if ( ERRORLEVEL>1 )
00385 fprintf(ERROR,
00386 "\n Warning: (opt->scale)[%i] too large - decreased to %e\n",
00387 i,large);
00388 (opt->scale)[i] = large;
00389 };
00390 };
00391 };
00392 if ( itlin_code == 0 )
00393 {
00394 if ( opt->i_max <= 0 )
00395 { opt->i_max = MIN(10,n);
00396 if ( ERRORLEVEL > 1 )
00397 fprintf(ERROR,
00398 " Warning: opt->i_max not set; is reset to i_max=%i\n",opt->i_max);
00399 }
00400 else if ( opt->i_max > n )
00401 { opt->i_max = MIN(10,n);
00402 if ( ERRORLEVEL > 1 )
00403 fprintf(ERROR,
00404 " Warning: opt->i_max is greater than n; is reset to i_max=%i\n",opt->i_max);
00405 };
00406 }
00407 else if ( itlin_code == 1 )
00408 {
00409 if ( opt->i_max <= 0 )
00410 { opt->i_max = 10;
00411 if ( ERRORLEVEL > 1 )
00412 fprintf(ERROR,
00413 " Warning: opt->i_max not set; is reset to i_max=%i\n",opt->i_max);
00414 };
00415 if ( opt->rho <= 0.0 )
00416 { opt->rho = 1.0;
00417 if ( ERRORLEVEL > 1 )
00418 fprintf(ERROR,
00419 " Warning: opt->rho not set; is reset to rho=%e\n",opt->rho);
00420 };
00421 };
00422 if ( MONITORLEVEL==0 ) return 0;
00423 fprintf(MONITOR,"\n Problem dimension: n = %i\n",n);
00424 if ( itlin_code == 0 )
00425 fprintf(MONITOR,"\n Prescribed residuum threshold: %e\n",opt->tol);
00426 else if ( itlin_code == 1 )
00427 fprintf(MONITOR,"\n Prescribed relative precision: %e\n",opt->tol);
00428 fprintf(MONITOR," The maximum permitted number of iteration steps is: %i\n",
00429 opt->maxiter);
00430 if ( itlin_code == GMRES || itlin_code == GBIT )
00431 fprintf(MONITOR," The maximum number of iterations before restart is: %i\n",
00432 opt->i_max);
00433 if ( itlin_code == 1 )
00434 {
00435 fprintf(MONITOR," The safety factor is rho = %6e\n",opt->rho);
00436 }
00437 else if ( itlin_code == 2 )
00438 {
00439 if ( opt->convcheck == Absolute )
00440 fprintf(MONITOR," The absolute error will be checked for termination\n");
00441 else if ( opt->convcheck == Relative )
00442 fprintf(MONITOR," The relative error will be checked for termination\n");
00443 };
00444 return 0;
00445 }