solve_sparse_crs_bicgstab Module Subroutine

    module subroutine solve_sparse_crs_bicgstab(self, A, b, x, status)
        implicit none
        class(type_solver_sparse_crs_bicgstab), intent(inout) :: self
        class(abst_matrix), intent(in) :: A
        real(real64), intent(inout) :: b(:)
        real(real64), intent(inout) :: x(:)
        integer(int32), intent(inout) :: status

        real(real64) :: rho, rho_old, alpha, beta, omega
        real(real64) :: resid
        integer(int32) :: iter, iN, vector_size

        select type (matrix => A)
        type is (type_crs)

            ! 1:Initialize
            rho = 1.0d0
            rho_old = 1.0d0
            alpha = 1.0d0
            beta = 1.0d0
            omega = 1.0d0

            do iN = 1, self%size
                self%p(iN) = 0.0d0
                self%s(iN) = 0.0d0
                self%phat(iN) = 0.0d0
                self%shat(iN) = 0.0d0

                ! 2: Set an initial value x0
                self%x(iN) = 0.0d0
            end do

            ! 3: r0 = b-Ax0
            call gemv(1.0d0, matrix, self%x(:), 0.0d0, self%r(:))
            ! self%r(:) = matrix * self%x(:)
            do iN = 1, self%size
                self%r(iN) = b(iN) - self%r(iN)
            end do
            ! 4: Create preconditioned matrix
            call self%Create_Preconditioner(matrix)

            ! 5: ^r0 = r0, (r*0, r0)!=0
            do iN = 1, self%size
                self%r0(iN) = self%r(iN)
            end do

            do iter = 1, self%max_iterations
                ! 7: (^r0, rk)
                rho = dot(self%r(:), self%r0(:))
                ! 8: rho check
                if (rho == 0.0d0) then
                    status = 0
                    do iN = 1, self%size
                        x(iN) = self%x(iN)
                    end do
                    return
                end if

                if (iter == 1) then
                    ! 10: p0 = r0
                    do iN = 1, self%size
                        self%p(iN) = self%r(iN)
                    end do
                else
                    ! 12: beta = (rho / rho_old) * (alpha_k / omega_k)
                    beta = (rho / rho_old) * (alpha / omega)
                    ! 13: p_k = r_k + beta_k(p_(k-1) - omega_k * Av)
                    do iN = 1, self%size
                        self%p(iN) = self%r(iN) + beta * (self%p(iN) - omega * self%v(iN))
                    end do
                end if
                ! 15: phat = M^-1 * p
                call self%Apply_Preconditioner(self%p(:), self%phat(:))
                ! 16: v = A * phat
                call gemv(1.0d0, matrix, self%phat(:), 0.0d0, self%v(:))
                ! call SpMV(self%CRS_A, self%phat, self%v)
                ! 17: alpha_k = rho / (^r0, v)
                alpha = rho / dot(self%r0(:), self%v(:))
                ! 18: s = r_k - alpha_k * v
                do iN = 1, self%size
                    self%s(iN) = self%r(iN) - alpha * self%v(iN)
                end do

                ! 19: shat = M^-1 * s
                call self%Apply_Preconditioner(self%s(:), self%shat(:))
                ! 20: t = A * shat
                call gemv(1.0d0, matrix, self%shat(:), 0.0d0, self%t(:))

                ! 21: omega_k = (t,s)/(t,t)
                omega = dot(self%t(:), self%s(:)) / dot(self%t(:), self%t(:))

                ! 22: omega breakdown check
                if (omega == 0.0d0) then
                    status = -1
                    return
                end if

                ! 23: x(i) = x(i-1) + alpha * M^-1 p(i-1) + omega * M^-1 s(i)
                ! 24: r(i) = s(i-1) - omega * AM^-1 s(i-1)
                do iN = 1, self%size
                    self%x(iN) = self%x(iN) + alpha * self%phat(iN) + omega * self%shat(iN)
                    self%r(iN) = self%s(iN) - omega * self%t(iN)
                end do

                ! 25: ||r_k+1||_2
                resid = norm_2(self%r(:))
                if (resid < self%tolerance) then
                    status = 0
                    do iN = 1, self%size
                        x(iN) = self%x(iN)
                    end do
                    return
                end if

                rho_old = rho

                if (was_interrupted()) stop
            end do
            status = -2

        end select
    end subroutine solve_sparse_crs_bicgstab