373,375d372
<   if (*normres <= tol)
<     return 0;
< 
377,424d373
<   for (ii = 0; ii < itmax; ii++)
<   {
<     rho = scalarproduct (rt, r, ntotal);
<     if (fabs (rho) < rhotol)
<       break;
< 
<     /* compute direction vector p */
<     if (ii == 0)
<     {
< #pragma omp parallel for
<       for (int j = 0; j < ntotal; j++)
< 	p[j] = r[j];
<     }
<     else
<     {
<       beta = (rho / rho1) * (alpha / omega);
< #pragma omp parallel for
<       for (int j = 0; j < ntotal; j++)
< 	p[j] = r[j] + beta * (p[j] - omega * vv[j]);
<     }
< 
<     /* compute direction adjusting vector ph and scalar alpha */
< #pragma omp parallel for
<     for (int j = 0; j < ntotal; j++)
<       ph.d0[j] = 0;
<     for (int j = 0; j < NRELAX; j++)	/* solves JFD*ph = p by relaxation*/
<       relax (ph.d0, nvar, n1, n2, n3, p, ncols, cols, JFD);
< 
<     J_times_dv (nvar, n1, n2, n3, ph, vv, u);	/* vv=J*ph*/
<     alpha = rho / scalarproduct (rt, vv, ntotal);
< #pragma omp parallel for
<     for (int j = 0; j < ntotal; j++)
<       s[j] = r[j] - alpha * vv[j];
< 
<     /* early check of tolerance */
<     *normres = norm2 (s, ntotal);
<     if (*normres <= tol)
<     {
< #pragma omp parallel for
<       for (int j = 0; j < ntotal; j++)
< 	dv.d0[j] += alpha * ph.d0[j];
<       if (output == 1) {
< 	printf ("bicgstab: %5d  %10.3e  %10.3e  %10.3e  %10.3e\n",
< 		ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
<         fflush(stdout);
<       }
<       break;
<     }
426,454c375,453
<     /* compute stabilizer vector sh and scalar omega */
< #pragma omp parallel for
<     for (int j = 0; j < ntotal; j++)
<       sh.d0[j] = 0;
<     for (int j = 0; j < NRELAX; j++)	/* solves JFD*sh = s by relaxation*/
<       relax (sh.d0, nvar, n1, n2, n3, s, ncols, cols, JFD);
< 
<     J_times_dv (nvar, n1, n2, n3, sh, t, u);	/* t=J*sh*/
<     omega = scalarproduct (t, s, ntotal) / scalarproduct (t, t, ntotal);
< 
<     /* compute new solution approximation */
< #pragma omp parallel for
<     for (int j = 0; j < ntotal; j++)
<     {
<       dv.d0[j] += alpha * ph.d0[j] + omega * sh.d0[j];
<       r[j] = s[j] - omega * t[j];
<     }
<     /* are we done? */
<     *normres = norm2 (r, ntotal);
<     if (output == 1) {
<       printf ("bicgstab: %5d  %10.3e  %10.3e  %10.3e  %10.3e\n",
< 	      ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
<       fflush(stdout);
<     }
<     if (*normres <= tol)
<       break;
<     rho1 = rho;
<     if (fabs (omega) < omegatol)
<       break;
---
>   if (*normres > tol) {
>     for (ii = 0; ii < itmax; ii++)
>       {
>         rho = scalarproduct (rt, r, ntotal);
>         if (fabs (rho) < rhotol)
>           break;
> 
>         /* compute direction vector p */
>         if (ii == 0)
>           {
> #pragma omp parallel for
>             for (int j = 0; j < ntotal; j++)
>               p[j] = r[j];
>           }
>         else
>           {
>             beta = (rho / rho1) * (alpha / omega);
> #pragma omp parallel for
>             for (int j = 0; j < ntotal; j++)
>               p[j] = r[j] + beta * (p[j] - omega * vv[j]);
>           }
> 
>         /* compute direction adjusting vector ph and scalar alpha */
> #pragma omp parallel for
>         for (int j = 0; j < ntotal; j++)
>           ph.d0[j] = 0;
>         for (int j = 0; j < NRELAX; j++)	/* solves JFD*ph = p by relaxation*/
>           relax (ph.d0, nvar, n1, n2, n3, p, ncols, cols, JFD);
> 
>         J_times_dv (nvar, n1, n2, n3, ph, vv, u);	/* vv=J*ph*/
>         alpha = rho / scalarproduct (rt, vv, ntotal);
> #pragma omp parallel for
>         for (int j = 0; j < ntotal; j++)
>           s[j] = r[j] - alpha * vv[j];
> 
>         /* early check of tolerance */
>         *normres = norm2 (s, ntotal);
>         if (*normres <= tol)
>           {
> #pragma omp parallel for
>             for (int j = 0; j < ntotal; j++)
>               dv.d0[j] += alpha * ph.d0[j];
>             if (output == 1) {
>               printf ("bicgstab: %5d  %10.3e  %10.3e  %10.3e  %10.3e\n",
>                       ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
>               fflush(stdout);
>             }
>             break;
>           }
> 
>         /* compute stabilizer vector sh and scalar omega */
> #pragma omp parallel for
>         for (int j = 0; j < ntotal; j++)
>           sh.d0[j] = 0;
>         for (int j = 0; j < NRELAX; j++)	/* solves JFD*sh = s by relaxation*/
>           relax (sh.d0, nvar, n1, n2, n3, s, ncols, cols, JFD);
> 
>         J_times_dv (nvar, n1, n2, n3, sh, t, u);	/* t=J*sh*/
>         omega = scalarproduct (t, s, ntotal) / scalarproduct (t, t, ntotal);
> 
>         /* compute new solution approximation */
> #pragma omp parallel for
>         for (int j = 0; j < ntotal; j++)
>           {
>             dv.d0[j] += alpha * ph.d0[j] + omega * sh.d0[j];
>             r[j] = s[j] - omega * t[j];
>           }
>         /* are we done? */
>         *normres = norm2 (r, ntotal);
>         if (output == 1) {
>           printf ("bicgstab: %5d  %10.3e  %10.3e  %10.3e  %10.3e\n",
>                   ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
>           fflush(stdout);
>         }
>         if (*normres <= tol)
>           break;
>         rho1 = rho;
>         if (fabs (omega) < omegatol)
>           break;
455a455
>       }
475a476,478
> 
>   if (*normres <= tol)
>     return 0;
