/*
GMRES - RESiduum minimization based iterative linear solver
* Written by L. Weimann
* Purpose Iterative solution of large linear systems
* Category ???. - Linear systems
* Keywords large linear system, iterative solver
* Version 1.0
* Revision June 2006
* Latest Change June 2006
* Library NewtonLib
* Code C, Double Precision
* Environment Standard C environment on PC's,
workstations and hosts.
* Copyright (c) Konrad-Zuse-Zentrum fuer
Informationstechnik Berlin (ZIB)
Takustrasse 7, D-14195 Berlin-Dahlem
phone : + 49/30/84185-0
fax : + 49/30/84185-125
* Contact Bodo Erdmann
ZIB, Division Scientific Computing,
Department Numerical Analysis and Modelling
phone : + 49/30/84185-185
fax : + 49/30/84185-107
e-mail: erdmann@zib.de
* References:
P. Deuflhard:
Newton Methods for Nonlinear Problems. -
Affine Invariance and Adaptive Algorithms.
Series Computational Mathematics 35, Springer (2004)
---------------------------------------------------------------
* Licence
You may use or modify this code for your own non commercial
purposes for an unlimited time.
In any case you should not deliver this code without a special
permission of ZIB.
In case you intend to use the code commercially, we oblige you
to sign an according licence agreement with ZIB.
* Warranty
This code has been tested up to a certain level. Defects and
weaknesses, which may be included in the code, do not establish
any warranties by ZIB. ZIB does not take over any liabilities
which may follow from acquisition or application of this code.
* Software status
This code is under care of ZIB and belongs to ZIB software class 2.
------------------------------------------------------------
* Parameters description
======================
The calling interface looks as follows:
extern void gmres(int n, double *y, MATVEC *matvec,
PRECON *preconr, PRECON *preconl, double *b,
struct ITLIN_OPT *opt, struct ITLIN_INFO *info)
The structures used within the parameter list are defined
as follows:
---
struct ITLIN_OPT
{
double tol, rho;
int i_max, maxiter;
TERM_CHECK termcheck; /* GMRES only * /
CONV_CHECK convcheck; /* PCG only * /
LOGICAL rescale; /* GBIT only * /
PRINT_LEVEL errorlevel, monitorlevel, datalevel;
FILE *errorfile, *monitorfile, *datafile,
*iterfile, *resfile, *miscfile;
double *scale;
};
where the applicable types used within this structure are defined by
typedef enum {None=0, Minimum=1, Verbose=2, Debug=3} PRINT_LEVEL;
typedef enum { CheckOnRestart=0, CheckEachIter=1 } TERM_CHECK ;
---
struct ITLIN_INFO
{
double precision, normdx, residuum;
int iter, rcode, subcode, nomatvec, noprecon, noprecl, noprecr;
};
---
A detailed description of the parameters follows:
int n :
The number of equations and unknown variables of the linear system.
double *y :
A pointer to an array of double values of size n.
The array must contain on input an initial guess of the linear system
solution, which is used as the start-vector of the iteration.
On output, the pointed array contains an approximate solution vector y*,
which fits the preconditioned residuum reduction condition
||B_l*r^i||/||B_l*r^0|| <= opt->tol,
where ||w|| denotes the Euclidian norm of w, and B_l denotes the matrix
of the left preconditioner.
void *matvec(int n, double *y, double *z);
A pointer to the matrix times vector multiplication user routine.
This routine is required - no default routine is supplied.
The parameters are:
int n input Number of vector components.
double *y input Vector of unknowns, of size n .
double *z output Vector which holds the matrix-vector product A*y.
void *preconr(int n, double *z, double *w);
A pointer to the right preconditioner user routine, which computes
w=B_r*z, where B_r should be an approximation of the inverse of the
matrix A. If a null pointer is supplied, then a dummy preconditioner
routine will be used which realizes the preconditioner matrix
B_r=identity.
int n input Number of vector components.
double *z input Preconditioned iterate, of size n .
double *w output Vector which holds the matrix-vector product B_r*z,
i.e. the original iterate.
void *preconl(int n, double *z, double *w);
A pointer to the left preconditioner user routine, which computes
w=B_l*z, where B_l should be an approximation of the inverse of the
matrix A. If a null pointer is supplied, then a dummy preconditioner
routine will be used which realizes the preconditioner matrix
B_l=identity.
int n input Number of vector components.
double *z input Residual vector, of size n .
double *w output Vector which holds the matrix-vector product B_l*z,
i.e. the preconditioned residuum.
double *b :
A pointer to an array of double values of size n.
The pointed array must hold the right hand side of the linear system
to solve.
struct ITLIN_OPT *opt:
A pointer to an options structure. The pointed fields of the structure
contain input options to GMRES.
opt->tol is of type double and must contain the error threshold
which the (left preconditioned) residuum norm reduction quotient must fit.
opt->maxiter is of type int and must contain the maximum number of allowed
iterations. if a zero or negative value is supplied, then the maximum
iteration count will be set to 100.
opt->i_max is of type int and must contain the maximum number of
iterations before a restart occurs. The main portion of used memory
by GMRES depends on i_max, such that n*i_max elements of double
storage will be used. If a nonpositive value is supplied, then i_max
is set to 10.
opt->termcheck is of type TERM_CHECK.
If set to CheckOnRestart, then the residuum norm reduction quotient
will be only checked when a restart occurs - such saving additional
computation effort which is necessary on each intermediate iteration
step to compute the quantity ||B_l*A*B_r*y^i||: One preconr call,
one matvec call, and one preconl call.
If set to CheckEachIter, then additional computations on each iteration
will be done to compute the (left preconditioned) residuum. This roughly
doubles the number of calls to the routines matvec, preconr and preconl.
For additional info, refer to the description of the parameter *y above.
opt->errorlevel is of type PRINT_LEVEL. If it is set to level None,
then no error message will be printed if an error condition occurs.
If it is set to level Minimum or any higher level, then error messages
will be printed, if appropriate.
opt->monitorlevel is of type PRINT_LEVEL. If it is set to level None,
then no monitor output will be printed.
If it is set to level Minimum, a few infomation will be printed.
If set to level Verbose, then some infomation about each iteration
step, fitting into a single line, will be printed. This only applies,
whenopt->termcheck=CheckEachIter is set, otherwise information will
be only printed out when a restart occurs. The higher level Debug
is reserved for future additional information output.
opt->datalevel is of type PRINT_LEVEL. If it is set to level None,
then no data output will be printed.
If it is set to level Minimum, then the values of the initial iteration
vector x and the final vector x will be printed.
If set to level Verbose, then the iteration vector x will be printed for
each step. The higher level Debug is reserved for future additional
information output.
opt->errorfile is of type pointer to FILE, as defined by the
header file. It must be set either to a NULL pointer, stderr, stdout,
or to another file pointer which has been initialized by a fopen call.
If it is set to NULL, opt->errorfile will be set to stdout. The error
messages will be printed to opt->errorfile.
opt->monitorfile is of type pointer to FILE, as defined by the
header file. It must be set either to a NULL pointer, stderr, stdout,
or to another file pointer which has been initialized by a fopen call.
If it is set to NULL, opt->monitorfile will be set to stdout. The monitor
output will be printed to opt->monitorfile.
opt->datafile is of type pointer to FILE, as defined by the
header file. It must be set either to a NULL pointer, stderr, stdout,
or to another file pointer which has been initialized by a fopen call.
If it is set to NULL, a file named "gmres.data" will be opened by a
fopen call and opt->datafile will be set to the filepointer which the
fopen returns. The data output will be printed to opt->datafile.
opt->iterfile is of type pointer to FILE, as defined by the
header file. It must be set either to a NULL pointer or to file pointer
which has been initialized by a fopen call. The iteration number and
the iteration vector will be written out to the associated file, for
each iteration step. If opt->iterfile is set to NULL, no such
data will be written out.
opt->resfile is of type pointer to FILE, as defined by the
header file. It must be set either to a NULL pointer or to file pointer
which has been initialized by a fopen call. The iteration number and
the residuum vector will be written out to the associated file, for
each iteration step. If opt->resfile is set to NULL, no such
data will be written out.
opt->miscfile is of type pointer to FILE, as defined by the
header file. It must be set either to a NULL pointer or to file pointer
which has been initialized by a fopen call. The iteration number, an
identification number of the calling code (0 for GMRES), the norm of the
(left preconditioned) residuum, followed by three zero floating point
values as placeholders, will be written out, for each iteration step,
if opt->termcheck=CheckEachIter is set. Otherwise, if
opt->termcheck=CheckOnRestart is set, data will be only written out
when a restart occurs.
If opt->miscfile is set to NULL, no such data will be written out.
Note: The output to the files opt->iterfile, opt->resfile and
opt->miscfile is written as a single for each iteration step.
Such, the data in these files are suitable as input to the
graphics utility GNUPLOT.
struct ITLIN_INFO *info:
A pointer to an info structure. The pointed fields of the structure
are set output info of GMRES.
info->precision is of type double and is set to the achieved residual
reduction ||r^i||/||r^0||.
info->iter is set to number of iteration steps done.
info->nomatvec is set to the number of done calls to the matrix times
vector multiplication user routine matvec.
info->noprecr is set to the number of done calls to the right
preconditioner user routine preconr or the dummy preconditioner routine,
if the user didn't supply a right preconditioner routine.
info->noprecl is set to the number of done calls to the left
preconditioner user routine preconl or the dummy preconditioner routine,
if the user didn't supply a left preconditioner routine.
info->rcode is set to the return-code of GMRES. A return-code 0
means that GMRES has terminated sucessfully. For the meaning of a
nonzero return-code, see the error messages list below.
The following error conditions may occur: (returned via info->rcode)
--------------------------------------------------------------------
-999 routine zibnum_fwalloc failed to allocate double memory via malloc.
-997 routine zibnum_pfwalloc failed to allocate double pointer memory
via malloc.
-995 Internal i/o control block could not be allocated via malloc.
-994 Internally used data structure could not be allocated via malloc.
-989 Default data-output file could not be opened via fopen call.
-99 NULL pointer supplied for matvec - the matrix times vector routine
must be defined!
2 Maximum number of iterations (as set by opt->maxiter) exceeded.
20 Nonpositive input for dimensional parameter n.
21 Nonpositive value for opt->tol supplied.
Summary of changes:
-------------------
Version Date Changes
1.0 2006/11/30 Initial Release.
*/
#include
#include
#include
#include "itlin.h"
double gmres_norm2(int n, double *v);
int gmres_qrfact(int n, int lda, double *a, double *q, int ijob);
void gmres_qrsolv(int n, int lda, double *a, double *q, double *b);
#define MAXITER_DEFAULT 100
extern struct ITLIN_IO *itlin_ioctl;
extern void gmres(int n, double *y, MATVEC *matvec,
PRECON *preconr, PRECON *preconl, double *b,
struct ITLIN_OPT *opt, struct ITLIN_INFO *info)
{
int i, j, l, m, m1, m11, l1, iter=0, fail,ijob,
nomatvec=0, nopreconl=0,
nopreconr=0, i_max, max_iter = opt->maxiter;
double **V, *r, *y0, *w, *vi, *vip1, *vl, *Hi, *h, *z, *q;
double s, normvip1, beta, beta0, beta1, eta, tol=opt->tol;
LOGICAL stop_iter, io_allocated=False;
TERM_CHECK termcheck = opt->termcheck;
struct ITLIN_DATA *data=malloc(sizeof(struct ITLIN_DATA));
if (!itlin_ioctl) itlin_ioctl=malloc(sizeof(struct ITLIN_IO));
if (!itlin_ioctl)
{ fprintf(stderr,"\n could not allocate output controlblock\n");
RCODE=-995; return; }
else
io_allocated = True;
if (!data)
{ fprintf(stderr,"\n could not allocate struct data\n");
RCODE=-994; return; };
data->codeid = GMRES;
data->normdx = 0.0;
data->tau = 0.0;
data->t = 0.0;
data->mode = Initial;
ERRORLEVEL = opt->errorlevel;
MONITORLEVEL = opt->monitorlevel;
DATALEVEL = opt->datalevel;
ERROR = opt->errorfile;
MONITOR = opt->monitorfile;
DATA = opt->datafile;
FITER = opt->iterfile;
FRES = opt->resfile;
FMISC = opt->miscfile;
if ( !ERROR && ERRORLEVEL>0 ) ERROR = stdout;
if ( !MONITOR && MONITORLEVEL>0 ) MONITOR = stdout;
if ( !DATA && DATALEVEL>0 )
{ DATA=fopen("gmres.data","w");
if (!DATA && ERRORLEVEL>0)
{ fprintf(ERROR,"\n fopen of file gmres.data failed\n");
RCODE=-989; return;
};
};
opt->errorfile = ERROR;
opt->monitorfile = MONITOR;
opt->datafile = DATA;
if ( MONITORLEVEL > 0 ) fprintf(MONITOR,"\n GMRES - Version 1.0\n");
if ( max_iter <= 0 ) max_iter = MAXITER_DEFAULT;
RCODE = itlin_parcheck_and_print(n,matvec,opt,0);
if ( RCODE !=0 )
{ if (io_allocated) {free(itlin_ioctl); itlin_ioctl=NULL;};
if (data) free(data);
return;
};
i_max = opt->i_max;
if (!preconr) preconr = &itlin_noprecon;
if (!preconl) preconl = &itlin_noprecon;
RCODE = zibnum_pfwalloc(i_max+2,&V,"V"); if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc(i_max+1,&h,"h"); if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc(i_max+1,&z,"z"); if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc((i_max+1)*(i_max+1),&Hi,"Hi");
if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc(n,&r,"r"); if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc(n,&y0,"y0"); if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc(n,&w,"w"); if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc(2*(i_max+1),&q,"q"); if ( RCODE !=0 ) return;
for (i=0;i<=i_max+1;i++)
{ RCODE = zibnum_fwalloc(n,&V[i],"V[i]");
if ( RCODE !=0 ) return;
};
/* Initialization */
if ( MONITORLEVEL>1 )
fprintf(MONITOR,"\n\n iter i norm(res) eta\n\n");
data->res = r;
restart:
if ( iter > 0 && MONITORLEVEL > 1 && termcheck==CheckEachIter )
fprintf(MONITOR," > RESTART\n");
for (j=0;jtol; i++)
{
if ( i==1 || termcheck==CheckEachIter )
{
data->residuum = beta1; itlin_dataout(iter,n,y,data);
data->mode = Intermediate;
if ( MONITORLEVEL>1 )
fprintf(MONITOR," %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
};
/* I. Orthogonalization */
vi = V[i-1];
preconr(n,vi,V[i]); nopreconr++;
matvec(n,V[i],w); nomatvec++;
preconl(n,w,r); nopreconl++;
for (l=0;l*1 )
{ for (l=0;l**1 )
fprintf(MONITOR," %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
if ( eta > tol && iter < max_iter ) goto restart;
else RCODE = ( eta <= tol ? 0 : 2 );
if ( MONITORLEVEL>1 && termcheck==CheckOnRestart )
fprintf(MONITOR," %6i %2i %10.3e %10.3e\n",iter,i,beta1,eta);
data->mode = ( RCODE == 0 ? Solution : Final );
data->residuum = beta1; itlin_dataout(iter,n,y,data);
if (info->res) for (j=0;jres)[j]=r[j];
errorexit:
if ( ERRORLEVEL > 0 && RCODE != 0 )
{
switch ( RCODE )
{
case 2:
fprintf(ERROR,"\n Error - Maximum allowed number of iterations exceeded\n");
break;
default :
fprintf(ERROR,"\n Error, code=%i\n",RCODE);
};
};
for (i=0;i<=i_max+1;i++) free(V[i]);
free(V); free(h); free(z); free(Hi); free(r); free(y0);
info->precision = eta;
info->iter = iter;
info->nomatvec = nomatvec;
info->noprecl = nopreconl;
info->noprecr = nopreconr;
}
double gmres_norm2(int n, double *v)
{ int i;
double rval = 0.0;
for (i=0;i=fabs(t1) )
{ t=t1/t2; s=-1.0/sqrt(1.0+t*t); c=-s*t; }
else
{ t=t2/t1; c=1.0/sqrt(1.0+t*t); s=-c*t; };
q[i]=c; q[i+1]=s; a[j1]=c*t1-s*t2;
if( a[j1]==0.0 ) info = k;
};
}
else
{ for (k=0;k=fabs(t1) )
{ t=t1/t2; s=-1.0/sqrt(1.0+t*t); c=-s*t; }
else
{ t=t2/t1; c=1.0/sqrt(1.0+t*t); s=-c*t; };
i=2*n-2; q[i]=c; q[i+1]=s; a[j1]=c*t1-s*t2;
if( a[j1]==0.0 ) info = n-1;
};
return info;
}
void gmres_qrsolv(int n, int lda, double *a, double *q, double *b)
/*
Translation of Fortran routine DHELS from file DGMRES of
library=slatec(slap) package.
*/
{ int i, k, kb;
double t, t1, t2, c, s;
for (k=0;k*