From 2ed8b1ac0fce2ae2367fd13a651a51af32c2c429 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 31 Jul 2007 09:51:43 +0000 Subject: [PATCH] Unpacked smoothed and raw aggregation assembly routines. --- mld_daggrmat_asb.F90 | 873 +------------------------------------ mld_daggrmat_raw_asb.F90 | 241 +++++++++++ mld_daggrmat_smth_asb.F90 | 709 ++++++++++++++++++++++++++++++ mld_prec_mod.f90 | 46 ++ mld_zaggrmat_asb.F90 | 876 +------------------------------------- mld_zaggrmat_raw_asb.F90 | 242 +++++++++++ mld_zaggrmat_smth_asb.F90 | 709 ++++++++++++++++++++++++++++++ 7 files changed, 1955 insertions(+), 1741 deletions(-) create mode 100644 mld_daggrmat_raw_asb.F90 create mode 100644 mld_daggrmat_smth_asb.F90 create mode 100644 mld_zaggrmat_raw_asb.F90 create mode 100644 mld_zaggrmat_smth_asb.F90 diff --git a/mld_daggrmat_asb.F90 b/mld_daggrmat_asb.F90 index 000cb6e0..13ee41e4 100644 --- a/mld_daggrmat_asb.F90 +++ b/mld_daggrmat_asb.F90 @@ -51,21 +51,20 @@ subroutine mld_daggrmat_asb(a,desc_a,ac,desc_ac,p,info) integer ::ictxt,np,me, err_act, icomm character(len=20) :: name, ch_err - name='mld_daggrmat_asb' + name='mld_aggrmat_asb' if(psb_get_errstatus().ne.0) return info=0 call psb_erractionsave(err_act) ictxt = psb_cd_get_context(desc_a) icomm = psb_cd_get_mpic(desc_a) - ictxt=psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) select case (p%iprcparm(aggr_kind_)) case (no_smooth_) - call raw_aggregate(info) + call mld_aggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info) if(info /= 0) then call psb_errpush(4010,name,a_err='raw_aggregate') @@ -75,7 +74,7 @@ subroutine mld_daggrmat_asb(a,desc_a,ac,desc_ac,p,info) case(tent_prol,biz_prol_) if (aggr_dump) call psb_csprt(70+me,a,head='% Input matrix') - call smooth_aggregate(info) + call mld_aggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info) if(info /= 0) then call psb_errpush(4010,name,a_err='smooth_aggregate') @@ -99,870 +98,4 @@ subroutine mld_daggrmat_asb(a,desc_a,ac,desc_ac,p,info) end if return -contains - - subroutine raw_aggregate(info) - use psb_base_mod - use mld_prec_type -#ifdef MPI_MOD - use mpi -#endif - implicit none -#ifdef MPI_H - include 'mpif.h' -#endif - - integer, intent(out) :: info - type(psb_dspmat_type) :: b - integer, pointer :: nzbr(:), idisp(:) - integer :: ictxt, nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& - & naggr, np, me, nzt,jl,nzl,nlr,& - & icomm,naggrm1, i, j, k, err_act - - name='raw_aggregate' - if(psb_get_errstatus().ne.0) return - info=0 - call psb_erractionsave(err_act) - - call psb_nullify_sp(b) - - ictxt = psb_cd_get_context(desc_a) - icomm = psb_cd_get_mpic(desc_a) - call psb_info(ictxt, me, np) - nglob = psb_cd_get_global_rows(desc_a) - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - - naggr = p%nlaggr(me+1) - ntaggr = sum(p%nlaggr) - allocate(nzbr(np), idisp(np),stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - naggrm1=sum(p%nlaggr(1:me)) - - if (p%iprcparm(coarse_mat_) == repl_mat_) then - do i=1, nrow - p%mlia(i) = p%mlia(i) + naggrm1 - end do - call psb_halo(p%mlia,desc_a,info) - end if - - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_halo') - goto 9999 - end if - - nzt = psb_sp_get_nnzeros(a) - - call psb_sp_all(b,nzt,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spall') - goto 9999 - end if - - call psb_sp_setifld(psb_dupl_ovwrt_,psb_dupl_,b,info) - call psb_sp_setifld(psb_upd_dflt_,psb_upd_,b,info) - b%fida = 'COO' - b%m=a%m - b%k=a%k - call psb_csdp(a,b,info) - if(info /= 0) then - info=4010 - ch_err='psb_csdp' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - nzt = psb_sp_get_nnzeros(b) - - j = 0 - do i=1, nzt - if ((1<=b%ia2(i)).and.(b%ia2(i)<=nrow)) then - j = j + 1 - b%aspk(j) = b%aspk(i) - b%ia1(j) = p%mlia(b%ia1(i)) - b%ia2(j) = p%mlia(b%ia2(i)) - end if - enddo - b%infoa(psb_nnz_)=j - call psb_fixcoo(b,info) - - nzt = psb_sp_get_nnzeros(b) - - call psb_sp_reall(b,nzt,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spreall') - goto 9999 - end if - b%m = naggr - b%k = naggr - - if (p%iprcparm(coarse_mat_) == repl_mat_) then - - call psb_cdrep(ntaggr,ictxt,desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdrep') - goto 9999 - end if - - nzbr(:) = 0 - nzbr(me+1) = nzt - call psb_sum(ictxt,nzbr(1:np)) - nzac = sum(nzbr) - call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spall') - goto 9999 - end if - - do ip=1,np - idisp(ip) = sum(nzbr(1:ip-1)) - enddo - ndx = nzbr(me+1) - - call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,& - & mpi_double_precision,icomm,info) - call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& - & mpi_integer,icomm,info) - call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& - & mpi_integer,icomm,info) - if(info /= 0) then - info=-1 - call psb_errpush(info,name) - goto 9999 - end if - - ac%m = ntaggr - ac%k = ntaggr - ac%infoa(psb_nnz_) = nzac - ac%fida='COO' - ac%descra='G' - call psb_fixcoo(ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_free') - goto 9999 - end if - - else if (p%iprcparm(coarse_mat_) == distr_mat_) then - - call psb_cdall(ictxt,desc_ac,info,nl=naggr) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdall') - goto 9999 - end if - call psb_cdasb(desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdasb') - goto 9999 - end if - - call psb_sp_clone(b,ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spclone') - goto 9999 - end if - call psb_sp_free(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_free') - goto 9999 - end if - - else - - write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(coarse_mat_) - end if - - call psb_ipcoo2csr(ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='ipcoo2csr') - goto 9999 - end if - - deallocate(nzbr,idisp) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine raw_aggregate - - - - subroutine smooth_aggregate(info) - use psb_base_mod - use mld_prec_type -#ifdef MPI_MOD - use mpi -#endif - implicit none -#ifdef MPI_H - include 'mpif.h' -#endif - - integer, intent(out) :: info - - type(psb_dspmat_type) :: b - integer, pointer :: nzbr(:), idisp(:) - integer :: ictxt, nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& - & naggr, np, me, & - & icomm, naggrm1,naggrp1,i,j,err_act,k,nzl - type(psb_dspmat_type), pointer :: am1,am2 - type(psb_dspmat_type) :: am3,am4 - logical :: ml_global_nmb - - logical, parameter :: test_dump=.false.,debug=.false. - integer, parameter :: ncmax=16 - real(kind(1.d0)) :: omega, anorm, tmp, dg - character(len=20) :: name - - - name='smooth_aggregate' - if(psb_get_errstatus().ne.0) return - info=0 - call psb_erractionsave(err_act) - - ictxt = psb_cd_get_context(desc_a) - icomm = psb_cd_get_mpic(desc_a) - - call psb_info(ictxt, me, np) - - call psb_nullify_sp(b) - call psb_nullify_sp(am3) - call psb_nullify_sp(am4) - - am2 => p%av(sm_pr_t_) - am1 => p%av(sm_pr_) - - nglob = psb_cd_get_global_rows(desc_a) - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - - naggr = p%nlaggr(me+1) - ntaggr = sum(p%nlaggr) - - allocate(nzbr(np), idisp(np),stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - - naggrm1 = sum(p%nlaggr(1:me)) - naggrp1 = sum(p%nlaggr(1:me+1)) - - ml_global_nmb = ( (p%iprcparm(aggr_kind_) == tent_prol).or.& - & ( (p%iprcparm(aggr_kind_) == biz_prol_).and.& - & (p%iprcparm(coarse_mat_) == repl_mat_)) ) - - - if (ml_global_nmb) then - p%mlia(1:nrow) = p%mlia(1:nrow) + naggrm1 - call psb_halo(p%mlia,desc_a,info) - - if(info /= 0) then - call psb_errpush(4010,name,a_err='f90_pshalo') - goto 9999 - end if - end if - - if (aggr_dump) then - open(30+me) - write(30+me,*) '% Aggregation map' - do i=1,ncol - write(30+me,*) i,p%mlia(i) - end do - close(30+me) - end if - - ! naggr: number of local aggregates - ! nrow: local rows. - ! - allocate(p%dorig(nrow),stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),& - & a_err='real(kind(1.d0))') - goto 9999 - end if - - ! Get diagonal D - call psb_sp_getdiag(a,p%dorig,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_getdiag') - goto 9999 - end if - - do i=1,size(p%dorig) - if (p%dorig(i) /= dzero) then - p%dorig(i) = done / p%dorig(i) - else - p%dorig(i) = done - end if - end do - - ! where (p%dorig /= dzero) - ! p%dorig = done / p%dorig - ! elsewhere - ! p%dorig = done - ! end where - - - ! 1. Allocate Ptilde in sparse matrix form - am4%fida='COO' - am4%m=ncol - if (ml_global_nmb) then - am4%k=ntaggr - call psb_sp_all(ncol,ntaggr,am4,ncol,info) - else - am4%k=naggr - call psb_sp_all(ncol,naggr,am4,ncol,info) - endif - if(info /= 0) then - call psb_errpush(4010,name,a_err='spall') - goto 9999 - end if - - if (ml_global_nmb) then - do i=1,ncol - am4%aspk(i) = done - am4%ia1(i) = i - am4%ia2(i) = p%mlia(i) - end do - am4%infoa(psb_nnz_) = ncol - else - do i=1,nrow - am4%aspk(i) = done - am4%ia1(i) = i - am4%ia2(i) = p%mlia(i) - end do - am4%infoa(psb_nnz_) = nrow - endif - - - - - call psb_ipcoo2csr(am4,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='ipcoo2csr') - goto 9999 - end if - - call psb_sp_clone(a,am3,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spclone') - goto 9999 - end if - - ! - ! WARNING: the cycles below assume that AM3 does have - ! its diagonal elements stored explicitly!!! - ! Should we switch to something safer? - ! - call psb_sp_scal(am3,p%dorig,info) - if(info /= 0) goto 9999 - - if (p%iprcparm(aggr_eig_) == max_norm_) then - - if (p%iprcparm(aggr_kind_) == biz_prol_) then - - ! - ! This only works with CSR. - ! - anorm = dzero - dg = done - do i=1,am3%m - tmp = dzero - do j=am3%ia2(i),am3%ia2(i+1)-1 - if (am3%ia1(j) <= am3%m) then - tmp = tmp + dabs(am3%aspk(j)) - endif - if (am3%ia1(j) == i ) then - dg = dabs(am3%aspk(j)) - end if - end do - anorm = max(anorm,tmp/dg) - enddo - - call psb_amx(ictxt,anorm) - else - anorm = psb_spnrmi(am3,desc_a,info) - endif - omega = 4.d0/(3.d0*anorm) - p%dprcparm(aggr_damp_) = omega - - else if (p%iprcparm(aggr_eig_) == user_choice_) then - - omega = p%dprcparm(aggr_damp_) - - else if (p%iprcparm(aggr_eig_) /= user_choice_) then - write(0,*) me,'Error: invalid choice for OMEGA in blaggrmat?? ',& - & p%iprcparm(aggr_eig_) - end if - - - if (am3%fida=='CSR') then - do i=1,am3%m - do j=am3%ia2(i),am3%ia2(i+1)-1 - if (am3%ia1(j) == i) then - am3%aspk(j) = done - omega*am3%aspk(j) - else - am3%aspk(j) = - omega*am3%aspk(j) - end if - end do - end do - else if (am3%fida=='COO') then - do j=1,am3%infoa(psb_nnz_) - if (am3%ia1(j) /= am3%ia2(j)) then - am3%aspk(j) = - omega*am3%aspk(j) - else - am3%aspk(j) = done - omega*am3%aspk(j) - endif - end do - call psb_ipcoo2csr(am3,info) - else - write(0,*) 'Missing implementation of I sum' - call psb_errpush(4010,name) - goto 9999 - end if - - if (test_dump) then - open(30+me) - write(30+me,*) 'OMEGA: ',omega - do i=1,size(p%dorig) - write(30+me,*) p%dorig(i) - end do - close(30+me) - end if - - if (test_dump) call & - & psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob) - if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,& - & ivc=desc_a%loc_to_glob) - if (debug) write(0,*) me,'Done gather, going for SYMBMM 1' - ! - ! Symbmm90 does the allocation for its result. - ! - ! am1 = (i-wDA)Ptilde - ! Doing it this way means to consider diag(Ai) - ! - ! - call psb_symbmm(am3,am4,am1,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='symbmm 1') - goto 9999 - end if - - call psb_numbmm(am3,am4,am1) - - if (debug) write(0,*) me,'Done NUMBMM 1' - - call psb_sp_free(am4,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_free') - goto 9999 - end if - - if (ml_global_nmb) then - ! - ! Now we have to gather the halo of am1, and add it to itself - ! to multiply it by A, - ! - call psb_sphalo(am1,desc_a,am4,info,clcnv=.false.) - - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sphalo') - goto 9999 - end if - - call psb_rwextd(ncol,am1,info,b=am4) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_rwextd') - goto 9999 - end if - - call psb_sp_free(am4,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - else - - call psb_rwextd(ncol,am1,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='rwextd') - goto 9999 - end if - endif - - if (test_dump) & - & call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob) - - call psb_symbmm(a,am1,am3,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='symbmm 2') - goto 9999 - end if - - call psb_numbmm(a,am1,am3) - if (debug) write(0,*) me,'Done NUMBMM 2' - - if (p%iprcparm(aggr_kind_) == tent_prol) then - call psb_transp(am1,am2,fmt='COO') - nzl = am2%infoa(psb_nnz_) - i=0 - ! - ! Now we have to fix this. The only rows of B that are correct - ! are those corresponding to "local" aggregates, i.e. indices in p%mlia(:) - ! - do k=1, nzl - if ((naggrm1 < am2%ia1(k)) .and.(am2%ia1(k) <= naggrp1)) then - i = i+1 - am2%aspk(i) = am2%aspk(k) - am2%ia1(i) = am2%ia1(k) - am2%ia2(i) = am2%ia2(k) - end if - end do - - am2%infoa(psb_nnz_) = i - call psb_ipcoo2csr(am2,info) - else - call psb_transp(am1,am2) - endif - if (debug) write(0,*) me,'starting sphalo/ rwxtd' - - if (p%iprcparm(aggr_kind_) == tent_prol) then - ! am2 = ((i-wDA)Ptilde)^T - call psb_sphalo(am3,desc_a,am4,info,clcnv=.false.) - - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sphalo') - goto 9999 - end if - call psb_rwextd(ncol,am3,info,b=am4) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_rwextd') - goto 9999 - end if - call psb_sp_free(am4,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - else if (p%iprcparm(aggr_kind_) == biz_prol_) then - - call psb_rwextd(ncol,am3,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_rwextd') - goto 9999 - end if - endif - - if (debug) write(0,*) me,'starting symbmm 3' - call psb_symbmm(am2,am3,b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='symbmm 3') - goto 9999 - end if - - if (debug) write(0,*) me,'starting numbmm 3' - call psb_numbmm(am2,am3,b) - if (debug) write(0,*) me,'Done NUMBMM 3' - -!!$ if (aggr_dump) call csprt(50+me,am1,head='% Operator PTrans.') - call psb_sp_free(am3,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - call psb_ipcsr2coo(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='ipcsr2coo') - goto 9999 - end if - - call psb_fixcoo(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='fixcoo') - goto 9999 - end if - - - if (test_dump) call psb_csprt(80+me,b,head='% Smoothed aggregate AC.') - - select case(p%iprcparm(aggr_kind_)) - - case(tent_prol) - - select case(p%iprcparm(coarse_mat_)) - - case(distr_mat_) - - call psb_sp_clone(b,ac,info) - if(info /= 0) goto 9999 - nzac = ac%infoa(psb_nnz_) - nzl = ac%infoa(psb_nnz_) - - call psb_cdall(ictxt,desc_ac,info,nl=p%nlaggr(me+1)) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdall') - goto 9999 - end if - - - call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdins') - goto 9999 - end if - - if (debug) write(0,*) me,'Created aux descr. distr.' - call psb_cdasb(desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdasb') - goto 9999 - end if - - - if (debug) write(0,*) me,'Asmbld aux descr. distr.' - - call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I') - if(info /= 0) then - call psb_errpush(4010,name,a_err='psglob_to_loc') - goto 9999 - end if - - - call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I') - if(info /= 0) then - call psb_errpush(4010,name,a_err='psglob_to_loc') - goto 9999 - end if - - - ac%m=desc_ac%matrix_data(psb_n_row_) - ac%k=desc_ac%matrix_data(psb_n_col_) - ac%fida='COO' - ac%descra='G' - - call psb_sp_free(b,info) - if (info == 0) deallocate(nzbr,idisp,stat=info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - if (np>1) then - nzl = psb_sp_get_nnzeros(am1) - call psb_glob_to_loc(am1%ia1(1:nzl),desc_ac,info,'I') - - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_glob_to_loc') - goto 9999 - end if - endif - am1%k=desc_ac%matrix_data(psb_n_col_) - - if (np>1) then - call psb_ipcsr2coo(am2,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_ipcsr2coo') - goto 9999 - end if - - nzl = am2%infoa(psb_nnz_) - call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I') - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_glob_to_loc') - goto 9999 - end if - - call psb_ipcoo2csr(am2,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_ipcoo2csr') - goto 9999 - end if - end if - am2%m=desc_ac%matrix_data(psb_n_col_) - - if (debug) write(0,*) me,'Done ac ' - case(repl_mat_) - ! - ! - call psb_cdrep(ntaggr,ictxt,desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdrep') - goto 9999 - end if - - nzbr(:) = 0 - nzbr(me+1) = b%infoa(psb_nnz_) - - call psb_sum(ictxt,nzbr(1:np)) - nzac = sum(nzbr) - call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) - if(info /= 0) goto 9999 - - - do ip=1,np - idisp(ip) = sum(nzbr(1:ip-1)) - enddo - ndx = nzbr(me+1) - - call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,& - & mpi_double_precision,icomm,info) - call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& - & mpi_integer,icomm,info) - call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& - & mpi_integer,icomm,info) - if(info /= 0) goto 9999 - - - ac%m = ntaggr - ac%k = ntaggr - ac%infoa(psb_nnz_) = nzac - ac%fida='COO' - ac%descra='G' - call psb_fixcoo(ac,info) - if(info /= 0) goto 9999 - call psb_sp_free(b,info) - if(info /= 0) goto 9999 - if (me==0) then - if (test_dump) call psb_csprt(80+me,ac,head='% Smoothed aggregate AC.') - endif - - deallocate(nzbr,idisp) - - case default - write(0,*) 'Inconsistent input in smooth_new_aggregate' - end select - - - case(biz_prol_) - - select case(p%iprcparm(coarse_mat_)) - - case(distr_mat_) - - call psb_sp_clone(b,ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spclone') - goto 9999 - end if - call psb_cdall(ictxt,desc_ac,info,nl=naggr) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdall') - goto 9999 - end if - - call psb_cdasb(desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdasb') - goto 9999 - end if - - call psb_sp_free(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_free') - goto 9999 - end if - - - case(repl_mat_) - ! - ! - - call psb_cdrep(ntaggr,ictxt,desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdrep') - goto 9999 - end if - - nzbr(:) = 0 - nzbr(me+1) = b%infoa(psb_nnz_) - call psb_sum(ictxt,nzbr(1:np)) - nzac = sum(nzbr) - call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_all') - goto 9999 - end if - - do ip=1,np - idisp(ip) = sum(nzbr(1:ip-1)) - enddo - ndx = nzbr(me+1) - - call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,& - & mpi_double_precision,icomm,info) - call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& - & mpi_integer,icomm,info) - call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& - & mpi_integer,icomm,info) - if(info /= 0) then - info=-1 - call psb_errpush(info,name) - goto 9999 - end if - - - ac%m = ntaggr - ac%k = ntaggr - ac%infoa(psb_nnz_) = nzac - ac%fida='COO' - ac%descra='G' - call psb_fixcoo(ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_fixcoo') - goto 9999 - end if - call psb_sp_free(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - end select - deallocate(nzbr,idisp) - - end select - - call psb_ipcoo2csr(ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_ipcoo2csr') - goto 9999 - end if - - if (debug) write(0,*) me,'Done smooth_aggregate ' - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name) - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - - - end subroutine smooth_aggregate - end subroutine mld_daggrmat_asb diff --git a/mld_daggrmat_raw_asb.F90 b/mld_daggrmat_raw_asb.F90 new file mode 100644 index 00000000..de919f89 --- /dev/null +++ b/mld_daggrmat_raw_asb.F90 @@ -0,0 +1,241 @@ +!!$ +!!$ +!!$ MD2P4 +!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS +!!$ 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 Second 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 MD2P4 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 MD2P4 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 mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info) + use psb_base_mod + use mld_prec_mod, mld_protect_name => mld_daggrmat_raw_asb + +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + type(psb_dspmat_type), intent(in), target :: a + type(mld_dbaseprc_type), intent(inout), target :: p + type(psb_dspmat_type), intent(inout), target :: ac + type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_ac + integer, intent(out) :: info + + logical, parameter :: aggr_dump=.false. + integer ::ictxt,np,me, err_act, icomm + character(len=20) :: name, ch_err + + type(psb_dspmat_type) :: b + integer, pointer :: nzbr(:), idisp(:) + integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& + & naggr, nzt,jl,nzl,nlr, naggrm1, i, j, k + + name='mld_aggrmat_raw_asb' + if(psb_get_errstatus().ne.0) return + info=0 + call psb_erractionsave(err_act) + + call psb_nullify_sp(b) + + ictxt = psb_cd_get_context(desc_a) + icomm = psb_cd_get_mpic(desc_a) + call psb_info(ictxt, me, np) + nglob = psb_cd_get_global_rows(desc_a) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + naggr = p%nlaggr(me+1) + ntaggr = sum(p%nlaggr) + allocate(nzbr(np), idisp(np),stat=info) + if (info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + naggrm1=sum(p%nlaggr(1:me)) + + if (p%iprcparm(coarse_mat_) == repl_mat_) then + do i=1, nrow + p%mlia(i) = p%mlia(i) + naggrm1 + end do + call psb_halo(p%mlia,desc_a,info) + end if + + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_halo') + goto 9999 + end if + + nzt = psb_sp_get_nnzeros(a) + + call psb_sp_all(b,nzt,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spall') + goto 9999 + end if + + call psb_sp_setifld(psb_dupl_ovwrt_,psb_dupl_,b,info) + call psb_sp_setifld(psb_upd_dflt_,psb_upd_,b,info) + b%fida = 'COO' + b%m=a%m + b%k=a%k + call psb_csdp(a,b,info) + if(info /= 0) then + info=4010 + ch_err='psb_csdp' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + nzt = psb_sp_get_nnzeros(b) + + j = 0 + do i=1, nzt + if ((1<=b%ia2(i)).and.(b%ia2(i)<=nrow)) then + j = j + 1 + b%aspk(j) = b%aspk(i) + b%ia1(j) = p%mlia(b%ia1(i)) + b%ia2(j) = p%mlia(b%ia2(i)) + end if + enddo + b%infoa(psb_nnz_)=j + call psb_fixcoo(b,info) + + nzt = psb_sp_get_nnzeros(b) + + call psb_sp_reall(b,nzt,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spreall') + goto 9999 + end if + b%m = naggr + b%k = naggr + + if (p%iprcparm(coarse_mat_) == repl_mat_) then + + call psb_cdrep(ntaggr,ictxt,desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdrep') + goto 9999 + end if + + nzbr(:) = 0 + nzbr(me+1) = nzt + call psb_sum(ictxt,nzbr(1:np)) + nzac = sum(nzbr) + call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spall') + goto 9999 + end if + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + + call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,& + & mpi_double_precision,icomm,info) + call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& + & mpi_integer,icomm,info) + call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& + & mpi_integer,icomm,info) + if(info /= 0) then + info=-1 + call psb_errpush(info,name) + goto 9999 + end if + + ac%m = ntaggr + ac%k = ntaggr + ac%infoa(psb_nnz_) = nzac + ac%fida='COO' + ac%descra='G' + call psb_fixcoo(ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_free') + goto 9999 + end if + + else if (p%iprcparm(coarse_mat_) == distr_mat_) then + + call psb_cdall(ictxt,desc_ac,info,nl=naggr) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdall') + goto 9999 + end if + call psb_cdasb(desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdasb') + goto 9999 + end if + + call psb_sp_clone(b,ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spclone') + goto 9999 + end if + call psb_sp_free(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_free') + goto 9999 + end if + + else + + write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(coarse_mat_) + end if + + call psb_ipcoo2csr(ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='ipcoo2csr') + goto 9999 + end if + + deallocate(nzbr,idisp) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_daggrmat_raw_asb diff --git a/mld_daggrmat_smth_asb.F90 b/mld_daggrmat_smth_asb.F90 new file mode 100644 index 00000000..caedd91b --- /dev/null +++ b/mld_daggrmat_smth_asb.F90 @@ -0,0 +1,709 @@ +!!$ +!!$ +!!$ MD2P4 +!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS +!!$ 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 Second 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 MD2P4 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 MD2P4 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 mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info) + use psb_base_mod + use mld_prec_mod, mld_protect_name => mld_daggrmat_smth_asb + +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + type(psb_dspmat_type), intent(in), target :: a + type(mld_dbaseprc_type), intent(inout), target :: p + type(psb_dspmat_type), intent(inout), target :: ac + type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_ac + integer, intent(out) :: info + + type(psb_dspmat_type) :: b + integer, pointer :: nzbr(:), idisp(:) + integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& + & naggr, nzt,jl,nzl,nlr,naggrm1,naggrp1, i, j, k + logical, parameter :: aggr_dump=.false. + integer ::ictxt,np,me, err_act, icomm + character(len=20) :: name, ch_err + type(psb_dspmat_type), pointer :: am1,am2 + type(psb_dspmat_type) :: am3,am4 + logical :: ml_global_nmb + + logical, parameter :: test_dump=.false.,debug=.false. + integer, parameter :: ncmax=16 + real(kind(1.d0)) :: omega, anorm, tmp, dg + + name='mld_aggrmat_smth_asb' + if(psb_get_errstatus().ne.0) return + info=0 + call psb_erractionsave(err_act) + + ictxt = psb_cd_get_context(desc_a) + icomm = psb_cd_get_mpic(desc_a) + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + + call psb_nullify_sp(b) + call psb_nullify_sp(am3) + call psb_nullify_sp(am4) + + am2 => p%av(sm_pr_t_) + am1 => p%av(sm_pr_) + + nglob = psb_cd_get_global_rows(desc_a) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + naggr = p%nlaggr(me+1) + ntaggr = sum(p%nlaggr) + + allocate(nzbr(np), idisp(np),stat=info) + if (info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + + naggrm1 = sum(p%nlaggr(1:me)) + naggrp1 = sum(p%nlaggr(1:me+1)) + + ml_global_nmb = ( (p%iprcparm(aggr_kind_) == tent_prol).or.& + & ( (p%iprcparm(aggr_kind_) == biz_prol_).and.& + & (p%iprcparm(coarse_mat_) == repl_mat_)) ) + + + if (ml_global_nmb) then + p%mlia(1:nrow) = p%mlia(1:nrow) + naggrm1 + call psb_halo(p%mlia,desc_a,info) + + if(info /= 0) then + call psb_errpush(4010,name,a_err='f90_pshalo') + goto 9999 + end if + end if + + if (aggr_dump) then + open(30+me) + write(30+me,*) '% Aggregation map' + do i=1,ncol + write(30+me,*) i,p%mlia(i) + end do + close(30+me) + end if + + ! naggr: number of local aggregates + ! nrow: local rows. + ! + allocate(p%dorig(nrow),stat=info) + if (info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),& + & a_err='real(kind(1.d0))') + goto 9999 + end if + + ! Get diagonal D + call psb_sp_getdiag(a,p%dorig,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_getdiag') + goto 9999 + end if + + do i=1,size(p%dorig) + if (p%dorig(i) /= dzero) then + p%dorig(i) = done / p%dorig(i) + else + p%dorig(i) = done + end if + end do + + ! where (p%dorig /= dzero) + ! p%dorig = done / p%dorig + ! elsewhere + ! p%dorig = done + ! end where + + + ! 1. Allocate Ptilde in sparse matrix form + am4%fida='COO' + am4%m=ncol + if (ml_global_nmb) then + am4%k=ntaggr + call psb_sp_all(ncol,ntaggr,am4,ncol,info) + else + am4%k=naggr + call psb_sp_all(ncol,naggr,am4,ncol,info) + endif + if(info /= 0) then + call psb_errpush(4010,name,a_err='spall') + goto 9999 + end if + + if (ml_global_nmb) then + do i=1,ncol + am4%aspk(i) = done + am4%ia1(i) = i + am4%ia2(i) = p%mlia(i) + end do + am4%infoa(psb_nnz_) = ncol + else + do i=1,nrow + am4%aspk(i) = done + am4%ia1(i) = i + am4%ia2(i) = p%mlia(i) + end do + am4%infoa(psb_nnz_) = nrow + endif + + + + + call psb_ipcoo2csr(am4,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='ipcoo2csr') + goto 9999 + end if + + call psb_sp_clone(a,am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spclone') + goto 9999 + end if + + ! + ! WARNING: the cycles below assume that AM3 does have + ! its diagonal elements stored explicitly!!! + ! Should we switch to something safer? + ! + call psb_sp_scal(am3,p%dorig,info) + if(info /= 0) goto 9999 + + if (p%iprcparm(aggr_eig_) == max_norm_) then + + if (p%iprcparm(aggr_kind_) == biz_prol_) then + + ! + ! This only works with CSR. + ! + anorm = dzero + dg = done + do i=1,am3%m + tmp = dzero + do j=am3%ia2(i),am3%ia2(i+1)-1 + if (am3%ia1(j) <= am3%m) then + tmp = tmp + dabs(am3%aspk(j)) + endif + if (am3%ia1(j) == i ) then + dg = dabs(am3%aspk(j)) + end if + end do + anorm = max(anorm,tmp/dg) + enddo + + call psb_amx(ictxt,anorm) + else + anorm = psb_spnrmi(am3,desc_a,info) + endif + omega = 4.d0/(3.d0*anorm) + p%dprcparm(aggr_damp_) = omega + + else if (p%iprcparm(aggr_eig_) == user_choice_) then + + omega = p%dprcparm(aggr_damp_) + + else if (p%iprcparm(aggr_eig_) /= user_choice_) then + write(0,*) me,'Error: invalid choice for OMEGA in blaggrmat?? ',& + & p%iprcparm(aggr_eig_) + end if + + + if (am3%fida=='CSR') then + do i=1,am3%m + do j=am3%ia2(i),am3%ia2(i+1)-1 + if (am3%ia1(j) == i) then + am3%aspk(j) = done - omega*am3%aspk(j) + else + am3%aspk(j) = - omega*am3%aspk(j) + end if + end do + end do + else if (am3%fida=='COO') then + do j=1,am3%infoa(psb_nnz_) + if (am3%ia1(j) /= am3%ia2(j)) then + am3%aspk(j) = - omega*am3%aspk(j) + else + am3%aspk(j) = done - omega*am3%aspk(j) + endif + end do + call psb_ipcoo2csr(am3,info) + else + write(0,*) 'Missing implementation of I sum' + call psb_errpush(4010,name) + goto 9999 + end if + + if (test_dump) then + open(30+me) + write(30+me,*) 'OMEGA: ',omega + do i=1,size(p%dorig) + write(30+me,*) p%dorig(i) + end do + close(30+me) + end if + + if (test_dump) call & + & psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob) + if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,& + & ivc=desc_a%loc_to_glob) + if (debug) write(0,*) me,'Done gather, going for SYMBMM 1' + ! + ! Symbmm90 does the allocation for its result. + ! + ! am1 = (i-wDA)Ptilde + ! Doing it this way means to consider diag(Ai) + ! + ! + call psb_symbmm(am3,am4,am1,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm 1') + goto 9999 + end if + + call psb_numbmm(am3,am4,am1) + + if (debug) write(0,*) me,'Done NUMBMM 1' + + call psb_sp_free(am4,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_free') + goto 9999 + end if + + if (ml_global_nmb) then + ! + ! Now we have to gather the halo of am1, and add it to itself + ! to multiply it by A, + ! + call psb_sphalo(am1,desc_a,am4,info,clcnv=.false.) + + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sphalo') + goto 9999 + end if + + call psb_rwextd(ncol,am1,info,b=am4) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_rwextd') + goto 9999 + end if + + call psb_sp_free(am4,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + else + + call psb_rwextd(ncol,am1,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='rwextd') + goto 9999 + end if + endif + + if (test_dump) & + & call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob) + + call psb_symbmm(a,am1,am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm 2') + goto 9999 + end if + + call psb_numbmm(a,am1,am3) + if (debug) write(0,*) me,'Done NUMBMM 2' + + if (p%iprcparm(aggr_kind_) == tent_prol) then + call psb_transp(am1,am2,fmt='COO') + nzl = am2%infoa(psb_nnz_) + i=0 + ! + ! Now we have to fix this. The only rows of B that are correct + ! are those corresponding to "local" aggregates, i.e. indices in p%mlia(:) + ! + do k=1, nzl + if ((naggrm1 < am2%ia1(k)) .and.(am2%ia1(k) <= naggrp1)) then + i = i+1 + am2%aspk(i) = am2%aspk(k) + am2%ia1(i) = am2%ia1(k) + am2%ia2(i) = am2%ia2(k) + end if + end do + + am2%infoa(psb_nnz_) = i + call psb_ipcoo2csr(am2,info) + else + call psb_transp(am1,am2) + endif + if (debug) write(0,*) me,'starting sphalo/ rwxtd' + + if (p%iprcparm(aggr_kind_) == tent_prol) then + ! am2 = ((i-wDA)Ptilde)^T + call psb_sphalo(am3,desc_a,am4,info,clcnv=.false.) + + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sphalo') + goto 9999 + end if + call psb_rwextd(ncol,am3,info,b=am4) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_rwextd') + goto 9999 + end if + call psb_sp_free(am4,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + else if (p%iprcparm(aggr_kind_) == biz_prol_) then + + call psb_rwextd(ncol,am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_rwextd') + goto 9999 + end if + endif + + if (debug) write(0,*) me,'starting symbmm 3' + call psb_symbmm(am2,am3,b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm 3') + goto 9999 + end if + + if (debug) write(0,*) me,'starting numbmm 3' + call psb_numbmm(am2,am3,b) + if (debug) write(0,*) me,'Done NUMBMM 3' + +!!$ if (aggr_dump) call csprt(50+me,am1,head='% Operator PTrans.') + call psb_sp_free(am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + call psb_ipcsr2coo(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='ipcsr2coo') + goto 9999 + end if + + call psb_fixcoo(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='fixcoo') + goto 9999 + end if + + + if (test_dump) call psb_csprt(80+me,b,head='% Smoothed aggregate AC.') + + select case(p%iprcparm(aggr_kind_)) + + case(tent_prol) + + select case(p%iprcparm(coarse_mat_)) + + case(distr_mat_) + + call psb_sp_clone(b,ac,info) + if(info /= 0) goto 9999 + nzac = ac%infoa(psb_nnz_) + nzl = ac%infoa(psb_nnz_) + + call psb_cdall(ictxt,desc_ac,info,nl=p%nlaggr(me+1)) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdall') + goto 9999 + end if + + + call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdins') + goto 9999 + end if + + if (debug) write(0,*) me,'Created aux descr. distr.' + call psb_cdasb(desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdasb') + goto 9999 + end if + + + if (debug) write(0,*) me,'Asmbld aux descr. distr.' + + call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I') + if(info /= 0) then + call psb_errpush(4010,name,a_err='psglob_to_loc') + goto 9999 + end if + + + call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I') + if(info /= 0) then + call psb_errpush(4010,name,a_err='psglob_to_loc') + goto 9999 + end if + + + ac%m=desc_ac%matrix_data(psb_n_row_) + ac%k=desc_ac%matrix_data(psb_n_col_) + ac%fida='COO' + ac%descra='G' + + call psb_sp_free(b,info) + if (info == 0) deallocate(nzbr,idisp,stat=info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + if (np>1) then + nzl = psb_sp_get_nnzeros(am1) + call psb_glob_to_loc(am1%ia1(1:nzl),desc_ac,info,'I') + + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_glob_to_loc') + goto 9999 + end if + endif + am1%k=desc_ac%matrix_data(psb_n_col_) + + if (np>1) then + call psb_ipcsr2coo(am2,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcsr2coo') + goto 9999 + end if + + nzl = am2%infoa(psb_nnz_) + call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I') + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_glob_to_loc') + goto 9999 + end if + + call psb_ipcoo2csr(am2,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcoo2csr') + goto 9999 + end if + end if + am2%m=desc_ac%matrix_data(psb_n_col_) + + if (debug) write(0,*) me,'Done ac ' + case(repl_mat_) + ! + ! + call psb_cdrep(ntaggr,ictxt,desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdrep') + goto 9999 + end if + + nzbr(:) = 0 + nzbr(me+1) = b%infoa(psb_nnz_) + + call psb_sum(ictxt,nzbr(1:np)) + nzac = sum(nzbr) + call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) + if(info /= 0) goto 9999 + + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + + call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,& + & mpi_double_precision,icomm,info) + call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& + & mpi_integer,icomm,info) + call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& + & mpi_integer,icomm,info) + if(info /= 0) goto 9999 + + + ac%m = ntaggr + ac%k = ntaggr + ac%infoa(psb_nnz_) = nzac + ac%fida='COO' + ac%descra='G' + call psb_fixcoo(ac,info) + if(info /= 0) goto 9999 + call psb_sp_free(b,info) + if(info /= 0) goto 9999 + if (me==0) then + if (test_dump) call psb_csprt(80+me,ac,head='% Smoothed aggregate AC.') + endif + + deallocate(nzbr,idisp) + + case default + write(0,*) 'Inconsistent input in smooth_new_aggregate' + end select + + + case(biz_prol_) + + select case(p%iprcparm(coarse_mat_)) + + case(distr_mat_) + + call psb_sp_clone(b,ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spclone') + goto 9999 + end if + call psb_cdall(ictxt,desc_ac,info,nl=naggr) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdall') + goto 9999 + end if + + call psb_cdasb(desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdasb') + goto 9999 + end if + + call psb_sp_free(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_free') + goto 9999 + end if + + + case(repl_mat_) + ! + ! + + call psb_cdrep(ntaggr,ictxt,desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdrep') + goto 9999 + end if + + nzbr(:) = 0 + nzbr(me+1) = b%infoa(psb_nnz_) + call psb_sum(ictxt,nzbr(1:np)) + nzac = sum(nzbr) + call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_all') + goto 9999 + end if + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + + call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,ac%aspk,nzbr,idisp,& + & mpi_double_precision,icomm,info) + call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& + & mpi_integer,icomm,info) + call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& + & mpi_integer,icomm,info) + if(info /= 0) then + info=-1 + call psb_errpush(info,name) + goto 9999 + end if + + + ac%m = ntaggr + ac%k = ntaggr + ac%infoa(psb_nnz_) = nzac + ac%fida='COO' + ac%descra='G' + call psb_fixcoo(ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_fixcoo') + goto 9999 + end if + call psb_sp_free(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + end select + deallocate(nzbr,idisp) + + end select + + call psb_ipcoo2csr(ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcoo2csr') + goto 9999 + end if + + if (debug) write(0,*) me,'Done smooth_aggregate ' + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + + + +end subroutine mld_daggrmat_smth_asb diff --git a/mld_prec_mod.f90 b/mld_prec_mod.f90 index 3c8d38e3..04282419 100644 --- a/mld_prec_mod.f90 +++ b/mld_prec_mod.f90 @@ -531,4 +531,50 @@ module mld_prec_mod end subroutine mld_zaggrmat_asb end interface + interface mld_aggrmat_raw_asb + subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info) + use psb_base_mod + use mld_prec_type + type(psb_dspmat_type), intent(in), target :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), intent(inout),target :: ac + type(psb_desc_type), intent(inout) :: desc_ac + type(mld_dbaseprc_type), intent(inout), target :: p + integer, intent(out) :: info + end subroutine mld_daggrmat_raw_asb + subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info) + use psb_base_mod + use mld_prec_type + type(psb_zspmat_type), intent(in), target :: a + type(mld_zbaseprc_type), intent(inout),target :: p + type(psb_zspmat_type), intent(inout),target :: ac + type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_ac + integer, intent(out) :: info + end subroutine mld_zaggrmat_raw_asb + end interface + + interface mld_aggrmat_smth_asb + subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info) + use psb_base_mod + use mld_prec_type + type(psb_dspmat_type), intent(in), target :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), intent(inout),target :: ac + type(psb_desc_type), intent(inout) :: desc_ac + type(mld_dbaseprc_type), intent(inout), target :: p + integer, intent(out) :: info + end subroutine mld_daggrmat_smth_asb + subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info) + use psb_base_mod + use mld_prec_type + type(psb_zspmat_type), intent(in), target :: a + type(mld_zbaseprc_type), intent(inout),target :: p + type(psb_zspmat_type), intent(inout),target :: ac + type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_ac + integer, intent(out) :: info + end subroutine mld_zaggrmat_smth_asb + end interface + end module mld_prec_mod diff --git a/mld_zaggrmat_asb.F90 b/mld_zaggrmat_asb.F90 index 631397d4..43428c89 100644 --- a/mld_zaggrmat_asb.F90 +++ b/mld_zaggrmat_asb.F90 @@ -51,21 +51,20 @@ subroutine mld_zaggrmat_asb(a,desc_a,ac,desc_ac,p,info) integer ::ictxt,np,me, err_act,icomm character(len=20) :: name, ch_err - name='mld_zaggrmat_asb' + name='mld_aggrmat_asb' if(psb_get_errstatus().ne.0) return info=0 call psb_erractionsave(err_act) ictxt = psb_cd_get_context(desc_a) icomm = psb_cd_get_mpic(desc_a) - ictxt=psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) select case (p%iprcparm(aggr_kind_)) case (no_smooth_) - call raw_aggregate(info) + call mld_aggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info) if(info /= 0) then call psb_errpush(4010,name,a_err='raw_aggregate') @@ -74,13 +73,14 @@ subroutine mld_zaggrmat_asb(a,desc_a,ac,desc_ac,p,info) if (aggr_dump) call psb_csprt(90+me,ac,head='% Raw aggregate.') case(tent_prol,biz_prol_) - - call smooth_aggregate(info) + if (aggr_dump) call psb_csprt(70+me,a,head='% Input matrix') + call mld_aggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info) if(info /= 0) then call psb_errpush(4010,name,a_err='smooth_aggregate') goto 9999 end if + if (aggr_dump) call psb_csprt(90+me,ac,head='% Smooth aggregate.') case default call psb_errpush(4010,name,a_err=name) goto 9999 @@ -98,870 +98,4 @@ subroutine mld_zaggrmat_asb(a,desc_a,ac,desc_ac,p,info) end if return -contains - - subroutine raw_aggregate(info) - use psb_base_mod - use mld_prec_type -#ifdef MPI_MOD - use mpi -#endif - implicit none -#ifdef MPI_H - include 'mpif.h' -#endif - - integer, intent(out) :: info - type(psb_zspmat_type) :: b - integer, pointer :: nzbr(:), idisp(:) - integer :: ictxt, nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& - & naggr, np, me, nzt,jl,nzl,nlr,& - & icomm,naggrm1, i, j, k, err_act - - name='raw_aggregate' - if(psb_get_errstatus().ne.0) return - info=0 - call psb_erractionsave(err_act) - - call psb_nullify_sp(b) - - ictxt = psb_cd_get_context(desc_a) - icomm = psb_cd_get_mpic(desc_a) - call psb_info(ictxt, me, np) - nglob = psb_cd_get_global_rows(desc_a) - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - - naggr = p%nlaggr(me+1) - ntaggr = sum(p%nlaggr) - allocate(nzbr(np), idisp(np),stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - naggrm1=sum(p%nlaggr(1:me)) - - if (p%iprcparm(coarse_mat_) == repl_mat_) then - do i=1, nrow - p%mlia(i) = p%mlia(i) + naggrm1 - end do - call psb_halo(p%mlia,desc_a,info) - end if - - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_halo') - goto 9999 - end if - - nzt = psb_sp_get_nnzeros(a) - - call psb_sp_all(b,nzt,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spall') - goto 9999 - end if - - call psb_sp_setifld(psb_dupl_ovwrt_,psb_dupl_,b,info) - call psb_sp_setifld(psb_upd_dflt_,psb_upd_,b,info) - b%fida = 'COO' - b%m=a%m - b%k=a%k - call psb_csdp(a,b,info) - if(info /= 0) then - info=4010 - ch_err='psb_csdp' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - nzt = psb_sp_get_nnzeros(b) - - j = 0 - do i=1, nzt - if ((1<=b%ia2(i)).and.(b%ia2(i)<=nrow)) then - j = j + 1 - b%aspk(j) = b%aspk(i) - b%ia1(j) = p%mlia(b%ia1(i)) - b%ia2(j) = p%mlia(b%ia2(i)) - end if - enddo - b%infoa(psb_nnz_)=j - call psb_fixcoo(b,info) - - nzt = psb_sp_get_nnzeros(b) - - call psb_sp_reall(b,nzt,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spreall') - goto 9999 - end if - b%m = naggr - b%k = naggr - - if (p%iprcparm(coarse_mat_) == repl_mat_) then - - call psb_cdrep(ntaggr,ictxt,desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdrep') - goto 9999 - end if - - nzbr(:) = 0 - nzbr(me+1) = nzt - call psb_sum(ictxt,nzbr(1:np)) - nzac = sum(nzbr) - call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spall') - goto 9999 - end if - - do ip=1,np - idisp(ip) = sum(nzbr(1:ip-1)) - enddo - ndx = nzbr(me+1) - - call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,& - & mpi_double_complex,icomm,info) - call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& - & mpi_integer,icomm,info) - call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& - & mpi_integer,icomm,info) - if(info /= 0) then - info=-1 - call psb_errpush(info,name) - goto 9999 - end if - - ac%m = ntaggr - ac%k = ntaggr - ac%infoa(psb_nnz_) = nzac - ac%fida='COO' - ac%descra='G' - call psb_fixcoo(ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_free') - goto 9999 - end if - - else if (p%iprcparm(coarse_mat_) == distr_mat_) then - - call psb_cdall(ictxt,desc_ac,info,nl=naggr) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdall') - goto 9999 - end if - call psb_cdasb(desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdasb') - goto 9999 - end if - - call psb_sp_clone(b,ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spclone') - goto 9999 - end if - call psb_sp_free(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_free') - goto 9999 - end if - - else - - write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(coarse_mat_) - end if - - call psb_ipcoo2csr(ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='ipcoo2csr') - goto 9999 - end if - - deallocate(nzbr,idisp) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine raw_aggregate - - - - subroutine smooth_aggregate(info) - use psb_base_mod - use mld_prec_type -#ifdef MPI_MOD - use mpi -#endif - implicit none -#ifdef MPI_H - include 'mpif.h' -#endif - - integer, intent(out) :: info - - type(psb_zspmat_type) :: b - integer, pointer :: nzbr(:), idisp(:) - integer :: ictxt, nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& - & naggr, np, me, & - & icomm, naggrm1,naggrp1,i,j,err_act,k,nzl - type(psb_zspmat_type), pointer :: am1,am2 - type(psb_zspmat_type) :: am3,am4 - logical :: ml_global_nmb - - logical, parameter :: test_dump=.false., debug=.false. - integer, parameter :: ncmax=16 - real(kind(1.d0)) :: omega, anorm, tmp, dg - character(len=20) :: name - - - name='smooth_aggregate' - if(psb_get_errstatus().ne.0) return - info=0 - call psb_erractionsave(err_act) - - ictxt = psb_cd_get_context(desc_a) - icomm = psb_cd_get_mpic(desc_a) - - call psb_info(ictxt, me, np) - - call psb_nullify_sp(b) - call psb_nullify_sp(am3) - call psb_nullify_sp(am4) - - am2 => p%av(sm_pr_t_) - am1 => p%av(sm_pr_) - - nglob = psb_cd_get_global_rows(desc_a) - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - - naggr = p%nlaggr(me+1) - ntaggr = sum(p%nlaggr) - - allocate(nzbr(np), idisp(np),stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - - naggrm1 = sum(p%nlaggr(1:me)) - naggrp1 = sum(p%nlaggr(1:me+1)) - - ml_global_nmb = ( (p%iprcparm(aggr_kind_) == tent_prol).or.& - & ( (p%iprcparm(aggr_kind_) == biz_prol_).and.& - & (p%iprcparm(coarse_mat_) == repl_mat_)) ) - - - if (ml_global_nmb) then - p%mlia(1:nrow) = p%mlia(1:nrow) + naggrm1 - call psb_halo(p%mlia,desc_a,info) - - if(info /= 0) then - call psb_errpush(4010,name,a_err='f90_pshalo') - goto 9999 - end if - end if - - if (aggr_dump) then - open(30+me) - write(30+me,*) '% Aggregation map' - do i=1,ncol - write(30+me,*) i,p%mlia(i) - end do - close(30+me) - end if - - ! naggr: number of local aggregates - ! nrow: local rows. - ! - allocate(p%dorig(nrow),stat=info) - if (info /= 0) then - info=4025 - call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),& - & a_err='real(kind(1.d0))') - goto 9999 - end if - - ! Get diagonal D - call psb_sp_getdiag(a,p%dorig,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_getdiag') - goto 9999 - end if - - do i=1,size(p%dorig) - if (p%dorig(i) /= zzero) then - p%dorig(i) = zone / p%dorig(i) - else - p%dorig(i) = zone - end if - end do - - ! where (p%dorig /= dzero) - ! p%dorig = done / p%dorig - ! elsewhere - ! p%dorig = done - ! end where - - - ! 1. Allocate Ptilde in sparse matrix form - am4%fida='COO' - am4%m=ncol - if (ml_global_nmb) then - am4%k=ntaggr - call psb_sp_all(ncol,ntaggr,am4,ncol,info) - else - am4%k=naggr - call psb_sp_all(ncol,naggr,am4,ncol,info) - endif - if(info /= 0) then - call psb_errpush(4010,name,a_err='spall') - goto 9999 - end if - - if (ml_global_nmb) then - do i=1,ncol - am4%aspk(i) = zone - am4%ia1(i) = i - am4%ia2(i) = p%mlia(i) - end do - am4%infoa(psb_nnz_) = ncol - else - do i=1,nrow - am4%aspk(i) = zone - am4%ia1(i) = i - am4%ia2(i) = p%mlia(i) - end do - am4%infoa(psb_nnz_) = nrow - endif - - - - - call psb_ipcoo2csr(am4,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='ipcoo2csr') - goto 9999 - end if - - call psb_sp_clone(a,am3,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spclone') - goto 9999 - end if - - ! - ! WARNING: the cycles below assume that AM3 does have - ! its diagonal elements stored explicitly!!! - ! Should we switch to something safer? - ! - call psb_sp_scal(am3,p%dorig,info) - if(info /= 0) goto 9999 - - if (p%iprcparm(aggr_eig_) == max_norm_) then - - if (p%iprcparm(aggr_kind_) == biz_prol_) then - - ! - ! This only works with CSR. - ! - anorm = dzero - dg = done - do i=1,am3%m - tmp = dzero - do j=am3%ia2(i),am3%ia2(i+1)-1 - if (am3%ia1(j) <= am3%m) then - tmp = tmp + abs(am3%aspk(j)) - endif - if (am3%ia1(j) == i ) then - dg = abs(am3%aspk(j)) - end if - end do - anorm = max(anorm,tmp/dg) - enddo - - call psb_amx(ictxt,anorm) - else - anorm = psb_spnrmi(am3,desc_a,info) - endif - omega = 4.d0/(3.d0*anorm) - p%dprcparm(aggr_damp_) = omega - - else if (p%iprcparm(aggr_eig_) == user_choice_) then - - omega = p%dprcparm(aggr_damp_) - - else if (p%iprcparm(aggr_eig_) /= user_choice_) then - write(0,*) me,'Error: invalid choice for OMEGA in blaggrmat?? ',& - & p%iprcparm(aggr_eig_) - end if - - - if (am3%fida=='CSR') then - do i=1,am3%m - do j=am3%ia2(i),am3%ia2(i+1)-1 - if (am3%ia1(j) == i) then - am3%aspk(j) = done - omega*am3%aspk(j) - else - am3%aspk(j) = - omega*am3%aspk(j) - end if - end do - end do - else if (am3%fida=='COO') then - do j=1,am3%infoa(psb_nnz_) - if (am3%ia1(j) /= am3%ia2(j)) then - am3%aspk(j) = - omega*am3%aspk(j) - else - am3%aspk(j) = done - omega*am3%aspk(j) - endif - end do - call psb_ipcoo2csr(am3,info) - else - write(0,*) 'Missing implementation of I sum' - call psb_errpush(4010,name) - goto 9999 - end if - - if (test_dump) then - open(30+me) - write(30+me,*) 'OMEGA: ',omega - do i=1,size(p%dorig) - write(30+me,*) p%dorig(i) - end do - close(30+me) - end if - - if (test_dump) call & - & psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob) - if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,& - & ivc=desc_a%loc_to_glob) - if (debug) write(0,*) me,'Done gather, going for SYMBMM 1' - ! - ! Symbmm90 does the allocation for its result. - ! - ! am1 = (i-wDA)Ptilde - ! Doing it this way means to consider diag(Ai) - ! - ! - call psb_symbmm(am3,am4,am1,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='symbmm 1') - goto 9999 - end if - - call psb_numbmm(am3,am4,am1) - - if (debug) write(0,*) me,'Done NUMBMM 1' - - call psb_sp_free(am4,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_free') - goto 9999 - end if - - if (ml_global_nmb) then - ! - ! Now we have to gather the halo of am1, and add it to itself - ! to multiply it by A, - ! - call psb_sphalo(am1,desc_a,am4,info,clcnv=.false.) - - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sphalo') - goto 9999 - end if - - call psb_rwextd(ncol,am1,info,b=am4) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_rwextd') - goto 9999 - end if - - call psb_sp_free(am4,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - else - - call psb_rwextd(ncol,am1,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='rwextd') - goto 9999 - end if - endif - - if (test_dump) & - & call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob) - - call psb_symbmm(a,am1,am3,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='symbmm 2') - goto 9999 - end if - - call psb_numbmm(a,am1,am3) - if (debug) write(0,*) me,'Done NUMBMM 2' - - if (p%iprcparm(aggr_kind_) == tent_prol) then - call psb_transc(am1,am2,fmt='COO') - nzl = am2%infoa(psb_nnz_) - i=0 - ! - ! Now we have to fix this. The only rows of B that are correct - ! are those corresponding to "local" aggregates, i.e. indices in p%mlia(:) - ! - do k=1, nzl - if ((naggrm1 < am2%ia1(k)) .and.(am2%ia1(k) <= naggrp1)) then - i = i+1 - am2%aspk(i) = am2%aspk(k) - am2%ia1(i) = am2%ia1(k) - am2%ia2(i) = am2%ia2(k) - end if - end do - - am2%infoa(psb_nnz_) = i - call psb_ipcoo2csr(am2,info) - else - call psb_transc(am1,am2) - endif - if (debug) write(0,*) me,'starting sphalo/ rwxtd' - - if (p%iprcparm(aggr_kind_) == tent_prol) then - ! am2 = ((i-wDA)Ptilde)^T - call psb_sphalo(am3,desc_a,am4,info,clcnv=.false.) - - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sphalo') - goto 9999 - end if - call psb_rwextd(ncol,am3,info,b=am4) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_rwextd') - goto 9999 - end if - call psb_sp_free(am4,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - else if (p%iprcparm(aggr_kind_) == biz_prol_) then - - call psb_rwextd(ncol,am3,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_rwextd') - goto 9999 - end if - endif - - if (debug) write(0,*) me,'starting symbmm 3' - call psb_symbmm(am2,am3,b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='symbmm 3') - goto 9999 - end if - - if (debug) write(0,*) me,'starting numbmm 3' - call psb_numbmm(am2,am3,b) - if (debug) write(0,*) me,'Done NUMBMM 3' - -!!$ if (aggr_dump) call csprt(50+me,am1,head='% Operator PTrans.') - call psb_sp_free(am3,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - call psb_ipcsr2coo(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='ipcsr2coo') - goto 9999 - end if - - call psb_fixcoo(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='fixcoo') - goto 9999 - end if - - - if (test_dump) call psb_csprt(80+me,b,head='% Smoothed aggregate AC.') - - select case(p%iprcparm(aggr_kind_)) - - case(tent_prol) - - select case(p%iprcparm(coarse_mat_)) - - case(distr_mat_) - - call psb_sp_clone(b,ac,info) - if(info /= 0) goto 9999 - nzac = ac%infoa(psb_nnz_) - nzl = ac%infoa(psb_nnz_) - - call psb_cdall(ictxt,desc_ac,info,nl=p%nlaggr(me+1)) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdall') - goto 9999 - end if - - - call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdins') - goto 9999 - end if - - if (debug) write(0,*) me,'Created aux descr. distr.' - call psb_cdasb(desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdasb') - goto 9999 - end if - - - if (debug) write(0,*) me,'Asmbld aux descr. distr.' - - call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I') - if(info /= 0) then - call psb_errpush(4010,name,a_err='psglob_to_loc') - goto 9999 - end if - - - call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I') - if(info /= 0) then - call psb_errpush(4010,name,a_err='psglob_to_loc') - goto 9999 - end if - - - ac%m=desc_ac%matrix_data(psb_n_row_) - ac%k=desc_ac%matrix_data(psb_n_col_) - ac%fida='COO' - ac%descra='G' - - call psb_sp_free(b,info) - if (info == 0) deallocate(nzbr,idisp,stat=info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - if (np>1) then - nzl = psb_sp_get_nnzeros(am1) - call psb_glob_to_loc(am1%ia1(1:nzl),desc_ac,info,'I') - - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_glob_to_loc') - goto 9999 - end if - endif - am1%k=desc_ac%matrix_data(psb_n_col_) - - if (np>1) then - call psb_ipcsr2coo(am2,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_ipcsr2coo') - goto 9999 - end if - - nzl = am2%infoa(psb_nnz_) - call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I') - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_glob_to_loc') - goto 9999 - end if - - call psb_ipcoo2csr(am2,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_ipcoo2csr') - goto 9999 - end if - end if - am2%m=desc_ac%matrix_data(psb_n_col_) - - case(repl_mat_) - ! - ! - call psb_cdrep(ntaggr,ictxt,desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdrep') - goto 9999 - end if - - nzbr(:) = 0 - nzbr(me+1) = b%infoa(psb_nnz_) - - call psb_sum(ictxt,nzbr(1:np)) - nzac = sum(nzbr) - call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) - if(info /= 0) goto 9999 - - - do ip=1,np - idisp(ip) = sum(nzbr(1:ip-1)) - enddo - ndx = nzbr(me+1) - - call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,& - & mpi_double_complex,icomm,info) - call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& - & mpi_integer,icomm,info) - call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& - & mpi_integer,icomm,info) - if(info /= 0) goto 9999 - - - ac%m = ntaggr - ac%k = ntaggr - ac%infoa(psb_nnz_) = nzac - ac%fida='COO' - ac%descra='G' - call psb_fixcoo(ac,info) - if(info /= 0) goto 9999 - call psb_sp_free(b,info) - if(info /= 0) goto 9999 - if (me==0) then - if (test_dump) call psb_csprt(80+me,ac,head='% Smoothed aggregate AC.') - endif - - deallocate(nzbr,idisp) - - case default - write(0,*) 'Inconsistent input in smooth_new_aggregate' - end select - - - case(biz_prol_) - - select case(p%iprcparm(coarse_mat_)) - - case(distr_mat_) - - call psb_sp_clone(b,ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='spclone') - goto 9999 - end if - call psb_cdall(ictxt,desc_ac,info,nl=naggr) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdall') - goto 9999 - end if - call psb_cdasb(desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdasb') - goto 9999 - end if - call psb_sp_free(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='sp_free') - goto 9999 - end if - - - case(repl_mat_) - ! - ! - - call psb_cdrep(ntaggr,ictxt,desc_ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_cdrep') - goto 9999 - end if - - nzbr(:) = 0 - nzbr(me+1) = b%infoa(psb_nnz_) - call psb_sum(ictxt,nzbr(1:np)) - nzac = sum(nzbr) - call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_all') - goto 9999 - end if - - call psb_get_mpicomm(ictxt,icomm) - do ip=1,np - idisp(ip) = sum(nzbr(1:ip-1)) - enddo - ndx = nzbr(me+1) - - call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,& - & mpi_double_complex,icomm,info) - call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& - & mpi_integer,icomm,info) - call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& - & mpi_integer,icomm,info) - if(info /= 0) then - info=-1 - call psb_errpush(info,name) - goto 9999 - end if - - - ac%m = ntaggr - ac%k = ntaggr - ac%infoa(psb_nnz_) = nzac - ac%fida='COO' - ac%descra='G' - call psb_fixcoo(ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_fixcoo') - goto 9999 - end if - call psb_sp_free(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_sp_free') - goto 9999 - end if - - end select - deallocate(nzbr,idisp) - - end select - - call psb_ipcoo2csr(ac,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_ipcoo2csr') - goto 9999 - end if - - if (debug) write(0,*) me,'Done smooth_aggregate ' - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name) - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() - return - end if - return - - - end subroutine smooth_aggregate - - - end subroutine mld_zaggrmat_asb diff --git a/mld_zaggrmat_raw_asb.F90 b/mld_zaggrmat_raw_asb.F90 new file mode 100644 index 00000000..1e1c7e8a --- /dev/null +++ b/mld_zaggrmat_raw_asb.F90 @@ -0,0 +1,242 @@ +!!$ +!!$ +!!$ MD2P4 +!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS +!!$ 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 Second 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 MD2P4 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 MD2P4 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 mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info) + use psb_base_mod + use mld_prec_mod, mld_protect_name => mld_zaggrmat_raw_asb + +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + type(psb_zspmat_type), intent(in), target :: a + type(mld_zbaseprc_type), intent(inout),target :: p + type(psb_zspmat_type), intent(inout), target :: ac + type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_ac + integer, intent(out) :: info + + logical, parameter :: aggr_dump=.false. + integer ::ictxt,np,me, err_act,icomm + character(len=20) :: name, ch_err + + type(psb_zspmat_type) :: b + integer, pointer :: nzbr(:), idisp(:) + integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& + & naggr, nzt,jl,nzl,nlr, naggrm1, i, j, k + + name='mld_aggrmat_raw_asb' + if(psb_get_errstatus().ne.0) return + info=0 + call psb_erractionsave(err_act) + + call psb_nullify_sp(b) + + ictxt = psb_cd_get_context(desc_a) + icomm = psb_cd_get_mpic(desc_a) + call psb_info(ictxt, me, np) + nglob = psb_cd_get_global_rows(desc_a) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + naggr = p%nlaggr(me+1) + ntaggr = sum(p%nlaggr) + allocate(nzbr(np), idisp(np),stat=info) + if (info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + naggrm1=sum(p%nlaggr(1:me)) + + if (p%iprcparm(coarse_mat_) == repl_mat_) then + do i=1, nrow + p%mlia(i) = p%mlia(i) + naggrm1 + end do + call psb_halo(p%mlia,desc_a,info) + end if + + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_halo') + goto 9999 + end if + + nzt = psb_sp_get_nnzeros(a) + + call psb_sp_all(b,nzt,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spall') + goto 9999 + end if + + call psb_sp_setifld(psb_dupl_ovwrt_,psb_dupl_,b,info) + call psb_sp_setifld(psb_upd_dflt_,psb_upd_,b,info) + b%fida = 'COO' + b%m=a%m + b%k=a%k + call psb_csdp(a,b,info) + if(info /= 0) then + info=4010 + ch_err='psb_csdp' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + nzt = psb_sp_get_nnzeros(b) + + j = 0 + do i=1, nzt + if ((1<=b%ia2(i)).and.(b%ia2(i)<=nrow)) then + j = j + 1 + b%aspk(j) = b%aspk(i) + b%ia1(j) = p%mlia(b%ia1(i)) + b%ia2(j) = p%mlia(b%ia2(i)) + end if + enddo + b%infoa(psb_nnz_)=j + call psb_fixcoo(b,info) + + nzt = psb_sp_get_nnzeros(b) + + call psb_sp_reall(b,nzt,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spreall') + goto 9999 + end if + b%m = naggr + b%k = naggr + + if (p%iprcparm(coarse_mat_) == repl_mat_) then + + call psb_cdrep(ntaggr,ictxt,desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdrep') + goto 9999 + end if + + nzbr(:) = 0 + nzbr(me+1) = nzt + call psb_sum(ictxt,nzbr(1:np)) + nzac = sum(nzbr) + call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spall') + goto 9999 + end if + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + + call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,& + & mpi_double_complex,icomm,info) + call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& + & mpi_integer,icomm,info) + call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& + & mpi_integer,icomm,info) + if(info /= 0) then + info=-1 + call psb_errpush(info,name) + goto 9999 + end if + + ac%m = ntaggr + ac%k = ntaggr + ac%infoa(psb_nnz_) = nzac + ac%fida='COO' + ac%descra='G' + call psb_fixcoo(ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_free') + goto 9999 + end if + + else if (p%iprcparm(coarse_mat_) == distr_mat_) then + + call psb_cdall(ictxt,desc_ac,info,nl=naggr) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdall') + goto 9999 + end if + call psb_cdasb(desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdasb') + goto 9999 + end if + + call psb_sp_clone(b,ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spclone') + goto 9999 + end if + call psb_sp_free(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_free') + goto 9999 + end if + + else + + write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(coarse_mat_) + end if + + call psb_ipcoo2csr(ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='ipcoo2csr') + goto 9999 + end if + + deallocate(nzbr,idisp) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_zaggrmat_raw_asb diff --git a/mld_zaggrmat_smth_asb.F90 b/mld_zaggrmat_smth_asb.F90 new file mode 100644 index 00000000..69c5bbc2 --- /dev/null +++ b/mld_zaggrmat_smth_asb.F90 @@ -0,0 +1,709 @@ +!!$ +!!$ +!!$ MD2P4 +!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS +!!$ 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 Second 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 MD2P4 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 MD2P4 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 mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info) + use psb_base_mod + use mld_prec_mod, mld_protect_name => mld_zaggrmat_smth_asb + +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + type(psb_zspmat_type), intent(in), target :: a + type(mld_zbaseprc_type), intent(inout),target :: p + type(psb_zspmat_type), intent(inout), target :: ac + type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_ac + integer, intent(out) :: info + + type(psb_zspmat_type) :: b + integer, pointer :: nzbr(:), idisp(:) + integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& + & naggr, nzt,jl,nzl,nlr,naggrm1,naggrp1, i, j, k + logical, parameter :: aggr_dump=.false. + integer ::ictxt,np,me, err_act, icomm + character(len=20) :: name, ch_err + type(psb_zspmat_type), pointer :: am1,am2 + type(psb_zspmat_type) :: am3,am4 + logical :: ml_global_nmb + + logical, parameter :: test_dump=.false., debug=.false. + integer, parameter :: ncmax=16 + real(kind(1.d0)) :: omega, anorm, tmp, dg + + name='mld_aggrmat_smth_asb' + if(psb_get_errstatus().ne.0) return + info=0 + call psb_erractionsave(err_act) + + ictxt = psb_cd_get_context(desc_a) + icomm = psb_cd_get_mpic(desc_a) + ictxt=psb_cd_get_context(desc_a) + + call psb_info(ictxt, me, np) + + + call psb_nullify_sp(b) + call psb_nullify_sp(am3) + call psb_nullify_sp(am4) + + am2 => p%av(sm_pr_t_) + am1 => p%av(sm_pr_) + + nglob = psb_cd_get_global_rows(desc_a) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + naggr = p%nlaggr(me+1) + ntaggr = sum(p%nlaggr) + + allocate(nzbr(np), idisp(np),stat=info) + if (info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + + naggrm1 = sum(p%nlaggr(1:me)) + naggrp1 = sum(p%nlaggr(1:me+1)) + + ml_global_nmb = ( (p%iprcparm(aggr_kind_) == tent_prol).or.& + & ( (p%iprcparm(aggr_kind_) == biz_prol_).and.& + & (p%iprcparm(coarse_mat_) == repl_mat_)) ) + + + if (ml_global_nmb) then + p%mlia(1:nrow) = p%mlia(1:nrow) + naggrm1 + call psb_halo(p%mlia,desc_a,info) + + if(info /= 0) then + call psb_errpush(4010,name,a_err='f90_pshalo') + goto 9999 + end if + end if + + if (aggr_dump) then + open(30+me) + write(30+me,*) '% Aggregation map' + do i=1,ncol + write(30+me,*) i,p%mlia(i) + end do + close(30+me) + end if + + ! naggr: number of local aggregates + ! nrow: local rows. + ! + allocate(p%dorig(nrow),stat=info) + if (info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),& + & a_err='real(kind(1.d0))') + goto 9999 + end if + + ! Get diagonal D + call psb_sp_getdiag(a,p%dorig,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_getdiag') + goto 9999 + end if + + do i=1,size(p%dorig) + if (p%dorig(i) /= zzero) then + p%dorig(i) = zone / p%dorig(i) + else + p%dorig(i) = zone + end if + end do + + ! where (p%dorig /= zzero) + ! p%dorig = zone / p%dorig + ! elsewhere + ! p%dorig = zone + ! end where + + + ! 1. Allocate Ptilde in sparse matrix form + am4%fida='COO' + am4%m=ncol + if (ml_global_nmb) then + am4%k=ntaggr + call psb_sp_all(ncol,ntaggr,am4,ncol,info) + else + am4%k=naggr + call psb_sp_all(ncol,naggr,am4,ncol,info) + endif + if(info /= 0) then + call psb_errpush(4010,name,a_err='spall') + goto 9999 + end if + + if (ml_global_nmb) then + do i=1,ncol + am4%aspk(i) = zone + am4%ia1(i) = i + am4%ia2(i) = p%mlia(i) + end do + am4%infoa(psb_nnz_) = ncol + else + do i=1,nrow + am4%aspk(i) = zone + am4%ia1(i) = i + am4%ia2(i) = p%mlia(i) + end do + am4%infoa(psb_nnz_) = nrow + endif + + + + + call psb_ipcoo2csr(am4,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='ipcoo2csr') + goto 9999 + end if + + call psb_sp_clone(a,am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spclone') + goto 9999 + end if + + ! + ! WARNING: the cycles below assume that AM3 does have + ! its diagonal elements stored explicitly!!! + ! Should we switch to something safer? + ! + call psb_sp_scal(am3,p%dorig,info) + if(info /= 0) goto 9999 + + if (p%iprcparm(aggr_eig_) == max_norm_) then + + if (p%iprcparm(aggr_kind_) == biz_prol_) then + + ! + ! This only works with CSR. + ! + anorm = dzero + dg = done + do i=1,am3%m + tmp = dzero + do j=am3%ia2(i),am3%ia2(i+1)-1 + if (am3%ia1(j) <= am3%m) then + tmp = tmp + abs(am3%aspk(j)) + endif + if (am3%ia1(j) == i ) then + dg = abs(am3%aspk(j)) + end if + end do + anorm = max(anorm,tmp/dg) + enddo + + call psb_amx(ictxt,anorm) + else + anorm = psb_spnrmi(am3,desc_a,info) + endif + omega = 4.d0/(3.d0*anorm) + p%dprcparm(aggr_damp_) = omega + + else if (p%iprcparm(aggr_eig_) == user_choice_) then + + omega = p%dprcparm(aggr_damp_) + + else if (p%iprcparm(aggr_eig_) /= user_choice_) then + write(0,*) me,'Error: invalid choice for OMEGA in blaggrmat?? ',& + & p%iprcparm(aggr_eig_) + end if + + + if (am3%fida=='CSR') then + do i=1,am3%m + do j=am3%ia2(i),am3%ia2(i+1)-1 + if (am3%ia1(j) == i) then + am3%aspk(j) = zone - omega*am3%aspk(j) + else + am3%aspk(j) = - omega*am3%aspk(j) + end if + end do + end do + else if (am3%fida=='COO') then + do j=1,am3%infoa(psb_nnz_) + if (am3%ia1(j) /= am3%ia2(j)) then + am3%aspk(j) = - omega*am3%aspk(j) + else + am3%aspk(j) = zone - omega*am3%aspk(j) + endif + end do + call psb_ipcoo2csr(am3,info) + else + write(0,*) 'Missing implementation of I sum' + call psb_errpush(4010,name) + goto 9999 + end if + + if (test_dump) then + open(30+me) + write(30+me,*) 'OMEGA: ',omega + do i=1,size(p%dorig) + write(30+me,*) p%dorig(i) + end do + close(30+me) + end if + + if (test_dump) call & + & psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob) + if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,& + & ivc=desc_a%loc_to_glob) + if (debug) write(0,*) me,'Done gather, going for SYMBMM 1' + ! + ! Symbmm90 does the allocation for its result. + ! + ! am1 = (i-wDA)Ptilde + ! Doing it this way means to consider diag(Ai) + ! + ! + call psb_symbmm(am3,am4,am1,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm 1') + goto 9999 + end if + + call psb_numbmm(am3,am4,am1) + + if (debug) write(0,*) me,'Done NUMBMM 1' + + call psb_sp_free(am4,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_free') + goto 9999 + end if + + if (ml_global_nmb) then + ! + ! Now we have to gather the halo of am1, and add it to itself + ! to multiply it by A, + ! + call psb_sphalo(am1,desc_a,am4,info,clcnv=.false.) + + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sphalo') + goto 9999 + end if + + call psb_rwextd(ncol,am1,info,b=am4) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_rwextd') + goto 9999 + end if + + call psb_sp_free(am4,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + else + + call psb_rwextd(ncol,am1,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='rwextd') + goto 9999 + end if + endif + + if (test_dump) & + & call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob) + + call psb_symbmm(a,am1,am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm 2') + goto 9999 + end if + + call psb_numbmm(a,am1,am3) + if (debug) write(0,*) me,'Done NUMBMM 2' + + if (p%iprcparm(aggr_kind_) == tent_prol) then + call psb_transp(am1,am2,fmt='COO') + nzl = am2%infoa(psb_nnz_) + i=0 + ! + ! Now we have to fix this. The only rows of B that are correct + ! are those corresponding to "local" aggregates, i.e. indices in p%mlia(:) + ! + do k=1, nzl + if ((naggrm1 < am2%ia1(k)) .and.(am2%ia1(k) <= naggrp1)) then + i = i+1 + am2%aspk(i) = am2%aspk(k) + am2%ia1(i) = am2%ia1(k) + am2%ia2(i) = am2%ia2(k) + end if + end do + + am2%infoa(psb_nnz_) = i + call psb_ipcoo2csr(am2,info) + else + call psb_transp(am1,am2) + endif + if (debug) write(0,*) me,'starting sphalo/ rwxtd' + + if (p%iprcparm(aggr_kind_) == tent_prol) then + ! am2 = ((i-wDA)Ptilde)^T + call psb_sphalo(am3,desc_a,am4,info,clcnv=.false.) + + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sphalo') + goto 9999 + end if + call psb_rwextd(ncol,am3,info,b=am4) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_rwextd') + goto 9999 + end if + call psb_sp_free(am4,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + else if (p%iprcparm(aggr_kind_) == biz_prol_) then + + call psb_rwextd(ncol,am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_rwextd') + goto 9999 + end if + endif + + if (debug) write(0,*) me,'starting symbmm 3' + call psb_symbmm(am2,am3,b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm 3') + goto 9999 + end if + + if (debug) write(0,*) me,'starting numbmm 3' + call psb_numbmm(am2,am3,b) + if (debug) write(0,*) me,'Done NUMBMM 3' + +!!$ if (aggr_dump) call csprt(50+me,am1,head='% Operator PTrans.') + call psb_sp_free(am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + call psb_ipcsr2coo(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='ipcsr2coo') + goto 9999 + end if + + call psb_fixcoo(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='fixcoo') + goto 9999 + end if + + + if (test_dump) call psb_csprt(80+me,b,head='% Smoothed aggregate AC.') + + select case(p%iprcparm(aggr_kind_)) + + case(tent_prol) + + select case(p%iprcparm(coarse_mat_)) + + case(distr_mat_) + + call psb_sp_clone(b,ac,info) + if(info /= 0) goto 9999 + nzac = ac%infoa(psb_nnz_) + nzl = ac%infoa(psb_nnz_) + + call psb_cdall(ictxt,desc_ac,info,nl=p%nlaggr(me+1)) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdall') + goto 9999 + end if + + + call psb_cdins(nzl,ac%ia1,ac%ia2,desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdins') + goto 9999 + end if + + if (debug) write(0,*) me,'Created aux descr. distr.' + call psb_cdasb(desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdasb') + goto 9999 + end if + + + if (debug) write(0,*) me,'Asmbld aux descr. distr.' + + call psb_glob_to_loc(ac%ia1(1:nzl),desc_ac,info,iact='I') + if(info /= 0) then + call psb_errpush(4010,name,a_err='psglob_to_loc') + goto 9999 + end if + + + call psb_glob_to_loc(ac%ia2(1:nzl),desc_ac,info,iact='I') + if(info /= 0) then + call psb_errpush(4010,name,a_err='psglob_to_loc') + goto 9999 + end if + + + ac%m=desc_ac%matrix_data(psb_n_row_) + ac%k=desc_ac%matrix_data(psb_n_col_) + ac%fida='COO' + ac%descra='G' + + call psb_sp_free(b,info) + if (info == 0) deallocate(nzbr,idisp,stat=info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + if (np>1) then + nzl = psb_sp_get_nnzeros(am1) + call psb_glob_to_loc(am1%ia1(1:nzl),desc_ac,info,'I') + + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_glob_to_loc') + goto 9999 + end if + endif + am1%k=desc_ac%matrix_data(psb_n_col_) + + if (np>1) then + call psb_ipcsr2coo(am2,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcsr2coo') + goto 9999 + end if + + nzl = am2%infoa(psb_nnz_) + call psb_glob_to_loc(am2%ia1(1:nzl),desc_ac,info,'I') + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_glob_to_loc') + goto 9999 + end if + + call psb_ipcoo2csr(am2,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcoo2csr') + goto 9999 + end if + end if + am2%m=desc_ac%matrix_data(psb_n_col_) + + if (debug) write(0,*) me,'Done ac ' + case(repl_mat_) + ! + ! + call psb_cdrep(ntaggr,ictxt,desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdrep') + goto 9999 + end if + + nzbr(:) = 0 + nzbr(me+1) = b%infoa(psb_nnz_) + + call psb_sum(ictxt,nzbr(1:np)) + nzac = sum(nzbr) + call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) + if(info /= 0) goto 9999 + + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + + call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,& + & mpi_double_complex,icomm,info) + call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& + & mpi_integer,icomm,info) + call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& + & mpi_integer,icomm,info) + if(info /= 0) goto 9999 + + + ac%m = ntaggr + ac%k = ntaggr + ac%infoa(psb_nnz_) = nzac + ac%fida='COO' + ac%descra='G' + call psb_fixcoo(ac,info) + if(info /= 0) goto 9999 + call psb_sp_free(b,info) + if(info /= 0) goto 9999 + if (me==0) then + if (test_dump) call psb_csprt(80+me,ac,head='% Smoothed aggregate AC.') + endif + + deallocate(nzbr,idisp) + + case default + write(0,*) 'Inconsistent input in smooth_new_aggregate' + end select + + + case(biz_prol_) + + select case(p%iprcparm(coarse_mat_)) + + case(distr_mat_) + + call psb_sp_clone(b,ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='spclone') + goto 9999 + end if + call psb_cdall(ictxt,desc_ac,info,nl=naggr) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdall') + goto 9999 + end if + + call psb_cdasb(desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdasb') + goto 9999 + end if + + call psb_sp_free(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='sp_free') + goto 9999 + end if + + + case(repl_mat_) + ! + ! + + call psb_cdrep(ntaggr,ictxt,desc_ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdrep') + goto 9999 + end if + + nzbr(:) = 0 + nzbr(me+1) = b%infoa(psb_nnz_) + call psb_sum(ictxt,nzbr(1:np)) + nzac = sum(nzbr) + call psb_sp_all(ntaggr,ntaggr,ac,nzac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_all') + goto 9999 + end if + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + + call mpi_allgatherv(b%aspk,ndx,mpi_double_complex,ac%aspk,nzbr,idisp,& + & mpi_double_complex,icomm,info) + call mpi_allgatherv(b%ia1,ndx,mpi_integer,ac%ia1,nzbr,idisp,& + & mpi_integer,icomm,info) + call mpi_allgatherv(b%ia2,ndx,mpi_integer,ac%ia2,nzbr,idisp,& + & mpi_integer,icomm,info) + if(info /= 0) then + info=-1 + call psb_errpush(info,name) + goto 9999 + end if + + + ac%m = ntaggr + ac%k = ntaggr + ac%infoa(psb_nnz_) = nzac + ac%fida='COO' + ac%descra='G' + call psb_fixcoo(ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_fixcoo') + goto 9999 + end if + call psb_sp_free(b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + end select + deallocate(nzbr,idisp) + + end select + + call psb_ipcoo2csr(ac,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcoo2csr') + goto 9999 + end if + + if (debug) write(0,*) me,'Done smooth_aggregate ' + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + + + +end subroutine mld_zaggrmat_smth_asb