/*
PCG - Preconditioned Conjugate Gradient -
Energy norm 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 pcg(int n, double *y,
MATVEC *matvec, PRECON *precon, 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 { Absolute=0, Relative=1 } CONV_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 either the absolute error condition epsilon_i <= opt->tol,
or the relative error condition sqrt(epsilon_i)/ ||y^i||_A <= opt->tol
where ||y||_A = sqrt() denotes the A-Matrix energy norm
(__ denotes the ordinary scalar product of u and v), and epsilon_i
is an estimate of the quantity square(||y^solution-y^i||_A).
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 *precon(int n, double *z, double *w);
A pointer to the preconditioner user routine, which computes w=B*z,
where B 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=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*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 pcg.
opt->tol is of type double and must contain the error threshold
which the energy norm error (absolute or relative) of the final iterate
y^i 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->convcheck is of type CONV_CHECK.
If set to Absolute, then the absolute error will be checked against the
prescribed tolerance opt->tol.
If set to Relative, then the relative error will be checked against the
prescribed tolerance.
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. 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 "pcg.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 (2 for PCG), the preconditioner
(B-matrix) energy norm of the residuum, the A-matrix energy norm of the
correction, the quantity epsilon_i of the error estimate, and the quantity
gamma_i^2 (used to compute epsilon_i) will be written out, for
each iteration step. 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 PCG.
info->precision is of type double and is set to the absolute error
estimate sqrt(epsilon_i), if opt->convcheck=Absolute is set, or to the
relative error estimate sqrt(epsilon_i)/||y^i||_A,
if opt->convcheck=Relative is set.
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->noprecon is set to the number of done calls to the preconditioner
user routine precon or the dummy preconditioner routine, if the user
didn't supply a preconditioner routine.
info->rcode is set to the return-code of PCG. A return-code 0
means that PCG 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.
-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.
5 The matrix A is not positive definite, and such, the energy norm
||y||_A is not properly defined. PCG cannot proceed the iteration.
6 The preconditioner B is not positive definite, and such, the energy
norm ||r||_B is not properly defined. PCG cannot proceed the
iteration.
20 Nonpositive input for dimensional parameter n.
21 Nonpositive value for opt->tol supplied.
Summary of changes:
-------------------
Version Date Changes
0.1 2006/06/13 Initial Prerelease.
*/
#include
#include
#include
#include "itlin.h"
#define GAMMA_SIZE 50 /* max. number of gamma_i values to be stored */
#define M_NEW 10 /* recompute m every N_NEW iterations */
#define M_INITIAL 5 /* initial m */
#define M_MAX 25 /* maximum permitted m */
#define MAXITER_DEFAULT 100
extern struct ITLIN_IO *itlin_ioctl;
double pcg_sprod(int n, double *v1, double *v2);
extern void pcg(int n, double *y, MATVEC *matvec, PRECON *precon, double *b,
struct ITLIN_OPT *opt, struct ITLIN_INFO *info)
{
int i, j, l, m=0, igamma=0, nomatvec=0, noprecon=0,
gamma_halfsize = GAMMA_SIZE/2,
max_iter = opt->maxiter;
CONV_CHECK convcheck;
double gamma[GAMMA_SIZE];
double *r, *rbar, *p, *w;
double log001 = log(0.01);
double gamma_sum=1.0e10, sigmai, sigmaip1, alpha, alphainv, beta,
normpsq, normysq, delta, tol=opt->tol;
LOGICAL io_allocated=False, stop_iter=False;
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 = PCG;
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("pcg.data","w");
if (!DATA && ERRORLEVEL>0)
{ fprintf(ERROR,"\n fopen of file pcg.data failed\n");
RCODE=-989; return;
};
};
opt->errorfile = ERROR;
opt->monitorfile = MONITOR;
opt->datafile = DATA;
if ( max_iter <= 0 ) max_iter = MAXITER_DEFAULT;
if ( MONITORLEVEL > 0 ) fprintf(MONITOR,"\n PCG - Version 0.1\n");
RCODE = itlin_parcheck_and_print(n,matvec,opt,2);
if ( RCODE !=0 )
{ if (io_allocated) {free(itlin_ioctl); itlin_ioctl=NULL;};
if (data) free(data);
return;
};
convcheck = opt->convcheck;
if (!precon) precon = &itlin_noprecon;
RCODE = 0;
RCODE = zibnum_fwalloc(n,&r,"r"); if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc(n,&rbar,"rbar"); if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc(n,&p,"p"); if ( RCODE !=0 ) return;
RCODE = zibnum_fwalloc(n,&w,"w"); if ( RCODE !=0 ) return;
data->res = r;
if ( MONITORLEVEL > 1 )
fprintf(MONITOR,"\n Iter norm(res) Norm(dy) gamma_sum gamma m\n\n");
/* Initialization */
matvec(n,y,r); nomatvec++; for (j=0;jresiduum = sqrt(sigmai);
data->normdx = alphainv*sqrt(normpsq);
itlin_dataout(i,n,y,data);
data->mode = Intermediate;
for (j=0;j= GAMMA_SIZE )
{ for (l=0;l M_MAX ) m = M_MAX; if ( m < 0 ) m=M_INITIAL;
};
/* recompute completely gamma_sum each time, otherwise leading
bits may be truncated by computing differences */
gamma_sum = 0.0;
for (l=MAX(0,igamma-m+1);l<=igamma;l++) gamma_sum += gamma[l];
data->tau = gamma_sum; data->t = gamma[igamma];
for (j=0;j 1 )
fprintf(MONITOR," %6i %10.3e %10.3e %10.3e %10.3e %3i\n",i,
data->residuum,data->normdx,gamma_sum,gamma[igamma],m);
for (j=0;jmode = ( RCODE==0 ? Solution : Final );
data->residuum = sigmai;
itlin_dataout(i,n,y,data);
if ( ERRORLEVEL > 0 && RCODE != 0 )
{
switch ( RCODE )
{
case 2:
fprintf(ERROR,"\n Error - Maximum allowed number of iterations exceeded\n");
break;
case 5:
fprintf(ERROR,"\n Error - Matrix not positive definite\n");
break;
case 6:
fprintf(ERROR,"\n Error - Preconditioner not positive definite\n");
break;
default :
fprintf(ERROR,"\n Error, code=%i\n",RCODE);
};
};
free(r); free(rbar); free(p); free(w);
if (io_allocated) {free(itlin_ioctl); itlin_ioctl=NULL;};
free(data);
info->precision = ( convcheck==Absolute ? sqrt(gamma_sum) : delta );
info->iter = i;
info->nomatvec = nomatvec;
info->noprecon = noprecon;
}
double pcg_sprod(int n, double *v1, double *v2)
{ int i;
double rval = 0.0;
for (i=0;i__