Added TRANSPOSE paths to Block Jacobi and AS preconditioner application.

mlprec/mld_das_aply.f90

@@ -103,15 +103,6 @@ subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
 
   call psb_info(ictxt, me, np)
 
-  select case(trans)
-  case('N','n')
-  case('T','t','C','c')
-  case default
-    info=40
-    int_err(1)=6
-    ch_err(2:2)=trans
-    goto 9999
-  end select
 
   select case(prec%iprcparm(mld_prec_type_))
 
@@ -178,6 +169,9 @@ subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
       tx(1:nrow_d)     = x(1:nrow_d) 
       tx(nrow_d+1:isz) = dzero
 
+      select case(toupper(trans))
+      case('N')
+
         !
         ! Get the overlap entries of tx (tx==x)
         ! 
@@ -211,7 +205,7 @@ subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
         ! block-Jacobi solver can be applied at the coarsest level of a multilevel
         ! preconditioner). The resulting vector is ty.
         !
-      call mld_bjac_aply(done,prec,tx,dzero,ty,prec%desc_data,trans,aux,info)
+        call mld_bjac_aply(done,prec,tx,dzero,ty,prec%desc_data,toupper(trans),aux,info)
         if(info /= 0) then
           info=4010
           ch_err='mld_bjac_aply'
@@ -256,6 +250,118 @@ subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
           goto 9999
         end select
 
+      case('T','C')
+
+
+        !
+        ! With transpose, we have to do it here
+        ! 
+        
+        select case (prec%iprcparm(mld_sub_prol_)) 
+
+        case(psb_none_)
+          ! 
+          ! Do nothing
+
+        case(psb_sum_) 
+          !
+          ! Transpose of sum is halo
+          !
+          call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_halo'
+            goto 9999
+          end if
+
+        case(psb_avg_) 
+          !
+          ! Tricky one: first we have to scale the overlap entries,
+          ! which we can do by assignind mode=0, i.e. no communication
+          ! (hence only scaling), then we do the halo
+          !
+          call psb_ovrl(tx,prec%desc_data,info,&
+               & update=psb_avg_,work=aux,mode=0)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_ovrl'
+            goto 9999
+          end if
+          call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_halo'
+            goto 9999
+          end if
+
+        case default
+          call psb_errpush(4001,name,a_err='Invalid mld_sub_prol_')
+          goto 9999
+        end select
+
+
+        !
+        ! If required, reorder tx according to the row/column permutation of the
+        ! local extended matrix, stored into the permutation vector prec%perm
+        !
+        if (prec%iprcparm(mld_sub_ren_)>0) then 
+          call psb_gelp('n',prec%perm,tx,info)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_gelp'
+            goto 9999
+          end if
+        endif
+
+        !
+        ! Apply to tx the block-Jacobi preconditioner/solver (multiple sweeps of the
+        ! block-Jacobi solver can be applied at the coarsest level of a multilevel
+        ! preconditioner). The resulting vector is ty.
+        !
+        call mld_bjac_aply(done,prec,tx,dzero,ty,prec%desc_data,toupper(trans),aux,info)
+        if(info /= 0) then
+          info=4010
+          ch_err='mld_bjac_aply'
+          goto 9999
+        end if
+
+        !
+        ! Apply to ty the inverse permutation of prec%perm
+        !
+        if (prec%iprcparm(mld_sub_ren_)>0) then 
+          call psb_gelp('n',prec%invperm,ty,info)
+          if(info /= 0) then
+            info=4010
+            ch_err='psb_gelp'
+            goto 9999
+          end if
+        endif
+
+        !
+        ! With transpose, we have to do it here
+        ! 
+        if (prec%iprcparm(mld_sub_restr_) == psb_halo_) then 
+          call psb_ovrl(ty,prec%desc_data,info,&
+               & update=psb_sum_,work=aux)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_ovrl'
+            goto 9999
+          end if
+        else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then 
+          call psb_errpush(4001,name,a_err='Invalid mld_sub_restr_')
+          goto 9999
+        end if
+
+      case default
+        info=40
+        int_err(1)=6
+        ch_err(2:2)=trans
+        goto 9999
+      end select
+
+
+
       !
       ! Compute y = beta*y + alpha*ty (ty==K^(-1)*tx)
       !

mlprec/mld_dbjac_aply.f90

@@ -227,9 +227,9 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
 
       case('T','C')
         call psb_spsm(done,prec%av(mld_u_pr_),x,dzero,ww,desc_data,info,&
-             & trans=trans,unit='L',diag=prec%d,choice=psb_none_,work=aux)
+             & trans=toupper(trans),unit='L',diag=prec%d,choice=psb_none_,work=aux)
         if (info == 0) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,&
