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
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;
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; };
339 {
DATA=fopen(
"gmres.data",
"w");
341 { fprintf(
ERROR,
"\n fopen of file gmres.data failed\n");
354 free(itlin_ioctl); itlin_ioctl=NULL;
356 if (data) free(data);
366 if (
RCODE !=0 )
return;
371 for (i=0;i<=i_max+1;i++)
373 if (
RCODE !=0 )
return;
377 fprintf(
MONITOR,
"\n\n iter i norm(res) eta\n\n");
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];
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++)
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];
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);
440 {
RCODE = 1; fprintf(
ERROR,
"\n gmres_qrfact failed with code=%i\n",fail);
445 for (l=0;l<i+1;l++) z[l]=0.0; z[0]=beta;
451 for (l=0;l<i;l++) s+= V[l][j]*z[l];
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];
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 );
472 fprintf(
MONITOR,
" %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
482 fprintf(
ERROR,
"\n Error - Maximum allowed number of iterations exceeded\n");
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);
500 free(itlin_ioctl); itlin_ioctl=NULL;
502 if (data) free(data);
509 for (i=0;i<n;i++) rval += v[i]*v[i];
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;
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];