diff -N -r -c HLib-1.3/Library/krylov.c HLib-1.3new/Library/krylov.c
*** HLib-1.3/Library/krylov.c	2004-12-12 17:42:32.000000000 +0100
--- HLib-1.3new/Library/krylov.c	2005-09-08 18:41:42.121728736 +0200
***************
*** 153,158 ****
--- 153,159 ----
  
    gm->nmax = nmax;
    gm->kmax = kmax+1;
+   gm->eps = 1e-8;
  
    return gm;
  }
***************
*** 179,185 ****
    double *g = gm->g;
    double *r = gm->r;
    double *tmp = gm->tmp;
!   double xx, yy, invnorm, rhsnorm=1.0;
    int kmax = gm->kmax;
    int kmax1 = kmax+1;
    int k, k1;
--- 180,186 ----
    double *g = gm->g;
    double *r = gm->r;
    double *tmp = gm->tmp;
!   double xx, yy, invnorm, newnorm, rhsnorm;
    int kmax = gm->kmax;
    int kmax1 = kmax+1;
    int k, k1;
***************
*** 203,218 ****
      }
      
      r[0] = dnrm2_(&n, V, eins_);
-     if(stop == HLIB_STOP_RELATIVE && iter==0) rhsnorm = r[0];
      invnorm = 1.0 / r[0];
      dscal_(&n, &invnorm, V, eins_);
  
      if(gm->disp)
        gm->disp(0, fabs(r[0]), 0);
      
      k = 0;
      k1 = 1;
      do {
        /* Multiply last search direction by P A or A */
        if(gm->prcd) {
  	gm->eval(V+k*n, tmp, gm->data);
--- 204,222 ----
      }
      
      r[0] = dnrm2_(&n, V, eins_);
      invnorm = 1.0 / r[0];
      dscal_(&n, &invnorm, V, eins_);
  
+     rhsnorm = (stop == HLIB_STOP_RELATIVE && iter == 0 ? r[0] : 1.0);
+ 
      if(gm->disp)
        gm->disp(0, fabs(r[0]), 0);
      
      k = 0;
      k1 = 1;
      do {
+       iter++;
+ 
        /* Multiply last search direction by P A or A */
        if(gm->prcd) {
  	gm->eval(V+k*n, tmp, gm->data);
***************
*** 220,225 ****
--- 224,232 ----
        }
        else
  	gm->eval(V+k*n, V+k1*n, gm->data);
+ 
+       /* Compute norm of new search direction */
+       newnorm = dnrm2_(&n, V+k1*n, eins_);
        
        /* Compute scalar products with old directions */
        dgemv_("Transposed", &n, &k1,
***************
*** 240,245 ****
--- 247,258 ----
  	H[k1+i*kmax1] = 0.0;
        H[k1+k*kmax1] = dnrm2_(&n, V+k1*n, eins_);
        
+       /* Check for invariant subspace */
+       if(!(fabs(H[k1+k*kmax1]) >= gm->eps * newnorm)) {
+ 	k = k1;
+ 	break;
+       }
+ 
        /* Normalize new search direction */
        invnorm = 1.0 / H[k1+k*kmax1];
        dscal_(&n, &invnorm, V+k1*n, eins_);
***************
*** 268,286 ****
        r[k] = g[2*k] * xx;
        r[k1] = -g[2*k+1] * xx;
  
-       /* Check for invariant subspace 
-       if(fabs(H[k+k*kmax1]) < eps * fabs(r[k]))
- 	break;
-       */
- 
        k++;
        k1 = k + 1;
-       iter++;
  
        if(gm->disp)
  	gm->disp(1, fabs(r[k]), iter);
      }
!     while(k1 < kmax && iter < maxiter && !(fabs(r[k]) < eps*rhsnorm));
  
      dtrsv_("Upper triangular", "Not Transposed", "Non-unit diagonal",
  	   &k, H, &kmax1, r, eins_);
--- 281,295 ----
        r[k] = g[2*k] * xx;
        r[k1] = -g[2*k+1] * xx;
  
        k++;
        k1 = k + 1;
  
        if(gm->disp)
  	gm->disp(1, fabs(r[k]), iter);
      }
!     while(k1 < kmax &&
! 	  iter < maxiter &&
! 	  !(fabs(r[k]) < eps*rhsnorm));
  
      dtrsv_("Upper triangular", "Not Transposed", "Non-unit diagonal",
  	   &k, H, &kmax1, r, eins_);
***************
*** 291,297 ****
  	   deins_,
  	   x, eins_);
    }
!   while(iter < maxiter && !(fabs(r[k]) < eps*rhsnorm));
  
    if(gm->disp)
      gm->disp(2, fabs(r[k]), iter);
--- 300,308 ----
  	   deins_,
  	   x, eins_);
    }
!   while(iter < maxiter &&
! 	!(fabs(r[k]) < eps*rhsnorm) &&
! 	k1 == kmax);
  
    if(gm->disp)
      gm->disp(2, fabs(r[k]), iter);
diff -N -r -c HLib-1.3/Library/krylov.h HLib-1.3new/Library/krylov.h
*** HLib-1.3/Library/krylov.h	2004-12-12 17:42:32.000000000 +0100
--- HLib-1.3new/Library/krylov.h	2005-09-08 17:36:37.731286304 +0200
***************
*** 39,54 ****
  #include "sparsematrix.h"
  
  struct _conjgrad {
!   double *r;
!   double *p;
!   double *a;
    
    void (*eval)(const double *x, double *y, void *data);
    void (*prcd)(const double *x, double *y, void *data);
    void (*disp)(int mode, double res, int iter);
    void *data;
  
!   int nmax;
  };
  
  pconjgrad
--- 39,54 ----
  #include "sparsematrix.h"
  
  struct _conjgrad {
!   double *r;			/* Residual */
!   double *p;			/* Search direction */
!   double *a;			/* Matrix applied to p */
    
    void (*eval)(const double *x, double *y, void *data);
    void (*prcd)(const double *x, double *y, void *data);
    void (*disp)(int mode, double res, int iter);
    void *data;
  
!   int nmax;			/* Maximal dimension */
  };
  
  pconjgrad
***************
*** 62,80 ****
  	       double eps, int maxiter);
  
  struct _gmres {
!   double *V;
!   double *H;
!   double *g;
!   double *r;
!   double *tmp;
  
    void (*eval)(const double *x, double *y, void *data);
    void (*prcd)(const double *x, double *y, void *data);
    void (*disp)(int mode, double res, int iter);
    void *data;
  
!   int nmax;
!   int kmax;
  };
  
  pgmres
--- 62,81 ----
  	       double eps, int maxiter);
  
  struct _gmres {
!   double *V;			/* Krylov basis */
!   double *H;			/* Matrix in Krylov basis */
!   double *g;			/* Givens rotations */
!   double *r;			/* Right-hand side in Krylov basis */
!   double *tmp;			/* Auxiliary vector */
  
    void (*eval)(const double *x, double *y, void *data);
    void (*prcd)(const double *x, double *y, void *data);
    void (*disp)(int mode, double res, int iter);
    void *data;
  
!   int nmax;			/* Maximal dimension */
!   int kmax;			/* Maximal dimension of Krylov spaces */
!   double eps;			/* Tolerance for detection of invariant subspaces */
  };
  
  pgmres