-             & trans=trans,unit='U',choice=psb_none_,work=aux)
+             & trans=toupper(trans),unit='U',choice=psb_none_,work=aux)
       case default
         call psb_errpush(4001,name,a_err='Invalid TRANS in ILU subsolve')
         goto 9999
@@ -338,6 +338,9 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
       !
       ! Use ILU(k)/MILU(k)/ILU(k,t) on the blocks.
       !
+
+      select case(toupper(trans))
+      case('N')
         do i=1, prec%iprcparm(mld_smooth_sweeps_) 
           !
           ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
@@ -357,6 +360,32 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
                & trans='N',unit='U',choice=psb_none_,work=aux)
           if (info /=0) exit
         end do
+      case('T','C')
+        do i=1, prec%iprcparm(mld_smooth_sweeps_) 
+          !
+          ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
+          ! block diagonal part and the remaining part of the local matrix
+          ! and Y(j) is the approximate solution at sweep j.
+          !
+          ty(1:n_row) = x(1:n_row)
+          call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,&
+               &   prec%desc_data,info,work=aux,trans=toupper(trans))
+          if (info /=0) exit
+          call psb_spsm(done,prec%av(mld_u_pr_),ty,dzero,ww,&
+               & prec%desc_data,info,&
+               & trans=toupper(trans),unit='L',diag=prec%d,choice=psb_none_,work=aux)
+          if (info /=0) exit
+          call psb_spsm(done,prec%av(mld_l_pr_),ww,dzero,tx,&
+               & prec%desc_data,info,&
+               & trans=toupper(trans),unit='U',choice=psb_none_,work=aux)
+          if (info /=0) exit
+        end do
+
+      case default
+        call psb_errpush(4001,name,a_err='Invalid TRANS in ILU subsolve')
+        goto 9999
+      end select
+
 
     case(mld_sludist_) 
       !
@@ -371,6 +400,8 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
       ! Use the LU factorization from SuperLU.
       !
 
+      select case(toupper(trans))
+      case('N')
         do i=1, prec%iprcparm(mld_smooth_sweeps_)
           ! 
           ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
@@ -386,12 +417,35 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
           if(info /= 0) exit 
           tx(1:n_row) = ty(1:n_row)        
         end do
+      case('T','C')
+        do i=1, prec%iprcparm(mld_smooth_sweeps_)
+          ! 
+          ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
+          ! block diagonal part and the remaining part of the local matrix
+          ! and Y(j) is the approximate solution at sweep j.
+          !
+          ty(1:n_row) = x(1:n_row)
+          call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,&
+               &   prec%desc_data,info,work=aux,trans=toupper(trans))
+          if(info /= 0) exit
+          call mld_dslu_solve(1,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info)
+          if(info /= 0) exit 
+          tx(1:n_row) = ty(1:n_row)        
+        end do
+
+      case default
+        call psb_errpush(4001,name,a_err='Invalid TRANS in SLU subsolve')
+        goto 9999
+      end select
+
 
     case(mld_umf_)
       !
       ! Use the LU factorization from UMFPACK.
       !
 
+      select case(toupper(trans))
+      case('N')
         do i=1, prec%iprcparm(mld_smooth_sweeps_) 
           ! 
           ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
@@ -408,6 +462,28 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
           if (info /= 0) exit
           tx(1:n_row) = ww(1:n_row)        
         end do
+      case('T','C')
+        do i=1, prec%iprcparm(mld_smooth_sweeps_) 
+          ! 
+          ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
+          ! block diagonal part and the remaining part of the local matrix
+          ! and Y(j) is the approximate solution at sweep j.
+          !
+          ty(1:n_row) = x(1:n_row)
+          call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,&
+               &   prec%desc_data,info,work=aux,trans=toupper(trans))
+          if (info /= 0) exit
+
+          call mld_dumf_solve(1,n_row,ww,ty,n_row,&
+               & prec%iprcparm(mld_umf_numptr_),info)
+          if (info /= 0) exit
+          tx(1:n_row) = ww(1:n_row)        
+        end do
+
+      case default
+        call psb_errpush(4001,name,a_err='Invalid TRANS in UMF subsolve')
+        goto 9999
+      end select
 
     case default
       call psb_errpush(4001,name,a_err='Invalid mld_sub_solve_')

mlprec/mld_zas_aply.f90

@@ -103,15 +103,6 @@ subroutine mld_zas_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
 
   call psb_info(ictxt, me, np)
 
-  select case(trans)
-  case('N','n')
-  case('T','t','C','c')
-  case default
-    info=40
-    int_err(1)=6
-    ch_err(2:2)=trans
-    goto 9999
-  end select
 
   select case(prec%iprcparm(mld_prec_type_))
 
