Merged changes for new precs.

src/modules/psb_prec_type.f90

@@ -11,7 +11,8 @@ module psb_prec_type
   integer, parameter :: min_prec_=0, noprec_=0, diagsc_=1, bja_=2,&
        & asm_=3, ras_=5, ash_=4, rash_=6, ras2lv_=7, ras2lvm_=8,&
        & lv2mras_=9, lv2smth_=10, lv2lsm_=11, sl2sm_=12, superlu_=13,&
-       & new_loc_smth_=14, new_glb_smth_=15, max_prec_=15
+       & new_loc_smth_=14, new_glb_smth_=15, ag2lsm_=16,&
+       & msy2l_=18, msy2g_=19, max_prec_=19
   integer, parameter   :: nohalo_=0,                 halo_=4
   integer, parameter   :: none_=0,                   sum_=1
   integer, parameter   :: avg_=2,                    square_root_=3
@@ -35,12 +36,13 @@ module psb_prec_type
   integer, parameter :: ilu_fill_in_=8, jac_sweeps_=9, ml_type_=10
   integer, parameter :: smth_pos_=11, aggr_alg_=12, smth_kind_=13
   integer, parameter :: om_choice_=14, glb_smth_=15, coarse_mat_=16
+  integer, parameter :: umf_symptr_=17, umf_numptr_=18
   integer, parameter :: ifpsz=20
   ! Entries in dprcparm: ILU(E) epsilon, smoother omega
   integer, parameter :: fact_eps_=1, smooth_omega_=2
   integer, parameter :: dfpsz=4
   ! Factorization types: none, ILU(N), ILU(E), SuperLU
-  integer, parameter :: f_none_=0,f_ilu_n_=1, f_ilu_e_=2,f_slu_=3
+  integer, parameter :: f_none_=0,f_ilu_n_=1,f_ilu_e_=2,f_slu_=3,f_umf_=4
   ! Fields for sparse matrices ensembles: 
   integer, parameter :: l_pr_=1, u_pr_=2, bp_ilu_avsz=2
   integer, parameter :: ap_nd_=3, ac_=4, sm_pr_t_=5, sm_pr_=6
@@ -69,7 +71,8 @@ module psb_prec_type
 
   
   character(len=15), parameter, private :: &
-       &  smooth_names(1:2)=(/'Pre-smoothing ','Post-smoothing'/)
+       &  smooth_names(1:3)=(/'Pre-smoothing ','Post-smoothing',&
+       & 'Smooth both   '/)
   character(len=15), parameter, private :: &
        &  smooth_kinds(0:2)=(/'No  smoother  ','Omega smoother',&
        &           'Bizr. smoother'/)
@@ -86,8 +89,8 @@ module psb_prec_type
        &  ml_names(0:3)=(/'None          ','Additive      ','Multiplicative',&
        & 'New ML        '/)
   character(len=15), parameter, private :: &
-       &  fact_names(0:3)=(/'None          ','ILU(n)        ',&
-       &  'ILU(eps)      ','Sparse LU     '/)
+       &  fact_names(0:4)=(/'None          ','ILU(n)        ',&
+       &  'ILU(eps)      ','Sparse SuperLU','UMFPACK Sp. LU'/)
 
   interface psb_base_precfree
     module procedure psb_dbase_precfree
@@ -359,6 +362,10 @@ contains
       if (p%iprcparm(f_type_)==f_slu_) then 
         call fort_slu_free(p%iprcparm(slu_ptr_),info)
       end if
+      if (p%iprcparm(f_type_)==f_umf_) then 
+        call fort_umf_free(p%iprcparm(umf_symptr_),&
+             & p%iprcparm(umf_numptr_),info)
+      end if
       deallocate(p%iprcparm,stat=info)
     end if
     call psb_nullify_baseprec(p)

src/modules/psb_serial_mod.f90

@@ -62,10 +62,11 @@ module psb_serial_mod
   end interface
 
   interface psb_fixcoo
