292 double gmres_norm2(
int n,
double *v);
293 int gmres_qrfact(
int n,
int lda,
double *a,
double *q,
int ijob);
294 void gmres_qrsolv(
int n,
int lda,
double *a,
double *q,
double *b);
296 #define MAXITER_DEFAULT 100
298 extern struct ITLIN_IO *itlin_ioctl;
300 extern void gmres(
int n,
double *y, MATVEC *
matvec,
301 PRECON *preconr, PRECON *preconl,
double *b,
302 struct ITLIN_OPT *opt,
struct ITLIN_INFO *info)
304 int i, j, l,
m, m1, m11, l1, iter=0, fail,ijob,
305 nomatvec=0, nopreconl=0,
306 nopreconr=0, i_max, max_iter = opt->maxiter;
307 double **V, *r, *y0, *w, *vi, *vip1, *vl, *Hi, *h, *z, *q;
308 double s, normvip1, beta, beta0, beta1, eta, tol=opt->tol;
309 LOGICAL stop_iter, io_allocated=False;
310 TERM_CHECK termcheck = opt->termcheck;
311 struct ITLIN_DATA *data=malloc(
sizeof(
struct ITLIN_DATA));
313 if (!itlin_ioctl) itlin_ioctl=malloc(
sizeof(
struct ITLIN_IO));
315 { fprintf(stderr,
"\n could not allocate output controlblock\n");
316 RCODE=-995;
return; }
320 { fprintf(stderr,
"\n could not allocate struct data\n");
321 RCODE=-994;
return; };
322 data->codeid = GMRES;
326 data->mode = Initial;
327 ERRORLEVEL = opt->errorlevel;
328 MONITORLEVEL = opt->monitorlevel;
329 DATALEVEL = opt->datalevel;
330 ERROR = opt->errorfile;
331 MONITOR = opt->monitorfile;
332 DATA = opt->datafile;
333 FITER = opt->iterfile;
335 FMISC = opt->miscfile;
336 if ( !ERROR && ERRORLEVEL>0 ) ERROR = stdout;
337 if ( !MONITOR && MONITORLEVEL>0 ) MONITOR = stdout;
338 if ( !DATA && DATALEVEL>0 )
339 { DATA=fopen(
"gmres.data",
"w");
340 if (!DATA && ERRORLEVEL>0)
341 { fprintf(ERROR,
"\n fopen of file gmres.data failed\n");
345 opt->errorfile = ERROR;
346 opt->monitorfile = MONITOR;
347 opt->datafile = DATA;
348 if ( MONITORLEVEL > 0 ) fprintf(MONITOR,
"\n GMRES - Version 0.1\n");
349 if ( max_iter <= 0 ) max_iter = MAXITER_DEFAULT;
350 RCODE = itlin_parcheck_and_print(n,matvec,opt,0);
354 free(itlin_ioctl); itlin_ioctl=NULL;
356 if (data) free(data);
360 if (!preconr) preconr = &itlin_noprecon;
361 if (!preconl) preconl = &itlin_noprecon;
362 RCODE = zibnum_pfwalloc(i_max+2,&V,
"V");
if ( RCODE !=0 )
return;
363 RCODE = zibnum_fwalloc(i_max+1,&h,
"h");
if ( RCODE !=0 )
return;
364 RCODE = zibnum_fwalloc(i_max+1,&z,
"z");
if ( RCODE !=0 )
return;
365 RCODE = zibnum_fwalloc((i_max+1)*(i_max+1),&Hi,
"Hi");
366 if ( RCODE !=0 )
return;
367 RCODE = zibnum_fwalloc(n,&r,
"r");
if ( RCODE !=0 )
return;
368 RCODE = zibnum_fwalloc(n,&y0,
"y0");
if ( RCODE !=0 )
return;
369 RCODE = zibnum_fwalloc(n,&w,
"w");
if ( RCODE !=0 )
return;
370 RCODE = zibnum_fwalloc(2*(i_max+1),&q,
"q");
if ( RCODE !=0 )
return;
371 for (i=0;i<=i_max+1;i++)
372 { RCODE = zibnum_fwalloc(n,&V[i],
"V[i]");
373 if ( RCODE !=0 )
return;
376 if ( MONITORLEVEL>1 )
377 fprintf(MONITOR,
"\n\n iter i norm(res) eta\n\n");
380 if ( iter > 0 && MONITORLEVEL > 1 && termcheck==CheckEachIter )
381 fprintf(MONITOR,
" > RESTART\n");
382 for (j=0;j<n;j++) y0[j]=y[j];
383 preconr(n,y,V[0]); nopreconr++;
384 matvec(n,V[0],w); nomatvec++;
385 preconl(n,w,r); nopreconl++;
386 for (j=0;j<n;j++) r[j] = b[j]-r[j];
387 beta = gmres_norm2(n,r);
389 { beta0 = beta; eta = ( beta0==0.0 ? 0.0 : 1.0 ); }
394 for (j=0;j<n;j++) vi[j] = r[j]/beta;
396 for (i=1; i<=i_max && eta>tol; i++)
398 if ( i==1 || termcheck==CheckEachIter )
400 data->residuum = beta1; itlin_dataout(iter,n,y,data);
401 data->mode = Intermediate;
402 if ( MONITORLEVEL>1 )
403 fprintf(MONITOR,
" %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
408 preconr(n,vi,V[i]); nopreconr++;
409 matvec(n,V[i],w); nomatvec++;
410 preconl(n,w,r); nopreconl++;
413 for (j=0;j<n;j++) s += vl[j]*r[j]; h[l] = s;
417 { s = 0.0;
for (l=0;l<i;l++) s += V[l][j]*h[l];
422 normvip1 = gmres_norm2(n,vip1);
423 for (j=0;j<n;j++) vip1[j] /= normvip1;
427 {
for (l=0;l<i;l++) Hi[(i_max+1)*(i-1)+l] = h[l];
428 for (l=0;l<i-1;l++) Hi[(i_max+1)*l+i] = 0.0;
429 Hi[(i_max+1)*(i-1)+i] = normvip1;
433 Hi[0] = h[0]; Hi[1] = normvip1;
437 ijob = ( i==1 ? 1 : 2);
438 fail = gmres_qrfact(i,i_max+1,Hi,q,ijob);
440 { RCODE = 1; fprintf(ERROR,
"\n gmres_qrfact failed with code=%i\n",fail);
443 if ( termcheck==CheckEachIter || i==i_max )
445 for (l=0;l<i+1;l++) z[l]=0.0; z[0]=beta;
446 gmres_qrsolv(i,i_max+1,Hi,q,z);
451 for (l=0;l<i;l++) s+= V[l][j]*z[l];
454 if ( termcheck==CheckEachIter )
456 preconr(n,y,V[i+1]); nopreconr++;
457 matvec(n,V[i+1],w); nomatvec++;
458 preconl(n,w,r); nopreconl++;
459 for (j=0;j<n;j++) r[j] = b[j]-r[j];
460 beta1 = gmres_norm2(n,r);
467 if ( termcheck==CheckEachIter && MONITORLEVEL>1 )
468 fprintf(MONITOR,
" %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
469 if ( eta > tol && iter < max_iter )
goto restart;
470 else RCODE = ( eta <= tol ? 0 : 2 );
471 if ( MONITORLEVEL>1 && termcheck==CheckOnRestart )
472 fprintf(MONITOR,
" %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
473 data->mode = ( RCODE == 0 ? Solution : Final );
474 data->residuum = beta1; itlin_dataout(iter,n,y,data);
477 if ( ERRORLEVEL > 0 && RCODE != 0 )
482 fprintf(ERROR,
"\n Error - Maximum allowed number of iterations exceeded\n");
485 fprintf(ERROR,
"\n Error, code=%i\n",RCODE);
488 for (i=0;i<=i_max+1;i++) free(V[i]);
489 free(V); free(h); free(z); free(Hi); free(r); free(y0);
491 info->precision = eta;
493 info->nomatvec = nomatvec;
494 info->noprecl = nopreconl;
495 info->noprecr = nopreconr;
500 free(itlin_ioctl); itlin_ioctl=NULL;
502 if (data) free(data);
506 double gmres_norm2(
int n,
double *v)
509 for (i=0;i<n;i++) rval += v[i]*v[i];
513 int gmres_qrfact(
int n,
int lda,
double *a,
double *q,
int ijob)
518 {
int i, j, j1, k,info=0;
519 double t, t1, t2, c, s;
522 {
for (j=0;j<k-1;j++)
523 { i=2*j; c=q[i]; s=q[i+1];
524 j1=j+k*lda; t1=a[j1]; t2=a[j1+1];
525 a[j1]=c*t1-s*t2; a[j1+1]=s*t1+c*t2;
527 i=2*k; j1=k+lda*k; t1=a[j1]; t2=a[j1+1];
528 if ( t2==0.0 ) { c=1.0; s=0.0; }
529 else if ( fabs(t2)>=fabs(t1) )
530 { t=t1/t2; s=-1.0/sqrt(1.0+t*t); c=-s*t; }
532 { t=t2/t1; c=1.0/sqrt(1.0+t*t); s=-c*t; };
533 q[i]=c; q[i+1]=s; a[j1]=c*t1-s*t2;
534 if( a[j1]==0.0 ) info = k;
538 {
for (k=0;k<n-1;k++)
540 i=2*k; j1=k+lda*(n-1);
541 t1=a[j1]; t2=a[j1+1];
543 a[j1]=c*t1-s*t2; a[j1+1]=s*t1 + c*t2;
546 j1 = n-1+lda*(n-1); t1 = a[j1]; t2 = a[j1+1];
547 if ( t2==0.0 ) { c=1.0; s=0.0; }
548 else if ( fabs(t2)>=fabs(t1) )
549 { t=t1/t2; s=-1.0/sqrt(1.0+t*t); c=-s*t; }
551 { t=t2/t1; c=1.0/sqrt(1.0+t*t); s=-c*t; };
552 i=2*n-2; q[i]=c; q[i+1]=s; a[j1]=c*t1-s*t2;
553 if( a[j1]==0.0 ) info = n-1;
558 void gmres_qrsolv(
int n,
int lda,
double *a,
double *q,
double *b)
564 double t, t1, t2, c, s;
566 { i=2*k; c=q[i]; s=q[i+1];
568 b[k]=c*t1-s*t2; b[k+1]=s*t1+c*t2;
572 k=n-1-kb; b[k] /= a[k+lda*k];
574 for (i=0;i<k;i++) b[i] += t*a[i+k*lda];
static const double m
mass of electron in grams
void matvec(int n, double *x, double *b)
Used for gmres.