@@ -178,6 +169,8 @@ subroutine mld_zas_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
       tx(1:nrow_d)     = x(1:nrow_d) 
       tx(nrow_d+1:isz) = dzero
 
+      select case(toupper(trans))
+      case('N')
         !
         ! Get the overlap entries of tx (tx==x)
         ! 
@@ -211,7 +204,7 @@ subroutine mld_zas_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
         ! block-Jacobi solver can be applied at the coarsest level of a multilevel
         ! preconditioner). The resulting vector is ty.
         !
-      call mld_bjac_aply(zone,prec,tx,zzero,ty,prec%desc_data,trans,aux,info)
+        call mld_bjac_aply(zone,prec,tx,zzero,ty,prec%desc_data,toupper(trans),aux,info)
         if(info /= 0) then
           info=4010
           ch_err='mld_bjac_aply'
@@ -256,6 +249,115 @@ subroutine mld_zas_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
           goto 9999
         end select
 
+      case('T','C')
+
+        !
+        ! With transpose, we have to do it here
+        ! 
+
+        select case (prec%iprcparm(mld_sub_prol_)) 
+
+        case(psb_none_)
+          ! 
+          ! Do nothing
+
+        case(psb_sum_) 
+          !
+          ! Transpose of sum is halo
+          !
+          call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_halo'
+            goto 9999
+          end if
+
+        case(psb_avg_) 
+          !
+          ! Tricky one: first we have to scale the overlap entries,
+          ! which we can do by assignind mode=0, i.e. no communication
+          ! (hence only scaling), then we do the halo
+          !
+          call psb_ovrl(tx,prec%desc_data,info,&
+               & update=psb_avg_,work=aux,mode=0)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_ovrl'
+            goto 9999
+          end if
+          call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_halo'
+            goto 9999
+          end if
+
+        case default
+          call psb_errpush(4001,name,a_err='Invalid mld_sub_prol_')
+          goto 9999
+        end select
+
+
+        !
+        ! If required, reorder tx according to the row/column permutation of the
+        ! local extended matrix, stored into the permutation vector prec%perm
+        !
+        if (prec%iprcparm(mld_sub_ren_)>0) then 
+          call psb_gelp('n',prec%perm,tx,info)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_gelp'
+            goto 9999
+          end if
+        endif
+
+        !
+        ! Apply to tx the block-Jacobi preconditioner/solver (multiple sweeps of the
+        ! block-Jacobi solver can be applied at the coarsest level of a multilevel
+        ! preconditioner). The resulting vector is ty.
+        !
+        call mld_bjac_aply(zone,prec,tx,zzero,ty,prec%desc_data,toupper(trans),aux,info)
+        if(info /= 0) then
+          info=4010
+          ch_err='mld_bjac_aply'
+          goto 9999
+        end if
+
+        !
+        ! Apply to ty the inverse permutation of prec%perm
+        !
+        if (prec%iprcparm(mld_sub_ren_)>0) then 
+          call psb_gelp('n',prec%invperm,ty,info)
+          if(info /= 0) then
+            info=4010
+            ch_err='psb_gelp'
+            goto 9999
+          end if
+        endif
+
+        !
+        ! With transpose, we have to do it here
+        ! 
+        if (prec%iprcparm(mld_sub_restr_) == psb_halo_) then 
+          call psb_ovrl(ty,prec%desc_data,info,&
+               & update=psb_sum_,work=aux)
+          if(info /=0) then
+            info=4010
+            ch_err='psb_ovrl'
+            goto 9999
+          end if
+        else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then 
+          call psb_errpush(4001,name,a_err='Invalid mld_sub_restr_')
+          goto 9999
+        end if
+
+      case default
+        info=40
+        int_err(1)=6
+        ch_err(2:2)=trans
+        goto 9999
+      end select
+
       !
       ! Compute y = beta*y + alpha*ty (ty==K^(-1)*tx)
       !

mlprec/mld_zbjac_aply.f90

@@ -339,11 +339,15 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
 
     tx = zzero
     ty = zzero
+
     select case(prec%iprcparm(mld_sub_solve_)) 
     case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) 
       !
       ! Use ILU(k)/MILU(k)/ILU(k,t) on the blocks.
       !
+      select case(toupper(trans))
+      case('N')
+
         do i=1, prec%iprcparm(mld_smooth_sweeps_) 
           !
           ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
@@ -363,6 +367,31 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
                & trans='N',unit='U',choice=psb_none_,work=aux)
           if (info /=0) exit
         end do