-     subroutine psb_dfixcoo(a,info)
+     subroutine psb_dfixcoo(a,info,idir)
        use psb_spmat_type
        type(psb_dspmat_type), intent(inout) :: a
        integer, intent(out)         :: info
+       integer, intent(in), optional :: idir
      end subroutine psb_dfixcoo
   end interface
 
@@ -78,6 +79,15 @@ module psb_serial_mod
      end subroutine psb_dipcoo2csr
   end interface
 
+  interface psb_ipcoo2csc
+     subroutine psb_dipcoo2csc(a,info,clshr)
+       use psb_spmat_type
+       type(psb_dspmat_type), intent(inout) :: a
+       integer, intent(out)         :: info
+       logical, optional :: clshr
+     end subroutine psb_dipcoo2csc
+  end interface
+
   interface psb_ipcsr2coo
      subroutine psb_dipcsr2coo(a,info)
        use psb_spmat_type

src/prec/Makefile

@@ -9,7 +9,7 @@ F90OBJS= psb_dcsrsetup.o psb_dprec.o \
 	psb_dgenaggrmap.o  psb_dsplu.o $(MPFOBJS) 
 #dcoocp.o dcoocpadd.o dcoofact.o dcoolu.o dcooluadd.o\
 
-COBJS=fort_slu_impl.o
+COBJS=fort_slu_impl.o fort_umf_impl.o
 INCDIRS=-I. -I.. -I$(LIBDIR) 
 
 OBJS=$(F90OBJS) $(COBJS)

src/prec/psb_dprec.f90

@@ -255,7 +255,8 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info)
     end if
 
     if (prec%iprcparm(iren_)>0) then 
-      call psb_dgelp('N',n_row,1,prec%perm,tx,isz,ww,isz,info)
+!!$      call psb_dgelp('N',n_row,1,prec%perm,tx,isz,ww,isz,info)
+      info = -1
       if(info /=0) then
          info=4010
          ch_err='psb_dgelp'
@@ -266,7 +267,8 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info)
     call psb_dbjacaply(prec,tx,zero,ty,prec%desc_data,trans,aux,info)
 
     if (prec%iprcparm(iren_)>0) then 
-      call psb_dgelp('N',n_row,1,prec%invperm,ty,isz,ww,isz,info)
+!!$      call psb_dgelp('N',n_row,1,prec%invperm,ty,isz,ww,isz,info)
+      info = -1 
       if(info /=0) then
          info=4010
          ch_err='psb_dgelp'
@@ -449,6 +451,27 @@ subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info)
 
       if(info /=0) goto 9999
 
+      if (beta == 0.d0) then 
+        y(1:n_row) = ww(1:n_row)
+      else if (beta==1.d0) then 
+        y(1:n_row) = ww(1:n_row) + y(1:n_row) 
+      else if (beta==-1.d0) then 
+        y(1:n_row) = ww(1:n_row) - y(1:n_row) 
+      else 
+        y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) 
+      endif
+    case (f_umf_) 
+
+
+      select case(trans)
+      case('N','n')
+        call fort_umf_solve(0,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info)
+      case('T','t','C','c')
+        call fort_umf_solve(1,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info)
+      end select
+
+      if(info /=0) goto 9999
+
       if (beta == 0.d0) then 
         y(1:n_row) = ww(1:n_row)
       else if (beta==1.d0) then 
@@ -507,6 +530,20 @@ subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info)
         if(info /=0) goto 9999
         tx(1:n_row) = ty(1:n_row)        
       end do
+    case(f_umf_) 
+      do i=1, prec%iprcparm(jac_sweeps_) 
+        !   X(k+1) = M^-1*(b-N*X(k))
+        ty(1:n_row) = x(1:n_row)
+        call psb_spmm(-one,prec%av(ap_nd_),tx,one,ty,&
+             &   prec%desc_data,info,work=aux)
+        if(info /=0) goto 9999
+
+        call fort_umf_solve(0,n_row,ww,ty,n_row,&
+             & prec%iprcparm(umf_numptr_),info)
+        if(info /=0) goto 9999
+        tx(1:n_row) = ww(1:n_row)        
+      end do
+
     end select
     
     if (beta == 0.d0) then 

