#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "slu_ddefs.h"
#define epsmac 1.0e-16
extern double ddot_(int *, double [], int *, double [], int *);
extern double dnrm2_(int *, double [], int *);
extern void daxpy_(int *, double *, double [], int *, double [], int *);
extern double dcopy_(int *, double [], int *, double [], int *);
int fgmr(int n,
void (*matvec) (double, double[], double, double[]),
void (*psolve) (int, double[], double[]),
double *rhs, double *sol, double tol, int im, int *itmax, FILE * fits)
{
/*----------------------------------------------------------------------
| *** Preconditioned FGMRES ***
+-----------------------------------------------------------------------
| This is a simple version of the ARMS preconditioned FGMRES algorithm.
+-----------------------------------------------------------------------
| Y. S. Dec. 2000. -- Apr. 2008
+-----------------------------------------------------------------------
| on entry:
|----------
|
| rhs = real vector of length n containing the right hand side.
| sol = real vector of length n containing an initial guess to the
| solution on input.
| tol = tolerance for stopping iteration
| im = Krylov subspace dimension
| (itmax) = max number of iterations allowed.
| fits = NULL: no output
| != NULL: file handle to output " resid vs time and its"
|
| on return:
|----------
| fgmr int = 0 --> successful return.
| int = 1 --> convergence not achieved in itmax iterations.
| sol = contains an approximate solution (upon successful return).
| itmax = has changed. It now contains the number of steps required
| to converge --
+-----------------------------------------------------------------------
| internal work arrays:
|----------
| vv = work array of length [im+1][n] (used to store the Arnoldi
| basis)
| hh = work array of length [im][im+1] (Householder matrix)
| z = work array of length [im][n] to store preconditioned vectors
+-----------------------------------------------------------------------
| subroutines called :
| matvec - matrix-vector multiplication operation
| psolve - (right) preconditionning operation
| psolve can be a NULL pointer (GMRES without preconditioner)
+---------------------------------------------------------------------*/
int maxits = *itmax;
int i, i1, ii, j, k, k1, its, retval, i_1 = 1;
double **hh, *c, *s, *rs, t, t0;
double beta, eps1 = 0.0, gam, **vv, **z;
its = 0;
vv = (double **)SUPERLU_MALLOC((im + 1) * sizeof(double *));
for (i = 0; i <= im; i++)
vv[i] = doubleMalloc(n);
z = (double **)SUPERLU_MALLOC(im * sizeof(double *));
hh = (double **)SUPERLU_MALLOC(im * sizeof(double *));
for (i = 0; i < im; i++)
{
hh[i] = doubleMalloc(i + 2);
z[i] = doubleMalloc(n);
}
c = doubleMalloc(im);
s = doubleMalloc(im);
rs = doubleMalloc(im + 1);
/*---- outer loop starts here ----*/
do
{
/*---- compute initial residual vector ----*/
matvec(1.0, sol, 0.0, vv[0]);
for (j = 0; j < n; j++)
vv[0][j] = rhs[j] - vv[0][j]; /* vv[0]= initial residual */
beta = dnrm2_(&n, vv[0], &i_1);
/*---- print info if fits != null ----*/
if (fits != NULL && its == 0)
fprintf(fits, "%8d %10.2e\n", its, beta);
/*if ( beta < tol * dnrm2_(&n, rhs, &i_1) )*/
if ( !(beta >= tol * dnrm2_(&n, rhs, &i_1)) )
break;
t = 1.0 / beta;
/*---- normalize: vv[0] = vv[0] / beta ----*/
for (j = 0; j < n; j++)
vv[0][j] = vv[0][j] * t;
if (its == 0)
eps1 = tol * beta;
/*---- initialize 1-st term of rhs of hessenberg system ----*/
rs[0] = beta;
for (i = 0; i < im; i++)
{
its++;
i1 = i + 1;
/*------------------------------------------------------------
| (Right) Preconditioning Operation z_{j} = M^{-1} v_{j}
+-----------------------------------------------------------*/
if (psolve)
psolve(n, z[i], vv[i]);
else
dcopy_(&n, vv[i], &i_1, z[i], &i_1);
/*---- matvec operation w = A z_{j} = A M^{-1} v_{j} ----*/
matvec(1.0, z[i], 0.0, vv[i1]);
/*------------------------------------------------------------
| modified gram - schmidt...
| h_{i,j} = (w,v_{i})
| w = w - h_{i,j} v_{i}
+------------------------------------------------------------*/
t0 = dnrm2_(&n, vv[i1], &i_1);
for (j = 0; j <= i; j++)
{
double negt;
t = ddot_(&n, vv[j], &i_1, vv[i1], &i_1);
hh[i][j] = t;
negt = -t;
daxpy_(&n, &negt, vv[j], &i_1, vv[i1], &i_1);
}
/*---- h_{j+1,j} = ||w||_{2} ----*/
t = dnrm2_(&n, vv[i1], &i_1);
while (t < 0.5 * t0)
{
t0 = t;
for (j = 0; j <= i; j++)
{
double negt;
t = ddot_(&n, vv[j], &i_1, vv[i1], &i_1);
hh[i][j] += t;
negt = -t;
daxpy_(&n, &negt, vv[j], &i_1, vv[i1], &i_1);
}
t = dnrm2_(&n, vv[i1], &i_1);
}
hh[i][i1] = t;
if (t != 0.0)
{
/*---- v_{j+1} = w / h_{j+1,j} ----*/
t = 1.0 / t;
for (k = 0; k < n; k++)
vv[i1][k] = vv[i1][k] * t;
}
/*---------------------------------------------------
| done with modified gram schimdt and arnoldi step
| now update factorization of hh
+--------------------------------------------------*/
/*--------------------------------------------------------
| perform previous transformations on i-th column of h
+-------------------------------------------------------*/
for (k = 1; k <= i; k++)
{
k1 = k - 1;
t = hh[i][k1];
hh[i][k1] = c[k1] * t + s[k1] * hh[i][k];
hh[i][k] = -s[k1] * t + c[k1] * hh[i][k];
}
gam = sqrt(pow(hh[i][i], 2) + pow(hh[i][i1], 2));
/*---------------------------------------------------
| if gamma is zero then any small value will do
| affect only residual estimate
+--------------------------------------------------*/
/* if (gam == 0.0) gam = epsmac; */
/*---- get next plane rotation ---*/
if (gam > 0.0)
{
c[i] = hh[i][i] / gam;
s[i] = hh[i][i1] / gam;
}
else
{
c[i] = 1.0;
s[i] = 0.0;
}
rs[i1] = -s[i] * rs[i];
rs[i] = c[i] * rs[i];
/*----------------------------------------------------
| determine residual norm and test for convergence
+---------------------------------------------------*/
hh[i][i] = c[i] * hh[i][i] + s[i] * hh[i][i1];
beta = fabs(rs[i1]);
if (fits != NULL)
fprintf(fits, "%8d %10.2e\n", its, beta);
if (beta <= eps1 || its >= maxits)
break;
}
if (i == im) i--;
/*---- now compute solution. 1st, solve upper triangular system ----*/
rs[i] = rs[i] / hh[i][i];
for (ii = 1; ii <= i; ii++)
{
k = i - ii;
k1 = k + 1;
t = rs[k];
for (j = k1; j <= i; j++)
t = t - hh[j][k] * rs[j];
rs[k] = t / hh[k][k];
}
/*---- linear combination of v[i]'s to get sol. ----*/
for (j = 0; j <= i; j++)
{
t = rs[j];
for (k = 0; k < n; k++)
sol[k] += t * z[j][k];
}
/* calculate the residual and output */
matvec(1.0, sol, 0.0, vv[0]);
for (j = 0; j < n; j++)
vv[0][j] = rhs[j] - vv[0][j]; /* vv[0]= initial residual */
/*---- print info if fits != null ----*/
beta = dnrm2_(&n, vv[0], &i_1);
/*---- restart outer loop if needed ----*/
/*if (beta >= eps1 / tol)*/
if ( !(beta < eps1 / tol) )
{
its = maxits + 10;
break;
}
if (beta <= eps1)
break;
} while(its < maxits);
retval = (its >= maxits);
for (i = 0; i <= im; i++)
SUPERLU_FREE(vv[i]);
SUPERLU_FREE(vv);
for (i = 0; i < im; i++)
{
SUPERLU_FREE(hh[i]);
SUPERLU_FREE(z[i]);
}
SUPERLU_FREE(hh);
SUPERLU_FREE(z);
SUPERLU_FREE(c);
SUPERLU_FREE(s);
SUPERLU_FREE(rs);
*itmax = its;
return retval;
} /*----end of fgmr ----*/
