diff --git a/src/modules/psb_methd_mod.f90 b/src/modules/psb_methd_mod.f90 index 5817c2f4..e51f95a5 100644 --- a/src/modules/psb_methd_mod.f90 +++ b/src/modules/psb_methd_mod.f90 @@ -63,7 +63,6 @@ Module psb_methd_mod use psb_serial_mod use psb_descriptor_type Use psb_prec_type -!!$ parameters Type(psb_dspmat_type), Intent(in) :: a Type(psb_desc_type), Intent(in) :: desc_a type(psb_dprec_type), intent(in) :: prec diff --git a/src/modules/psb_psblas_mod.f90 b/src/modules/psb_psblas_mod.f90 index 86ca12ca..bdd89b5f 100644 --- a/src/modules/psb_psblas_mod.f90 +++ b/src/modules/psb_psblas_mod.f90 @@ -229,5 +229,22 @@ module psb_psblas_mod end subroutine psb_dspsv end interface - + + interface psb_gelp + subroutine psb_dgelp(trans,iperm,x,desc_a,info) + use psb_descriptor_type + type(psb_desc_type), intent(in) :: desc_a + real(kind(1.d0)), intent(inout) :: x(:,:) + integer, intent(inout) :: iperm(:),info + character, intent(in) :: trans + end subroutine psb_dgelp + subroutine psb_dgelpv(trans,iperm,x,desc_a,info) + use psb_descriptor_type + type(psb_desc_type), intent(in) :: desc_a + real(kind(1.d0)), intent(inout) :: x(:) + integer, intent(inout) :: iperm(:),info + character, intent(in) :: trans + end subroutine psb_dgelpv + end interface + end module psb_psblas_mod diff --git a/src/prec/psb_dbldaggrmat.f90 b/src/prec/psb_dbldaggrmat.f90 index b5d75292..55eb2342 100644 --- a/src/prec/psb_dbldaggrmat.f90 +++ b/src/prec/psb_dbldaggrmat.f90 @@ -185,7 +185,6 @@ contains enddo end if - call psb_fixcoo(b,info) if(info /= 0) then call psb_errpush(4010,name,a_err='fixcoo') diff --git a/src/prec/psb_dprec.f90 b/src/prec/psb_dprec.f90 index 26da3a4c..2bf23007 100644 --- a/src/prec/psb_dprec.f90 +++ b/src/prec/psb_dprec.f90 @@ -1,3 +1,31 @@ +<<<<<<< variant A + /home/buttari/NUMERICAL/psblas-2.0/src/prec: + total used in directory 1000 available 4183416 + drwxr-xr-x 2 buttari users 4096 Oct 14 11:26 CVS + -rw-r--r-- 1 buttari users 9723 Sep 7 11:16 fort_slu_impl.c + -rw-r--r-- 1 buttari users 3616 Oct 13 17:19 fort_slu_impl.o + -rw-r--r-- 1 buttari users 20777 Sep 7 11:16 gps.f + -rw-r--r-- 1 buttari users 47396 Oct 13 17:19 gps.o + -rw-r--r-- 1 buttari users 623 Sep 12 17:38 Makefile + -rw-r--r-- 1 buttari users 23481 Oct 13 20:40 psb_dbldaggrmat.f90 + -rw-r--r-- 1 buttari users 158084 Oct 13 20:40 psb_dbldaggrmat.o + -rw-r--r-- 1 buttari users 19821 Sep 27 11:35 psb_dcslu.f90 + -rw-r--r-- 1 buttari users 109004 Oct 13 17:19 psb_dcslu.o + -rw-r--r-- 1 buttari users 6503 Sep 12 14:28 psb_dcsrsetup.f90 + -rw-r--r-- 1 buttari users 15096 Oct 13 17:19 psb_dcsrsetup.o + -rw-r--r-- 1 buttari users 6868 Oct 14 10:38 psb_dgenaggrmap.f90 + -rw-r--r-- 1 buttari users 37432 Oct 14 10:39 psb_dgenaggrmap.o + -rw-r--r-- 1 buttari users 17221 Oct 14 10:38 psb_dprecbld.f90 + -rw-r--r-- 1 buttari users 111376 Oct 14 10:39 psb_dprecbld.o + -rw-r--r-- 1 buttari users 27689 Oct 13 12:40 psb_dprec.f90 + -rw-r--r-- 1 buttari users 1669 Oct 7 16:25 psb_dprecfree.f90 + -rw-r--r-- 1 buttari users 5632 Oct 13 17:19 psb_dprecfree.o + -rw-r--r-- 1 buttari users 200168 Oct 13 17:19 psb_dprec.o + -rw-r--r-- 1 buttari users 5353 Sep 30 14:29 psb_dprecset.f90 + -rw-r--r-- 1 buttari users 72964 Oct 13 17:19 psb_dprecset.o + -rw-r--r-- 1 buttari users 11328 Oct 5 17:28 psb_dsplu.f90 + -rw-r--r-- 1 buttari users 44108 Oct 13 17:19 psb_dsplu.o +>>>>>>> variant B subroutine psb_dprecaply(prec,x,y,desc_data,info,trans, work) use psb_serial_mod @@ -970,3 +998,977 @@ subroutine psb_dprecaply1(prec,x,desc_data,info,trans) return end subroutine psb_dprecaply1 +####### Ancestor +subroutine psb_dprecaply(prec,x,y,desc_data,info,trans, work) + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type),intent(in) :: desc_data + type(psb_dprec_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(kind(0.d0)), optional, target :: work(:) + + ! Local variables + character ::trans_ + real(kind(1.d0)), pointer :: work_(:) + integer :: icontxt,nprow,npcol,me,mycol,err_act, int_err(5) + logical,parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + external mpi_wtime + character(len=20) :: name, ch_err + + interface + subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbase_prec), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dbaseprcaply + end interface + + interface + subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbase_prec), intent(in) :: baseprecv(:) + real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(inout) :: x(:), y(:) + character :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dmlprcaply + end interface + + name='psb_dprecaply' + info = 0 + call psb_erractionsave(err_act) + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + + if (present(trans)) then + trans_=trans + else + trans_='N' + end if + + if (present(work)) then + work_ => work + else + allocate(work_(4*desc_data%matrix_data(psb_n_col_))) + end if + + if (.not.(associated(prec%baseprecv))) then + write(0,*) 'Inconsistent preconditioner: neither SMTH nor BASE?' + end if + if (size(prec%baseprecv) >1) then + call psb_dmlprcaply(prec%baseprecv,x,zero,y,desc_data,trans_,work_,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_dmlprcaply') + goto 9999 + end if + + else if (size(prec%baseprecv) == 1) then + call psb_dbaseprcaply(prec%baseprecv(1),x,zero,y,desc_data,trans_, work_,info) + else + write(0,*) 'Inconsistent preconditioner: size of baseprecv???' + endif + + if (present(work)) then + else + deallocate(work_) + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + +end subroutine psb_dprecaply + + + +subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) + ! + ! Compute Y <- beta*Y + K^-1 X + ! where K is a a basic preconditioner stored in prec + ! + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbase_prec), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + + ! Local variables + integer :: n_row,n_col, int_err(5) + real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:) + character ::diagl, diagu + integer :: icontxt,nprow,npcol,me,mycol,i, isz, nrg, err_act + real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime + logical,parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + external mpi_wtime + character(len=20) :: name, ch_err + + interface + subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type), intent(in) :: desc_data + type(psb_dbase_prec), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dbjacaply + end interface + + name='psb_dbaseprcaply' + info = 0 + call psb_erractionsave(err_act) + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + + diagl='U' + diagu='U' + + select case(trans) + case('N','n') + case('T','t','C','c') + case default + info=40 + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + select case(prec%iprcparm(p_type_)) + + case(noprec_) + + n_row=desc_data%matrix_data(psb_n_row_) + if (beta==zero) then + y(1:n_row) = x(1:n_row) + else if (beta==one) then + y(1:n_row) = x(1:n_row) + y(1:n_row) + else if (beta==-one) then + y(1:n_row) = x(1:n_row) - y(1:n_row) + else + y(1:n_row) = x(1:n_row) + beta*y(1:n_row) + end if + + case(diagsc_) + + n_row=desc_data%matrix_data(psb_n_row_) + if (beta==zero) then + y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + else if (beta==one) then + y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + y(1:n_row) + else if (beta==-one) then + y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - y(1:n_row) + else + y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + beta*y(1:n_row) + end if + + case(bja_) + + call psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) + if(info.ne.0) then + info=4010 + ch_err='psb_bjacaply' + goto 9999 + end if + + case(asm_,ras_,ash_,rash_) + + ! Note: currently trans is unused. + n_row=prec%desc_data%matrix_data(psb_n_row_) + n_col=prec%desc_data%matrix_data(psb_n_col_) + + isz=max(n_row,N_COL) + if ((6*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + aux => work(3*isz+1:) + else if ((4*isz) <= size(work)) then + aux => work(1:) + allocate(ww(isz),tx(isz),ty(isz)) + else if ((3*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + allocate(aux(4*isz)) + else + allocate(ww(isz),tx(isz),ty(isz),& + &aux(4*isz)) + endif + + if (debug) write(0,*)' vdiag: ',prec%d(:) + if (debug) write(0,*) 'Bi-CGSTAB with Additive Schwarz prec' + + tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_)) + tx(desc_data%matrix_data(psb_n_row_)+1:isz) = zero + + if (prec%iprcparm(restr_)==psb_halo_) then + call psb_halo(tx,prec%desc_data,info,work=aux) + if(info /=0) then + info=4010 + ch_err='psb_halo' + goto 9999 + end if + else if (prec%iprcparm(restr_) /= psb_none_) then + write(0,*) 'Problem in PRCAPLY: Unknown value for restriction ',& + &prec%iprcparm(restr_) + end if + + if (prec%iprcparm(iren_)>0) then + call psb_dgelp('N',n_row,1,prec%perm,tx,isz,ww,isz,info) + if(info /=0) then + info=4010 + ch_err='psb_dgelp' + goto 9999 + end if + endif + + call psb_dbjacaply(prec,tx,zero,ty,prec%desc_data,trans,aux,info) + + if (prec%iprcparm(iren_)>0) then + call psb_dgelp('N',n_row,1,prec%invperm,ty,isz,ww,isz,info) + if(info /=0) then + info=4010 + ch_err='psb_dgelp' + goto 9999 + end if + endif + + select case (prec%iprcparm(prol_)) + + case(psb_none_) + ! Would work anyway, but since it's supposed to do nothing... + ! call f90_psovrl(ty,prec%desc_data,update_type=prec%a_restrict) + + case(psb_sum_,psb_avg_) + call psb_ovrl(ty,prec%desc_data,info,& + & update_type=prec%iprcparm(prol_),work=aux) + if(info /=0) then + info=4010 + ch_err='psb_ovrl' + goto 9999 + end if + + case default + write(0,*) 'Problem in PRCAPLY: Unknown value for prolongation ',& + & prec%iprcparm(prol_) + end select + + if (beta == zero) then + y(1:desc_data%matrix_data(psb_n_row_)) = ty(1:desc_data%matrix_data(psb_n_row_)) + else if (beta == one) then + y(1:desc_data%matrix_data(psb_n_row_)) = y(1:desc_data%matrix_data(psb_n_row_)) + & + & ty(1:desc_data%matrix_data(psb_n_row_)) + else if (beta == -one) then + y(1:desc_data%matrix_data(psb_n_row_)) = -y(1:desc_data%matrix_data(psb_n_row_)) + & + & ty(1:desc_data%matrix_data(psb_n_row_)) + else + y(1:desc_data%matrix_data(psb_n_row_)) = beta*y(1:desc_data%matrix_data(psb_n_row_)) + & + & ty(1:desc_data%matrix_data(psb_n_row_)) + end if + + + if ((6*isz) <= size(work)) then + else if ((4*isz) <= size(work)) then + deallocate(ww,tx,ty) + else if ((3*isz) <= size(work)) then + deallocate(aux) + else + deallocate(ww,aux,tx,ty) + endif + + case default + write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(p_type_),':',& + & min_prec_,noprec_,diagsc_,bja_,& + & asm_,ras_,ash_,rash_ + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name,i_err=int_err,a_err=ch_err) + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + +end subroutine psb_dbaseprcaply + + + +subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) + ! + ! Compute Y <- beta*Y + K^-1 X + ! where K is a a Block Jacobi preconditioner stored in prec + ! Note that desc_data may or may not be the same as prec%desc_data, + ! but since both are INTENT(IN) this should be legal. + ! + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type), intent(in) :: desc_data + type(psb_dbase_prec), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + + ! Local variables + integer :: n_row,n_col + real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:),tb(:) + character ::diagl, diagu + integer :: icontxt,nprow,npcol,me,mycol,i, isz, nrg, err_act, int_err(5) + real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime + logical,parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + external mpi_wtime + character(len=20) :: name, ch_err + + name='psb_dbjacaply' + info = 0 + call psb_erractionsave(err_act) + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + + diagl='U' + diagu='U' + + select case(trans) + case('N','n') + case('T','t','C','c') + case default + call psb_errpush(40,name) + goto 9999 + end select + + + n_row=desc_data%matrix_data(psb_n_row_) + n_col=desc_data%matrix_data(psb_n_col_) + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col)) + endif + else + allocate(ww(n_col),aux(4*n_col)) + endif + + + if (prec%iprcparm(jac_sweeps_) == 1) then + + + select case(prec%iprcparm(f_type_)) + case(f_ilu_n_,f_ilu_e_) + + select case(trans) + case('N','n') + + call psb_spsm(one,prec%av(l_pr_),x,zero,ww,desc_data,info,& + & trans='N',unit=diagl,choice=psb_none_,work=aux) + if(info /=0) goto 9999 + ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) + call psb_spsm(one,prec%av(u_pr_),ww,beta,y,desc_data,info,& + & trans='N',unit=diagu,choice=psb_none_, work=aux) + if(info /=0) goto 9999 + + case('T','t','C','c') + call psb_spsm(one,prec%av(u_pr_),x,zero,ww,desc_data,info,& + & trans=trans,unit=diagu,choice=psb_none_, work=aux) + if(info /=0) goto 9999 + ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) + call psb_spsm(one,prec%av(l_pr_),ww,beta,y,desc_data,info,& + & trans=trans,unit=diagl,choice=psb_none_,work=aux) + if(info /=0) goto 9999 + + end select + + case(f_slu_) + + ww(1:n_row) = x(1:n_row) + + select case(trans) + case('N','n') + call fort_slu_solve(0,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info) + case('T','t','C','c') + call fort_slu_solve(1,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info) + end select + + if(info /=0) goto 9999 + + if (beta == 0.d0) then + y(1:n_row) = ww(1:n_row) + else if (beta==1.d0) then + y(1:n_row) = ww(1:n_row) + y(1:n_row) + else if (beta==-1.d0) then + y(1:n_row) = ww(1:n_row) - y(1:n_row) + else + y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) + endif + + case default + write(0,*) 'Unknown factorization type in bjac_prcaply',prec%iprcparm(f_type_) + end select + if (debug) write(0,*)' Y: ',y(:) + + else if (prec%iprcparm(jac_sweeps_) > 1) then + + ! Note: we have to add TRANS to this one !!!!!!!!! + + if (size(prec%av) < ap_nd_) then + info = 4011 + goto 9999 + endif + + allocate(tx(n_col),ty(n_col)) + tx = zero + ty = zero + select case(prec%iprcparm(f_type_)) + case(f_ilu_n_,f_ilu_e_) + do i=1, prec%iprcparm(jac_sweeps_) + ! X(k+1) = M^-1*(b-N*X(k)) + ty(1:n_row) = x(1:n_row) + call psb_spmm(-one,prec%av(ap_nd_),tx,one,ty,& + & prec%desc_data,info,work=aux) + if(info /=0) goto 9999 + call psb_spsm(one,prec%av(l_pr_),ty,zero,ww,& + & prec%desc_data,info,& + & trans='N',unit='U',choice=psb_none_,work=aux) + if(info /=0) goto 9999 + ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) + call psb_spsm(one,prec%av(u_pr_),ww,zero,tx,& + & prec%desc_data,info,& + & trans='N',unit='U',choice=psb_none_,work=aux) + if(info /=0) goto 9999 + end do + + case(f_slu_) + do i=1, prec%iprcparm(jac_sweeps_) + ! X(k+1) = M^-1*(b-N*X(k)) + ty(1:n_row) = x(1:n_row) + call psb_spmm(-one,prec%av(ap_nd_),tx,one,ty,& + & prec%desc_data,info,work=aux) + if(info /=0) goto 9999 + + call fort_slu_solve(0,n_row,1,ty,n_row,prec%iprcparm(slu_ptr_),info) + if(info /=0) goto 9999 + tx(1:n_row) = ty(1:n_row) + end do + end select + + if (beta == 0.d0) then + y(1:n_row) = tx(1:n_row) + else if (beta==1.d0) then + y(1:n_row) = tx(1:n_row) + y(1:n_row) + else if (beta==-1.d0) then + y(1:n_row) = tx(1:n_row) - y(1:n_row) + else + y(1:n_row) = tx(1:n_row) + beta*y(1:n_row) + endif + + deallocate(tx,ty) + + + else + + goto 9999 + + endif + + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name,i_err=int_err,a_err=ch_err) + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + +end subroutine psb_dbjacaply + + +subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) + ! + ! Compute Y <- beta*Y + K^-1 X + ! where K is a multilevel (actually 2-level) preconditioner stored in prec + ! + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_blacs_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbase_prec), intent(in) :: baseprecv(:) + real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(inout) :: x(:), y(:) + character :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + + + ! Local variables + integer :: n_row,n_col + real(kind(1.d0)), allocatable :: tx(:),ty(:),t2l(:),w2l(:),& + & x2l(:),b2l(:),tz(:),tty(:) + character ::diagl, diagu + integer :: icontxt,nprow,npcol,me,mycol,i, isz, nrg,nr2l,err_act, iptype, int_err(5) + real(kind(1.d0)) :: omega + real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime + logical, parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + integer :: ismth + external mpi_wtime + character(len=20) :: name, ch_err + + interface + subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbase_prec), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dbaseprcaply + end interface + + name='psb_dmlprcaply' + info = 0 + call psb_erractionsave(err_act) + + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + + omega=baseprecv(2)%dprcparm(smooth_omega_) + ismth=baseprecv(2)%iprcparm(smth_kind_) + + select case(baseprecv(2)%iprcparm(ml_type_)) + case(no_ml_) + ! Should not really get here. + write(0,*) 'Smooth preconditioner with no multilevel in MLPRCAPLY????' + + case(add_ml_prec_) + + ! + ! Additive multilevel + ! + t1 = mpi_wtime() + n_row = desc_data%matrix_data(psb_n_row_) + n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + call psb_dbaseprcaply(baseprecv(1),x,beta,y,desc_data,trans,work,info) + if(info /=0) goto 9999 + + nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) + allocate(t2l(nr2l),w2l(nr2l)) + t2l(:) = zero + w2l(:) = zero + + if (ismth /= no_smth_) then + ! + ! Smoothed aggregation + ! + allocate(tx(max(n_row,n_col)),ty(max(n_row,n_col)),& + & tz(max(n_row,n_col))) + tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_)) + tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + ty(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + tz(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + + + if (baseprecv(2)%iprcparm(glb_smth_) >0) then + call psb_halo(tx,desc_data,info,work=work) + if(info /=0) goto 9999 + else + tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + end if + + call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) + if(info /=0) goto 9999 + + else + ! + ! Raw aggregation, may take shortcut + ! + do i=1,desc_data%matrix_data(psb_n_row_) + t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + x(i) + end do + + end if + + if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then + call gsum2d(icontxt,'All',t2l(1:nrg)) + else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(2)%iprcparm(coarse_mat_) + endif + + w2l=t2l + call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + + + if (ismth /= no_smth_) then + + call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) + if(info /=0) goto 9999 + ! + ! Finally add back into Y. + ! + call psb_axpby(one,ty,one,y,desc_data,info) + if(info /=0) goto 9999 + deallocate(tx,ty,tz) + + else + + do i=1, desc_data%matrix_data(psb_n_row_) + y(i) = y(i) + t2l(baseprecv(2)%mlia(i)) + enddo + + end if + + if (debug) write(0,*)' Y2: ',Y(:) + + deallocate(t2l,w2l) + + case(mult_ml_prec_) + + ! + ! Multiplicative multilevel + ! Pre/post smoothing versions. + + select case(baseprecv(2)%iprcparm(smth_pos_)) + + case(post_smooth_) + + + t1 = mpi_wtime() + n_row = desc_data%matrix_data(psb_n_row_) + n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) + allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col)) + t2l(:) = zero + w2l(:) = zero + + ! + ! Need temp copies to handle Y<- betaY + K^-1 X + ! One of the temp copies is not strictly needed when beta==zero + ! + + if (debug) write(0,*)' mult_ml_apply omega ',omega + if (debug) write(0,*)' mult_ml_apply X: ',X(:) + call psb_axpby(one,x,zero,tx,desc_data,info) + if(info /=0) goto 9999 + + if (ismth /= no_smth_) then + ! + ! Smoothed aggregation + ! + allocate(tz(max(n_row,n_col))) + + if (baseprecv(2)%iprcparm(glb_smth_) >0) then + call psb_halo(tx,desc_data,info,work=work) + if(info /=0) goto 9999 + else + tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + end if + + call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) + if(info /=0) goto 9999 + + else + ! + ! Raw aggregation, may take shortcut + ! + do i=1,desc_data%matrix_data(psb_n_row_) + t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) + end do + end if + + if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then + call gsum2d(icontxt,'All',t2l(1:nrg)) + else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(2)%iprcparm(coarse_mat_) + endif + + t6 = mpi_wtime() + w2l=t2l + call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + if(info /=0) goto 9999 + + if (ismth /= no_smth_) then + if (ismth == smth_omg_) & + & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) + call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) + if(info /=0) goto 9999 + ! + ! Finally add back into Y. + ! + deallocate(tz) + else + ty(:) = zero + do i=1, desc_data%matrix_data(psb_n_row_) + ty(i) = ty(i) + t2l(baseprecv(2)%mlia(i)) + enddo + + end if + deallocate(t2l,w2l) + + call psb_spmm(-one,baseprecv(2)%aorig,ty,one,tx,desc_data,info,work=work) + if(info /=0) goto 9999 + + call psb_dbaseprcaply(baseprecv(1),tx,one,ty,desc_data,trans,work,info) + if(info /=0) goto 9999 + + call psb_axpby(one,ty,beta,y,desc_data,info) + if(info /=0) goto 9999 + + deallocate(tx,ty) + + + + case(pre_smooth_) + + t1 = mpi_wtime() + n_row = desc_data%matrix_data(psb_n_row_) + n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) + allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col)) + t2l(:) = zero + w2l(:) = zero + + ! + ! Need temp copies to handle Y<- betaY + K^-1 X + ! One of the temp copies is not strictly needed when beta==zero + ! + call psb_axpby(one,x,zero,tx,desc_data,info) + call psb_axpby(one,y,zero,ty,desc_data,info) + if(info /=0) goto 9999 + + call psb_dbaseprcaply(baseprecv(1),x,zero,tty,desc_data,trans,work,info) + if(info /=0) goto 9999 + + call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) + if(info /=0) goto 9999 + + if (ismth /= no_smth_) then + allocate(tz(max(n_row,n_col))) + + if (baseprecv(2)%iprcparm(glb_smth_) >0) then + call psb_halo(tx,desc_data,info,work=work) + if(info /=0) goto 9999 + else + tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + end if + + call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) + if(info /=0) goto 9999 + + else + ! + ! Raw aggregation, may take shortcuts + ! + do i=1,desc_data%matrix_data(psb_n_row_) + t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) + end do + end if + + if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then + call gsum2d(icontxt,'All',t2l(1:nrg)) + else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(2)%iprcparm(coarse_mat_) + endif + + t6 = mpi_wtime() + w2l=t2l + call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + if(info /=0) goto 9999 + + if (ismth /= no_smth_) then + + if (ismth == smth_omg_) & + & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) + call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) + if(info /=0) goto 9999 + + call psb_axpby(one,ty,one,tty,desc_data,info) + if(info /=0) goto 9999 + + deallocate(tz) + else + + do i=1, desc_data%matrix_data(psb_n_row_) + tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i)) + enddo + + end if + + call psb_axpby(one,tty,beta,y,desc_data,info) + if(info /=0) goto 9999 + + deallocate(t2l,w2l,tx,ty,tty) + + + + + case default + + write(0,*) 'Unknown value for ml_smooth_pos',baseprecv(2)%iprcparm(smth_pos_) + + end select + + case default + write(0,*) me, 'Wrong mltype into PRCAPLY ',& + & baseprecv(2)%iprcparm(ml_type_) + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + +end subroutine psb_dmlprcaply + + +subroutine psb_dprecaply1(prec,x,desc_data,info,trans) + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type),intent(in) :: desc_data + type(psb_dprec_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + logical,parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + + + ! Local variables + character :: trans_ + integer :: icontxt,nprow,npcol,me,mycol,i, isz, err_act, int_err(5) + real(kind(1.d0)), pointer :: WW(:), w1(:) + character(len=20) :: name, ch_err + name='psb_dprec1' + info = 0 + call psb_erractionsave(err_act) + + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + if (present(trans)) then + trans_=trans + else + trans_='N' + end if + + allocate(ww(size(x)),w1(size(x))) + call psb_dprecaply(prec,x,ww,desc_data,info,trans_,w1) + if(info /=0) goto 9999 + x(:) = ww(:) + deallocate(ww,W1) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return +end subroutine psb_dprecaply1 + +======= end diff --git a/src/prec/psb_dprecbld.f90 b/src/prec/psb_dprecbld.f90 index 2ac3c4ef..7da64648 100644 --- a/src/prec/psb_dprecbld.f90 +++ b/src/prec/psb_dprecbld.f90 @@ -11,7 +11,7 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd) Implicit None integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target :: a + type(psb_dspmat_type), target :: a type(psb_dprec_type),intent(inout) :: p type(psb_desc_type), intent(in) :: desc_a character, intent(in), optional :: upd @@ -125,7 +125,7 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd) call psb_errpush(info,name) goto 9999 end if - call psb_dgelp('n',a%Pl,p%baseprecv(1)%d,desc_a,info) + call psb_gelp('n',a%pl,p%baseprecv(1)%d,desc_a,info) if(info /= 0) then info=4010 ch_err='psb_dgelp' @@ -239,13 +239,16 @@ 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) - if(info /= 0) then - info=4010 - ch_err='psb_mlprec_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if + call blacs_barrier(icontxt,'All') ! to be removed + write(0,'(i2," Calling mlprecbld")')me + call blacs_barrier(icontxt,'All') ! to be removed + call psb_mlprec_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 + end if endif @@ -496,6 +499,8 @@ subroutine psb_mlprec_bld(a,desc_a,p,info) integer, intent(out) :: info end subroutine psb_dbldaggrmat end interface + + integer :: icontxt, nprow, npcol, me, mycol name='psb_mlprec_bld' info=0 diff --git a/src/serial/psb_dneigh.f90 b/src/serial/psb_dneigh.f90 index 4101e999..f55dbe46 100644 --- a/src/serial/psb_dneigh.f90 +++ b/src/serial/psb_dneigh.f90 @@ -16,23 +16,9 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev) integer, pointer :: neigh(:) ! the neighbours integer, optional :: lev ! level of neighbours to find - interface psb_spgtrow - subroutine psb_dspgtrow(irw,a,b,info,append,iren,lrw) - use psb_spmat_type - type(psb_dspmat_type), intent(in) :: a - integer, intent(in) :: irw - type(psb_dspmat_type), intent(inout) :: b - logical, intent(in), optional :: append - integer, intent(in), target, optional :: iren(:) - integer, intent(in), optional :: lrw - integer, intent(out) :: info - end subroutine psb_dspgtrow - end interface - - type(psb_dspmat_type) :: atmp integer, pointer :: tmpn(:) integer :: level, dim, i, j, k, r, c, brow,& - & elem_pt, ii, n1, col_idx, ne, err_act + & elem_pt, ii, n1, col_idx, ne, err_act, nn, nidx character(len=20) :: name, ch_err name='psb_dneigh' @@ -52,105 +38,245 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev) level=1 end if - if ((a%fida /= 'CSR')) then - - call psb_spall(atmp,1,info) - if(info.ne.0) then - call psb_errpush(4010,name) - goto 9999 - end if - call psb_spgtrow(idx,a,atmp,info) - if(info.ne.0) then - call psb_errpush(4010,name) - goto 9999 - end if - dim=atmp%infoa(psb_nnz_) - allocate(tmpn(dim)) - tmpn(1:dim)=atmp%ia2(1:dim) - - if(level.eq.2) then - do i=1,dim - if ((1<=tmpn(i)).and.(tmpn(i)<=a%m).and.(tmpn(i).ne.idx)) then - call psb_spgtrow(tmpn(i),a,atmp,info,append=.true.) - if(info.ne.0) then - call psb_errpush(4010,name) - goto 9999 - end if - end if - end do - end if + call psb_dneigh1l(a,idx,neigh,n) + + if(level.eq.2) then + allocate(tmpn(1)) + n1=n + do i=1,n1 + nidx=neigh(i) + if((nidx.ne.idx).and.(nidx.gt.0).and.(nidx.le.a%m)) then + call psb_dneigh1l(a,nidx,tmpn,nn) + if((n+nn).gt.size(neigh)) call psb_realloc(n+nn,neigh,info) + neigh(n+1:n+nn)=tmpn(1:nn) + n=n+nn + end if + end do + 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 - dim=atmp%infoa(psb_nnz_) - if(dim.gt.size(neigh)) & - & call psb_realloc(dim,neigh,info) - if(info.ne.0) then - call psb_errpush(4010,name) - goto 9999 - end if - call psb_spfree(atmp,info) - if(info.ne.0) then - call psb_errpush(4010,name) - goto 9999 - end if - call psb_nullify_sp(atmp) - deallocate(tmpn) - - else if(a%fida.eq.'CSR') then +contains - dim=0 - if(level.eq.1) dim=(a%ia2(idx+1)-a%ia2(idx)) - if(dim >size(neigh)) call psb_realloc(dim,neigh,info) - if(info.ne.izero) then - info=4010 - ch_err='psrealloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif + subroutine psb_dneigh1l(a,idx,neigh,n) + + use psb_realloc_mod + use psb_const_mod + use psb_spmat_type + implicit none + + type(psb_dspmat_type), intent(in) :: a ! the sparse matrix + integer, intent(in) :: idx ! the index whose neighbours we want to find + integer, intent(out) :: n ! the number of neighbours and the info + integer, pointer :: neigh(:) ! the neighbours + + + select case(a%fida(1:3)) + case('CSR') + call csr_dneigh1l(a,idx,neigh,n) + case('COO') + call coo_dneigh1l(a,idx,neigh,n) + case('JAD') + call jad_dneigh1l(a,idx,neigh,n) + end select + + end subroutine psb_dneigh1l + + subroutine csr_dneigh1l(a,idx,neigh,n) + + use psb_realloc_mod + use psb_const_mod + use psb_spmat_type + implicit none + + + type(psb_dspmat_type), intent(in) :: a ! the sparse matrix + integer, intent(in) :: idx ! the index whose neighbours we want to find + integer, intent(out) :: n ! the number of neighbours and the info + integer, pointer :: neigh(:) ! the neighbours + + integer :: dim, i, iidx + + if(a%pl(1).ne.0) then + iidx=a%pl(idx) + else + iidx=idx + end if + + dim=a%ia2(iidx+1)-a%ia2(iidx) + if(dim.gt.size(neigh)) call psb_realloc(dim,neigh,info) + n=0 - if(level.eq.1) then - do i=a%ia2(idx), a%ia2(idx+1)-1 - n=n+1 - neigh(n)=a%ia1(i) - end do + do i=a%ia2(iidx), a%ia2(iidx+1)-1 + n=n+1 + neigh(n)=a%ia1(i) + end do + + end subroutine csr_dneigh1l + + + subroutine coo_dneigh1l(a,idx,neigh,n) + + use psb_realloc_mod + use psb_const_mod + use psb_spmat_type + implicit none + + + type(psb_dspmat_type), intent(in) :: a ! the sparse matrix + integer, intent(in) :: idx ! the index whose neighbours we want to find + integer, intent(out) :: n ! the number of neighbours and the info + integer, pointer :: neigh(:) ! the neighbours + + integer :: dim, i, iidx, ip, nza + if(a%pl(1).ne.0) then + iidx=a%pl(idx) else + iidx=idx + end if + + nza=a%infoa(psb_nnz_) - do i=a%ia2(idx), a%ia2(idx+1)-1 - - j=a%ia1(i) - if ((1<=j).and.(j<=a%m).and.(j.ne.idx)) then - - dim=dim+ a%ia2(j+1)-a%ia2(j) - if(dim.gt.size(neigh)) call psb_realloc(dim,neigh,info) - if(info.ne.izero) then - info=4010 - ch_err='psrealloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - - do k=a%ia2(j), a%ia2(j+1)-1 - n=n+1 - neigh(n)=a%ia1(k) + if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then + call ibsrch(ip,iidx,nza,a%ia1) + if (ip /= -1) then + ! bring ip backward to the beginning of the row + do + if (ip < 2) exit + if (a%ia1(ip-1) == iidx) then + ip = ip -1 + else + exit + end if end do - end if - end do + end if + + dim=0 + do + if(a%ia1(ip).eq.iidx) then + dim=dim+1 + if(dim.gt.size(neigh)) call psb_realloc(dim*3/2,neigh,info) + neigh(dim)=a%ia2(ip) + ip=ip+1 + else + exit + end if + end do + + else + dim=0 + do i=1,nza + if(a%ia1(i).eq.iidx) then + dim=dim+1 + if(dim.gt.size(neigh)) call psb_realloc(dim*3/2,neigh,info) + neigh(dim)=a%ia2(ip) + end if + end do end if - end if + n=dim + + end subroutine coo_dneigh1l + - 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 + + subroutine jad_dneigh1l(a,idx,neigh,n) + + use psb_realloc_mod + use psb_const_mod + use psb_spmat_type + implicit none + + + type(psb_dspmat_type), intent(in) :: a ! the sparse matrix + integer, intent(in) :: idx ! the index whose neighbours we want to find + integer, intent(out) :: n ! the number of neighbours and the info + integer, pointer :: neigh(:) ! the neighbours + + integer :: dim, i, iidx, ip, nza + integer, pointer :: ia1(:), ia2(:), ia3(:),& + & ja(:), ka(:) + integer :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ng + + if(a%pl(1).ne.0) then + iidx=a%pl(idx) + else + iidx=idx + end if + + nza=a%infoa(psb_nnz_) + + png = a%ia2(1) ! points to the number of blocks + pia = a%ia2(2) ! points to the beginning of ia(3,png) + pja = a%ia2(3) ! points to the beginning of ja(:) + + ng = a%ia2(png) ! the number of blocks + ja => a%ia2(pja:) ! the array containing the pointers to ka and aspk + ka => a%ia1(:) ! the array containing the column indices + ia1 => a%ia2(pia:pja-1:3) ! the array containing the first row index of each block + ia2 => a%ia2(pia+1:pja-1:3) ! the array containing a pointer to the pos. in ja to the first jad column + ia3 => a%ia2(pia+2:pja-1:3) ! the array containing a pointer to the pos. in ja to the first csr column + + + i=0 + dim=0 + blkfnd: do + i=i+1 + if(ia1(i).eq.iidx) then + blk=i + dim=dim+ia3(i)-ia2(i) + ipx = ia1(i) ! the first row index of the block + rb = iidx-ipx ! the row offset within the block + row = ia3(i)+rb + dim = dim+ja(row+1)-ja(row) + exit blkfnd + else if(ia1(i).gt.iidx) then + blk=i-1 + dim=dim+ia3(i-1)-ia2(i-1) + ipx = ia1(i-1) ! the first row index of the block + rb = iidx-ipx ! the row offset within the block + row = ia3(i-1)+rb + dim = dim+ja(row+1)-ja(row) + exit blkfnd + end if + end do blkfnd + + if(dim.gt.size(neigh)) call psb_realloc(dim,neigh,info) + + ipx = ia1(blk) ! the first row index of the block + k_pt= ia2(blk) ! the pointer to the beginning of a column in ja + rb = iidx-ipx ! the row offset within the block + npg = ja(k_pt+1)-ja(k_pt) ! the number of rows in the block + + k=0 + do col = ia2(blk), ia3(blk)-1 + k=k+1 + neigh(k) = ka(ja(col)+rb) + end do + + ! extract second part of the row from the csr tail + row=ia3(blk)+rb + do j=ja(row), ja(row+1)-1 + k=k+1 + neigh(k) = ka(j) + end do + + n=k + + end subroutine jad_dneigh1l end subroutine psb_dneigh + diff --git a/src/serial/psb_dspgtdiag.f90 b/src/serial/psb_dspgtdiag.f90 index 1515e293..bc2b1cbd 100644 --- a/src/serial/psb_dspgtdiag.f90 +++ b/src/serial/psb_dspgtdiag.f90 @@ -68,6 +68,8 @@ subroutine psb_dspgtdiag(a,d,info) else if (a%fida == 'JAD') then rng=min(a%m,a%k) + nrb=16 + write(0,*)'in spgtdiag' do i=1, rng, nrb irb=min(i+nrb-1,rng) call psb_spgtrow(i,a,tmpa,info,lrw=irb) @@ -86,6 +88,7 @@ subroutine psb_dspgtdiag(a,d,info) enddo end do + write(0,*)'leaving spgtdiag' end if diff --git a/src/serial/psb_dspgtrow.f90 b/src/serial/psb_dspgtrow.f90 index 3d039847..450e9eec 100644 --- a/src/serial/psb_dspgtrow.f90 +++ b/src/serial/psb_dspgtrow.f90 @@ -429,21 +429,38 @@ contains rb = indices(i)-ipx ! the row offset within the block npg = ja(k_pt+1)-ja(k_pt) ! the number of rows in the block - do col = ia2(blk), ia3(blk)-1 - k=k+1 - b%aspk(nzb+k) = a%aspk(ja(col)+rb) - b%ia1(nzb+k) = irw+i-1 - b%ia2(nzb+k) = ka(ja(col)+rb) - end do - + if(associated(iren))then + do col = ia2(blk), ia3(blk)-1 + k=k+1 + b%aspk(nzb+k) = a%aspk(ja(col)+rb) + b%ia1(nzb+k) = iren(irw+i-1) + b%ia2(nzb+k) = iren(ka(ja(col)+rb)) + end do + else + do col = ia2(blk), ia3(blk)-1 + k=k+1 + b%aspk(nzb+k) = a%aspk(ja(col)+rb) + b%ia1(nzb+k) = irw+i-1 + b%ia2(nzb+k) = ka(ja(col)+rb) + end do + end if ! extract second part of the row from the csr tail row=ia3(blk)+rb - do j=ja(row), ja(row+1)-1 - k=k+1 - b%aspk(nzb+k) = a%aspk(j) - b%ia1(nzb+k) = irw+i-1 - b%ia2(nzb+k) = ka(j) - end do + if(associated(iren))then + do j=ja(row), ja(row+1)-1 + k=k+1 + b%aspk(nzb+k) = a%aspk(j) + b%ia1(nzb+k) = iren(irw+i-1) + b%ia2(nzb+k) = iren(ka(j)) + end do + else + do j=ja(row), ja(row+1)-1 + k=k+1 + b%aspk(nzb+k) = a%aspk(j) + b%ia1(nzb+k) = irw+i-1 + b%ia2(nzb+k) = ka(j) + end do + end if end do b%infoa(psb_nnz_) = nzb+k diff --git a/test/Fileread/RUNS/rtst.inp b/test/Fileread/RUNS/rtst.inp index 41615778..023616a4 100644 --- a/test/Fileread/RUNS/rtst.inp +++ b/test/Fileread/RUNS/rtst.inp @@ -1,13 +1,13 @@ 11 Number of inputs -kivap001.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ +bcsstk35.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ NONE -BICGSTAB +CGS ILU !!!! Actually, it's IPREC below. Should take this line out. CSR 2 IPART: Partition method 1 ISTOPC 00800 ITMAX --1 ITRACE +6 ITRACE 7 IPREC 0:NONE 1:DIAGSC 2:ILU 3: AS 4: RAS 5,6: variants 1 ML 1.d-6 EPS diff --git a/test/Fileread/df_sample.f90 b/test/Fileread/df_sample.f90 index 6becb2fe..4be5269b 100644 --- a/test/Fileread/df_sample.f90 +++ b/test/Fileread/df_sample.f90 @@ -64,8 +64,8 @@ program df_sample integer :: internal, m,ii,nnzero real(kind(1.d0)) :: mpi_wtime, t1, t2, tprec, r_amax, b_amax,& &scale,resmx,resmxp - integer :: nrhs, nrow, n_row, dim, nv - integer, pointer :: ivg(:), ipv(:) + integer :: nrhs, nrow, n_row, dim, nv, ne + integer, pointer :: ivg(:), ipv(:), neigh(:) external mpi_wtime @@ -260,6 +260,9 @@ program df_sample else if (cmethd.eq.'CGS') then call psb_cgs(a,pre,b_col,x_col,eps,desc_a,info,& & itmax,iter,err,itrace) + else if (cmethd.eq.'CG') then + call psb_cg(a,pre,b_col,x_col,eps,desc_a,info,& + & itmax,iter,err,itrace) else if (cmethd.eq.'BICGSTABL') then call psb_bicgstabl(a,pre,b_col,x_col,eps,desc_a,info,& & itmax,iter,err,ierr,itrace,ml) @@ -297,12 +300,12 @@ program df_sample if (amroot) then write(0,'(" ")') write(0,'("Saving x on file")') - write(20,*) 'matrix: ',mtrx_file - write(20,*) 'computed solution on ',nprow,' processors.' - write(20,*) 'iterations to convergence: ',iter - write(20,*) 'error indicator (infinity norm) on exit:', & - & ' ||r||/(||a||||x||+||b||) = ',err - write(20,*) 'max residual = ',resmx, resmxp +!!$ write(20,*) 'matrix: ',mtrx_file +!!$ write(20,*) 'computed solution on ',nprow,' processors.' +!!$ write(20,*) 'iterations to convergence: ',iter +!!$ write(20,*) 'error indicator (infinity norm) on exit:', & +!!$ & ' ||r||/(||a||||x||+||b||) = ',err +!!$ write(20,*) 'max residual = ',resmx, resmxp do i=1,m_problem write(20,998) i,x_col_glob(i),r_col_glob(i),b_col_glob(i) enddo diff --git a/test/Fileread/getp.f90 b/test/Fileread/getp.f90 index 2e992021..1f7451a1 100644 --- a/test/Fileread/getp.f90 +++ b/test/Fileread/getp.f90 @@ -104,6 +104,7 @@ CONTAINS write(*,'("Preconditioner : ",i)')iprec if(iprec.gt.2) write(*,'("Overlapping levels : ",i)')novr write(*,'("Iterative method : ",a)')cmethd + write(*,'("Storage format : ",a3)')afmt(1:3) write(*,'(" ")') else CALL PR_USAGE(0) diff --git a/test/Fileread/mmio.f90 b/test/Fileread/mmio.f90 index 06241570..b93266ff 100644 --- a/test/Fileread/mmio.f90 +++ b/test/Fileread/mmio.f90 @@ -149,7 +149,6 @@ contains a%ia1(1:nzr) = ia2_loc(1:nzr) tmp(1:nzr) = ia1_loc(1:nzr) else - write(0,*) 'After DESYM: ',nzr,ia2_loc(1:10) do i=1,nzr a%aspk(i) = as_loc(i) a%ia1(i) = ia2_loc(i)