src/prec/psb_dprecbld.f90

@@ -29,6 +29,35 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd)
      end subroutine psb_dcslu
   end interface
 
+  interface 
+    subroutine psb_splu_bld(a,desc_a,p,info)
+      use psb_serial_mod
+      use psb_descriptor_type
+      use psb_prec_type
+      use psb_const_mod
+      implicit none 
+      
+      type(psb_dspmat_type), intent(in)   :: a
+      type(psb_desc_type), intent(in)     :: desc_a
+      type(psb_dbase_prec), intent(inout) :: p
+      integer, intent(out)                :: info
+    end subroutine psb_splu_bld
+  end interface
+  interface 
+    subroutine psb_umf_bld(a,desc_a,p,info)
+      use psb_serial_mod
+      use psb_descriptor_type
+      use psb_prec_type
+      use psb_const_mod
+      implicit none 
+      
+      type(psb_dspmat_type), intent(in)   :: a
+      type(psb_desc_type), intent(in)     :: desc_a
+      type(psb_dbase_prec), intent(inout) :: p
+      integer, intent(out)                :: info
+    end subroutine psb_umf_bld
+  end interface
+
   ! Local scalars
   Integer      :: err, nnzero, n_row, n_col,I,j,icontxt,&
        & me,mycol,nprow,npcol,mglob,lw, mtype, nrg, nzg, err_act
@@ -198,6 +227,17 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd)
            goto 9999
         end if
 
+     case(f_umf_)
+        p%baseprecv(1)%av => null()
+        if(debug) write(0,*)me,': calling umf_bld'
+        call psb_umf_bld(a,desc_a,p%baseprecv(1),info)
+        if(info /= 0) then
+           info=4010
+           ch_err='umf_bld'
+           call psb_errpush(info,name,a_err=ch_err)
+           goto 9999
+        end if
+
      case(f_none_) 
         write(0,*) 'Fact=None in PRECBLD Bja/ASM??'
 
@@ -416,9 +456,7 @@ subroutine psb_splu_bld(a,desc_a,p,info)
      write(0,*) me, 'SPLUBLD: Done slu_Factor',info,p%iprcparm(slu_ptr_)
      call blacs_barrier(icontxt,'All')
   endif
-  if (info /= 0) then 
-     write(0,*) 'From fort_slu_factor: ',info
-  endif
+
   call psb_spfree(blck,info)
   call psb_spfree(atmp,info)
   if(info /= 0) then
@@ -443,6 +481,184 @@ end subroutine psb_splu_bld
 
 
 
