From 51d966a10929ef3374ccca319f73a25f53c0e8d3 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 7 Mar 2006 15:10:34 +0000 Subject: [PATCH] Second tranche of renaming, especially in the preconditioner section. Repackaging the preconditioners in a more sensible way. WARNING: this is not complete yet, but we need a checkpoint to fall back to if needed. --- Makefile | 5 +- src/modules/psb_prec_mod.f90 | 6 +- src/modules/psb_prec_type.f90 | 11 +- src/prec/Makefile | 8 +- .../{psb_dcsrsetup.f90 => psb_dasmatbld.f90} | 12 +- src/prec/psb_dcslu.f90 | 761 ++++---------- src/prec/{psb_dsplu.f90 => psb_dilu_bld.f90} | 14 +- src/prec/psb_dmlprc_bld.f90 | 241 +++++ src/prec/psb_dprecbld.f90 | 978 +++--------------- src/prec/psb_dslu_bld.f90 | 213 ++++ src/prec/psb_dsp_renum.f90 | 467 +++++++++ src/prec/psb_dumf_bld.f90 | 215 ++++ 12 files changed, 1513 insertions(+), 1418 deletions(-) rename src/prec/{psb_dcsrsetup.f90 => psb_dasmatbld.f90} (96%) rename src/prec/{psb_dsplu.f90 => psb_dilu_bld.f90} (98%) create mode 100644 src/prec/psb_dmlprc_bld.f90 create mode 100644 src/prec/psb_dslu_bld.f90 create mode 100644 src/prec/psb_dsp_renum.f90 create mode 100644 src/prec/psb_dumf_bld.f90 diff --git a/Makefile b/Makefile index 0b22b253..e5f7ed84 100644 --- a/Makefile +++ b/Makefile @@ -10,8 +10,9 @@ lib: clean: (cd src; make clean) -veryclean: - (cd src; make veryclean) +cleanlib: (cd lib; /bin/rm -f *.a *$(.mod) *$(.fh)) +veryclean: cleanlib + (cd src; make veryclean) .PHONY: lib diff --git a/src/modules/psb_prec_mod.f90 b/src/modules/psb_prec_mod.f90 index d3ab97e4..c4a2184f 100644 --- a/src/modules/psb_prec_mod.f90 +++ b/src/modules/psb_prec_mod.f90 @@ -111,8 +111,8 @@ end interface end subroutine psb_dcslu end interface - interface psb_csrsetup - Subroutine psb_dcsrsetup(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) + interface psb_asmatbld + Subroutine psb_dasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) use psb_serial_mod Use psb_descriptor_type Use psb_prec_type @@ -124,7 +124,7 @@ end interface Character, Intent(in) :: upd integer, intent(out) :: info character(len=5), optional :: outfmt - end Subroutine psb_dcsrsetup + end Subroutine psb_dasmatbld end interface interface psb_prcaply diff --git a/src/modules/psb_prec_type.f90 b/src/modules/psb_prec_type.f90 index 6ae2fd00..329ef7ce 100644 --- a/src/modules/psb_prec_type.f90 +++ b/src/modules/psb_prec_type.f90 @@ -66,13 +66,15 @@ 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 + !Renumbering. SEE BELOW + integer, parameter :: renum_none_=0, renum_glb_=1, renum_gps_=2 !! 2 ints for 64 bit versions integer, parameter :: slu_ptr_=17, umf_symptr_=17, umf_numptr_=19 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 + ! Factorization types: none, ILU(N), ILU(E), SuperLU, UMFPACK 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 @@ -291,6 +293,13 @@ contains is_legal_n_ovr = (ip >=0) return end function is_legal_n_ovr + function is_legal_renum(ip) + integer, intent(in) :: ip + logical :: is_legal_renum + ! For the time being we are disabling renumbering options. + is_legal_renum = (ip ==0) + return + end function is_legal_renum function is_legal_jac_sweeps(ip) integer, intent(in) :: ip logical :: is_legal_jac_sweeps diff --git a/src/prec/Makefile b/src/prec/Makefile index 6615df35..704f0cb8 100644 --- a/src/prec/Makefile +++ b/src/prec/Makefile @@ -4,10 +4,10 @@ include ../../Make.inc LIBDIR=../../lib/ MPFOBJS=psb_dcslu.o psb_dbldaggrmat.o -F90OBJS= psb_dcsrsetup.o psb_dprec.o \ - psb_dprecbld.o gps.o psb_dprecfree.o psb_dprecset.o \ - psb_dgenaggrmap.o psb_dsplu.o $(MPFOBJS) -#dcoocp.o dcoocpadd.o dcoofact.o dcoolu.o dcooluadd.o\ +F90OBJS=psb_dasmatbld.o psb_dslu_bld.o psb_dumf_bld.o psb_dilu_bld.o\ + psb_dmlprc_bld.o psb_dsp_renum.o\ + psb_dprec.o psb_dprecbld.o gps.o psb_dprecfree.o psb_dprecset.o \ + psb_dgenaggrmap.o $(MPFOBJS) COBJS=fort_slu_impl.o fort_umf_impl.o INCDIRS=-I. -I.. -I$(LIBDIR) diff --git a/src/prec/psb_dcsrsetup.f90 b/src/prec/psb_dasmatbld.f90 similarity index 96% rename from src/prec/psb_dcsrsetup.f90 rename to src/prec/psb_dasmatbld.f90 index 32b60889..e25298fe 100644 --- a/src/prec/psb_dcsrsetup.f90 +++ b/src/prec/psb_dasmatbld.f90 @@ -50,7 +50,7 @@ !* * !* * !***************************************************************************** -Subroutine psb_dcsrsetup(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) +Subroutine psb_dasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) use psb_serial_mod use psb_descriptor_type @@ -82,12 +82,12 @@ Subroutine psb_dcsrsetup(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) & tot_recv, ircode, n_row, nztot,nhalo, nrow_a,err_act Logical,Parameter :: debug=.false., debugprt=.false. character(len=20) :: name, ch_err - name='psb_dcsrsetup' + name='psb_dasmatbld' if(psb_get_errstatus().ne.0) return info=0 call psb_erractionsave(err_act) - If(debug) Write(0,*)'IN DCSRSETUP ', upd + If(debug) Write(0,*)'IN DASMATBLD ', upd icontxt=desc_data%matrix_data(psb_ctxt_) tot_recv=0 @@ -161,7 +161,7 @@ Subroutine psb_dcsrsetup(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - if (debug) write(0,*) 'Early return from dcsrsetup: P>=3 N_OVR=0' + if (debug) write(0,*) 'Early return from asmatbld: P>=3 N_OVR=0' endif return endif @@ -169,7 +169,7 @@ Subroutine psb_dcsrsetup(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) call blacs_get(icontxt,10,icomm ) Call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) - If(debug)Write(0,*)'BEGIN dcsrsetup',me,upd,novr + If(debug)Write(0,*)'BEGIN dasmatbld',me,upd,novr t1 = mpi_wtime() If (upd == 'F') Then @@ -236,5 +236,5 @@ Subroutine psb_dcsrsetup(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) end if Return -End Subroutine psb_dcsrsetup +End Subroutine psb_dasmatbld diff --git a/src/prec/psb_dcslu.f90 b/src/prec/psb_dcslu.f90 index 68d75ebc..7e85d23a 100644 --- a/src/prec/psb_dcslu.f90 +++ b/src/prec/psb_dcslu.f90 @@ -35,9 +35,8 @@ !!$ !***************************************************************************** !* * -!* This is where the action takes place. As you may notice, the only * -!* piece that's really enabled is that for CSR. This is to be fixed. * -!* CSRSETUP does the setup: building the prec descriptor plus retrieving * +!* This is where the action takes place. * +!* ASMATBLD does the setup: building the prec descriptor plus retrieving * !* matrix rows if needed * !* * !* * @@ -45,7 +44,6 @@ !* * !* some open code does the renumbering * !* * -!* DSPLU90 does the actual factorization. * !* * !* * !* * @@ -80,23 +78,24 @@ subroutine psb_dcslu(a,desc_a,p,upd,info) external mpi_wtime logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false. integer istpb, istpe, ifctb, ifcte, err_act, irank, icomm, nztota, nztotb,& - & nztmp, nzl, nnr, ir, mglob, mtype, n_row, nrow_a,n_col, nhalo,lovr, ind, iind, pi,nr,ns + & nztmp, nzl, nnr, ir, mglob, mtype, n_row, nrow_a,n_col, nhalo,lovr, & + & ind, iind, pi,nr,ns integer ::icontxt,nprow,npcol,me,mycol character(len=20) :: name, ch_err - interface - subroutine psb_dsplu(a,l,u,d,info,blck) - use psb_spmat_type - integer, intent(out) :: info - type(psb_dspmat_type),intent(in) :: a - type(psb_dspmat_type),intent(inout) :: l,u - type(psb_dspmat_type),intent(in), optional, target :: blck - real(kind(1.d0)), intent(inout) :: d(:) - end subroutine psb_dsplu + interface psb_ilu_bld + subroutine psb_dilu_bld(a,l,u,d,info,blck) + use psb_spmat_type + integer, intent(out) :: info + type(psb_dspmat_type),intent(in) :: a + type(psb_dspmat_type),intent(inout) :: l,u + type(psb_dspmat_type),intent(in), optional, target :: blck + real(kind(1.d0)), intent(inout) :: d(:) + end subroutine psb_dilu_bld end interface - interface psb_csrsetup - Subroutine psb_dcsrsetup(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) + interface psb_asmatbld + Subroutine psb_dasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) use psb_serial_mod Use psb_descriptor_type Use psb_prec_type @@ -108,8 +107,24 @@ subroutine psb_dcslu(a,desc_a,p,upd,info) Character, Intent(in) :: upd integer, intent(out) :: info character(len=5), optional :: outfmt - end Subroutine psb_dcsrsetup - end interface + end Subroutine psb_dasmatbld + end interface + + interface psb_sp_renum + subroutine psb_dsp_renum(a,desc_a,blck,p,atmp,info) + use psb_prec_type + use psb_descriptor_type + use psb_spmat_type + implicit none + + ! .. array Arguments .. + type(psb_dspmat_type), intent(in) :: a,blck + type(psb_dspmat_type), intent(inout) :: atmp + type(psb_dbase_prec), intent(inout) :: p + type(psb_desc_type), intent(in) :: desc_a + integer, intent(out) :: info + end subroutine psb_dsp_renum + end interface if(psb_get_errstatus().ne.0) return info=0 @@ -121,20 +136,20 @@ subroutine psb_dcslu(a,desc_a,p,upd,info) m = a%m if (m < 0) then - info = 10 - int_err(1) = 1 - int_err(2) = m - call psb_errpush(info,name,i_err=int_err) - goto 9999 + info = 10 + int_err(1) = 1 + int_err(2) = m + call psb_errpush(info,name,i_err=int_err) + goto 9999 endif trans = 'N' unitd = 'U' if (p%iprcparm(n_ovr_) < 0) then - info = 11 - int_err(1) = 1 - int_err(2) = p%iprcparm(n_ovr_) - call psb_errpush(info,name,i_err=int_err) - goto 9999 + info = 11 + int_err(1) = 1 + int_err(2) = p%iprcparm(n_ovr_) + call psb_errpush(info,name,i_err=int_err) + goto 9999 endif ! call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) @@ -144,45 +159,45 @@ subroutine psb_dcslu(a,desc_a,p,upd,info) call psb_nullify_sp(blck) t1= mpi_wtime() - if(debug) write(0,*)me,': calling psb_csrsetup',p%iprcparm(p_type_),p%iprcparm(n_ovr_) - call psb_csrsetup(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& + if(debug) write(0,*)me,': calling psb_asmatbld',p%iprcparm(p_type_),p%iprcparm(n_ovr_) + call psb_asmatbld(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& & blck,desc_a,upd,p%desc_data,info) if(info/=0) then - info=4010 - ch_err='psb_csrsetup' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + info=4010 + ch_err='psb_asmatbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 end if t2= mpi_wtime() - if(debug) write(0,*)me,': out of psb_csrsetup' + if(debug) write(0,*)me,': out of psb_asmatbld' if (associated(p%av)) then - if (size(p%av) < bp_ilu_avsz) then - do k=1,size(p%av) - call psb_spfree(p%av(k),info) - end do - deallocate(p%av) - p%av => null() - endif + if (size(p%av) < bp_ilu_avsz) then + do k=1,size(p%av) + call psb_spfree(p%av(k),info) + end do + deallocate(p%av) + p%av => null() + endif endif if (.not.associated(p%av)) then - allocate(p%av(bp_ilu_avsz),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if + allocate(p%av(bp_ilu_avsz),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if endif do k=1,size(p%av) - call psb_nullify_sp(p%av(k)) + call psb_nullify_sp(p%av(k)) end do nrow_a = desc_a%matrix_data(psb_n_row_) call psb_spinfo(psb_nztotreq_,a,nztota,info) if(info/=0) then - info=4010 - ch_err='psb_spinfo' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + info=4010 + ch_err='psb_spinfo' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 end if n_col = desc_a%matrix_data(psb_n_col_) nhalo = n_col-nrow_a @@ -195,571 +210,165 @@ subroutine psb_dcslu(a,desc_a,p,upd,info) call psb_spall(n_row,n_row,p%av(l_pr_),nztota+lovr,info) call psb_spall(n_row,n_row,p%av(u_pr_),nztota+lovr,info) if(info/=0) then - info=4010 - ch_err='psb_spall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + info=4010 + ch_err='psb_spall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 end if if (associated(p%d)) then - if (size(p%d) < n_row) then - deallocate(p%d) - endif + if (size(p%d) < n_row) then + deallocate(p%d) + endif endif if (.not.associated(p%d)) then - allocate(p%d(n_row),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if + allocate(p%d(n_row),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if endif if (debug) then - write(0,*) me,'Done psb_csrsetup' - call blacs_barrier(icontxt,'All') + write(0,*) me,'Done psb_asmatbld' + call blacs_barrier(icontxt,'All') endif if (p%iprcparm(iren_) > 0) then - ! - ! Here we allocate a full copy to hold local A and received BLK - ! - - call psb_spinfo(psb_nztotreq_,a,nztota,info) - call psb_spinfo(psb_nztotreq_,blck,nztotb,info) - call psb_spall(atmp,nztota+nztotb,info) - if(info/=0) then - info=4011 - call psb_errpush(info,name) - goto 9999 - end if - - -! write(0,*) 'DCSLU ',nztota,nztotb,a%m - - - call apply_renum(info) - if(info/=0) then - info=4010 - ch_err='apply_ernum' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - t3 = mpi_wtime() - if (debugprt) then - call blacs_barrier(icontxt,'All') - open(40+me) - call psb_csprt(40+me,atmp,head='% Local matrix') - close(40+me) - endif - if (debug) write(0,*) me,' Factoring rows ',& - &atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1 - - ! - ! Ok, factor the matrix. - ! - t5 = mpi_wtime() - blck%m=0 - call psb_dsplu(atmp,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) - if(info/=0) then - info=4010 - ch_err='psb_dsplu' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - 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 - - - else if (p%iprcparm(iren_) == 0) then - t3 = mpi_wtime() - ! This is where we have mo renumbering, thus no need - ! for ATMP - - if (debugprt) then - open(40+me) - call blacs_barrier(icontxt,'All') - call psb_csprt(40+me,a,iv=p%desc_data%loc_to_glob,& - & head='% Local matrix') - if (p%iprcparm(p_type_)==asm_) then - call psb_csprt(40+me,blck,iv=p%desc_data%loc_to_glob,& - & irs=a%m,head='% Received rows') - endif - close(40+me) - endif - - t5= mpi_wtime() - if (debug) write(0,*) me,' Going for dsplu' - call psb_dsplu(a,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) - if(info/=0) then - info=4010 - ch_err='psb_dsplu' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (debug) write(0,*) me,' Done dsplu' - endif - - - if (debugprt) then - ! - ! Print out the factors on file. - ! - open(80+me) - - call psb_csprt(80+me,p%av(l_pr_),head='% Local L factor') - write(80+me,*) '% Diagonal: ',p%av(l_pr_)%m - do i=1,p%av(l_pr_)%m - write(80+me,*) i,i,p%d(i) - enddo - call psb_csprt(80+me,p%av(u_pr_),head='% Local U factor') - - close(80+me) - endif - - -! ierr = MPE_Log_event( ifcte, 0, "st SIMPLE" ) - t6 = mpi_wtime() -! -! write(0,'(i3,1x,a,3(1x,g18.9))') me,'renum/factor time',t3-t2,t6-t5 -! if (me==0) write(0,'(a,3(1x,g18.9))') 'renum/factor time',t3-t2,t6-t5 - - call psb_spfree(blck,info) - if(info/=0) then - info=4010 - ch_err='psb_spfree' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - if (debug) write(0,*) me,'End of cslu' + ! + ! Here we allocate a full copy to hold local A and received BLK + ! - call psb_erractionrestore(err_act) - return + call psb_spinfo(psb_nztotreq_,a,nztota,info) + call psb_spinfo(psb_nztotreq_,blck,nztotb,info) + call psb_spall(atmp,nztota+nztotb,info) + if(info/=0) then + info=4011 + call psb_errpush(info,name) + goto 9999 + end if -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error() - return - end if - return -contains + ! write(0,*) 'DCSLU ',nztota,nztotb,a%m - subroutine apply_renum(info) + call psb_sp_renum(a,desc_a,blck,p,atmp,info) - integer, intent(out) :: info - character(len=20) :: name, ch_err + if(info/=0) then + info=4010 + ch_err='psb_sp_renum' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if - if(psb_get_errstatus().ne.0) return - info=0 - name='apply_renum' - call psb_erractionsave(err_act) + t3 = mpi_wtime() + if (debugprt) then + call blacs_barrier(icontxt,'All') + open(40+me) + call psb_csprt(40+me,atmp,head='% Local matrix') + close(40+me) + endif + if (debug) write(0,*) me,' Factoring rows ',& + &atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1 -!!!!!!!!!!!!!!!! CHANGE FOR NON-CSR A ! - ! Renumbering type: - ! 1. Global column indices - ! (2. GPS band reduction disabled for the time being) - - if (p%iprcparm(iren_)==1) then - atmp%m = a%m + blck%m - atmp%k = a%k - atmp%fida='CSR' - atmp%descra = 'GUN' - - ! This is the renumbering coherent with global indices.. - mglob = desc_a%matrix_data(psb_m_) - ! - ! Remember: we have switched IA1=COLS and IA2=ROWS - ! Now identify the set of distinct local column indices - ! - - nnr = p%desc_data%matrix_data(psb_n_row_) - allocate(p%perm(nnr),p%invperm(nnr),itmp2(nnr),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - do k=1,nnr - itmp2(k) = p%desc_data%loc_to_glob(k) - enddo - ! - ! We want: NEW(I) = OLD(PERM(I)) - ! - call isrx(nnr,itmp2,p%perm) - - do k=1, nnr - p%invperm(p%perm(k)) = k - enddo - t3 = mpi_wtime() - - ! Build ATMP with new numbering. - - allocate(itmp(max(8,atmp%m+2,nztmp+2)),rtmp(atmp%m),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - j = 1 - atmp%ia2(1) = 1 - do i=1, atmp%m - ir = p%perm(i) - - if (ir <= a%m ) then - - nzl = a%ia2(ir+1) - a%ia2(ir) - if (nzl > size(rtmp)) then - call psb_realloc(nzl,rtmp,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - jj = a%ia2(ir) - k=0 - do kk=1, nzl - if (a%ia1(jj+kk-1)<=atmp%m) then - k = k + 1 - rtmp(k) = a%aspk(jj+kk-1) - atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) - endif - enddo - call isrx(k,atmp%ia1(j:j+k-1),itmp2) - do kk=1,k - atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) - enddo - - else if (ir <= atmp%m ) then - - ir = ir - a%m - nzl = blck%ia2(ir+1) - blck%ia2(ir) - if (nzl > size(rtmp)) then - call psb_realloc(nzl,rtmp,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - jj = blck%ia2(ir) - k=0 - do kk=1, nzl - if (blck%ia1(jj+kk-1)<=atmp%m) then - k = k + 1 - rtmp(k) = blck%aspk(jj+kk-1) - atmp%ia1(j+k-1) = p%invperm(blck%ia1(jj+kk-1)) - endif - enddo - call isrx(k,atmp%ia1(j:j+k-1),itmp2) - do kk=1,k - atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) - enddo - - else - write(0,*) 'Row index error 1 :',i,ir - endif - - j = j + k - atmp%ia2(i+1) = j - - enddo - - t4 = mpi_wtime() - - - deallocate(itmp,itmp2,rtmp) - - else if (p%iprcparm(iren_)==2) then - - atmp%m = a%m + blck%m - atmp%k = a%k - atmp%fida='CSR' - atmp%descra = 'GUN' - do i=1, a%m - atmp%ia2(i) = a%ia2(i) - do j= a%ia2(i), a%ia2(i+1)-1 - atmp%ia1(j) = a%ia1(j) - enddo - enddo - atmp%ia2(a%m+1) = a%ia2(a%m+1) - - if (blck%m>0) then - do i=1, blck%m - atmp%ia2(a%m+i) = nztota+blck%ia2(i) - do j= blck%ia2(i), blck%ia2(i+1)-1 - atmp%ia1(nztota+j) = blck%ia1(j) - enddo - enddo - atmp%ia2(atmp%m+1) = nztota+blck%ia2(blck%m+1) - endif - nztmp = atmp%ia2(atmp%m+1) - 1 - - - ! This is a renumbering with Gibbs-Poole-Stockmeyer - ! band reduction. Switched off for now. To be fixed, - ! gps_reduction should get p%perm. - -! write(0,*) me,' Renumbering: realloc perms',atmp%m - call psb_realloc(atmp%m,p%perm,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_realloc(atmp%m,p%invperm,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - allocate(itmp(max(8,atmp%m+2,nztmp+2)),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - itmp(1:8) = 0 -! write(0,*) me,' Renumbering: Calling Metis' -! call blacs_barrier(icontxt,'All') - -! write(0,*) size(p%av(u_pr_)%pl),size(p%av(l_pr_)%pr) - call gps_reduction(atmp%m,atmp%ia2,atmp%ia1,p%perm,p%invperm,info) - if(info/=0) then - info=4010 - ch_err='gps_reduction' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - -! write(0,*) me,' Renumbering: Done GPS' - call blacs_barrier(icontxt,'All') - do i=1, atmp%m - if (p%perm(i) /= i) then - write(0,*) me,' permutation is not identity ' - exit - endif - enddo - - - - do k=1, nnr - p%invperm(p%perm(k)) = k - enddo - t3 = mpi_wtime() - - ! Build ATMP with new numbering. - - allocate(itmp2(max(8,atmp%m+2,nztmp+2)),rtmp(atmp%m),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - j = 1 - atmp%ia2(1) = 1 - do i=1, atmp%m - ir = p%perm(i) - - if (ir <= a%m ) then - - nzl = a%ia2(ir+1) - a%ia2(ir) - if (nzl > size(rtmp)) then - call psb_realloc(nzl,rtmp,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - jj = a%ia2(ir) - k=0 - do kk=1, nzl - if (a%ia1(jj+kk-1)<=atmp%m) then - k = k + 1 - rtmp(k) = a%aspk(jj+kk-1) - atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) - endif - enddo - call isrx(k,atmp%ia1(j:j+k-1),itmp2) - do kk=1,k - atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) - enddo - - else if (ir <= atmp%m ) then - - ir = ir - a%m - nzl = blck%ia2(ir+1) - blck%ia2(ir) - if (nzl > size(rtmp)) then - call psb_realloc(nzl,rtmp,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - jj = blck%ia2(ir) - k=0 - do kk=1, nzl - if (blck%ia1(jj+kk-1)<=atmp%m) then - k = k + 1 - rtmp(k) = blck%aspk(jj+kk-1) - atmp%ia1(j+k-1) = p%invperm(blck%ia1(jj+kk-1)) - endif - enddo - call isrx(k,atmp%ia1(j:j+k-1),itmp2) - do kk=1,k - atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) - enddo - - else - write(0,*) 'Row index error 1 :',i,ir - endif - - j = j + k - atmp%ia2(i+1) = j - - enddo - - t4 = mpi_wtime() - - - - deallocate(itmp,itmp2,rtmp) - + ! Ok, factor the matrix. + ! + t5 = mpi_wtime() + blck%m=0 + call psb_ilu_bld(atmp,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) + if(info/=0) then + info=4010 + ch_err='psb_ilu_bld' + 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 + 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 - return - end subroutine apply_renum + else if (p%iprcparm(iren_) == 0) then + t3 = mpi_wtime() + ! This is where we have mo renumbering, thus no need + ! for ATMP + + if (debugprt) then + open(40+me) + call blacs_barrier(icontxt,'All') + call psb_csprt(40+me,a,iv=p%desc_data%loc_to_glob,& + & head='% Local matrix') + if (p%iprcparm(p_type_)==asm_) then + call psb_csprt(40+me,blck,iv=p%desc_data%loc_to_glob,& + & irs=a%m,head='% Received rows') + endif + close(40+me) + endif - subroutine gps_reduction(m,ia,ja,perm,iperm,info) - integer i,j,dgConn,Npnt,m - integer n,idpth,ideg,ibw2,ipf2 - integer,dimension(:) :: perm,iperm,ia,ja - integer, intent(out) :: info + t5= mpi_wtime() + if (debug) write(0,*) me,' Going for dilu_bld' + call psb_ilu_bld(a,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) + if(info/=0) then + info=4010 + ch_err='psb_ilu_bld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (debug) write(0,*) me,' Done dilu_bld' + endif - integer,dimension(:,:),allocatable::NDstk - integer,dimension(:),allocatable::iOld,renum,ndeg,lvl,lvls1,lvls2,ccstor - !--- Per la common area. - common /gra/ n,iDpth,iDeg + if (debugprt) then + ! + ! Print out the factors on file. + ! + open(80+me) - character(len=20) :: name, ch_err + call psb_csprt(80+me,p%av(l_pr_),head='% Local L factor') + write(80+me,*) '% Diagonal: ',p%av(l_pr_)%m + do i=1,p%av(l_pr_)%m + write(80+me,*) i,i,p%d(i) + enddo + call psb_csprt(80+me,p%av(u_pr_),head='% Local U factor') - if(psb_get_errstatus().ne.0) return - info=0 - name='gps_reduction' - call psb_erractionsave(err_act) + close(80+me) + endif - !--- Calcolo il massimo grado di connettivita'. - npnt = m - write(6,*) ' GPS su ',npnt - dgConn=0 - do i=1,m - dgconn = max(dgconn,(ia(i+1)-ia(i))) - enddo - !--- Il max valore di connettivita' e "dgConn" - - !--- Valori della common - n=Npnt !--- Numero di righe - iDeg=dgConn !--- Massima connettivita' -! iDpth= !--- Numero di livelli non serve settarlo - - allocate(NDstk(Npnt,dgConn),stat=info) - if (info/=0) then - info=4000 - call psb_errpush(info,name) - goto 9999 - else - write(0,*) 'gps_reduction first alloc OK' - endif - allocate(iOld(Npnt),renum(Npnt+1),ndeg(Npnt),lvl(Npnt),lvls1(Npnt),& - &lvls2(Npnt),ccstor(Npnt),stat=info) - if (info/=0) then - info=4000 - call psb_errpush(info,name) - goto 9999 - else - write(0,*) 'gps_reduction 2nd alloc OK' - endif + ! ierr = MPE_Log_event( ifcte, 0, "st SIMPLE" ) + t6 = mpi_wtime() + ! + ! write(0,'(i3,1x,a,3(1x,g18.9))') me,'renum/factor time',t3-t2,t6-t5 + ! if (me==0) write(0,'(a,3(1x,g18.9))') 'renum/factor time',t3-t2,t6-t5 - !--- Prepariamo il grafo della matrice - Ndstk(:,:)=0 - do i=1,Npnt - k=0 - do j = ia(i),ia(i+1) - 1 - if ((1<=ja(j)).and.( ja( j ) /= i ).and.(ja(j)<=npnt)) then - k = k+1 - Ndstk(i,k)=ja(j) - endif - enddo - ndeg(i)=k - enddo + call psb_spfree(blck,info) + if(info/=0) then + info=4010 + ch_err='psb_spfree' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if - !--- Numerazione. - do i=1,Npnt - iOld(i)=i - enddo - write(0,*) 'gps_red : Preparation done' - !--- - !--- Chiamiamo funzione reduce. - call reduce(Ndstk,Npnt,iOld,renum,ndeg,lvl,lvls1, lvls2,ccstor,ibw2,ipf2) - write(0,*) 'gps_red : Done reduce' - !--- Permutazione - perm(1:Npnt)=renum(1:Npnt) - !--- Inversa permutazione - do i=1,Npnt - iperm(perm(i))=i - enddo - !--- Puliamo tutto. - deallocate(NDstk,iOld,renum,ndeg,lvl,lvls1,lvls2,ccstor) + if (debug) write(0,*) me,'End of cslu' - call psb_erractionrestore(err_act) - return + 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 + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() return - - end subroutine gps_reduction + end if + return end subroutine psb_dcslu diff --git a/src/prec/psb_dsplu.f90 b/src/prec/psb_dilu_bld.f90 similarity index 98% rename from src/prec/psb_dsplu.f90 rename to src/prec/psb_dilu_bld.f90 index cf976eb0..c1ac28e7 100644 --- a/src/prec/psb_dsplu.f90 +++ b/src/prec/psb_dilu_bld.f90 @@ -33,7 +33,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_dsplu(a,l,u,d,info,blck) +subroutine psb_dilu_bld(a,l,u,d,info,blck) ! ! This routine copies and factors "on the fly" from A and BLCK @@ -86,11 +86,11 @@ subroutine psb_dsplu(a,l,u,d,info,blck) blck_%m=0 endif - call psb_dspluint(m,a%m,a,blck_%m,blck_,& + call psb_dilu_bldint(m,a%m,a,blck_%m,blck_,& & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info) if(info.ne.0) then info=4010 - ch_err='psb_dspluint' + ch_err='psb_dilu_bldint' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -130,7 +130,7 @@ subroutine psb_dsplu(a,l,u,d,info,blck) return contains - subroutine psb_dspluint(m,ma,a,mb,b,& + subroutine psb_dilu_bldint(m,ma,a,mb,b,& & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info) implicit none @@ -148,7 +148,7 @@ contains integer :: int_err(5) character(len=20) :: name, ch_err - name='psb_dspluint' + name='psb_dilu_bldint' if(psb_get_errstatus().ne.0) return info=0 call psb_erractionsave(err_act) @@ -473,5 +473,5 @@ contains return end if return - end subroutine psb_dspluint -end subroutine psb_dsplu + end subroutine psb_dilu_bldint +end subroutine psb_dilu_bld diff --git a/src/prec/psb_dmlprc_bld.f90 b/src/prec/psb_dmlprc_bld.f90 new file mode 100644 index 00000000..1c646508 --- /dev/null +++ b/src/prec/psb_dmlprc_bld.f90 @@ -0,0 +1,241 @@ +!!$ +!!$ +!!$ MPcube: Multilevel Parallel Preconditioners Package +!!$ for +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Daniela Di Serafino II University of Naples +!!$ Pasqua D'Ambra ICAR-CNR +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MPCUBE group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine psb_dmlprc_bld(a,desc_a,p,info) + + use psb_serial_mod + use psb_tools_mod + use psb_descriptor_type + use psb_prec_type + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_dspmat_type), intent(in), target :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dbase_prec), intent(inout) :: p + integer, intent(out) :: info + + integer :: i, nrg, nzg, err_act,k + character(len=20) :: name, ch_err + + interface psb_ilu_bld + subroutine psb_dilu_bld(a,l,u,d,info,blck) + use psb_spmat_type + integer, intent(out) :: info + type(psb_dspmat_type),intent(in) :: a + type(psb_dspmat_type),intent(inout) :: l,u + type(psb_dspmat_type),intent(in), optional, target :: blck + real(kind(1.d0)), intent(inout) :: d(:) + end subroutine psb_dilu_bld + end interface + + interface psb_genaggrmap + subroutine psb_dgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info) + use psb_spmat_type + use psb_descriptor_type + implicit none + integer, intent(in) :: aggr_type + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer, pointer :: ilaggr(:),nlaggr(:) + integer, intent(out) :: info + end subroutine psb_dgenaggrmap + end interface + + interface psb_bldaggrmat + subroutine psb_dbldaggrmat(a,desc_a,p,info) + use psb_prec_type + use psb_descriptor_type + use psb_spmat_type + type(psb_dspmat_type), intent(in), target :: a + type(psb_dbase_prec), intent(inout) :: p + type(psb_desc_type), intent(in) :: desc_a + integer, intent(out) :: info + end subroutine psb_dbldaggrmat + end interface + + integer :: icontxt, nprow, npcol, me, mycol + + name='psb_mlprec_bld' + if(psb_get_errstatus().ne.0) return + info=0 + call psb_erractionsave(err_act) + + p%aorig => a + allocate(p%av(smth_avsz),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + do i=1, smth_avsz + call psb_nullify_sp(p%av(i)) + call psb_spall(0,0,p%av(i),1,info) + if(info /= 0) then + info=4010 + ch_err='psb_spall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + end do + + ! Currently this is ignored by gen_aggrmap, but it could be + ! changed in the future. Need to package nlaggr & mlia in a + ! private data structure? + + call psb_genaggrmap(p%iprcparm(aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info) + if(info /= 0) then + info=4010 + ch_err='psb_gen_aggrmap' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call psb_bldaggrmat(a,desc_a,p,info) + if(info /= 0) then + info=4010 + ch_err='psb_bld_aggrmat' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + nrg = p%av(ac_)%m + call psb_spinfo(psb_nztotreq_,p%av(ac_),nzg,info) + call psb_ipcoo2csr(p%av(ac_),info) + if(info /= 0) then + info=4011 + ch_err='psb_ipcoo2csr' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(p%d(nrg),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + select case(p%iprcparm(f_type_)) + case(f_ilu_n_,f_ilu_e_) + call psb_spreall(p%av(l_pr_),nzg,info) + call psb_spreall(p%av(u_pr_),nzg,info) + call psb_ilu_bld(p%av(ac_),p%av(l_pr_),p%av(u_pr_),p%d,info) + if(info /= 0) then + info=4011 + ch_err='psb_ilu_bld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(f_slu_) +!!$ call psb_spall(0,0,p%av(l_pr_),1,info) +!!$ call psb_spall(0,0,p%av(u_pr_),1,info) + call psb_ipcsr2coo(p%av(ac_),info) + if(info /= 0) then + info=4011 + ch_err='psb_ipcsr2coo' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + k=0 + do i=1,p%av(ac_)%infoa(psb_nnz_) + if (p%av(ac_)%ia2(i) <= p%av(ac_)%m) then + k = k + 1 + p%av(ac_)%aspk(k) = p%av(ac_)%aspk(i) + p%av(ac_)%ia1(k) = p%av(ac_)%ia1(i) + p%av(ac_)%ia2(k) = p%av(ac_)%ia2(i) + end if + end do + p%av(ac_)%infoa(psb_nnz_) = k + call psb_ipcoo2csr(p%av(ac_),info) + call psb_spinfo(psb_nztotreq_,p%av(ac_),nzg,info) + call fort_slu_factor(nrg,nzg,& + & p%av(ac_)%aspk,p%av(ac_)%ia2,p%av(ac_)%ia1,p%iprcparm(slu_ptr_),info) + if(info /= 0) then + info=4011 + ch_err='psb_fort_slu_factor' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(f_umf_) +!!$ call psb_spall(0,0,p%av(l_pr_),1,info) +!!$ call psb_spall(0,0,p%av(u_pr_),1,info) + call psb_ipcsr2coo(p%av(ac_),info) + if(info /= 0) then + info=4011 + ch_err='psb_ipcsr2coo' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + k=0 + do i=1,p%av(ac_)%infoa(psb_nnz_) + if (p%av(ac_)%ia2(i) <= p%av(ac_)%m) then + k = k + 1 + p%av(ac_)%aspk(k) = p%av(ac_)%aspk(i) + p%av(ac_)%ia1(k) = p%av(ac_)%ia1(i) + p%av(ac_)%ia2(k) = p%av(ac_)%ia2(i) + end if + end do + p%av(ac_)%infoa(psb_nnz_) = k + call psb_ipcoo2csc(p%av(ac_),info) + call psb_spinfo(psb_nztotreq_,p%av(ac_),nzg,info) + call fort_umf_factor(nrg,nzg,& + & p%av(ac_)%aspk,p%av(ac_)%ia1,p%av(ac_)%ia2,& + & p%iprcparm(umf_symptr_),p%iprcparm(umf_numptr_),info) + if(info /= 0) then + info=4011 + ch_err='psb_fort_umf_factor' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case default + write(0,*) 'Invalid fact type for multi level',(p%iprcparm(f_type_)) + end select + + 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_dmlprc_bld diff --git a/src/prec/psb_dprecbld.f90 b/src/prec/psb_dprecbld.f90 index 1db7ca8f..2a5828dc 100644 --- a/src/prec/psb_dprecbld.f90 +++ b/src/prec/psb_dprecbld.f90 @@ -53,45 +53,61 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd) character, intent(in), optional :: upd interface psb_cslu - subroutine psb_dcslu(a,desc_data,p,upd,info) - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbase_prec), intent(inout) :: p - character, intent(in) :: upd - end subroutine psb_dcslu + subroutine psb_dcslu(a,desc_data,p,upd,info) + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target :: a + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbase_prec), intent(inout) :: p + character, intent(in) :: upd + end subroutine psb_dcslu end interface - interface - subroutine psb_superlu_bld(a,desc_a,p,info) + interface psb_slu_bld + subroutine psb_dslu_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_superlu_bld + end subroutine psb_dslu_bld end interface - interface - subroutine psb_umf_bld(a,desc_a,p,info) + + interface psb_umf_bld + subroutine psb_dumf_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 subroutine psb_dumf_bld + end interface + + interface psb_mlprc_bld + subroutine psb_dmlprc_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_dmlprc_bld end interface ! Local scalars @@ -123,14 +139,14 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd) call blacs_gridinfo(icontxt, nprow, npcol, me, mycol) if (present(upd)) then - if (debug) write(0,*) 'UPD ', upd - if ((UPD.eq.'F').or.(UPD.eq.'T')) then - IUPD=UPD - else - IUPD='F' - endif + if (debug) write(0,*) 'UPD ', upd + if ((UPD.eq.'F').or.(UPD.eq.'T')) then + IUPD=UPD + else + IUPD='F' + endif else - IUPD='F' + IUPD='F' endif if (.not.associated(p%baseprecv)) then @@ -154,151 +170,147 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd) select case(p%baseprecv(1)%iprcparm(p_type_)) case (noprec_) - ! Do nothing. + ! Do nothing. case (diagsc_) - if (debug) write(0,*) 'Precond: Diagonal scaling' - ! diagonal scaling + if (debug) write(0,*) 'Precond: Diagonal scaling' + ! diagonal scaling - call psb_realloc(n_col,p%baseprecv(1)%d,info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='psb_realloc') - goto 9999 - end if + call psb_realloc(n_col,p%baseprecv(1)%d,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_realloc') + goto 9999 + end if + + call psb_csrws(p%baseprecv(1)%d,a,info,trans='N') + if(info /= 0) then + info=4010 + ch_err='psb_csrws' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if - call psb_csrws(p%baseprecv(1)%d,a,info,trans='N') - if(info /= 0) then + if (debug) write(ilout+me,*) 'VDIAG ',n_row + do i=1,n_row + if (p%baseprecv(1)%d(i).eq.0.0d0) then + p%baseprecv(1)%d(i)=1.d0 + else + p%baseprecv(1)%d(i) = 1.d0/p%baseprecv(1)%d(i) + endif + + if (debug) write(ilout+me,*) i,desc_a%loc_to_glob(i),& + & p%baseprecv(1)%d(i) + if (p%baseprecv(1)%d(i).lt.0.d0) then + write(0,*) me,'Negative RWS? ',i,p%baseprecv(1)%d(i) + endif + end do + if (a%pl(1) /= 0) then + allocate(work(n_row),stat=info) + if (info /= 0) then + info=4000 + call psb_errpush(info,name) + goto 9999 + end if + call psb_gelp('n',a%pl,p%baseprecv(1)%d,desc_a,info) + if(info /= 0) then info=4010 - ch_err='psb_csrws' + ch_err='psb_dgelp' call psb_errpush(info,name,a_err=ch_err) goto 9999 - end if - - if (debug) write(ilout+me,*) 'VDIAG ',n_row - do i=1,n_row - if (p%baseprecv(1)%d(i).eq.0.0d0) then - p%baseprecv(1)%d(i)=1.d0 - else - p%baseprecv(1)%d(i) = 1.d0/p%baseprecv(1)%d(i) - endif - - if (debug) write(ilout+me,*) i,desc_a%loc_to_glob(i),& - & p%baseprecv(1)%d(i) - if (p%baseprecv(1)%d(i).lt.0.d0) then - write(0,*) me,'Negative RWS? ',i,p%baseprecv(1)%d(i) - endif - end do - if (a%pl(1) /= 0) then - allocate(work(n_row),stat=info) - if (info /= 0) then - info=4000 - call psb_errpush(info,name) - goto 9999 - end if - call psb_gelp('n',a%pl,p%baseprecv(1)%d,desc_a,info) - if(info /= 0) then - info=4010 - ch_err='psb_dgelp' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if + end if - deallocate(work) - endif + deallocate(work) + endif - if (debug) then - allocate(gd(mglob),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if + if (debug) then + allocate(gd(mglob),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if - call psb_gather(gd, p%baseprecv(1)%d, desc_a, info, iroot=iroot) - if(info /= 0) then - info=4010 - ch_err='psb_dgatherm' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if + call psb_gather(gd, p%baseprecv(1)%d, desc_a, info, iroot=iroot) + if(info /= 0) then + info=4010 + ch_err='psb_dgatherm' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if - if (me.eq.iroot) then - write(iout+nprow,*) 'VDIAG CHECK ',mglob - do i=1,mglob - write(iout+nprow,*) i,gd(i) - enddo - endif - deallocate(gd) - endif - if (debug) write(*,*) 'Preconditioner DIAG computed OK' + if (me.eq.iroot) then + write(iout+nprow,*) 'VDIAG CHECK ',mglob + do i=1,mglob + write(iout+nprow,*) i,gd(i) + enddo + endif + deallocate(gd) + endif + if (debug) write(*,*) 'Preconditioner DIAG computed OK' case (bja_,asm_) - call psb_check_def(p%baseprecv(1)%iprcparm(n_ovr_),'overlap',& - & 0,is_legal_n_ovr) - call psb_check_def(p%baseprecv(1)%iprcparm(restr_),'restriction',& - & psb_halo_,is_legal_restrict) - call psb_check_def(p%baseprecv(1)%iprcparm(prol_),'prolongator',& - & psb_none_,is_legal_prolong) - - if ((p%baseprecv(1)%iprcparm(iren_)<0).or.(p%baseprecv(1)%iprcparm(iren_)>2)) then - write(0,*) 'Bad PREC%IRENUM value, defaulting to 0', & - & p%baseprecv(1)%iprcparm(iren_) - p%baseprecv(1)%iprcparm(iren_) = 0 - endif + call psb_check_def(p%baseprecv(1)%iprcparm(n_ovr_),'overlap',& + & 0,is_legal_n_ovr) + call psb_check_def(p%baseprecv(1)%iprcparm(restr_),'restriction',& + & psb_halo_,is_legal_restrict) + call psb_check_def(p%baseprecv(1)%iprcparm(prol_),'prolongator',& + & psb_none_,is_legal_prolong) + call psb_check_def(p%baseprecv(1)%iprcparm(iren_),'renumbering',& + & renum_none_,is_legal_renum) - if (debug) write(0,*)me, ': Calling PSB_DCSLU' + if (debug) write(0,*)me, ': Calling PSB_DCSLU' - select case(p%baseprecv(1)%iprcparm(f_type_)) + select case(p%baseprecv(1)%iprcparm(f_type_)) - case(f_ilu_n_,f_ilu_e_) - call psb_cslu(a,desc_a,p%baseprecv(1),iupd,info) - if(debug) write(0,*)me,': out of psb_dcslu' - if(info /= 0) then - info=4010 - ch_err='psb_dcslu' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if + case(f_ilu_n_,f_ilu_e_) + call psb_cslu(a,desc_a,p%baseprecv(1),iupd,info) + if(debug) write(0,*)me,': out of psb_dcslu' + if(info /= 0) then + info=4010 + ch_err='psb_dcslu' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if - case(f_slu_) - p%baseprecv(1)%av => null() - if(debug) write(0,*)me,': calling superlu_bld' - call psb_superlu_bld(a,desc_a,p%baseprecv(1),info) - if(info /= 0) then - info=4010 - ch_err='superlu_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if + case(f_slu_) + p%baseprecv(1)%av => null() + if(debug) write(0,*)me,': calling slu_bld' + call psb_slu_bld(a,desc_a,p%baseprecv(1),info) + if(info /= 0) then + info=4010 + ch_err='slu_bld' + call psb_errpush(info,name,a_err=ch_err) + 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_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??' + case(f_none_) + write(0,*) 'Fact=None in PRECBLD Bja/ASM??' - case default - write(0,*) 'Unknown factor type in precbld bja/asm: ',& - &p%baseprecv(1)%iprcparm(f_type_) - end select + case default + write(0,*) 'Unknown factor type in precbld bja/asm: ',& + &p%baseprecv(1)%iprcparm(f_type_) + end select end select if (size(p%baseprecv) >1) then - + if (.not.associated(p%baseprecv(2)%iprcparm)) then info = 2222 call psb_errpush(info,name) @@ -332,12 +344,12 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd) call psb_check_def(p%baseprecv(2)%iprcparm(jac_sweeps_),'Jacobi sweeps',& & 1,is_legal_jac_sweeps) - call psb_mlprec_bld(a,desc_a,p%baseprecv(2),info) + call psb_mlprc_bld(a,desc_a,p%baseprecv(2),info) if(info /= 0) then - info=4010 - ch_err='psb_mlprec_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + info=4010 + ch_err='psb_mlprc_bld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 end if endif @@ -348,682 +360,10 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd) 9999 continue call psb_erractionrestore(err_act) if (err_act.eq.act_abort) then - call psb_error() - return + call psb_error() + return end if return end subroutine psb_dprecbld - -!!$ -!!$ -!!$ MPcube: Multilevel Parallel Preconditioners Package -!!$ for -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ Daniela Di Serafino II University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MPCUBE group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_superlu_bld(a,desc_a,p,info) - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - use psb_tools_mod - 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 - - if(psb_get_errstatus().ne.0) return - info=0 - name='psb_superlu_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, 'SPLUBLD: 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, 'SPLUBLD: 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, 'SPLUBLD: Done csrsetup',info,nzb,blck%fida - call blacs_barrier(icontxt,'All') - endif - if (nzb > 0 ) then - if (size(atmp%aspk) 0 ) then - if (size(atmp%aspk) a - allocate(p%av(smth_avsz),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - do i=1, smth_avsz - call psb_nullify_sp(p%av(i)) - call psb_spall(0,0,p%av(i),1,info) - if(info /= 0) then - info=4010 - ch_err='psb_spall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - end do - - ! Currently this is ignored by gen_aggrmap, but it could be - ! changed in the future. Need to package nlaggr & mlia in a - ! private data structure? - - call psb_genaggrmap(p%iprcparm(aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info) - if(info /= 0) then - info=4010 - ch_err='psb_gen_aggrmap' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_bldaggrmat(a,desc_a,p,info) - if(info /= 0) then - info=4010 - ch_err='psb_bld_aggrmat' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - nrg = p%av(ac_)%m - call psb_spinfo(psb_nztotreq_,p%av(ac_),nzg,info) - call psb_ipcoo2csr(p%av(ac_),info) - if(info /= 0) then - info=4011 - ch_err='psb_ipcoo2csr' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - allocate(p%d(nrg),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - select case(p%iprcparm(f_type_)) - case(f_ilu_n_,f_ilu_e_) - call psb_spreall(p%av(l_pr_),nzg,info) - call psb_spreall(p%av(u_pr_),nzg,info) - call psb_splu(p%av(ac_),p%av(l_pr_),p%av(u_pr_),p%d,info) - if(info /= 0) then - info=4011 - ch_err='psb_splu' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case(f_slu_) -!!$ call psb_spall(0,0,p%av(l_pr_),1,info) -!!$ call psb_spall(0,0,p%av(u_pr_),1,info) - call psb_ipcsr2coo(p%av(ac_),info) - if(info /= 0) then - info=4011 - ch_err='psb_ipcsr2coo' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - k=0 - do i=1,p%av(ac_)%infoa(psb_nnz_) - if (p%av(ac_)%ia2(i) <= p%av(ac_)%m) then - k = k + 1 - p%av(ac_)%aspk(k) = p%av(ac_)%aspk(i) - p%av(ac_)%ia1(k) = p%av(ac_)%ia1(i) - p%av(ac_)%ia2(k) = p%av(ac_)%ia2(i) - end if - end do - p%av(ac_)%infoa(psb_nnz_) = k - call psb_ipcoo2csr(p%av(ac_),info) - call psb_spinfo(psb_nztotreq_,p%av(ac_),nzg,info) - call fort_slu_factor(nrg,nzg,& - & p%av(ac_)%aspk,p%av(ac_)%ia2,p%av(ac_)%ia1,p%iprcparm(slu_ptr_),info) - if(info /= 0) then - info=4011 - ch_err='psb_fort_slu_factor' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case(f_umf_) -!!$ call psb_spall(0,0,p%av(l_pr_),1,info) -!!$ call psb_spall(0,0,p%av(u_pr_),1,info) - call psb_ipcsr2coo(p%av(ac_),info) - if(info /= 0) then - info=4011 - ch_err='psb_ipcsr2coo' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - k=0 - do i=1,p%av(ac_)%infoa(psb_nnz_) - if (p%av(ac_)%ia2(i) <= p%av(ac_)%m) then - k = k + 1 - p%av(ac_)%aspk(k) = p%av(ac_)%aspk(i) - p%av(ac_)%ia1(k) = p%av(ac_)%ia1(i) - p%av(ac_)%ia2(k) = p%av(ac_)%ia2(i) - end if - end do - p%av(ac_)%infoa(psb_nnz_) = k - call psb_ipcoo2csc(p%av(ac_),info) - call psb_spinfo(psb_nztotreq_,p%av(ac_),nzg,info) - call fort_umf_factor(nrg,nzg,& - & p%av(ac_)%aspk,p%av(ac_)%ia1,p%av(ac_)%ia2,& - & p%iprcparm(umf_symptr_),p%iprcparm(umf_numptr_),info) - if(info /= 0) then - info=4011 - ch_err='psb_fort_umf_factor' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case default - write(0,*) 'Invalid fact type for multi level',(p%iprcparm(f_type_)) - end select - - 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_mlprec_bld diff --git a/src/prec/psb_dslu_bld.f90 b/src/prec/psb_dslu_bld.f90 new file mode 100644 index 00000000..df97215e --- /dev/null +++ b/src/prec/psb_dslu_bld.f90 @@ -0,0 +1,213 @@ +!!$ +!!$ +!!$ MPcube: Multilevel Parallel Preconditioners Package +!!$ for +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Daniela Di Serafino II University of Naples +!!$ Pasqua D'Ambra ICAR-CNR +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MPCUBE group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine psb_dslu_bld(a,desc_a,p,info) + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_tools_mod + 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_asmatbld + Subroutine psb_dasmatbld(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_dasmatbld + end interface + + if(psb_get_errstatus().ne.0) return + info=0 + name='psb_slu_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, 'SPLUBLD: 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, 'SPLUBLD: Done csdp',info,nza,atmp%m,atmp%k + call blacs_barrier(icontxt,'All') + endif + call psb_asmatbld(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_asmatbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + nzb = blck%infoa(psb_nnz_) + if (Debug) then + write(0,*) me, 'SPLUBLD: Done asmatbld',info,nzb,blck%fida + call blacs_barrier(icontxt,'All') + endif + if (nzb > 0 ) then + if (size(atmp%aspk) size(rtmp)) then + call psb_realloc(nzl,rtmp,info) + if(info/=0) then + info=4010 + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + endif + jj = a%ia2(ir) + k=0 + do kk=1, nzl + if (a%ia1(jj+kk-1)<=atmp%m) then + k = k + 1 + rtmp(k) = a%aspk(jj+kk-1) + atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) + endif + enddo + call isrx(k,atmp%ia1(j:j+k-1),itmp2) + do kk=1,k + atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) + enddo + + else if (ir <= atmp%m ) then + + ir = ir - a%m + nzl = blck%ia2(ir+1) - blck%ia2(ir) + if (nzl > size(rtmp)) then + call psb_realloc(nzl,rtmp,info) + if(info/=0) then + info=4010 + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + endif + jj = blck%ia2(ir) + k=0 + do kk=1, nzl + if (blck%ia1(jj+kk-1)<=atmp%m) then + k = k + 1 + rtmp(k) = blck%aspk(jj+kk-1) + atmp%ia1(j+k-1) = p%invperm(blck%ia1(jj+kk-1)) + endif + enddo + call isrx(k,atmp%ia1(j:j+k-1),itmp2) + do kk=1,k + atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) + enddo + + else + write(0,*) 'Row index error 1 :',i,ir + endif + + j = j + k + atmp%ia2(i+1) = j + + enddo + + t4 = mpi_wtime() + + + deallocate(itmp,itmp2,rtmp) + + else if (p%iprcparm(iren_)==renum_gps_) then + + atmp%m = a%m + blck%m + atmp%k = a%k + atmp%fida='CSR' + atmp%descra = 'GUN' + do i=1, a%m + atmp%ia2(i) = a%ia2(i) + do j= a%ia2(i), a%ia2(i+1)-1 + atmp%ia1(j) = a%ia1(j) + enddo + enddo + atmp%ia2(a%m+1) = a%ia2(a%m+1) + + if (blck%m>0) then + do i=1, blck%m + atmp%ia2(a%m+i) = nztota+blck%ia2(i) + do j= blck%ia2(i), blck%ia2(i+1)-1 + atmp%ia1(nztota+j) = blck%ia1(j) + enddo + enddo + atmp%ia2(atmp%m+1) = nztota+blck%ia2(blck%m+1) + endif + nztmp = atmp%ia2(atmp%m+1) - 1 + + + ! This is a renumbering with Gibbs-Poole-Stockmeyer + ! band reduction. Switched off for now. To be fixed, + ! gps_reduction should get p%perm. + + ! write(0,*) me,' Renumbering: realloc perms',atmp%m + call psb_realloc(atmp%m,p%perm,info) + if(info/=0) then + info=4010 + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call psb_realloc(atmp%m,p%invperm,info) + if(info/=0) then + info=4010 + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(itmp(max(8,atmp%m+2,nztmp+2)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + itmp(1:8) = 0 + ! write(0,*) me,' Renumbering: Calling Metis' + ! call blacs_barrier(icontxt,'All') + + ! write(0,*) size(p%av(u_pr_)%pl),size(p%av(l_pr_)%pr) + call gps_reduction(atmp%m,atmp%ia2,atmp%ia1,p%perm,p%invperm,info) + if(info/=0) then + info=4010 + ch_err='gps_reduction' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! write(0,*) me,' Renumbering: Done GPS' + call blacs_barrier(icontxt,'All') + do i=1, atmp%m + if (p%perm(i) /= i) then + write(0,*) me,' permutation is not identity ' + exit + endif + enddo + + + + do k=1, nnr + p%invperm(p%perm(k)) = k + enddo + t3 = mpi_wtime() + + ! Build ATMP with new numbering. + + allocate(itmp2(max(8,atmp%m+2,nztmp+2)),rtmp(atmp%m),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + j = 1 + atmp%ia2(1) = 1 + do i=1, atmp%m + ir = p%perm(i) + + if (ir <= a%m ) then + + nzl = a%ia2(ir+1) - a%ia2(ir) + if (nzl > size(rtmp)) then + call psb_realloc(nzl,rtmp,info) + if(info/=0) then + info=4010 + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + endif + jj = a%ia2(ir) + k=0 + do kk=1, nzl + if (a%ia1(jj+kk-1)<=atmp%m) then + k = k + 1 + rtmp(k) = a%aspk(jj+kk-1) + atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) + endif + enddo + call isrx(k,atmp%ia1(j:j+k-1),itmp2) + do kk=1,k + atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) + enddo + + else if (ir <= atmp%m ) then + + ir = ir - a%m + nzl = blck%ia2(ir+1) - blck%ia2(ir) + if (nzl > size(rtmp)) then + call psb_realloc(nzl,rtmp,info) + if(info/=0) then + info=4010 + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + endif + jj = blck%ia2(ir) + k=0 + do kk=1, nzl + if (blck%ia1(jj+kk-1)<=atmp%m) then + k = k + 1 + rtmp(k) = blck%aspk(jj+kk-1) + atmp%ia1(j+k-1) = p%invperm(blck%ia1(jj+kk-1)) + endif + enddo + call isrx(k,atmp%ia1(j:j+k-1),itmp2) + do kk=1,k + atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) + enddo + + else + write(0,*) 'Row index error 1 :',i,ir + endif + + j = j + k + atmp%ia2(i+1) = j + + enddo + + t4 = mpi_wtime() + + + + deallocate(itmp,itmp2,rtmp) + + 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 + +contains + + + subroutine gps_reduction(m,ia,ja,perm,iperm,info) + integer i,j,dgConn,Npnt,m + integer n,idpth,ideg,ibw2,ipf2 + integer,dimension(:) :: perm,iperm,ia,ja + integer, intent(out) :: info + + integer,dimension(:,:),allocatable::NDstk + integer,dimension(:),allocatable::iOld,renum,ndeg,lvl,lvls1,lvls2,ccstor + !--- Per la common area. + + common /gra/ n,iDpth,iDeg + + character(len=20) :: name, ch_err + + if(psb_get_errstatus().ne.0) return + info=0 + name='gps_reduction' + call psb_erractionsave(err_act) + + + !--- Calcolo il massimo grado di connettivita'. + npnt = m + write(6,*) ' GPS su ',npnt + dgConn=0 + do i=1,m + dgconn = max(dgconn,(ia(i+1)-ia(i))) + enddo + !--- Il max valore di connettivita' e "dgConn" + + !--- Valori della common + n=Npnt !--- Numero di righe + iDeg=dgConn !--- Massima connettivita' + ! iDpth= !--- Numero di livelli non serve settarlo + + allocate(NDstk(Npnt,dgConn),stat=info) + if (info/=0) then + info=4000 + call psb_errpush(info,name) + goto 9999 + else + write(0,*) 'gps_reduction first alloc OK' + endif + allocate(iOld(Npnt),renum(Npnt+1),ndeg(Npnt),lvl(Npnt),lvls1(Npnt),& + &lvls2(Npnt),ccstor(Npnt),stat=info) + if (info/=0) then + info=4000 + call psb_errpush(info,name) + goto 9999 + else + write(0,*) 'gps_reduction 2nd alloc OK' + endif + + !--- Prepariamo il grafo della matrice + Ndstk(:,:)=0 + do i=1,Npnt + k=0 + do j = ia(i),ia(i+1) - 1 + if ((1<=ja(j)).and.( ja( j ) /= i ).and.(ja(j)<=npnt)) then + k = k+1 + Ndstk(i,k)=ja(j) + endif + enddo + ndeg(i)=k + enddo + + !--- Numerazione. + do i=1,Npnt + iOld(i)=i + enddo + write(0,*) 'gps_red : Preparation done' + !--- + !--- Chiamiamo funzione reduce. + call reduce(Ndstk,Npnt,iOld,renum,ndeg,lvl,lvls1, lvls2,ccstor,ibw2,ipf2) + write(0,*) 'gps_red : Done reduce' + !--- Permutazione + perm(1:Npnt)=renum(1:Npnt) + !--- Inversa permutazione + do i=1,Npnt + iperm(perm(i))=i + enddo + !--- Puliamo tutto. + deallocate(NDstk,iOld,renum,ndeg,lvl,lvls1,lvls2,ccstor) + + 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 gps_reduction + +end subroutine psb_dsp_renum diff --git a/src/prec/psb_dumf_bld.f90 b/src/prec/psb_dumf_bld.f90 new file mode 100644 index 00000000..dd1e3f2f --- /dev/null +++ b/src/prec/psb_dumf_bld.f90 @@ -0,0 +1,215 @@ +!!$ +!!$ +!!$ MPcube: Multilevel Parallel Preconditioners Package +!!$ for +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Daniela Di Serafino II University of Naples +!!$ Pasqua D'Ambra ICAR-CNR +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MPCUBE group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ + +subroutine psb_dumf_bld(a,desc_a,p,info) + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_tools_mod + 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_asmatbld + Subroutine psb_dasmatbld(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_dasmatbld + 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_asmatbld(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_asmatbld' + 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 asmatbld',info,nzb,blck%fida + call blacs_barrier(icontxt,'All') + endif + if (nzb > 0 ) then + if (size(atmp%aspk)