+      case('T','C')
+        do i=1, prec%iprcparm(mld_smooth_sweeps_) 
+          !
+          ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
+          ! block diagonal part and the remaining part of the local matrix
+          ! and Y(j) is the approximate solution at sweep j.
+          !
+          ty(1:n_row) = x(1:n_row)
+          call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
+               &   prec%desc_data,info,work=aux,trans=toupper(trans))
+          if (info /=0) exit
+          call psb_spsm(zone,prec%av(mld_u_pr_),ty,zzero,ww,&
+               & prec%desc_data,info,&
+               & trans=toupper(trans),unit='L',diag=prec%d,choice=psb_none_,work=aux)
+          if (info /=0) exit
+          call psb_spsm(zone,prec%av(mld_l_pr_),ww,zzero,tx,&
+               & prec%desc_data,info,&
+               & trans=toupper(trans),unit='U',choice=psb_none_,work=aux)
+          if (info /=0) exit
+        end do
+      case default
+        call psb_errpush(4001,name,a_err='Invalid TRANS in ILU subsolve')
+        goto 9999
+      end select
+
 
     case(mld_sludist_) 
       !
@@ -377,6 +406,8 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
       ! Use the LU factorization from SuperLU.
       !
 
+      select case(toupper(trans))
+      case('N')
         do i=1, prec%iprcparm(mld_smooth_sweeps_)
           ! 
           ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
@@ -392,12 +423,51 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
           if (info /= 0) exit
           tx(1:n_row) = ty(1:n_row)        
         end do
+      case('T')
+        do i=1, prec%iprcparm(mld_smooth_sweeps_)
+          ! 
+          ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
+          ! block diagonal part and the remaining part of the local matrix
+          ! and Y(j) is the approximate solution at sweep j.
+          !
+          ty(1:n_row) = x(1:n_row)
+          call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
+               &   prec%desc_data,info,work=aux,trans=toupper(trans))
+          if (info /= 0) exit
+
+          call mld_zslu_solve(1,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info)
+          if (info /= 0) exit
+          tx(1:n_row) = ty(1:n_row)        
+        end do
+      case('C')
+        do i=1, prec%iprcparm(mld_smooth_sweeps_)
+          ! 
+          ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
+          ! block diagonal part and the remaining part of the local matrix
+          ! and Y(j) is the approximate solution at sweep j.
+          !
+          ty(1:n_row) = x(1:n_row)
+          call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
+               &   prec%desc_data,info,work=aux,trans=toupper(trans))
+          if (info /= 0) exit
+
+          call mld_zslu_solve(2,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info)
+          if (info /= 0) exit
+          tx(1:n_row) = ty(1:n_row)        
+        end do
+
+      case default
+        call psb_errpush(4001,name,a_err='Invalid TRANS in SLU subsolve')
+        goto 9999
+      end select
 
     case(mld_umf_)
       !
       ! Use the LU factorization from UMFPACK.
       !
 
+      select case(toupper(trans))
+      case('N')
         do i=1, prec%iprcparm(mld_smooth_sweeps_) 
           ! 
           ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
@@ -414,6 +484,45 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
           if (info /= 0) exit
           tx(1:n_row) = ww(1:n_row)        
         end do
+      case('T')
+        do i=1, prec%iprcparm(mld_smooth_sweeps_) 
+          ! 
+          ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
+          ! block diagonal part and the remaining part of the local matrix
+          ! and Y(j) is the approximate solution at sweep j.
+          !
+          ty(1:n_row) = x(1:n_row)
+          call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
+               &   prec%desc_data,info,work=aux,trans=toupper(trans))
+          if (info /= 0) exit
+          
+          call mld_zumf_solve(1,n_row,ww,ty,n_row,&
+               & prec%iprcparm(mld_umf_numptr_),info)
+          if (info /= 0) exit
+          tx(1:n_row) = ww(1:n_row)        
+        end do
+      case('C')
+        do i=1, prec%iprcparm(mld_smooth_sweeps_) 
+          ! 
+          ! Compute Y(k+1) = D^(-1)*(X-ND*Y(k)), where D and ND are the
+          ! block diagonal part and the remaining part of the local matrix
+          ! and Y(j) is the approximate solution at sweep j.
+          !
+          ty(1:n_row) = x(1:n_row)
+          call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
+               &   prec%desc_data,info,work=aux,trans=toupper(trans))
+          if (info /= 0) exit
+          
+          call mld_zumf_solve(2,n_row,ww,ty,n_row,&
+               & prec%iprcparm(mld_umf_numptr_),info)
+          if (info /= 0) exit
+          tx(1:n_row) = ww(1:n_row)        
+        end do
+
+      case default
+        call psb_errpush(4001,name,a_err='Invalid TRANS in UMF subsolve')
+        goto 9999
+      end select
 
     case default
       call psb_errpush(4001,name,a_err='Invalid mld_sub_solve_')