+subroutine psb_umf_bld(a,desc_a,p,info)
+  use psb_serial_mod
+  use psb_descriptor_type
+  use psb_prec_type
+  use psb_const_mod
+  implicit none 
+
+  type(psb_dspmat_type), intent(in)   :: a
+  type(psb_desc_type), intent(in)     :: desc_a
+  type(psb_dbase_prec), intent(inout) :: p
+  integer, intent(out)                :: info
+
+
+  type(psb_dspmat_type)    :: blck, atmp
+  character(len=5)         :: fmt
+  character                :: upd='F'
+  integer                  :: i,j,nza,nzb,nzt,icontxt, me,mycol,nprow,npcol,err_act
+  logical, parameter :: debug=.false.
+  character(len=20)   :: name, ch_err
+
+  interface psb_csrsetup
+    Subroutine psb_dcsrsetup(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
+      use psb_serial_mod
+      Use psb_descriptor_type
+      Use psb_prec_type
+      integer, intent(in)                  :: ptype,novr
+      Type(psb_dspmat_type), Intent(in)    ::  a
+      Type(psb_dspmat_type), Intent(inout) ::  blk
+      Type(psb_desc_type), Intent(inout)   :: desc_p
+      Type(psb_desc_type), Intent(in)      :: desc_data 
+      Character, Intent(in)                :: upd
+      integer, intent(out)                 :: info
+      character(len=5), optional           :: outfmt
+    end Subroutine psb_dcsrsetup
+ end interface
+
+  info=0
+  name='psb_umf_bld'
+  call psb_erractionsave(err_act)
+
+  icontxt = desc_A%matrix_data(psb_ctxt_)
+  call blacs_gridinfo(icontxt, nprow, npcol, me, mycol)
+
+
+  fmt = 'COO'
+  call psb_nullify_sp(blck)    
+  call psb_nullify_sp(atmp)    
+  atmp%fida='COO'
+  if (Debug) then 
+     write(0,*) me, 'UMFBLD: Calling  csdp'
+     call blacs_barrier(icontxt,'All')
+  endif
+
+  call psb_dcsdp(a,atmp,info)
+  if(info /= 0) then
+     info=4010
+     ch_err='psb_dcsdp'
+     call psb_errpush(info,name,a_err=ch_err)
+     goto 9999
+  end if
+  nza = atmp%infoa(psb_nnz_)
+  if (Debug) then 
+     write(0,*) me, 'UMFBLD: Done csdp',info,nza,atmp%m,atmp%k
+     call blacs_barrier(icontxt,'All')
+  endif
+  call psb_csrsetup(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,&
+       & blck,desc_a,upd,p%desc_data,info,outfmt=fmt)  
+  if(info /= 0) then
+     info=4010
+     ch_err='psb_csrsetup'
+     call psb_errpush(info,name,a_err=ch_err)
+     goto 9999
+  end if
+
+  nzb = blck%infoa(psb_nnz_)
+  if (Debug) then 
+     write(0,*) me, 'UMFBLD: Done csrsetup',info,nzb,blck%fida
+     call blacs_barrier(icontxt,'All')
+  endif
+  if (nzb > 0 ) then 
+     if (size(atmp%aspk)<nza+nzb) then 
+        call psb_spreall(atmp,nza+nzb,info)
+        if(info /= 0) then
+           info=4010
+           ch_err='psb_spreall'
+           call psb_errpush(info,name,a_err=ch_err)
+           goto 9999
+        end if
+     endif
+     if (Debug) then 
+        write(0,*) me, 'UMFBLD: Done realloc',info,nza+nzb,atmp%fida
+        call blacs_barrier(icontxt,'All')
+     endif
+
+     do j=1,nzb
+        atmp%aspk(nza+j) = blck%aspk(j)
+        atmp%ia1(nza+j)  = blck%ia1(j)
+        atmp%ia2(nza+j)  = blck%ia2(j)
+     end do
+     atmp%infoa(psb_nnz_) = nza+nzb
+     atmp%m = atmp%m + blck%m
+     atmp%k = max(a%k,blck%k)
+  else
+     atmp%infoa(psb_nnz_) = nza
+     atmp%m = a%m 
+     atmp%k = a%k
+  endif
+
+  i=0
+  do j=1, atmp%infoa(psb_nnz_) 
+     if (atmp%ia2(j) <= atmp%m) then 
+        i = i + 1
+        atmp%aspk(i) = atmp%aspk(j)
+        atmp%ia1(i) = atmp%ia1(j)
+        atmp%ia2(i) = atmp%ia2(j)
+     endif
+  enddo
+  atmp%infoa(psb_nnz_) = i
+
+
+  call psb_ipcoo2csc(atmp,info)
+  if(info /= 0) then
+     info=4010
+     ch_err='psb_ipcoo2csc'
+     call psb_errpush(info,name,a_err=ch_err)
+     goto 9999
+  end if
+  call psb_spinfo(psb_nztotreq_,atmp,nzt,info)
+  if(info /= 0) then
+     info=4010
+     ch_err='psb_spinfo'
+     call psb_errpush(info,name,a_err=ch_err)
+     goto 9999
+  end if
+  if (Debug) then 
+     write(0,*) me,'Calling fort_slu_factor ',nzt,atmp%m,&
+          & atmp%k,p%desc_data%matrix_data(psb_n_row_)
+     call blacs_barrier(icontxt,'All')
+  endif
+
+  call fort_umf_factor(atmp%m,nzt,&
+       & atmp%aspk,atmp%ia2,atmp%ia1,&
+       & p%iprcparm(umf_symptr_),p%iprcparm(umf_numptr_),info)
+  if(info /= 0) then
+     info=4010
+     ch_err='fort_umf_fact'
+     call psb_errpush(info,name,a_err=ch_err)
+     goto 9999
+  end if
+
+  if (Debug) then 
+     write(0,*) me, 'UMFBLD: Done umf_Factor',info,p%iprcparm(umf_numptr_)
+     call blacs_barrier(icontxt,'All')
+  endif
+  call psb_spfree(blck,info)
+  call psb_spfree(atmp,info)
+  if(info /= 0) then
+     info=4010
+     ch_err='psb_spfree'
+     call psb_errpush(info,name,a_err=ch_err)
+     goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_erractionrestore(err_act)
+  if (err_act.eq.act_abort) then
+     call psb_error()
+     return
+  end if
+  return
+
+end subroutine psb_umf_bld
+
+
+
 
 subroutine psb_mlprec_bld(a,desc_a,p,info)
 

src/serial/Makefile

@@ -6,7 +6,7 @@ FOBJS = psb_cest.o psb_dcoins.o psb_dcsdp.o psb_dcsmm.o psb_dcsmv.o \
 	psb_dfixcoo.o psb_dipcoo2csr.o psb_dipcsr2coo.o psb_dneigh.o \
 	psb_dnumbmm.o psb_drwextd.o psb_dspgtdiag.o psb_dspgtrow.o \
 	psb_dspinfo.o psb_dspscal.o psb_dsymbmm.o psb_dtransp.o \
-	string_impl.o
+	string_impl.o psb_dipcoo2csc.o 
 
 
 INCDIRS = -I ../../lib -I .

src/serial/psb_dfixcoo.f90

@@ -2,7 +2,7 @@
 ! Subroutine: 
 ! Parameters:
 
-Subroutine psb_dfixcoo(A,INFO)
+Subroutine psb_dfixcoo(A,INFO,idir)
   use psb_spmat_type
   use psb_const_mod
   implicit none
@@ -10,10 +10,11 @@ Subroutine psb_dfixcoo(A,INFO)
   !....Parameters...
   Type(psb_dspmat_type), intent(inout) :: A
   Integer, intent(out)                 :: info
+  integer, intent(in), optional :: idir
 
   integer, allocatable :: iaux(:)
   !locals
-  Integer              :: nza, nzl,iret
+  Integer              :: nza, nzl,iret,idir_
   integer              :: i,j, irw, icl
   logical, parameter   :: debug=.false.
 
@@ -24,51 +25,103 @@ Subroutine psb_dfixcoo(A,INFO)
     info = -1
     return
   end if
+  if (present(idir)) then 
+    idir_ = idir
+  else
+    idir_ = 0
+  endif
 
   nza = a%infoa(psb_nnz_)
   if (nza < 2) return
-  
+
   allocate(iaux(nza+2),stat=info) 
   if (info /= 0) return
 
-  call mrgsrt(nza,a%ia1,iaux,iret)
-  if (iret.eq.0) call reordvn(nza,a%aspk,a%ia1,a%ia2,iaux)
-  i    = 1
-  j    = i
-  do while (i.le.nza)
-    do while ((a%ia1(j).eq.a%ia1(i)))
-      j = j+1
+  select case(idir_) 
+
+  case(0) !  Row major order
+
+    call mrgsrt(nza,a%ia1,iaux,iret)
+    if (iret.eq.0) call reordvn(nza,a%aspk,a%ia1,a%ia2,iaux)
+    i    = 1
+    j    = i
+    do while (i.le.nza)
+      do while ((a%ia1(j).eq.a%ia1(i)))
+        j = j+1
+        if (j > nza) exit
+      enddo
+      nzl = j - i
+      call mrgsrt(nzl,a%ia2(i),iaux,iret)
+      if (iret.eq.0) &
+           & call reordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux)
+      i = j
+    enddo
+
+    i = 1
+    irw = a%ia1(i)
+    icl = a%ia2(i)
+    j = 1
+    do 
+      j = j + 1
       if (j > nza) exit
