From 02fb43ba82ad5f3a352d6f8575f7b426faca0767 Mon Sep 17 00:00:00 2001
From: gabrielequatrana <gabriele.quatrana@outlook.com>
Date: Thu, 28 Mar 2024 13:36:02 +0100
Subject: [PATCH] Fixed convergence

---
 base/psblas/psb_dqrfact.f90 |  8 --------
 krylov/psb_dbgmres.f90      | 39 ++++++++++++-------------------------
 2 files changed, 12 insertions(+), 35 deletions(-)

diff --git a/base/psblas/psb_dqrfact.f90 b/base/psblas/psb_dqrfact.f90
index 2b7f2039..0de4ed16 100644
--- a/base/psblas/psb_dqrfact.f90
+++ b/base/psblas/psb_dqrfact.f90
@@ -69,14 +69,6 @@ function psb_dqrfact(x, desc_a, info) result(res)
    if (me == psb_root_) then
       call qr_temp%bld(temp)
       res = qr_temp%qr_fact(info)
-
-      ! TODO Check sulla diagonale di R
-      do i=1,n
-         if (res(i,i) == dzero) then
-            write(*,*) 'DIAGONAL 0'
-         end if
-      end do
-
       temp = qr_temp%get_vect()
       call psb_bcast(ctxt,res)
    else
diff --git a/krylov/psb_dbgmres.f90 b/krylov/psb_dbgmres.f90
index e55bd034..da2e7640 100644
--- a/krylov/psb_dbgmres.f90
+++ b/krylov/psb_dbgmres.f90
@@ -67,7 +67,6 @@ subroutine psb_dbgmres_multivect(a, prec, b, x, eps, desc_a, info, itmax, iter,
    type(psb_d_multivect_type)                :: w, r0, rm
 
    real(psb_dpk_)                            :: t1, t2
-   logical                                   :: solution = .false.
 
    real(psb_dpk_)                            :: rti, rti1
    integer(psb_ipk_)                         :: litmax, naux, itrace_, n_row, n_col, nrhs, nrep
@@ -225,8 +224,9 @@ subroutine psb_dbgmres_multivect(a, prec, b, x, eps, desc_a, info, itmax, iter,
 
    ! BGMRES algorithm
 
-   ! TODO Con tanti ITRS e tanti NRHS si ottengono NaN
-   ! TODO Deflazione e restart dopo aver trovato una colonna, difficile...
+   ! TODO Con tanti ITRS e tanti NRHS si ottengono NaN, deflazione e restart dopo aver trovato una colonna, difficile...
+
+   ! TODO Provare a compilare su GPU remota (Vedere REC)
 
    ! STEP 1: Compute R(0) = B - A*X(0)
 
@@ -369,23 +369,8 @@ subroutine psb_dbgmres_multivect(a, prec, b, x, eps, desc_a, info, itmax, iter,
          end do
       end if
 
-      ! Check convergence
-      if (maxval(errnum) <= eps*maxval(errden)) then
-
-         ! Compute result and exit
-         if (istop_ == 1) then
-
-            !Compute X(j) = X(0) + VT(j)*Y(j)
-            call psb_geaxpby(done,psb_geprod(vt(:,1:j*nrhs),y(1:j*nrhs,:),desc_a,info,global=.false.),done,x,desc_a,info)
-            if (info /= psb_success_) then
-               info=psb_err_from_subroutine_non_
-               call psb_errpush(info,name)
-               goto 9999
-            end if
-
-            solution = .true.
-
-         end if
+      ! TODO Check convergence (max o media)
+      if (all(errnum.le.(eps*errden))) then
 
          ! Exit algorithm
          exit outer
@@ -398,13 +383,11 @@ subroutine psb_dbgmres_multivect(a, prec, b, x, eps, desc_a, info, itmax, iter,
    end do outer
 
    ! STEP 10: X(m) = X(0) + VT(m)*Y(m)
-   if (.not.solution) then
-      call psb_geaxpby(done,psb_geprod(vt(:,1:nrep*nrhs),y,desc_a,info,global=.false.),done,x,desc_a,info)
-      if (info /= psb_success_) then
-         info=psb_err_from_subroutine_non_
-         call psb_errpush(info,name)
-         goto 9999
-      end if
+   call psb_geaxpby(done,psb_geprod(vt(:,1:j*nrhs),y(1:j*nrhs,:),desc_a,info,global=.false.),done,x,desc_a,info)
+   if (info /= psb_success_) then
+      info=psb_err_from_subroutine_non_
+      call psb_errpush(info,name)
+      goto 9999
    end if
 
    ! END algorithm
@@ -460,6 +443,8 @@ contains
       beta_e1 = dzero
       beta_e1(1:nrhs,1:nrhs) = beta
 
+      ! TODO DGELS ha anche i residui (con i residui fai come MATLAB e poi si prova se esce uguale)
+
       ! Compute min Frobenius norm
       call dgels('N',m_h,n_h,nrhs,h_temp(1:m_h,1:n_h),m_h,beta_e1,m_h,work,lwork,info)