diff --git a/base/modules/psb_serial_mod.f90 b/base/modules/psb_serial_mod.f90 index 88acca98..e653786a 100644 --- a/base/modules/psb_serial_mod.f90 +++ b/base/modules/psb_serial_mod.f90 @@ -43,6 +43,20 @@ module psb_serial_mod use psb_mat_mod interface psb_symbmm + subroutine psb_ssymbmm(a,b,c,info) + use psb_mat_mod + implicit none + type(psb_s_sparse_mat), intent(in) :: a,b + type(psb_s_sparse_mat), intent(out) :: c + integer, intent(out) :: info + end subroutine psb_ssymbmm + subroutine psb_sbase_symbmm(a,b,c,info) + use psb_mat_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a,b + type(psb_s_csr_sparse_mat), intent(out) :: c + integer, intent(out) :: info + end subroutine psb_sbase_symbmm subroutine psb_dsymbmm(a,b,c,info) use psb_mat_mod implicit none @@ -57,8 +71,49 @@ module psb_serial_mod type(psb_d_csr_sparse_mat), intent(out) :: c integer, intent(out) :: info end subroutine psb_dbase_symbmm + subroutine psb_csymbmm(a,b,c,info) + use psb_mat_mod + implicit none + type(psb_c_sparse_mat), intent(in) :: a,b + type(psb_c_sparse_mat), intent(out) :: c + integer, intent(out) :: info + end subroutine psb_csymbmm + subroutine psb_cbase_symbmm(a,b,c,info) + use psb_mat_mod + implicit none + class(psb_c_base_sparse_mat), intent(in) :: a,b + type(psb_c_csr_sparse_mat), intent(out) :: c + integer, intent(out) :: info + end subroutine psb_cbase_symbmm + subroutine psb_zsymbmm(a,b,c,info) + use psb_mat_mod + implicit none + type(psb_z_sparse_mat), intent(in) :: a,b + type(psb_z_sparse_mat), intent(out) :: c + integer, intent(out) :: info + end subroutine psb_zsymbmm + subroutine psb_zbase_symbmm(a,b,c,info) + use psb_mat_mod + implicit none + class(psb_z_base_sparse_mat), intent(in) :: a,b + type(psb_z_csr_sparse_mat), intent(out) :: c + integer, intent(out) :: info + end subroutine psb_zbase_symbmm end interface + interface psb_numbmm + subroutine psb_snumbmm(a,b,c) + use psb_mat_mod + implicit none + type(psb_s_sparse_mat), intent(in) :: a,b + type(psb_s_sparse_mat), intent(inout) :: c + end subroutine psb_snumbmm + subroutine psb_sbase_numbmm(a,b,c) + use psb_mat_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a,b + type(psb_s_csr_sparse_mat), intent(inout) :: c + end subroutine psb_sbase_numbmm subroutine psb_dnumbmm(a,b,c) use psb_mat_mod implicit none @@ -71,9 +126,51 @@ module psb_serial_mod class(psb_d_base_sparse_mat), intent(in) :: a,b type(psb_d_csr_sparse_mat), intent(inout) :: c end subroutine psb_dbase_numbmm + subroutine psb_cnumbmm(a,b,c) + use psb_mat_mod + implicit none + type(psb_c_sparse_mat), intent(in) :: a,b + type(psb_c_sparse_mat), intent(inout) :: c + end subroutine psb_cnumbmm + subroutine psb_cbase_numbmm(a,b,c) + use psb_mat_mod + implicit none + class(psb_c_base_sparse_mat), intent(in) :: a,b + type(psb_c_csr_sparse_mat), intent(inout) :: c + end subroutine psb_cbase_numbmm + subroutine psb_znumbmm(a,b,c) + use psb_mat_mod + implicit none + type(psb_z_sparse_mat), intent(in) :: a,b + type(psb_z_sparse_mat), intent(inout) :: c + end subroutine psb_znumbmm + subroutine psb_zbase_numbmm(a,b,c) + use psb_mat_mod + implicit none + class(psb_z_base_sparse_mat), intent(in) :: a,b + type(psb_z_csr_sparse_mat), intent(inout) :: c + end subroutine psb_zbase_numbmm end interface interface psb_rwextd + subroutine psb_srwextd(nr,a,info,b,rowscale) + use psb_mat_mod + implicit none + integer, intent(in) :: nr + type(psb_s_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + type(psb_s_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale + end subroutine psb_srwextd + subroutine psb_sbase_rwextd(nr,a,info,b,rowscale) + use psb_mat_mod + implicit none + integer, intent(in) :: nr + class(psb_s_base_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + class(psb_s_base_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale + end subroutine psb_sbase_rwextd subroutine psb_drwextd(nr,a,info,b,rowscale) use psb_mat_mod implicit none @@ -92,6 +189,42 @@ module psb_serial_mod class(psb_d_base_sparse_mat), intent(in), optional :: b logical,intent(in), optional :: rowscale end subroutine psb_dbase_rwextd + subroutine psb_crwextd(nr,a,info,b,rowscale) + use psb_mat_mod + implicit none + integer, intent(in) :: nr + type(psb_c_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + type(psb_c_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale + end subroutine psb_crwextd + subroutine psb_cbase_rwextd(nr,a,info,b,rowscale) + use psb_mat_mod + implicit none + integer, intent(in) :: nr + class(psb_c_base_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + class(psb_c_base_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale + end subroutine psb_cbase_rwextd + subroutine psb_zrwextd(nr,a,info,b,rowscale) + use psb_mat_mod + implicit none + integer, intent(in) :: nr + type(psb_z_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + type(psb_z_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale + end subroutine psb_zrwextd + subroutine psb_zbase_rwextd(nr,a,info,b,rowscale) + use psb_mat_mod + implicit none + integer, intent(in) :: nr + class(psb_z_base_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + class(psb_z_base_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale + end subroutine psb_zbase_rwextd end interface diff --git a/base/serial/Makefile b/base/serial/Makefile index 872287c9..ded3bb17 100644 --- a/base/serial/Makefile +++ b/base/serial/Makefile @@ -1,8 +1,10 @@ include ../../Make.inc -FOBJS = psb_lsame.o psb_dsymbmm.o psb_dnumbmm.o \ - psb_drwextd.o +FOBJS = psb_lsame.o \ + psb_ssymbmm.o psb_dsymbmm.o psb_csymbmm.o psb_zsymbmm.o \ + psb_snumbmm.o psb_dnumbmm.o psb_cnumbmm.o psb_znumbmm.o \ + psb_srwextd.o psb_drwextd.o psb_crwextd.o psb_zrwextd.o # psb_regen_mod.o psb_lsame.o psb_zspgetrow.o\ # psb_zcsmm.o psb_zcsmv.o psb_zspgtdiag.o psb_zspgtblk.o\ # psb_zcsnmi.o psb_zcsrws.o psb_zcssm.o psb_zcssv.o psb_zspcnv.o\ diff --git a/base/serial/psb_cnumbmm.f90 b/base/serial/psb_cnumbmm.f90 index 89105984..7029d00d 100644 --- a/base/serial/psb_cnumbmm.f90 +++ b/base/serial/psb_cnumbmm.f90 @@ -29,7 +29,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: psb_snumbmm.f90 +! File: psb_cnumbmm.f90 ! Subroutine: ! Arguments: ! @@ -41,41 +41,146 @@ ! subroutine psb_cnumbmm(a,b,c) - use psb_spmat_type + use psb_mat_mod + use psb_string_mod use psb_serial_mod, psb_protect_name => psb_cnumbmm - implicit none + implicit none - type(psb_cspmat_type) :: a,b,c + type(psb_c_sparse_mat), intent(in) :: a,b + type(psb_c_sparse_mat), intent(inout) :: c + integer :: info + integer :: err_act + character(len=*), parameter :: name='psb_numbmm' + + call psb_erractionsave(err_act) + info = 0 + + if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + select type(aa=>c%a) + type is (psb_c_csr_sparse_mat) + call psb_numbmm(a%a,b%a,aa) + class default + info = 1121 + call psb_errpush(info,name) + goto 9999 + end select + + call c%set_asb() + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_cnumbmm + +subroutine psb_cbase_numbmm(a,b,c) + use psb_mat_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_cbase_numbmm + implicit none + + class(psb_c_base_sparse_mat), intent(in) :: a,b + type(psb_c_csr_sparse_mat), intent(inout) :: c + integer, allocatable :: itemp(:) + integer :: nze, ma,na,mb,nb + character(len=20) :: name complex(psb_spk_), allocatable :: temp(:) - integer :: info - logical :: csra, csrb - - allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info) - if (info /= 0) then - return + integer :: info + integer :: err_act + name='psb_numbmm' + call psb_erractionsave(err_act) + info = 0 + + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb endif - call psb_realloc(size(c%ia1),c%aspk,info) + allocate(temp(max(ma,na,mb,nb)),stat=info) + if (info /= 0) then + info = 4000 + call psb_Errpush(info,name) + goto 9999 + endif + ! ! Note: we still have to test about possible performance hits. ! ! - csra = (psb_toupper(a%fida(1:3))=='CSR') - csrb = (psb_toupper(b%fida(1:3))=='CSR') - - if (csra.and.csrb) then - call cnumbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,& - & b%ia2,b%ia1,0,b%aspk,& - & c%ia2,c%ia1,0,c%aspk,temp) - else - call inner_numbmm(a,b,c,temp,info) + call psb_ensure_size(size(c%ja),c%val,info) + select type(a) + type is (psb_c_csr_sparse_mat) + select type(b) + type is (psb_c_csr_sparse_mat) + call csr_numbmm(a,b,c,temp,info) + class default + call gen_numbmm(a,b,c,temp,info) + end select + class default + call gen_numbmm(a,b,c,temp,info) + end select + + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 end if + + call c%set_asb() deallocate(temp) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if return contains - subroutine inner_numbmm(a,b,c,temp,info) - type(psb_cspmat_type) :: a,b,c + subroutine csr_numbmm(a,b,c,temp,info) + type(psb_c_csr_sparse_mat), intent(in) :: a,b + type(psb_c_csr_sparse_mat), intent(inout) :: c + complex(psb_spk_) :: temp(:) + integer, intent(out) :: info + integer :: nze, ma,na,mb,nb + + info = 0 + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + call cnumbmm(ma,na,nb,a%irp,a%ja,0,a%val,& + & b%irp,b%ja,0,b%val,& + & c%irp,c%ja,0,c%val,temp) + + + end subroutine csr_numbmm + + subroutine gen_numbmm(a,b,c,temp,info) + class(psb_c_base_sparse_mat), intent(in) :: a,b + type(psb_c_csr_sparse_mat), intent(inout) :: c integer :: info complex(psb_spk_) :: temp(:) integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) @@ -83,56 +188,60 @@ contains integer :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln complex(psb_spk_) :: ajj - - n = a%m - m = a%k - l = b%k + n = a%get_nrows() + m = a%get_ncols() + l = b%get_ncols() maxlmn = max(l,m,n) allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& & aval(maxlmn),bval(maxlmn), stat=info) - if (info /= 0) then + if (info /= 0) then + info = 4000 return endif do i = 1,maxlmn - temp(i) = czero + temp(i) = dzero end do minlm = min(l,m) minln = min(l,n) minmn = min(m,n) do i = 1,n - call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info) - + call a%csget(i,i,nazr,iarw,iacl,aval,info) do jj=1, nazr j=iacl(jj) ajj = aval(jj) if ((j<1).or.(j>m)) then write(0,*) ' NUMBMM: Problem with A ',i,jj,j,m + info = 1 + return + endif - call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info) + call b%csget(j,j,nbzr,ibrw,ibcl,bval,info) do k=1,nbzr if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn + info = 2 + return else temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k) endif enddo end do - do j = c%ia2(i),c%ia2(i+1)-1 - if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn)) then - write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn + do j = c%irp(i),c%irp(i+1)-1 + if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then + write(0,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn + info = 3 + return else - c%aspk(j) = temp(c%ia1(j)) - temp(c%ia1(j)) = czero + c%val(j) = temp(c%ja(j)) + temp(c%ja(j)) = dzero endif end do end do + + end subroutine gen_numbmm +end subroutine psb_cbase_numbmm - - end subroutine inner_numbmm - - -end subroutine psb_cnumbmm diff --git a/base/serial/psb_crwextd.f90 b/base/serial/psb_crwextd.f90 index 2ee918be..ee3ccbdd 100644 --- a/base/serial/psb_crwextd.f90 +++ b/base/serial/psb_crwextd.f90 @@ -39,126 +39,209 @@ ! ! subroutine psb_crwextd(nr,a,info,b,rowscale) - use psb_spmat_type use psb_error_mod use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_crwextd implicit none ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes) - integer, intent(in) :: nr - type(psb_cspmat_type), intent(inout) :: a - integer,intent(out) :: info - type(psb_cspmat_type), intent(in), optional :: b - logical,intent(in), optional :: rowscale + integer, intent(in) :: nr + type(psb_c_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + type(psb_c_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale integer :: i,j,ja,jb,err_act,nza,nzb character(len=20) :: name, ch_err + type(psb_c_coo_sparse_mat) :: actmp logical rowscale_ name='psb_crwextd' info = 0 call psb_erractionsave(err_act) + if (nr > a%get_nrows()) then + select type(aa=> a%a) + type is (psb_c_csr_sparse_mat) + if (present(b)) then + call psb_rwextd(nr,aa,info,b%a,rowscale) + else + call psb_rwextd(nr,aa,info,rowscale=rowscale) + end if + type is (psb_c_coo_sparse_mat) + if (present(b)) then + call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale) + else + call psb_rwextd(nr,aa,info,rowscale=rowscale) + end if + class default + call aa%mv_to_coo(actmp,info) + if (info == 0) then + if (present(b)) then + call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale) + else + call psb_rwextd(nr,actmp,info,rowscale=rowscale) + end if + end if + if (info == 0) call aa%mv_from_coo(actmp,info) + end select + end if + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_crwextd +subroutine psb_cbase_rwextd(nr,a,info,b,rowscale) + use psb_error_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_cbase_rwextd + implicit none + + ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes) + integer, intent(in) :: nr + class(psb_c_base_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + class(psb_c_base_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale + + integer :: i,j,ja,jb,err_act,nza,nzb, ma, mb, na, nb + character(len=20) :: name, ch_err + logical rowscale_ + + name='psb_cbase_rwextd' + info = 0 + call psb_erractionsave(err_act) + if (present(rowscale)) then rowscale_ = rowscale else rowscale_ = .true. end if - if (nr > a%m) then - if (psb_toupper(a%fida) == 'CSR') then - call psb_realloc(nr+1,a%ia2,info) - if (present(b)) then - nzb = psb_sp_get_nnzeros(b) - call psb_realloc(size(a%ia1)+nzb,a%ia1,info) - call psb_realloc(size(a%aspk)+nzb,a%aspk,info) - if (psb_toupper(b%fida)=='CSR') then - - do i=1, min(nr-a%m,b%m) - a%ia2(a%m+i+1) = a%ia2(a%m+i) + b%ia2(i+1) - b%ia2(i) - ja = a%ia2(a%m+i) - jb = b%ia2(i) - do - if (jb >= b%ia2(i+1)) exit - a%aspk(ja) = b%aspk(jb) - a%ia1(ja) = b%ia1(jb) - ja = ja + 1 - jb = jb + 1 - end do - end do - do j=i,nr-a%m - a%ia2(a%m+i+1) = a%ia2(a%m+i) + ma = a%get_nrows() + na = a%get_ncols() + + + select type(a) + type is (psb_c_csr_sparse_mat) + + call psb_ensure_size(nr+1,a%irp,info) + + if (present(b)) then + mb = b%get_nrows() + nb = b%get_ncols() + nzb = b%get_nzeros() + + select type (b) + type is (psb_c_csr_sparse_mat) + call psb_ensure_size(size(a%ja)+nzb,a%ja,info) + call psb_ensure_size(size(a%val)+nzb,a%val,info) + do i=1, min(nr-ma,mb) + a%irp(ma+i+1) = a%irp(ma+i) + b%irp(i+1) - b%irp(i) + ja = a%irp(ma+i) + jb = b%irp(i) + do + if (jb >= b%irp(i+1)) exit + a%val(ja) = b%val(jb) + a%ja(ja) = b%ja(jb) + ja = ja + 1 + jb = jb + 1 end do - else - write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' - endif - else - do i=a%m+2,nr+1 - a%ia2(i) = a%ia2(i-1) end do - end if - a%m = nr - a%k = max(a%k,b%k) + do j=i,nr-ma + a%irp(ma+i+1) = a%irp(ma+i) + end do + class default - else if (psb_toupper(a%fida) == 'COO') then + write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' + end select + call a%set_ncols(max(na,nb)) - if (present(b)) then - nza = psb_sp_get_nnzeros(a) - nzb = psb_sp_get_nnzeros(b) - call psb_sp_reall(a,nza+nzb,info) - if (psb_toupper(b%fida)=='COO') then - if (rowscale_) then - do j=1,nzb - if ((a%m + b%ia1(j)) <= nr) then - nza = nza + 1 - a%ia1(nza) = a%m + b%ia1(j) - a%ia2(nza) = b%ia2(j) - a%aspk(nza) = b%aspk(j) - end if - enddo - else - do j=1,nzb - if ((b%ia1(j)) <= nr) then - nza = nza + 1 - a%ia1(nza) = b%ia1(j) - a%ia2(nza) = b%ia2(j) - a%aspk(nza) = b%aspk(j) - endif - enddo - endif - a%infoa(psb_nnz_) = nza - else if(psb_toupper(b%fida)=='CSR') then - do i=1, min(nr-a%m,b%m) - do - jb = b%ia2(i) - if (jb >= b%ia2(i+1)) exit - nza = nza + 1 - a%aspk(nza) = b%aspk(jb) - a%ia1(nza) = a%m + i - a%ia2(nza) = b%ia1(jb) - jb = jb + 1 - end do - end do - a%infoa(psb_nnz_) = nza - else - write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' - endif - endif - a%m = nr - a%k = max(a%k,b%k) - else if (a%fida == 'JAD') then - info=135 - ch_err=a%fida(1:3) - call psb_errpush(info,name,a_err=ch_err) - goto 9999 else - info=136 - ch_err=a%fida(1:3) - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + + do i=ma+2,nr+1 + a%irp(i) = a%irp(i-1) + end do + end if - end if + call a%set_nrows(nr) + + + type is (psb_c_coo_sparse_mat) + nza = a%get_nzeros() + + if (present(b)) then + mb = b%get_nrows() + nb = b%get_ncols() + nzb = b%get_nzeros() + call a%reallocate(nza+nzb) + + select type(b) + type is (psb_c_coo_sparse_mat) + + if (rowscale_) then + do j=1,nzb + if ((ma + b%ia(j)) <= nr) then + nza = nza + 1 + a%ia(nza) = ma + b%ia(j) + a%ja(nza) = b%ja(j) + a%val(nza) = b%val(j) + end if + enddo + else + do j=1,nzb + if ((ma + b%ia(j)) <= nr) then + nza = nza + 1 + a%ia(nza) = b%ia(j) + a%ja(nza) = b%ja(j) + a%val(nza) = b%val(j) + end if + enddo + endif + call a%set_nzeros(nza) + + type is (psb_c_csr_sparse_mat) + + do i=1, min(nr-ma,mb) + do + jb = b%irp(i) + if (jb >= b%irp(i+1)) exit + nza = nza + 1 + a%val(nza) = b%val(jb) + a%ia(nza) = ma + i + a%ja(nza) = b%ja(jb) + jb = jb + 1 + end do + end do + call a%set_nzeros(nza) + + class default + write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' + + end select + + call a%set_ncols(max(na,nb)) + endif + + call a%set_nrows(nr) + + class default + info = 135 + ch_err=a%get_fmt() + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + + end select call psb_erractionrestore(err_act) return @@ -166,9 +249,9 @@ subroutine psb_crwextd(nr,a,info,b,rowscale) 9999 continue call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error() - return + call psb_error() + return end if return -end subroutine psb_crwextd +end subroutine psb_cbase_rwextd diff --git a/base/serial/psb_csymbmm.f90 b/base/serial/psb_csymbmm.f90 index 8716103c..e7b32411 100644 --- a/base/serial/psb_csymbmm.f90 +++ b/base/serial/psb_csymbmm.f90 @@ -40,62 +40,107 @@ ! subroutine psb_csymbmm(a,b,c,info) - use psb_spmat_type + use psb_mat_mod use psb_string_mod use psb_serial_mod, psb_protect_name => psb_csymbmm implicit none - type(psb_cspmat_type) :: a,b,c + type(psb_c_sparse_mat), intent(in) :: a,b + type(psb_c_sparse_mat), intent(out) :: c + integer, intent(out) :: info + type(psb_c_csr_sparse_mat), allocatable :: ccsr + integer :: err_act + character(len=*), parameter :: name='psb_symbmm' + call psb_erractionsave(err_act) + info = 0 + + if ((a%is_null()) .or.(b%is_null())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + allocate(ccsr, stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + end if + call psb_symbmm(a%a,b%a,ccsr,info) + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + call move_alloc(ccsr,c%a) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_csymbmm + +subroutine psb_cbase_symbmm(a,b,c,info) + use psb_mat_mod + use psb_serial_mod, psb_protect_name => psb_cbase_symbmm + implicit none + + class(psb_c_base_sparse_mat), intent(in) :: a,b + type(psb_c_csr_sparse_mat), intent(out) :: c + integer, intent(out) :: info integer, allocatable :: itemp(:) - integer :: nze,info - - interface - subroutine symbmm (n, m, l, ia, ja, diaga, & - & ib, jb, diagb, ic, jc, diagc, index) - integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,& - & diagc, index(*) - integer, allocatable :: ic(:),jc(:) - end subroutine symbmm - end interface - - character(len=20) :: name - integer :: err_act - logical :: csra, csrb + integer :: nze, ma,na,mb,nb + character(len=20) :: name + integer :: err_act name='psb_symbmm' call psb_erractionsave(err_act) + info = 0 - csra = (psb_toupper(a%fida(1:3))=='CSR') - csrb = (psb_toupper(b%fida(1:3))=='CSR') + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() - if (b%m /= a%k) then - write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k + + if ( mb /= na ) then + write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb endif - allocate(itemp(max(a%m,a%k,b%m,b%k)),stat=info) + allocate(itemp(max(ma,na,mb,nb)),stat=info) if (info /= 0) then - return + info = 4000 + call psb_Errpush(info,name) + goto 9999 endif - nze = max(a%m+1,2*a%m) - call psb_sp_reall(c,nze,info) ! ! Note: we need to test whether there is a performance impact ! in not using the original Douglas & Bank code. ! - if (csra.and.csrb) then - call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,& - & b%ia2,b%ia1,0,& - & c%ia2,c%ia1,0,itemp) - else - call inner_symbmm(a,b,c,itemp,info) - endif - call psb_realloc(size(c%ia1),c%aspk,info) - - c%pl(1) = 0 - c%pr(1) = 0 - c%m=a%m - c%k=b%k - c%fida='CSR' - c%descra='GUN' + select type(a) + type is (psb_c_csr_sparse_mat) + select type(b) + type is (psb_c_csr_sparse_mat) + call csr_symbmm(a,b,c,itemp,info) + class default + call gen_symbmm(a,b,c,itemp,info) + end select + class default + call gen_symbmm(a,b,c,itemp,info) + end select + + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_realloc(size(c%ja),c%val,info) deallocate(itemp) + call psb_erractionrestore(err_act) return @@ -108,82 +153,119 @@ subroutine psb_csymbmm(a,b,c,info) return contains - subroutine inner_symbmm(a,b,c,index,info) - type(psb_cspmat_type) :: a,b,c + + subroutine csr_symbmm(a,b,c,itemp,info) + type(psb_c_csr_sparse_mat), intent(in) :: a,b + type(psb_c_csr_sparse_mat), intent(out) :: c + integer :: itemp(:) + integer, intent(out) :: info + interface + subroutine symbmm (n, m, l, ia, ja, diaga, & + & ib, jb, diagb, ic, jc, diagc, index) + integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,& + & diagc, index(*) + integer, allocatable :: ic(:),jc(:) + end subroutine symbmm + end interface + integer :: nze, ma,na,mb,nb + + info = 0 + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + call symbmm(ma,na,nb,a%irp,a%ja,0,& + & b%irp,b%ja,0,& + & c%irp,c%ja,0,itemp) + + end subroutine csr_symbmm + subroutine gen_symbmm(a,b,c,index,info) + class(psb_c_base_sparse_mat), intent(in) :: a,b + type(psb_c_csr_sparse_mat), intent(out) :: c integer :: index(:),info integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) - complex(psb_spk_), allocatable :: aval(:),bval(:) integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn + integer :: nze, ma,na,mb,nb + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() - n = a%m - m = a%k - l = b%k + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + + n = ma + m = na + l = nb maxlmn = max(l,m,n) allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& - & aval(maxlmn),bval(maxlmn), stat=info) + & stat=info) if (info /= 0) then + info = 4000 return endif - - if (size(c%ia2) < n+1) then - - call psb_realloc(n+1,c%ia2,info) - endif do i=1,maxlmn index(i)=0 end do - c%ia2(1)=1 - minlm = min(l,m) - minmn = min(m,n) - - main: do i=1,n - istart=-1 - length=0 - call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info) - do jj=1, nazr - - j=iacl(jj) - - if ((j<1).or.(j>m)) then - write(0,*) ' SymbMM: Problem with A ',i,jj,j,m - endif - call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info) - do k=1,nbzr - if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then - write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn - else - if(index(ibcl(k)) == 0) then - index(ibcl(k))=istart - istart=ibcl(k) - length=length+1 - endif - endif - end do - end do + c%irp(1)=1 + minlm = min(l,m) + minmn = min(m,n) + + main: do i=1,n + istart=-1 + length=0 + call a%csget(i,i,nazr,iarw,iacl,info) + do jj=1, nazr + + j=iacl(jj) - c%ia2(i+1)=c%ia2(i)+length - - if (c%ia2(i+1) > size(c%ia1)) then - if (n > (2*i)) then - nze = max(c%ia2(i+1), c%ia2(i)*((n+i-1)/i)) + if ((j<1).or.(j>m)) then + write(0,*) ' SymbMM: Problem with A ',i,jj,j,m + info = 1 + return + endif + call b%csget(j,j,nbzr,ibrw,ibcl,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn + info=2 + return else - nze = max(c%ia2(i+1), nint((dble(c%ia2(i))*(dble(n)/i))) ) - endif - call psb_realloc(nze,c%ia1,info) - end if - do j= c%ia2(i),c%ia2(i+1)-1 - c%ia1(j)=istart - istart=index(istart) - index(c%ia1(j))=0 + if(index(ibcl(k)) == 0) then + index(ibcl(k))=istart + istart=ibcl(k) + length=length+1 + endif + endif end do - call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1)) - index(i) = 0 - end do main + end do - end subroutine inner_symbmm + c%irp(i+1)=c%irp(i)+length -end subroutine psb_csymbmm + if (c%irp(i+1) > size(c%ja)) then + if (n > (2*i)) then + nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i)) + else + nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) ) + endif + call psb_realloc(nze,c%ja,info) + end if + do j= c%irp(i),c%irp(i+1)-1 + c%ja(j)=istart + istart=index(istart) + index(c%ja(j))=0 + end do + call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1)) + index(i) = 0 + end do main + + end subroutine gen_symbmm + +end subroutine psb_cbase_symbmm diff --git a/base/serial/psb_dsymbmm.f90 b/base/serial/psb_dsymbmm.f90 index 211bed40..34952549 100644 --- a/base/serial/psb_dsymbmm.f90 +++ b/base/serial/psb_dsymbmm.f90 @@ -187,7 +187,6 @@ contains type(psb_d_csr_sparse_mat), intent(out) :: c integer :: index(:),info integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) - real(psb_dpk_), allocatable :: aval(:),bval(:) integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn integer :: nze, ma,na,mb,nb @@ -205,7 +204,7 @@ contains maxlmn = max(l,m,n) allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& - & aval(maxlmn),bval(maxlmn), stat=info) + & stat=info) if (info /= 0) then info = 4000 return @@ -222,7 +221,7 @@ contains main: do i=1,n istart=-1 length=0 - call a%csget(i,i,nazr,iarw,iacl,aval,info) + call a%csget(i,i,nazr,iarw,iacl,info) do jj=1, nazr j=iacl(jj) @@ -232,7 +231,7 @@ contains info = 1 return endif - call b%csget(j,j,nbzr,ibrw,ibcl,bval,info) + call b%csget(j,j,nbzr,ibrw,ibcl,info) do k=1,nbzr if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn diff --git a/base/serial/psb_snumbmm.f90 b/base/serial/psb_snumbmm.f90 index e240e5bc..eff2cb3c 100644 --- a/base/serial/psb_snumbmm.f90 +++ b/base/serial/psb_snumbmm.f90 @@ -41,72 +41,173 @@ ! subroutine psb_snumbmm(a,b,c) - use psb_spmat_type + use psb_mat_mod + use psb_string_mod use psb_serial_mod, psb_protect_name => psb_snumbmm - implicit none + implicit none - type(psb_sspmat_type) :: a,b,c - real(psb_spk_), allocatable :: temp(:) - integer :: info - logical :: csra, csrb - - allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info) - if (info /= 0) then - return + type(psb_s_sparse_mat), intent(in) :: a,b + type(psb_s_sparse_mat), intent(inout) :: c + integer :: info + integer :: err_act + character(len=*), parameter :: name='psb_numbmm' + + call psb_erractionsave(err_act) + info = 0 + + if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 endif - call psb_realloc(size(c%ia1),c%aspk,info) + + select type(aa=>c%a) + type is (psb_s_csr_sparse_mat) + call psb_numbmm(a%a,b%a,aa) + class default + info = 1121 + call psb_errpush(info,name) + goto 9999 + end select + + call c%set_asb() + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_snumbmm + +subroutine psb_sbase_numbmm(a,b,c) + use psb_mat_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_sbase_numbmm + implicit none + + class(psb_s_base_sparse_mat), intent(in) :: a,b + type(psb_s_csr_sparse_mat), intent(inout) :: c + integer, allocatable :: itemp(:) + integer :: nze, ma,na,mb,nb + character(len=20) :: name + real(psb_spk_), allocatable :: temp(:) + integer :: info + integer :: err_act + name='psb_numbmm' + call psb_erractionsave(err_act) + info = 0 + + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb + endif + allocate(temp(max(ma,na,mb,nb)),stat=info) + if (info /= 0) then + info = 4000 + call psb_Errpush(info,name) + goto 9999 + endif + ! ! Note: we still have to test about possible performance hits. ! ! - csra = (psb_toupper(a%fida(1:3))=='CSR') - csrb = (psb_toupper(b%fida(1:3))=='CSR') - - if (csra.and.csrb) then - call snumbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,& - & b%ia2,b%ia1,0,b%aspk,& - & c%ia2,c%ia1,0,c%aspk,temp) - else - call inner_numbmm(a,b,c,temp,info) - if (info /= 0) then - write(0,*) 'Error ',info,' from inner numbmm' - end if + call psb_ensure_size(size(c%ja),c%val,info) + select type(a) + type is (psb_s_csr_sparse_mat) + select type(b) + type is (psb_s_csr_sparse_mat) + call csr_numbmm(a,b,c,temp,info) + class default + call gen_numbmm(a,b,c,temp,info) + end select + class default + call gen_numbmm(a,b,c,temp,info) + end select + + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 end if - call psb_sp_setifld(psb_spmat_asb_,psb_state_,c,info) + + call c%set_asb() deallocate(temp) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if return contains - subroutine inner_numbmm(a,b,c,temp,info) - type(psb_sspmat_type) :: a,b,c + subroutine csr_numbmm(a,b,c,temp,info) + type(psb_s_csr_sparse_mat), intent(in) :: a,b + type(psb_s_csr_sparse_mat), intent(inout) :: c + real(psb_spk_) :: temp(:) + integer, intent(out) :: info + integer :: nze, ma,na,mb,nb + + info = 0 + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + call snumbmm(ma,na,nb,a%irp,a%ja,0,a%val,& + & b%irp,b%ja,0,b%val,& + & c%irp,c%ja,0,c%val,temp) + + + end subroutine csr_numbmm + + subroutine gen_numbmm(a,b,c,temp,info) + class(psb_s_base_sparse_mat), intent(in) :: a,b + type(psb_s_csr_sparse_mat), intent(inout) :: c integer :: info - real(psb_spk_) :: temp(:) + real(psb_spk_) :: temp(:) integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) real(psb_spk_), allocatable :: aval(:),bval(:) integer :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln real(psb_spk_) :: ajj - n = a%m - m = a%k - l = b%k + n = a%get_nrows() + m = a%get_ncols() + l = b%get_ncols() maxlmn = max(l,m,n) allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& & aval(maxlmn),bval(maxlmn), stat=info) - if (info /= 0) then + if (info /= 0) then + info = 4000 return endif do i = 1,maxlmn - temp(i) = szero + temp(i) = dzero end do minlm = min(l,m) minln = min(l,n) minmn = min(m,n) do i = 1,n - call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info) - + call a%csget(i,i,nazr,iarw,iacl,aval,info) do jj=1, nazr j=iacl(jj) ajj = aval(jj) @@ -116,7 +217,7 @@ contains return endif - call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info) + call b%csget(j,j,nbzr,ibrw,ibcl,bval,info) do k=1,nbzr if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn @@ -127,19 +228,19 @@ contains endif enddo end do - do j = c%ia2(i),c%ia2(i+1)-1 - if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn)) then - write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn + do j = c%irp(i),c%irp(i+1)-1 + if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then + write(0,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn info = 3 return else - c%aspk(j) = temp(c%ia1(j)) - temp(c%ia1(j)) = szero + c%val(j) = temp(c%ja(j)) + temp(c%ja(j)) = dzero endif end do end do + + end subroutine gen_numbmm - end subroutine inner_numbmm - -end subroutine psb_snumbmm +end subroutine psb_sbase_numbmm diff --git a/base/serial/psb_srwextd.f90 b/base/serial/psb_srwextd.f90 index 89a6d81e..953040cd 100644 --- a/base/serial/psb_srwextd.f90 +++ b/base/serial/psb_srwextd.f90 @@ -39,127 +39,209 @@ ! ! subroutine psb_srwextd(nr,a,info,b,rowscale) - use psb_spmat_type use psb_error_mod use psb_string_mod use psb_serial_mod, psb_protect_name => psb_srwextd implicit none ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes) - integer, intent(in) :: nr - type(psb_sspmat_type), intent(inout) :: a - integer,intent(out) :: info - type(psb_sspmat_type), intent(in), optional :: b - logical,intent(in), optional :: rowscale + integer, intent(in) :: nr + type(psb_s_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + type(psb_s_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale integer :: i,j,ja,jb,err_act,nza,nzb character(len=20) :: name, ch_err + type(psb_s_coo_sparse_mat) :: actmp logical rowscale_ name='psb_srwextd' info = 0 call psb_erractionsave(err_act) + if (nr > a%get_nrows()) then + select type(aa=> a%a) + type is (psb_s_csr_sparse_mat) + if (present(b)) then + call psb_rwextd(nr,aa,info,b%a,rowscale) + else + call psb_rwextd(nr,aa,info,rowscale=rowscale) + end if + type is (psb_s_coo_sparse_mat) + if (present(b)) then + call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale) + else + call psb_rwextd(nr,aa,info,rowscale=rowscale) + end if + class default + call aa%mv_to_coo(actmp,info) + if (info == 0) then + if (present(b)) then + call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale) + else + call psb_rwextd(nr,actmp,info,rowscale=rowscale) + end if + end if + if (info == 0) call aa%mv_from_coo(actmp,info) + end select + end if + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_srwextd +subroutine psb_sbase_rwextd(nr,a,info,b,rowscale) + use psb_error_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_sbase_rwextd + implicit none + + ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes) + integer, intent(in) :: nr + class(psb_s_base_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + class(psb_s_base_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale + + integer :: i,j,ja,jb,err_act,nza,nzb, ma, mb, na, nb + character(len=20) :: name, ch_err + logical rowscale_ + + name='psb_sbase_rwextd' + info = 0 + call psb_erractionsave(err_act) + if (present(rowscale)) then rowscale_ = rowscale else rowscale_ = .true. end if - if (nr > a%m) then - if (psb_toupper(a%fida) == 'CSR') then - call psb_ensure_size(nr+1,a%ia2,info) - if (present(b)) then - nzb = psb_sp_get_nnzeros(b) - call psb_ensure_size(size(a%ia1)+nzb,a%ia1,info) - call psb_ensure_size(size(a%aspk)+nzb,a%aspk,info) - if (psb_toupper(b%fida)=='CSR') then - - do i=1, min(nr-a%m,b%m) - a%ia2(a%m+i+1) = a%ia2(a%m+i) + b%ia2(i+1) - b%ia2(i) - ja = a%ia2(a%m+i) - jb = b%ia2(i) - do - if (jb >= b%ia2(i+1)) exit - a%aspk(ja) = b%aspk(jb) - a%ia1(ja) = b%ia1(jb) - ja = ja + 1 - jb = jb + 1 - end do - end do - do j=i,nr-a%m - a%ia2(a%m+i+1) = a%ia2(a%m+i) + ma = a%get_nrows() + na = a%get_ncols() + + + select type(a) + type is (psb_s_csr_sparse_mat) + + call psb_ensure_size(nr+1,a%irp,info) + + if (present(b)) then + mb = b%get_nrows() + nb = b%get_ncols() + nzb = b%get_nzeros() + + select type (b) + type is (psb_s_csr_sparse_mat) + call psb_ensure_size(size(a%ja)+nzb,a%ja,info) + call psb_ensure_size(size(a%val)+nzb,a%val,info) + do i=1, min(nr-ma,mb) + a%irp(ma+i+1) = a%irp(ma+i) + b%irp(i+1) - b%irp(i) + ja = a%irp(ma+i) + jb = b%irp(i) + do + if (jb >= b%irp(i+1)) exit + a%val(ja) = b%val(jb) + a%ja(ja) = b%ja(jb) + ja = ja + 1 + jb = jb + 1 end do - else - write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' - endif - else - do i=a%m+2,nr+1 - a%ia2(i) = a%ia2(i-1) end do - end if - a%m = nr - a%k = max(a%k,b%k) + do j=i,nr-ma + a%irp(ma+i+1) = a%irp(ma+i) + end do + class default - else if (psb_toupper(a%fida) == 'COO') then + write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' + end select + call a%set_ncols(max(na,nb)) - if (present(b)) then - nza = psb_sp_get_nnzeros(a) - nzb = psb_sp_get_nnzeros(b) - call psb_sp_reall(a,nza+nzb,info) - if (psb_toupper(b%fida)=='COO') then - if (rowscale_) then - do j=1,nzb - if ((a%m + b%ia1(j)) <= nr) then - nza = nza + 1 - a%ia1(nza) = a%m + b%ia1(j) - a%ia2(nza) = b%ia2(j) - a%aspk(nza) = b%aspk(j) - end if - enddo - else - do j=1,nzb - if ((b%ia1(j)) <= nr) then - nza = nza + 1 - a%ia1(nza) = b%ia1(j) - a%ia2(nza) = b%ia2(j) - a%aspk(nza) = b%aspk(j) - endif - enddo - endif - a%infoa(psb_nnz_) = nza - else if(psb_toupper(b%fida)=='CSR') then - do i=1, min(nr-a%m,b%m) - do - jb = b%ia2(i) - if (jb >= b%ia2(i+1)) exit - nza = nza + 1 - a%aspk(nza) = b%aspk(jb) - a%ia1(nza) = a%m + i - a%ia2(nza) = b%ia1(jb) - jb = jb + 1 - end do - end do - a%infoa(psb_nnz_) = nza - else - write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' - endif - endif - a%m = nr - a%k = max(a%k,b%k) - else if (a%fida == 'JAD') then - info=135 - ch_err=a%fida(1:3) - call psb_errpush(info,name,a_err=ch_err) - goto 9999 else - info=136 - ch_err=a%fida(1:3) - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + + do i=ma+2,nr+1 + a%irp(i) = a%irp(i-1) + end do + end if - end if + call a%set_nrows(nr) + + + type is (psb_s_coo_sparse_mat) + nza = a%get_nzeros() + + if (present(b)) then + mb = b%get_nrows() + nb = b%get_ncols() + nzb = b%get_nzeros() + call a%reallocate(nza+nzb) + + select type(b) + type is (psb_s_coo_sparse_mat) + + if (rowscale_) then + do j=1,nzb + if ((ma + b%ia(j)) <= nr) then + nza = nza + 1 + a%ia(nza) = ma + b%ia(j) + a%ja(nza) = b%ja(j) + a%val(nza) = b%val(j) + end if + enddo + else + do j=1,nzb + if ((ma + b%ia(j)) <= nr) then + nza = nza + 1 + a%ia(nza) = b%ia(j) + a%ja(nza) = b%ja(j) + a%val(nza) = b%val(j) + end if + enddo + endif + call a%set_nzeros(nza) + + type is (psb_s_csr_sparse_mat) + + do i=1, min(nr-ma,mb) + do + jb = b%irp(i) + if (jb >= b%irp(i+1)) exit + nza = nza + 1 + a%val(nza) = b%val(jb) + a%ia(nza) = ma + i + a%ja(nza) = b%ja(jb) + jb = jb + 1 + end do + end do + call a%set_nzeros(nza) + + class default + write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' + + end select + + call a%set_ncols(max(na,nb)) + endif + + call a%set_nrows(nr) + + class default + info = 135 + ch_err=a%get_fmt() + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + + end select call psb_erractionrestore(err_act) return @@ -167,9 +249,9 @@ subroutine psb_srwextd(nr,a,info,b,rowscale) 9999 continue call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error() - return + call psb_error() + return end if return -end subroutine psb_srwextd +end subroutine psb_sbase_rwextd diff --git a/base/serial/psb_ssymbmm.f90 b/base/serial/psb_ssymbmm.f90 index 0404b69b..4bda9892 100644 --- a/base/serial/psb_ssymbmm.f90 +++ b/base/serial/psb_ssymbmm.f90 @@ -40,62 +40,107 @@ ! subroutine psb_ssymbmm(a,b,c,info) - use psb_spmat_type + use psb_mat_mod use psb_string_mod use psb_serial_mod, psb_protect_name => psb_ssymbmm implicit none - type(psb_sspmat_type) :: a,b,c + type(psb_s_sparse_mat), intent(in) :: a,b + type(psb_s_sparse_mat), intent(out) :: c + integer, intent(out) :: info + type(psb_s_csr_sparse_mat), allocatable :: ccsr + integer :: err_act + character(len=*), parameter :: name='psb_symbmm' + call psb_erractionsave(err_act) + info = 0 + + if ((a%is_null()) .or.(b%is_null())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + allocate(ccsr, stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + end if + call psb_symbmm(a%a,b%a,ccsr,info) + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + call move_alloc(ccsr,c%a) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_ssymbmm + +subroutine psb_sbase_symbmm(a,b,c,info) + use psb_mat_mod + use psb_serial_mod, psb_protect_name => psb_sbase_symbmm + implicit none + + class(psb_s_base_sparse_mat), intent(in) :: a,b + type(psb_s_csr_sparse_mat), intent(out) :: c + integer, intent(out) :: info integer, allocatable :: itemp(:) - integer :: nze,info - - interface - subroutine symbmm (n, m, l, ia, ja, diaga, & - & ib, jb, diagb, ic, jc, diagc, index) - integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,& - & diagc, index(*) - integer, allocatable :: ic(:),jc(:) - end subroutine symbmm - end interface + integer :: nze, ma,na,mb,nb character(len=20) :: name integer :: err_act - logical :: csra, csrb name='psb_symbmm' call psb_erractionsave(err_act) + info = 0 + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() - csra = (psb_toupper(a%fida(1:3))=='CSR') - csrb = (psb_toupper(b%fida(1:3))=='CSR') - if (b%m /= a%k) then - write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k + if ( mb /= na ) then + write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb endif - allocate(itemp(max(a%m,a%k,b%m,b%k)),stat=info) + allocate(itemp(max(ma,na,mb,nb)),stat=info) if (info /= 0) then - return + info = 4000 + call psb_Errpush(info,name) + goto 9999 endif - nze = max(a%m+1,2*a%m) - call psb_sp_reall(c,nze,info) ! ! Note: we need to test whether there is a performance impact ! in not using the original Douglas & Bank code. ! - if (csra.and.csrb) then - call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,& - & b%ia2,b%ia1,0,& - & c%ia2,c%ia1,0,itemp) - else - call inner_symbmm(a,b,c,itemp,info) - endif - call psb_realloc(size(c%ia1),c%aspk,info) - - c%pl(1) = 0 - c%pr(1) = 0 - c%m=a%m - c%k=b%k - c%fida='CSR' - c%descra='GUN' - + select type(a) + type is (psb_s_csr_sparse_mat) + select type(b) + type is (psb_s_csr_sparse_mat) + call csr_symbmm(a,b,c,itemp,info) + class default + call gen_symbmm(a,b,c,itemp,info) + end select + class default + call gen_symbmm(a,b,c,itemp,info) + end select + + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_realloc(size(c%ja),c%val,info) deallocate(itemp) + call psb_erractionrestore(err_act) return @@ -109,86 +154,118 @@ subroutine psb_ssymbmm(a,b,c,info) contains - subroutine inner_symbmm(a,b,c,index,info) - type(psb_sspmat_type) :: a,b,c + subroutine csr_symbmm(a,b,c,itemp,info) + type(psb_s_csr_sparse_mat), intent(in) :: a,b + type(psb_s_csr_sparse_mat), intent(out) :: c + integer :: itemp(:) + integer, intent(out) :: info + interface + subroutine symbmm (n, m, l, ia, ja, diaga, & + & ib, jb, diagb, ic, jc, diagc, index) + integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,& + & diagc, index(*) + integer, allocatable :: ic(:),jc(:) + end subroutine symbmm + end interface + integer :: nze, ma,na,mb,nb + + info = 0 + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + call symbmm(ma,na,nb,a%irp,a%ja,0,& + & b%irp,b%ja,0,& + & c%irp,c%ja,0,itemp) + + end subroutine csr_symbmm + subroutine gen_symbmm(a,b,c,index,info) + class(psb_s_base_sparse_mat), intent(in) :: a,b + type(psb_s_csr_sparse_mat), intent(out) :: c integer :: index(:),info integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) - real(psb_spk_), allocatable :: aval(:),bval(:) integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn + integer :: nze, ma,na,mb,nb + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() - n = a%m - m = a%k - l = b%k + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + + n = ma + m = na + l = nb maxlmn = max(l,m,n) allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& - & aval(maxlmn),bval(maxlmn), stat=info) + & stat=info) if (info /= 0) then + info = 4000 return endif - - if (size(c%ia2) < n+1) then - - call psb_realloc(n+1,c%ia2,info) - endif do i=1,maxlmn index(i)=0 end do - c%ia2(1)=1 - minlm = min(l,m) - minmn = min(m,n) - - main: do i=1,n - istart=-1 - length=0 - call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info) - do jj=1, nazr - - j=iacl(jj) - - if ((j<1).or.(j>m)) then - write(0,*) ' SymbMM: Problem with A ',i,jj,j,m - info = 1 + c%irp(1)=1 + minlm = min(l,m) + minmn = min(m,n) + + main: do i=1,n + istart=-1 + length=0 + call a%csget(i,i,nazr,iarw,iacl,info) + do jj=1, nazr + + j=iacl(jj) + + if ((j<1).or.(j>m)) then + write(0,*) ' SymbMM: Problem with A ',i,jj,j,m + info = 1 + return + endif + call b%csget(j,j,nbzr,ibrw,ibcl,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn + info=2 return - endif - call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info) - do k=1,nbzr - if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then - write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn - info=2 - return - else - if(index(ibcl(k)) == 0) then - index(ibcl(k))=istart - istart=ibcl(k) - length=length+1 - endif + else + if(index(ibcl(k)) == 0) then + index(ibcl(k))=istart + istart=ibcl(k) + length=length+1 endif - end do + endif end do + end do - c%ia2(i+1)=c%ia2(i)+length - - if (c%ia2(i+1) > size(c%ia1)) then - if (n > (2*i)) then - nze = max(c%ia2(i+1), c%ia2(i)*((n+i-1)/i)) - else - nze = max(c%ia2(i+1), nint((dble(c%ia2(i))*(dble(n)/i))) ) - endif - call psb_realloc(nze,c%ia1,info) - end if - do j= c%ia2(i),c%ia2(i+1)-1 - c%ia1(j)=istart - istart=index(istart) - index(c%ia1(j))=0 - end do - call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1)) - index(i) = 0 - end do main + c%irp(i+1)=c%irp(i)+length - end subroutine inner_symbmm + if (c%irp(i+1) > size(c%ja)) then + if (n > (2*i)) then + nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i)) + else + nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) ) + endif + call psb_realloc(nze,c%ja,info) + end if + do j= c%irp(i),c%irp(i+1)-1 + c%ja(j)=istart + istart=index(istart) + index(c%ja(j))=0 + end do + call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1)) + index(i) = 0 + end do main -end subroutine psb_ssymbmm + end subroutine gen_symbmm + +end subroutine psb_sbase_symbmm diff --git a/base/serial/psb_znumbmm.f90 b/base/serial/psb_znumbmm.f90 index 25353619..909ff748 100644 --- a/base/serial/psb_znumbmm.f90 +++ b/base/serial/psb_znumbmm.f90 @@ -29,7 +29,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: psb_dnumbmm.f90 +! File: psb_znumbmm.f90 ! Subroutine: ! Arguments: ! @@ -41,41 +41,146 @@ ! subroutine psb_znumbmm(a,b,c) - use psb_spmat_type + use psb_mat_mod + use psb_string_mod use psb_serial_mod, psb_protect_name => psb_znumbmm - implicit none + implicit none - type(psb_zspmat_type) :: a,b,c + type(psb_z_sparse_mat), intent(in) :: a,b + type(psb_z_sparse_mat), intent(inout) :: c + integer :: info + integer :: err_act + character(len=*), parameter :: name='psb_numbmm' + + call psb_erractionsave(err_act) + info = 0 + + if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + select type(aa=>c%a) + type is (psb_z_csr_sparse_mat) + call psb_numbmm(a%a,b%a,aa) + class default + info = 1121 + call psb_errpush(info,name) + goto 9999 + end select + + call c%set_asb() + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_znumbmm + +subroutine psb_zbase_numbmm(a,b,c) + use psb_mat_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_zbase_numbmm + implicit none + + class(psb_z_base_sparse_mat), intent(in) :: a,b + type(psb_z_csr_sparse_mat), intent(inout) :: c + integer, allocatable :: itemp(:) + integer :: nze, ma,na,mb,nb + character(len=20) :: name complex(psb_dpk_), allocatable :: temp(:) - integer :: info - logical :: csra, csrb - - allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info) - if (info /= 0) then - return + integer :: info + integer :: err_act + name='psb_numbmm' + call psb_erractionsave(err_act) + info = 0 + + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb endif - call psb_realloc(size(c%ia1),c%aspk,info) + allocate(temp(max(ma,na,mb,nb)),stat=info) + if (info /= 0) then + info = 4000 + call psb_Errpush(info,name) + goto 9999 + endif + ! ! Note: we still have to test about possible performance hits. ! ! - csra = (psb_toupper(a%fida(1:3))=='CSR') - csrb = (psb_toupper(b%fida(1:3))=='CSR') - - if (csra.and.csrb) then - call znumbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,& - & b%ia2,b%ia1,0,b%aspk,& - & c%ia2,c%ia1,0,c%aspk,temp) - else - call inner_numbmm(a,b,c,temp,info) + call psb_ensure_size(size(c%ja),c%val,info) + select type(a) + type is (psb_z_csr_sparse_mat) + select type(b) + type is (psb_z_csr_sparse_mat) + call csr_numbmm(a,b,c,temp,info) + class default + call gen_numbmm(a,b,c,temp,info) + end select + class default + call gen_numbmm(a,b,c,temp,info) + end select + + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 end if + + call c%set_asb() deallocate(temp) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if return contains - subroutine inner_numbmm(a,b,c,temp,info) - type(psb_zspmat_type) :: a,b,c + subroutine csr_numbmm(a,b,c,temp,info) + type(psb_z_csr_sparse_mat), intent(in) :: a,b + type(psb_z_csr_sparse_mat), intent(inout) :: c + complex(psb_dpk_) :: temp(:) + integer, intent(out) :: info + integer :: nze, ma,na,mb,nb + + info = 0 + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + call cnumbmm(ma,na,nb,a%irp,a%ja,0,a%val,& + & b%irp,b%ja,0,b%val,& + & c%irp,c%ja,0,c%val,temp) + + + end subroutine csr_numbmm + + subroutine gen_numbmm(a,b,c,temp,info) + class(psb_z_base_sparse_mat), intent(in) :: a,b + type(psb_z_csr_sparse_mat), intent(inout) :: c integer :: info complex(psb_dpk_) :: temp(:) integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) @@ -83,14 +188,14 @@ contains integer :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln complex(psb_dpk_) :: ajj - - n = a%m - m = a%k - l = b%k + n = a%get_nrows() + m = a%get_ncols() + l = b%get_ncols() maxlmn = max(l,m,n) allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& & aval(maxlmn),bval(maxlmn), stat=info) - if (info /= 0) then + if (info /= 0) then + info = 4000 return endif @@ -102,37 +207,41 @@ contains minmn = min(m,n) do i = 1,n - call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info) - + call a%csget(i,i,nazr,iarw,iacl,aval,info) do jj=1, nazr j=iacl(jj) ajj = aval(jj) if ((j<1).or.(j>m)) then write(0,*) ' NUMBMM: Problem with A ',i,jj,j,m + info = 1 + return + endif - call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info) + call b%csget(j,j,nbzr,ibrw,ibcl,bval,info) do k=1,nbzr if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn + info = 2 + return else temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k) endif enddo end do - do j = c%ia2(i),c%ia2(i+1)-1 - if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn)) then - write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn + do j = c%irp(i),c%irp(i+1)-1 + if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then + write(0,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn + info = 3 + return else - c%aspk(j) = temp(c%ia1(j)) - temp(c%ia1(j)) = dzero + c%val(j) = temp(c%ja(j)) + temp(c%ja(j)) = dzero endif end do end do + + end subroutine gen_numbmm +end subroutine psb_zbase_numbmm - - end subroutine inner_numbmm - - -end subroutine psb_znumbmm diff --git a/base/serial/psb_zrwextd.f90 b/base/serial/psb_zrwextd.f90 index 3ea707d1..0b76c11b 100644 --- a/base/serial/psb_zrwextd.f90 +++ b/base/serial/psb_zrwextd.f90 @@ -39,126 +39,209 @@ ! ! subroutine psb_zrwextd(nr,a,info,b,rowscale) - use psb_spmat_type use psb_error_mod use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_zrwextd implicit none ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes) - integer, intent(in) :: nr - type(psb_zspmat_type), intent(inout) :: a - integer,intent(out) :: info - type(psb_zspmat_type), intent(in), optional :: b - logical,intent(in), optional :: rowscale + integer, intent(in) :: nr + type(psb_z_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + type(psb_z_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale integer :: i,j,ja,jb,err_act,nza,nzb character(len=20) :: name, ch_err + type(psb_z_coo_sparse_mat) :: actmp logical rowscale_ name='psb_zrwextd' info = 0 call psb_erractionsave(err_act) + if (nr > a%get_nrows()) then + select type(aa=> a%a) + type is (psb_z_csr_sparse_mat) + if (present(b)) then + call psb_rwextd(nr,aa,info,b%a,rowscale) + else + call psb_rwextd(nr,aa,info,rowscale=rowscale) + end if + type is (psb_z_coo_sparse_mat) + if (present(b)) then + call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale) + else + call psb_rwextd(nr,aa,info,rowscale=rowscale) + end if + class default + call aa%mv_to_coo(actmp,info) + if (info == 0) then + if (present(b)) then + call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale) + else + call psb_rwextd(nr,actmp,info,rowscale=rowscale) + end if + end if + if (info == 0) call aa%mv_from_coo(actmp,info) + end select + end if + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_zrwextd +subroutine psb_zbase_rwextd(nr,a,info,b,rowscale) + use psb_error_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_zbase_rwextd + implicit none + + ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes) + integer, intent(in) :: nr + class(psb_z_base_sparse_mat), intent(inout) :: a + integer,intent(out) :: info + class(psb_z_base_sparse_mat), intent(in), optional :: b + logical,intent(in), optional :: rowscale + + integer :: i,j,ja,jb,err_act,nza,nzb, ma, mb, na, nb + character(len=20) :: name, ch_err + logical rowscale_ + + name='psb_zbase_rwextd' + info = 0 + call psb_erractionsave(err_act) + if (present(rowscale)) then rowscale_ = rowscale else rowscale_ = .true. end if - if (nr > a%m) then - if (psb_toupper(a%fida) == 'CSR') then - call psb_realloc(nr+1,a%ia2,info) - if (present(b)) then - nzb = psb_sp_get_nnzeros(b) - call psb_realloc(size(a%ia1)+nzb,a%ia1,info) - call psb_realloc(size(a%aspk)+nzb,a%aspk,info) - if (psb_toupper(b%fida)=='CSR') then - - do i=1, min(nr-a%m,b%m) - a%ia2(a%m+i+1) = a%ia2(a%m+i) + b%ia2(i+1) - b%ia2(i) - ja = a%ia2(a%m+i) - jb = b%ia2(i) - do - if (jb >= b%ia2(i+1)) exit - a%aspk(ja) = b%aspk(jb) - a%ia1(ja) = b%ia1(jb) - ja = ja + 1 - jb = jb + 1 - end do - end do - do j=i,nr-a%m - a%ia2(a%m+i+1) = a%ia2(a%m+i) + ma = a%get_nrows() + na = a%get_ncols() + + + select type(a) + type is (psb_z_csr_sparse_mat) + + call psb_ensure_size(nr+1,a%irp,info) + + if (present(b)) then + mb = b%get_nrows() + nb = b%get_ncols() + nzb = b%get_nzeros() + + select type (b) + type is (psb_z_csr_sparse_mat) + call psb_ensure_size(size(a%ja)+nzb,a%ja,info) + call psb_ensure_size(size(a%val)+nzb,a%val,info) + do i=1, min(nr-ma,mb) + a%irp(ma+i+1) = a%irp(ma+i) + b%irp(i+1) - b%irp(i) + ja = a%irp(ma+i) + jb = b%irp(i) + do + if (jb >= b%irp(i+1)) exit + a%val(ja) = b%val(jb) + a%ja(ja) = b%ja(jb) + ja = ja + 1 + jb = jb + 1 end do - else - write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' - endif - else - do i=a%m+2,nr+1 - a%ia2(i) = a%ia2(i-1) end do - end if - a%m = nr - a%k = max(a%k,b%k) + do j=i,nr-ma + a%irp(ma+i+1) = a%irp(ma+i) + end do + class default - else if (psb_toupper(a%fida) == 'COO') then + write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' + end select + call a%set_ncols(max(na,nb)) - if (present(b)) then - nza = psb_sp_get_nnzeros(a) - nzb = psb_sp_get_nnzeros(b) - call psb_sp_reall(a,nza+nzb,info) - if (psb_toupper(b%fida)=='COO') then - if (rowscale_) then - do j=1,nzb - if ((a%m + b%ia1(j)) <= nr) then - nza = nza + 1 - a%ia1(nza) = a%m + b%ia1(j) - a%ia2(nza) = b%ia2(j) - a%aspk(nza) = b%aspk(j) - end if - enddo - else - do j=1,nzb - if ((b%ia1(j)) <= nr) then - nza = nza + 1 - a%ia1(nza) = b%ia1(j) - a%ia2(nza) = b%ia2(j) - a%aspk(nza) = b%aspk(j) - endif - enddo - endif - a%infoa(psb_nnz_) = nza - else if(psb_toupper(b%fida)=='CSR') then - do i=1, min(nr-a%m,b%m) - do - jb = b%ia2(i) - if (jb >= b%ia2(i+1)) exit - nza = nza + 1 - a%aspk(nza) = b%aspk(jb) - a%ia1(nza) = a%m + i - a%ia2(nza) = b%ia1(jb) - jb = jb + 1 - end do - end do - a%infoa(psb_nnz_) = nza - else - write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' - endif - endif - a%m = nr - a%k = max(a%k,b%k) - else if (a%fida == 'JAD') then - info=135 - ch_err=a%fida(1:3) - call psb_errpush(info,name,a_err=ch_err) - goto 9999 else - info=136 - ch_err=a%fida(1:3) - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + + do i=ma+2,nr+1 + a%irp(i) = a%irp(i-1) + end do + end if - end if + call a%set_nrows(nr) + + + type is (psb_z_coo_sparse_mat) + nza = a%get_nzeros() + + if (present(b)) then + mb = b%get_nrows() + nb = b%get_ncols() + nzb = b%get_nzeros() + call a%reallocate(nza+nzb) + + select type(b) + type is (psb_z_coo_sparse_mat) + + if (rowscale_) then + do j=1,nzb + if ((ma + b%ia(j)) <= nr) then + nza = nza + 1 + a%ia(nza) = ma + b%ia(j) + a%ja(nza) = b%ja(j) + a%val(nza) = b%val(j) + end if + enddo + else + do j=1,nzb + if ((ma + b%ia(j)) <= nr) then + nza = nza + 1 + a%ia(nza) = b%ia(j) + a%ja(nza) = b%ja(j) + a%val(nza) = b%val(j) + end if + enddo + endif + call a%set_nzeros(nza) + + type is (psb_z_csr_sparse_mat) + + do i=1, min(nr-ma,mb) + do + jb = b%irp(i) + if (jb >= b%irp(i+1)) exit + nza = nza + 1 + a%val(nza) = b%val(jb) + a%ia(nza) = ma + i + a%ja(nza) = b%ja(jb) + jb = jb + 1 + end do + end do + call a%set_nzeros(nza) + + class default + write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!' + + end select + + call a%set_ncols(max(na,nb)) + endif + + call a%set_nrows(nr) + + class default + info = 135 + ch_err=a%get_fmt() + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + + end select call psb_erractionrestore(err_act) return @@ -166,9 +249,9 @@ subroutine psb_zrwextd(nr,a,info,b,rowscale) 9999 continue call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then - call psb_error() - return + call psb_error() + return end if return -end subroutine psb_zrwextd +end subroutine psb_zbase_rwextd diff --git a/base/serial/psb_zsymbmm.f90 b/base/serial/psb_zsymbmm.f90 index b038b297..287a4eed 100644 --- a/base/serial/psb_zsymbmm.f90 +++ b/base/serial/psb_zsymbmm.f90 @@ -40,62 +40,107 @@ ! subroutine psb_zsymbmm(a,b,c,info) - use psb_spmat_type + use psb_mat_mod use psb_string_mod use psb_serial_mod, psb_protect_name => psb_zsymbmm implicit none - type(psb_zspmat_type) :: a,b,c + type(psb_z_sparse_mat), intent(in) :: a,b + type(psb_z_sparse_mat), intent(out) :: c + integer, intent(out) :: info + type(psb_z_csr_sparse_mat), allocatable :: ccsr + integer :: err_act + character(len=*), parameter :: name='psb_symbmm' + call psb_erractionsave(err_act) + info = 0 + + if ((a%is_null()) .or.(b%is_null())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + allocate(ccsr, stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + end if + call psb_symbmm(a%a,b%a,ccsr,info) + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + call move_alloc(ccsr,c%a) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_zsymbmm + +subroutine psb_zbase_symbmm(a,b,c,info) + use psb_mat_mod + use psb_serial_mod, psb_protect_name => psb_zbase_symbmm + implicit none + + class(psb_z_base_sparse_mat), intent(in) :: a,b + type(psb_z_csr_sparse_mat), intent(out) :: c + integer, intent(out) :: info integer, allocatable :: itemp(:) - integer :: nze,info - - interface - subroutine symbmm (n, m, l, ia, ja, diaga, & - & ib, jb, diagb, ic, jc, diagc, index) - integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,& - & diagc, index(*) - integer, allocatable :: ic(:),jc(:) - end subroutine symbmm - end interface - - character(len=20) :: name - integer :: err_act - logical :: csra, csrb + integer :: nze, ma,na,mb,nb + character(len=20) :: name + integer :: err_act name='psb_symbmm' call psb_erractionsave(err_act) + info = 0 - csra = (psb_toupper(a%fida(1:3))=='CSR') - csrb = (psb_toupper(b%fida(1:3))=='CSR') + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() - if (b%m /= a%k) then - write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k + + if ( mb /= na ) then + write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb endif - allocate(itemp(max(a%m,a%k,b%m,b%k)),stat=info) + allocate(itemp(max(ma,na,mb,nb)),stat=info) if (info /= 0) then - return + info = 4000 + call psb_Errpush(info,name) + goto 9999 endif - nze = max(a%m+1,2*a%m) - call psb_sp_reall(c,nze,info) ! ! Note: we need to test whether there is a performance impact ! in not using the original Douglas & Bank code. ! - if (csra.and.csrb) then - call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,& - & b%ia2,b%ia1,0,& - & c%ia2,c%ia1,0,itemp) - else - call inner_symbmm(a,b,c,itemp,info) - endif - call psb_realloc(size(c%ia1),c%aspk,info) - - c%pl(1) = 0 - c%pr(1) = 0 - c%m=a%m - c%k=b%k - c%fida='CSR' - c%descra='GUN' + select type(a) + type is (psb_z_csr_sparse_mat) + select type(b) + type is (psb_z_csr_sparse_mat) + call csr_symbmm(a,b,c,itemp,info) + class default + call gen_symbmm(a,b,c,itemp,info) + end select + class default + call gen_symbmm(a,b,c,itemp,info) + end select + + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_realloc(size(c%ja),c%val,info) deallocate(itemp) + call psb_erractionrestore(err_act) return @@ -108,82 +153,119 @@ subroutine psb_zsymbmm(a,b,c,info) return contains - subroutine inner_symbmm(a,b,c,index,info) - type(psb_zspmat_type) :: a,b,c + + subroutine csr_symbmm(a,b,c,itemp,info) + type(psb_z_csr_sparse_mat), intent(in) :: a,b + type(psb_z_csr_sparse_mat), intent(out) :: c + integer :: itemp(:) + integer, intent(out) :: info + interface + subroutine symbmm (n, m, l, ia, ja, diaga, & + & ib, jb, diagb, ic, jc, diagc, index) + integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,& + & diagc, index(*) + integer, allocatable :: ic(:),jc(:) + end subroutine symbmm + end interface + integer :: nze, ma,na,mb,nb + + info = 0 + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + call symbmm(ma,na,nb,a%irp,a%ja,0,& + & b%irp,b%ja,0,& + & c%irp,c%ja,0,itemp) + + end subroutine csr_symbmm + subroutine gen_symbmm(a,b,c,index,info) + class(psb_z_base_sparse_mat), intent(in) :: a,b + type(psb_z_csr_sparse_mat), intent(out) :: c integer :: index(:),info integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) - complex(psb_dpk_), allocatable :: aval(:),bval(:) integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn + integer :: nze, ma,na,mb,nb + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() - n = a%m - m = a%k - l = b%k + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + + n = ma + m = na + l = nb maxlmn = max(l,m,n) allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& - & aval(maxlmn),bval(maxlmn), stat=info) + & stat=info) if (info /= 0) then + info = 4000 return endif - - if (size(c%ia2) < n+1) then - - call psb_realloc(n+1,c%ia2,info) - endif do i=1,maxlmn index(i)=0 end do - c%ia2(1)=1 - minlm = min(l,m) - minmn = min(m,n) - - main: do i=1,n - istart=-1 - length=0 - call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info) - do jj=1, nazr - - j=iacl(jj) - - if ((j<1).or.(j>m)) then - write(0,*) ' SymbMM: Problem with A ',i,jj,j,m - endif - call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info) - do k=1,nbzr - if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then - write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn - else - if(index(ibcl(k)) == 0) then - index(ibcl(k))=istart - istart=ibcl(k) - length=length+1 - endif - endif - end do - end do + c%irp(1)=1 + minlm = min(l,m) + minmn = min(m,n) + + main: do i=1,n + istart=-1 + length=0 + call a%csget(i,i,nazr,iarw,iacl,info) + do jj=1, nazr + + j=iacl(jj) - c%ia2(i+1)=c%ia2(i)+length - - if (c%ia2(i+1) > size(c%ia1)) then - if (n > (2*i)) then - nze = max(c%ia2(i+1), c%ia2(i)*((n+i-1)/i)) + if ((j<1).or.(j>m)) then + write(0,*) ' SymbMM: Problem with A ',i,jj,j,m + info = 1 + return + endif + call b%csget(j,j,nbzr,ibrw,ibcl,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn + info=2 + return else - nze = max(c%ia2(i+1), nint((dble(c%ia2(i))*(dble(n)/i))) ) - endif - call psb_realloc(nze,c%ia1,info) - end if - do j= c%ia2(i),c%ia2(i+1)-1 - c%ia1(j)=istart - istart=index(istart) - index(c%ia1(j))=0 + if(index(ibcl(k)) == 0) then + index(ibcl(k))=istart + istart=ibcl(k) + length=length+1 + endif + endif end do - call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1)) - index(i) = 0 - end do main + end do - end subroutine inner_symbmm + c%irp(i+1)=c%irp(i)+length -end subroutine psb_zsymbmm + if (c%irp(i+1) > size(c%ja)) then + if (n > (2*i)) then + nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i)) + else + nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) ) + endif + call psb_realloc(nze,c%ja,info) + end if + do j= c%irp(i),c%irp(i+1)-1 + c%ja(j)=istart + istart=index(istart) + index(c%ja(j))=0 + end do + call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1)) + index(i) = 0 + end do main + + end subroutine gen_symbmm + +end subroutine psb_zbase_symbmm