+      if ((a%ia1(j) == irw).and.(a%ia2(j) == icl)) then 
+        a%aspk(i) = a%aspk(i) + a%aspk(j)
+      else
+        i = i+1
+        a%aspk(i) = a%aspk(j)
+        a%ia1(i) = a%ia1(j)
+        a%ia2(i) = a%ia2(j)
+        irw = a%ia1(i) 
+        icl = a%ia2(i) 
+      endif
     enddo
-    nzl = j - i
-    call mrgsrt(nzl,a%ia2(i),iaux,iret)
-    if (iret.eq.0) &
-         & call reordvn(nzl,a%aspk(i),a%ia1(i),a%ia2(i),iaux)
-    i = j
-  enddo
+    a%infoa(psb_nnz_) = i    
+    a%infoa(psb_srtd_) = psb_isrtdcoo_
 
-  i = 1
-  irw = a%ia1(i)
-  icl = a%ia2(i)
-  j = 1
-  do 
-    j = j + 1
-    if (j > nza) exit
-    if ((a%ia1(j) == irw).and.(a%ia2(j) == icl)) then 
-      a%aspk(i) = a%aspk(i) + a%aspk(j)
-    else
-      i = i+1
-      a%aspk(i) = a%aspk(j)
-      a%ia1(i) = a%ia1(j)
-      a%ia2(i) = a%ia2(j)
-      irw = a%ia1(i) 
-      icl = a%ia2(i) 
-    endif
-  enddo
-  a%infoa(psb_nnz_) = i    
-  a%infoa(psb_srtd_) = psb_isrtdcoo_
+    if(debug) write(0,*)'FIXCOO: end second loop'
+  case(1) !  Col major order
+
+    call mrgsrt(nza,a%ia2,iaux,iret)
+    if (iret.eq.0) call reordvn(nza,a%aspk,a%ia1,a%ia2,iaux)
+    i    = 1
+    j    = i
+    do while (i.le.nza)
+      do while ((a%ia2(j).eq.a%ia2(i)))
+        j = j+1
+        if (j > nza) exit
+      enddo
+      nzl = j - i
+      call mrgsrt(nzl,a%ia1(i:i+nzl-1),iaux,iret)
+      if (iret.eq.0) &
+           & call reordvn(nzl,a%aspk(i:i+nzl-1),a%ia1(i:i+nzl-1),a%ia2(i:i+nzl-1),iaux)
+      i = j
+    enddo
+
+    i = 1
+    irw = a%ia1(i)
+    icl = a%ia2(i)
+    j = 1
+    do 
+      j = j + 1
+      if (j > nza) exit
+      if ((a%ia1(j) == irw).and.(a%ia2(j) == icl)) then 
+        a%aspk(i) = a%aspk(i) + a%aspk(j)
+      else
+        i = i+1
+        a%aspk(i) = a%aspk(j)
+        a%ia1(i) = a%ia1(j)
+        a%ia2(i) = a%ia2(j)
+        irw = a%ia1(i) 
+        icl = a%ia2(i) 
+      endif
+    enddo
+    a%infoa(psb_nnz_) = i    
+    a%infoa(psb_srtd_) = psb_isrtdcoo_
 
-  if(debug) write(0,*)'FIXCOO: end second loop'
+    if(debug) write(0,*)'FIXCOO: end second loop'
+  case default
+    write(0,*) 'Fixcoo: unknown direction ',idir_
+  end select
 
   deallocate(iaux)
   return