From 4a5acb014cdc98d6866d99ff38160b1b31757bd3 Mon Sep 17 00:00:00 2001 From: sfilippone Date: Fri, 5 Jun 2026 16:33:25 +0200 Subject: [PATCH] Put mmio_impl into templates --- util/psb_c_mmio_impl.f90 | 171 ++++++++++---------- util/psb_d_mmio_impl.f90 | 53 +++++-- util/psb_i_mmio_impl.f90 | 328 +++++++++++++++++++++++++++++++++++++++ util/psb_s_mmio_impl.f90 | 124 +++++++++++---- util/psb_z_mmio_impl.f90 | 171 ++++++++++---------- 5 files changed, 644 insertions(+), 203 deletions(-) create mode 100644 util/psb_i_mmio_impl.f90 diff --git a/util/psb_c_mmio_impl.f90 b/util/psb_c_mmio_impl.f90 index 68690d06c..ad653462a 100644 --- a/util/psb_c_mmio_impl.f90 +++ b/util/psb_c_mmio_impl.f90 @@ -36,12 +36,13 @@ subroutine mm_cvet_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_spk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -52,6 +53,7 @@ subroutine mm_cvet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -76,15 +78,15 @@ subroutine mm_cvet_read(b, info, iunit, filename) read(line,fmt=*)nrow,ncol - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow),stat = ircode) if (ircode /= 0) goto 993 - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i) = cmplx(bre,bim,kind=psb_spk_) + do i=1,nrow + read(infile,fmt=*,end=902) b(i) end do + end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -109,12 +111,13 @@ subroutine mm_cvet2_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_spk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -125,6 +128,7 @@ subroutine mm_cvet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -152,15 +156,10 @@ subroutine mm_cvet2_read(b, info, iunit, filename) if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow,ncol),stat = ircode) if (ircode /= 0) goto 993 - do j=1, ncol - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i,j) = cmplx(bre,bim,kind=psb_spk_) - end do - end do + read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -189,8 +188,10 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -201,6 +202,7 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -209,17 +211,19 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename) outfile=6 endif endif - - write(outfile,'(a)') '%%MatrixMarket matrix array complex general' + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = size(b,2) - write(outfile,*) nrow,ncol + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' - write(outfile,fmt='(2(es26.18,1x))') ((b(i,j), i=1,nrow),j=1,ncol) + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -241,8 +245,10 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -253,6 +259,7 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -262,20 +269,20 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename) endif endif - write(outfile,'(a)') '%%MatrixMarket matrix array complex general' + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = 1 write(outfile,*) nrow,ncol - write(frmtv,'(a,i0,a)') '(',2*ncol,'(es26.18,1x))' + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' do i=1,size(b,1) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -298,7 +305,7 @@ subroutine mm_cvect_read(b, info, iunit, filename) complex(psb_spk_), allocatable :: bv(:) call mm_array_read(bv, info, iunit, filename) - call b%bld(bv) + if (info == 0) call b%bld(bv) end subroutine mm_cvect_read @@ -331,9 +338,10 @@ subroutine cmm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_c_coo_sparse_mat), allocatable :: acoo - real(psb_spk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -344,6 +352,7 @@ subroutine cmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -369,24 +378,27 @@ subroutine cmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = cone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -398,37 +410,34 @@ subroutine cmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) - - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = cone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -456,10 +465,11 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -470,6 +480,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -481,7 +492,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -492,6 +503,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) return end subroutine cmm_mat_write + subroutine lcmm_mat_read(a, info, iunit, filename) use psb_base_mod implicit none @@ -501,12 +513,13 @@ subroutine lcmm_mat_read(a, info, iunit, filename) character(len=*), optional, intent(in) :: filename character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character(1024) :: line - integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr - integer(psb_ipk_) :: ircode,infile + integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr + integer(psb_ipk_) :: ircode, infile type(psb_lc_coo_sparse_mat), allocatable :: acoo - real(psb_spk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -517,6 +530,7 @@ subroutine lcmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -542,24 +556,27 @@ subroutine lcmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = cone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -571,37 +588,34 @@ subroutine lcmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) - - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = cone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -629,10 +643,10 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -643,6 +657,7 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -654,7 +669,7 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return diff --git a/util/psb_d_mmio_impl.f90 b/util/psb_d_mmio_impl.f90 index 374eeed92..4fc917a61 100644 --- a/util/psb_d_mmio_impl.f90 +++ b/util/psb_d_mmio_impl.f90 @@ -39,8 +39,10 @@ subroutine mm_dvet_read(b, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -51,6 +53,7 @@ subroutine mm_dvet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -83,7 +86,7 @@ subroutine mm_dvet_read(b, info, iunit, filename) end do end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -111,8 +114,10 @@ subroutine mm_dvet2_read(b, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -123,6 +128,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -153,7 +159,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename) read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -182,8 +188,10 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -194,6 +202,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -202,7 +211,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) outfile=6 endif endif - + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' @@ -210,9 +219,11 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) ncol = size(b,2) write(outfile,*) nrow, ncol - write(outfile,fmt='(es26.18,1x)') ((b(i,j), i=1,nrow),j=1,ncol) + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' + + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -234,8 +245,10 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -246,6 +259,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -268,7 +282,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -324,9 +338,10 @@ subroutine dmm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_d_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -337,6 +352,7 @@ subroutine dmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -421,7 +437,7 @@ subroutine dmm_mat_read(a, info, iunit, filename) call a%mv_from(acoo) end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -449,10 +465,11 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -463,6 +480,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -474,7 +492,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -498,9 +516,10 @@ subroutine ldmm_mat_read(a, info, iunit, filename) integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr integer(psb_ipk_) :: ircode, infile type(psb_ld_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -511,6 +530,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -595,7 +615,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename) call a%mv_from(acoo) end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -623,10 +643,10 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -637,6 +657,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -648,7 +669,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -658,5 +679,3 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) write(psb_err_unit,*) 'Error while opening ',filename return end subroutine ldmm_mat_write - - diff --git a/util/psb_i_mmio_impl.f90 b/util/psb_i_mmio_impl.f90 new file mode 100644 index 000000000..64d02bcd7 --- /dev/null +++ b/util/psb_i_mmio_impl.f90 @@ -0,0 +1,328 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine mm_ivet_read(b, info, iunit, filename) + use psb_base_mod + implicit none + integer(psb_ipk_), allocatable, intent(out) :: b(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile + character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& + & line*1024 + logical :: opened + + info = psb_success_ + opened = .false. + if (present(filename)) then + if (filename == '-') then + infile=5 + else + if (present(iunit)) then + infile=iunit + else + infile=99 + endif + open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. + endif + else + if (present(iunit)) then + infile=iunit + else + infile=5 + endif + endif + + read(infile,fmt=*, end=902) mmheader, object, fmt, type, sym + + if ( (object /= 'matrix').or.(fmt /= 'array')) then + write(psb_err_unit,*) 'read_rhs: input file type not yet supported' + info = -3 + return + end if + + do + read(infile,fmt='(a)') line + if (line(1:1) /= '%') exit + end do + + read(line,fmt=*)nrow,ncol + + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then + allocate(b(nrow),stat = ircode) + if (ircode /= 0) goto 993 + do i=1,nrow + read(infile,fmt=*,end=902) b(i) + end do + + end if ! read right hand sides + if (opened) close(infile) + + return + ! open failed +901 write(psb_err_unit,*) 'mm_vet_read: could not open file ',& + & infile,' for input' + info = -1 + return + +902 write(psb_err_unit,*) 'mmv_vet_read: unexpected end of file ',infile,& + & ' during input' + info = -2 + return +993 write(psb_err_unit,*) 'mm_vet_read: memory allocation failure' + info = -3 + return +end subroutine mm_ivet_read + +subroutine mm_ivet2_read(b, info, iunit, filename) + use psb_base_mod + implicit none + integer(psb_ipk_), allocatable, intent(out) :: b(:,:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile + character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& + & line*1024 + logical :: opened + + info = psb_success_ + opened = .false. + if (present(filename)) then + if (filename == '-') then + infile=5 + else + if (present(iunit)) then + infile=iunit + else + infile=99 + endif + open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. + endif + else + if (present(iunit)) then + infile=iunit + else + infile=5 + endif + endif + + read(infile,fmt=*, end=902) mmheader, object, fmt, type, sym + + if ( (object /= 'matrix').or.(fmt /= 'array')) then + write(psb_err_unit,*) 'read_rhs: input file type not yet supported' + info = -3 + return + end if + + do + read(infile,fmt='(a)') line + if (line(1:1) /= '%') exit + end do + + read(line,fmt=*)nrow,ncol + + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then + allocate(b(nrow,ncol),stat = ircode) + if (ircode /= 0) goto 993 + read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) + + end if ! read right hand sides + if (opened) close(infile) + + return + ! open failed +901 write(psb_err_unit,*) 'mm_vet_read: could not open file ',& + & infile,' for input' + info = -1 + return + +902 write(psb_err_unit,*) 'mmv_vet_read: unexpected end of file ',infile,& + & ' during input' + info = -2 + return +993 write(psb_err_unit,*) 'mm_vet_read: memory allocation failure' + info = -3 + return +end subroutine mm_ivet2_read + +subroutine mm_ivet2_write(b, header, info, iunit, filename) + use psb_base_mod + implicit none + integer(psb_ipk_), intent(in) :: b(:,:) + character(len=*), intent(in) :: header + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile + + character(len=80) :: frmtv + logical :: opened + + info = psb_success_ + opened = .false. + if (present(filename)) then + if (filename == '-') then + outfile=6 + else + if (present(iunit)) then + outfile=iunit + else + outfile=99 + endif + open(outfile,file=filename, err=901, action='WRITE') + opened = .true. + endif + else + if (present(iunit)) then + outfile=iunit + else + outfile=6 + endif + endif + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' + write(outfile,'(a)') '% '//trim(header) + write(outfile,'(a)') '% ' + nrow = size(b,1) + ncol = size(b,2) + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' + + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) + + if (opened) close(outfile) + + return + ! open failed +901 write(psb_err_unit,*) 'mm_vet_write: could not open file ',& + & outfile,' for output' + info = -1 + return + +end subroutine mm_ivet2_write + +subroutine mm_ivet1_write(b, header, info, iunit, filename) + use psb_base_mod + implicit none + integer(psb_ipk_), intent(in) :: b(:) + character(len=*), intent(in) :: header + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile + + character(len=80) :: frmtv + logical :: opened + + info = psb_success_ + opened = .false. + if (present(filename)) then + if (filename == '-') then + outfile=6 + else + if (present(iunit)) then + outfile=iunit + else + outfile=99 + endif + open(outfile,file=filename, err=901, action='WRITE') + opened = .true. + endif + else + if (present(iunit)) then + outfile=iunit + else + outfile=6 + endif + endif + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' + write(outfile,'(a)') '% '//trim(header) + write(outfile,'(a)') '% ' + nrow = size(b,1) + ncol = 1 + write(outfile,*) nrow,ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' + + do i=1,size(b,1) + write(outfile,frmtv) b(i) + end do + + if (opened) close(outfile) + + return + ! open failed +901 write(psb_err_unit,*) 'mm_vet_write: could not open file ',& + & outfile,' for output' + info = -1 + return + +end subroutine mm_ivet1_write + +subroutine mm_ivect_read(b, info, iunit, filename) + use psb_base_mod + use psb_mmio_mod, psb_protect_name => mm_ivect_read + implicit none + type(psb_i_vect_type), intent(inout) :: b + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + ! + integer(psb_ipk_), allocatable :: bv(:) + + call mm_array_read(bv, info, iunit, filename) + if (info == 0) call b%bld(bv) + +end subroutine mm_ivect_read + +subroutine mm_ivect_write(b, header, info, iunit, filename) + use psb_base_mod + use psb_mmio_mod, psb_protect_name => mm_ivect_write + implicit none + type(psb_i_vect_type), intent(inout) :: b + character(len=*), intent(in) :: header + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + info = psb_success_ + if (.not.allocated(b%v)) return + call b%sync() + + call mm_array_write(b%v%v,header,info,iunit,filename) + +end subroutine mm_ivect_write + diff --git a/util/psb_s_mmio_impl.f90 b/util/psb_s_mmio_impl.f90 index 94e886849..0ffd9a93d 100644 --- a/util/psb_s_mmio_impl.f90 +++ b/util/psb_s_mmio_impl.f90 @@ -36,11 +36,13 @@ subroutine mm_svet_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -51,6 +53,7 @@ subroutine mm_svet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -83,8 +86,7 @@ subroutine mm_svet_read(b, info, iunit, filename) end do end if ! read right hand sides - - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -109,11 +111,13 @@ subroutine mm_svet2_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -124,6 +128,7 @@ subroutine mm_svet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -154,8 +159,7 @@ subroutine mm_svet2_read(b, info, iunit, filename) read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -184,8 +188,10 @@ subroutine mm_svet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -196,6 +202,7 @@ subroutine mm_svet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -204,17 +211,19 @@ subroutine mm_svet2_write(b, header, info, iunit, filename) outfile=6 endif endif - + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = size(b,2) - write(outfile,*) nrow,ncol + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' - write(outfile,fmt='(es26.18,1x)') ((b(i,j), i=1,nrow),j=1,ncol) + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -236,8 +245,10 @@ subroutine mm_svet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -248,6 +259,7 @@ subroutine mm_svet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -270,7 +282,7 @@ subroutine mm_svet1_write(b, header, info, iunit, filename) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -326,9 +338,10 @@ subroutine smm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_s_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -339,6 +352,7 @@ subroutine smm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -370,9 +384,14 @@ subroutine smm_mat_read(a, info, iunit, filename) read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + call acoo%set_nzeros(nnzero) else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so @@ -391,17 +410,34 @@ subroutine smm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then + call acoo%allocate(nrow,ncol,2*nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + nzr = nnzero + do i=1,nnzero + if (acoo%ia(i) /= acoo%ja(i)) then + nzr = nzr + 1 + acoo%ia(nzr) = acoo%ja(i) + acoo%ja(nzr) = acoo%ia(i) + end if + end do + call acoo%set_nzeros(nzr) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -429,10 +465,11 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -443,6 +480,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -454,7 +492,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -465,6 +503,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) return end subroutine smm_mat_write + subroutine lsmm_mat_read(a, info, iunit, filename) use psb_base_mod implicit none @@ -474,12 +513,13 @@ subroutine lsmm_mat_read(a, info, iunit, filename) character(len=*), optional, intent(in) :: filename character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character(1024) :: line - integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr - integer(psb_ipk_) :: ircode,infile + integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr + integer(psb_ipk_) :: ircode, infile type(psb_ls_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -490,6 +530,7 @@ subroutine lsmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -521,9 +562,14 @@ subroutine lsmm_mat_read(a, info, iunit, filename) read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + call acoo%set_nzeros(nnzero) else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so @@ -542,17 +588,34 @@ subroutine lsmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then + call acoo%allocate(nrow,ncol,2*nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + nzr = nnzero + do i=1,nnzero + if (acoo%ia(i) /= acoo%ja(i)) then + nzr = nzr + 1 + acoo%ia(nzr) = acoo%ja(i) + acoo%ja(nzr) = acoo%ia(i) + end if + end do + call acoo%set_nzeros(nzr) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -580,10 +643,10 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -594,6 +657,7 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -605,7 +669,7 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return diff --git a/util/psb_z_mmio_impl.f90 b/util/psb_z_mmio_impl.f90 index 948de2833..aa66f0edd 100644 --- a/util/psb_z_mmio_impl.f90 +++ b/util/psb_z_mmio_impl.f90 @@ -36,12 +36,13 @@ subroutine mm_zvet_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_dpk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -52,6 +53,7 @@ subroutine mm_zvet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -76,15 +78,15 @@ subroutine mm_zvet_read(b, info, iunit, filename) read(line,fmt=*)nrow,ncol - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow),stat = ircode) if (ircode /= 0) goto 993 - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i) = cmplx(bre,bim,kind=psb_dpk_) + do i=1,nrow + read(infile,fmt=*,end=902) b(i) end do + end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -109,12 +111,13 @@ subroutine mm_zvet2_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_dpk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -125,6 +128,7 @@ subroutine mm_zvet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -149,18 +153,13 @@ subroutine mm_zvet2_read(b, info, iunit, filename) read(line,fmt=*)nrow,ncol - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow,ncol),stat = ircode) if (ircode /= 0) goto 993 - do j=1, ncol - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i,j) = cmplx(bre,bim,kind=psb_dpk_) - end do - end do + read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -189,8 +188,10 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -201,6 +202,7 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -209,17 +211,19 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename) outfile=6 endif endif - - write(outfile,'(a)') '%%MatrixMarket matrix array complex general' + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = size(b,2) - write(outfile,*) nrow,ncol + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' - write(outfile,fmt='(2(es26.18,1x))') ((b(i,j), i=1,nrow),j=1,ncol) + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -241,8 +245,10 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -253,6 +259,7 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -269,13 +276,13 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename) ncol = 1 write(outfile,*) nrow,ncol - write(frmtv,'(a,i0,a)') '(',2*ncol,'(es26.18,1x))' + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' do i=1,size(b,1) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -331,9 +338,10 @@ subroutine zmm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_z_coo_sparse_mat), allocatable :: acoo - real(psb_dpk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -344,6 +352,7 @@ subroutine zmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -369,24 +378,27 @@ subroutine zmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = zone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -398,37 +410,34 @@ subroutine zmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = zone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -446,6 +455,7 @@ subroutine zmm_mat_read(a, info, iunit, filename) return end subroutine zmm_mat_read + subroutine zmm_mat_write(a,mtitle,info,iunit,filename) use psb_base_mod implicit none @@ -455,10 +465,11 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -469,6 +480,7 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -480,7 +492,7 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -501,12 +513,13 @@ subroutine lzmm_mat_read(a, info, iunit, filename) character(len=*), optional, intent(in) :: filename character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character(1024) :: line - integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr - integer(psb_ipk_) :: ircode,infile + integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr + integer(psb_ipk_) :: ircode, infile type(psb_lz_coo_sparse_mat), allocatable :: acoo - real(psb_dpk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -517,6 +530,7 @@ subroutine lzmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -542,24 +556,27 @@ subroutine lzmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = zone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -571,37 +588,34 @@ subroutine lzmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) - - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = zone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -619,6 +633,7 @@ subroutine lzmm_mat_read(a, info, iunit, filename) return end subroutine lzmm_mat_read + subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) use psb_base_mod implicit none @@ -628,10 +643,10 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -642,6 +657,7 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -653,7 +669,7 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -663,4 +679,3 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) write(psb_err_unit,*) 'Error while opening ',filename return end subroutine lzmm_mat_write -