DOUG 0.2

SpMtx_util.f90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 !!----------------------------------------------------------
00023 !!Some useful procedures for sparse matrixes
00024 !!----------------------------------------------------------
00025 Module SpMtx_util
00026   use DOUG_utils
00027   use RealKind
00028   use SpMtx_class
00029   use globals
00030   use DenseMtx_mod
00031 
00032   implicit none
00033 
00034 Contains
00037   subroutine SpMtx_findRowElems(A,rows,inds)
00038     type(SpMtx),intent(in) :: A
00039     integer,intent(in) :: rows(:) !< matrix row numbers
00040     integer,pointer :: inds(:) !< found matrix elements
00041 
00042     logical,dimension(:),pointer :: isin
00043     integer :: i,n,nz
00044 
00045     allocate(isin(A%nrows))
00046     isin=.false.
00047     n=size(rows)
00048     do i=1,n
00049       isin(rows(i))=.true.
00050     enddo
00051     ! count
00052     nz=0
00053     do i=1,A%nnz
00054       if (isin(A%indi(i))) then
00055         nz=nz+1
00056       endif
00057     enddo
00058     allocate(inds(nz))
00059     nz = 0
00060     do i=1,A%nnz
00061       if (isin(A%indi(i))) then
00062         nz=nz+1
00063         inds(nz)=i
00064       endif
00065     enddo
00066     deallocate(isin)
00067   end subroutine SpMtx_findRowElems
00068 
00071   subroutine SpMtx_findColumnElems(A,cols,inds)
00072     type(SpMtx),intent(in) :: A
00073     integer,intent(in) :: cols(:) !< matrix column numbers
00074     integer,pointer :: inds(:) !< found matrix elements
00075 
00076     logical,dimension(:),pointer :: isin
00077     integer :: i,n,nz
00078 
00079     allocate(isin(A%ncols))
00080     isin=.false.
00081     n=size(cols)
00082     do i=1,n
00083       isin(cols(i))=.true.
00084     enddo
00085     ! count
00086     nz=0
00087     do i=1,A%nnz
00088       if (isin(A%indj(i))) then
00089         nz=nz+1
00090       endif
00091     enddo
00092     allocate(inds(nz))
00093     nz = 0
00094     do i=1,A%nnz
00095       if (isin(A%indj(i))) then
00096         nz=nz+1
00097         inds(nz)=i
00098       endif
00099     enddo
00100     deallocate(isin)
00101   end subroutine SpMtx_findColumnElems
00102 
00105   subroutine SpMtx_findElems(A,x_inds,inds)
00106     type(SpMtx),intent(in) :: A
00107     integer,intent(in) :: x_inds(:) !< matrix row and column numbers
00108     integer,pointer :: inds(:) !< found matrix elements
00109 
00110     logical,dimension(:),pointer :: isin
00111     integer :: i,n,nz
00112 
00113     allocate(isin(max(A%nrows,A%ncols)))
00114     isin=.false.
00115     n=size(x_inds)
00116     do i=1,n
00117       isin(x_inds(i))=.true.
00118     enddo
00119     ! count
00120     nz=0
00121     do i=1,A%nnz
00122       if (isin(A%indi(i)).and.isin(A%indj(i))) then
00123         nz=nz+1
00124       endif
00125     enddo
00126     allocate(inds(nz))
00127     nz = 0
00128     do i=1,A%nnz
00129       if (isin(A%indi(i)).and.isin(A%indj(i))) then
00130         nz=nz+1
00131         inds(nz)=i
00132       endif
00133     enddo
00134     deallocate(isin)
00135   end subroutine SpMtx_findElems
00136 
00139   subroutine GetGivenRowsElements(A,nodes,indi,indj,val)
00140     Type(SpMtx),intent(in) :: A ! the fine level matrix
00141     integer,dimension(:),intent(in) :: nodes
00142     integer,dimension(:),pointer :: indi,indj,inds
00143     real(kind=rk),dimension(:),pointer :: val
00144     integer :: nz, ierr
00145 
00146     ! copy
00147     call SpMtx_findRowElems(A,nodes,inds)
00148     nz = size(inds)
00149     allocate(indi(nz),indj(nz),val(nz))
00150     indi = A%indi(inds)
00151     indj = A%indj(inds)
00152     val = A%val(inds)
00153     deallocate(inds)
00154   end subroutine GetGivenRowsElements
00155 
00158   subroutine GetGivenElements(A,nodes,indi,indj,val)
00159     Type(SpMtx),intent(in) :: A ! the fine level matrix
00160     integer,dimension(:),intent(in) :: nodes
00161     integer,dimension(:),pointer :: indi,indj,inds
00162     real(kind=rk),dimension(:),pointer :: val
00163     integer :: nz, ierr
00164 
00165     ! copy
00166     call SpMtx_findElems(A,nodes,inds)
00167     nz = size(inds)
00168     allocate(indi(nz),indj(nz),val(nz))
00169     indi = A%indi(inds)
00170     indj = A%indj(inds)
00171     val = A%val(inds)
00172     deallocate(inds)
00173   end subroutine GetGivenElements
00174 
00175   subroutine KeepGivenRowIndeces(A,inds,keepShape)
00176     implicit none
00177     Type(SpMtx),intent(inout) :: A ! the fine level matrix
00178     integer,dimension(:),intent(in) :: inds
00179     logical,intent(in) :: keepShape
00180     integer,dimension(:),pointer :: indi,indj
00181     real(kind=rk),dimension(:),pointer :: val
00182     logical,dimension(:),pointer :: isin
00183     integer :: i,n,nz
00184     allocate(isin(A%nrows))
00185     isin=.false.
00186     n=size(inds)
00187     do i=1,n
00188       isin(inds(i))=.true.
00189     enddo
00190     nz=0
00191     do i=1,A%nnz
00192       if (isin(A%indi(i))) then
00193         nz=nz+1
00194         A%indi(nz)=A%indi(i)
00195         A%indj(nz)=A%indj(i)
00196         A%val(nz)=A%val(i)
00197       endif
00198     enddo
00199     deallocate(isin)
00200     if (nz<A%nnz) then
00201       allocate(indi(nz))
00202       indi=A%indi(1:nz)
00203       deallocate(A%indi)
00204       allocate(A%indi(nz))
00205       A%indi=indi
00206       deallocate(indi)
00207       allocate(indj(nz))
00208       indj=A%indj(1:nz)
00209       deallocate(A%indj)
00210       allocate(A%indj(nz))
00211       A%indj=indj
00212       deallocate(indj)
00213       allocate(val(nz))
00214       val=A%val(1:nz)
00215       deallocate(A%val)
00216       allocate(A%val(nz))
00217       A%val=val
00218       deallocate(val)
00219       A%nnz=nz
00220       if (.not.keepShape) then
00221         A%nrows=maxval(A%indi)
00222       end if
00223     endif
00224 
00225     if (A%arrange_type/=D_SPMTX_ARRNG_NO) then
00226       deallocate(A%m_bound)
00227       A%arrange_type=D_SPMTX_ARRNG_NO
00228     end if
00229 
00230   end subroutine KeepGivenRowIndeces
00231 
00232   subroutine KeepGivenColumnIndeces(A,inds,keepShape)
00233     implicit none
00234     Type(SpMtx),intent(inout) :: A ! the fine level matrix
00235     integer,dimension(:),intent(in) :: inds
00236     logical,intent(in) :: keepShape
00237     integer,dimension(:),pointer :: indi,indj
00238     real(kind=rk),dimension(:),pointer :: val
00239     logical,dimension(:),pointer :: isin
00240     integer :: i,n,nz
00241     allocate(isin(A%ncols))
00242     isin=.false.
00243     n=size(inds)
00244     do i=1,n
00245       isin(inds(i))=.true.
00246     enddo
00247     nz=0
00248     do i=1,A%nnz
00249       if (isin(A%indj(i))) then
00250         nz=nz+1
00251         A%indi(nz)=A%indi(i)
00252         A%indj(nz)=A%indj(i)
00253         A%val(nz)=A%val(i)
00254       endif
00255     enddo
00256     deallocate(isin)
00257     if (nz<A%nnz) then
00258       allocate(indi(nz))
00259       indi=A%indi(1:nz)
00260       deallocate(A%indi)
00261       allocate(A%indi(nz))
00262       A%indi=indi
00263       deallocate(indi)
00264       allocate(indj(nz))
00265       indj=A%indj(1:nz)
00266       deallocate(A%indj)
00267       allocate(A%indj(nz))
00268       A%indj=indj
00269       deallocate(indj)
00270       allocate(val(nz))
00271       val=A%val(1:nz)
00272       deallocate(A%val)
00273       allocate(A%val(nz))
00274       A%val=val
00275       deallocate(val)
00276       A%nnz=nz
00277       if (.not.keepShape) then
00278         A%ncols=maxval(A%indj)
00279       end if
00280     endif
00281 
00282     if (A%arrange_type/=D_SPMTX_ARRNG_NO) then
00283       deallocate(A%m_bound)
00284       A%arrange_type=D_SPMTX_ARRNG_NO
00285     end if
00286   end subroutine KeepGivenColumnIndeces
00287 
00288     Function SpMtx_findElem(A, i, j) result(n)
00289         type(SpMtx), intent(in)  :: A
00290         integer, intent(in)         :: i, j
00291         integer                     :: n, k
00292         integer                     :: inds, inde
00293 
00294         if (A%arrange_type == D_SpMtx_ARRNG_ROWS) then
00295           if (i+1<=size(A%M_bound)) then
00296             inds = A%M_bound(i)
00297             inde = A%M_bound(i+1)-1
00298           else
00299             inds=1
00300             inde=0
00301           endif
00302         else if (A%arrange_type == D_SpMtx_ARRNG_COLS) then
00303           if (j+1<=size(A%M_bound)) then
00304             inds = A%M_bound(j)
00305             inde = A%M_bound(j+1)-1
00306           else
00307             inds=1
00308             inde=0
00309           endif
00310         else
00311           inds = 1
00312           inde = A%nnz
00313         endif
00314         n=0
00315         do k=inds,inde
00316           if ((A%indi(k)==i).AND.(A%indj(k)==j)) then
00317             n=k
00318             return
00319           end if
00320         end do
00321 
00322     End Function SpMtx_findElem
00323 
00324     Subroutine SpMtx_setVal(A, n, i, j, val)
00325         type(SpMtx), intent(inout)  :: A
00326         integer, intent(in)         :: n, i, j
00327         real(kind=rk), intent(in)   :: val
00328 
00329         A%indi(n) = i
00330         A%indj(n) = j
00331         A%val(n) = val
00332     End Subroutine
00333 
00334     Subroutine SpMtx_printMat(M)
00335 
00336       Implicit None
00337       type(SpMtx)      :: M
00338       integer          :: maxi, maxj, i, j, n
00339 
00340 !!$        maxi = maxval(M%indi)
00341 !!$        maxj = maxval(M%indj)
00342 
00343       maxi = M%nrows ! + ks
00344       maxj = M%ncols ! + ks
00345 
00346       do i = 1, maxi
00347          write (stream, "('[')", advance="no")
00348          do j = 1, maxj
00349             if (j /= 1) write (stream, "(',')", advance="no")
00350             n = SpMtx_findElem(M, i, j)
00351             if (n == 0) then
00352                write (stream, "(f6.3)", advance="no") 0.0
00353             else
00354                write (stream, "(f6.3)", advance="no") M%val(n)
00355             end if
00356          end do
00357 
00358          write (stream, "(']')", advance="yes")
00359 
00360       end do
00361 
00362     End Subroutine SpMtx_printMat
00363 
00364     subroutine SpMtx_printMat_in_arrays(n,indi,indj,val)
00365       Implicit None
00366       integer,dimension(:) :: indi,indj
00367       real(kind=rk),dimension(:) :: val
00368       integer          :: maxi, maxj, i, j, n, k
00369 
00370       do k=1,n
00371         write (stream,*)indi(k),indj(k),val(k)
00372       enddo
00373 return
00374       maxi = maxval(indi(:n))
00375       maxj = maxval(indj(:n))
00376       do i = 1, maxi
00377          write (stream, "('[')", advance="no")
00378          do j = 1, maxj
00379             if (j /= 1) write (stream, "(',')", advance="no")
00380        inn: do k=1,n
00381               if (indi(k)==i.and.indj(k)==j) exit inn
00382             enddo inn
00383             if (k>n) then
00384                write (stream, "(f6.3)", advance="no") 0.0
00385             else
00386                write (stream, "(f6.3)", advance="no") val(k)
00387             end if
00388          end do
00389          write (stream, "(']')", advance="yes")
00390       end do
00391     end Subroutine SpMtx_printMat_in_arrays
00392 
00393     subroutine SpMtx_printRaw(A,transp,startnz,endnz)
00394       implicit none
00395       type(SpMtx), intent(in)      :: A
00396       logical,optional             :: transp
00397       integer,optional             :: startnz,endnz
00398       integer                      :: i,i1,i2
00399       logical :: t
00400 
00401       if (present(transp).and.transp) then
00402         t=.true.
00403       else
00404         t=.false.
00405       endif
00406       if (present(startnz)) then
00407         i1=startnz
00408       else
00409         i1=1
00410       endif
00411       if (present(endnz)) then
00412         i2=endnz
00413       else
00414         i2=A%nnz
00415       endif
00416       write (stream,'(I0," x ",I0)') A%nrows, A%ncols
00417       write (stream,'(A5)',advance='no') "N"
00418       write (stream,'(A5)',advance='no') "indi"
00419       write (stream,'(A5)',advance='no') "indj"
00420       write (stream,'(A9)') "val"
00421       write (stream,*)"-----------------------------"
00422       do i=i1,i2
00423          write (stream,'(i5)',advance='no') i
00424          if (t) then
00425            write (stream,'(i5)',advance='no') A%indj(i)
00426            write (stream,'(i5)',advance='no') A%indi(i)
00427          else
00428            write (stream,'(i5)',advance='no') A%indi(i)
00429            write (stream,'(i5)',advance='no') A%indj(i)
00430          endif
00431          write (stream,'(e13.4)') A%val(i)
00432       end do
00433     end subroutine SpMtx_printRaw
00434 
00435 
00436     !-----------------------------
00437     ! Print info on sparse matirix
00438     !-----------------------------
00439     subroutine SpMtx_printInfo(S)
00440       implicit none
00441 
00442       type(SpMtx), intent(in) :: S
00443 
00444       write(stream,*)"-----------------------------"
00445       write(stream,*) 'Sparse matrix:'
00446       write(stream,*) '  nnz = ', S%nnz
00447       write(stream,*) 'nrows = ', S%nrows
00448       write(stream,*) 'ncols = ', S%ncols
00449       write(stream,'(a)',advance='no') ' arrange type: '
00450       if (S%arrange_type == D_SpMtx_ARRNG_NO) then
00451          write(stream,*) 'not arranged'
00452       else if (S%arrange_type == D_SpMtx_ARRNG_ROWS) then
00453          write(stream,*) 'arranged for rows'
00454       else if (S%arrange_type == D_SpMtx_ARRNG_COLS) then
00455          write(stream,*) 'arranged for columns'
00456       else
00457          call DOUG_abort('[SpMtx_printInfo] : Wrong matrix arrange type.',-1)
00458       end if
00459       write(stream,'(a)',advance='no') ' shape: '
00460       if (S%shape == D_SpMtx_SHAPE_UNDEF) then
00461          write(stream,*) 'undefined'
00462       else if (S%shape == D_SpMtx_SHAPE_SQUARE) then
00463          write(stream,*) 'square'
00464       else if (S%shape == D_SpMtx_SHAPE_MOREROWS) then
00465          write(stream,*) 'more rows'
00466       else if (S%shape == D_SpMtx_SHAPE_MORECOLS) then
00467          write(stream,*) 'more columns'
00468       end if
00469       write(stream,*) 'symmetric structure:   ',merge('Yes','No ', &
00470            S%symmstruct.eqv.(.true.))
00471       write(stream,*) 'symmetric numerically: ',merge('Yes','No ', &
00472            S%symmnumeric.eqv.(.true.))
00473       write(stream,*) 'number of blocks     = ', S%nblocks
00474       write(stream,*) 'inner freedoms bound = ', S%mtx_inner_bound
00475       write(stream,*) 'mtx_bbs = '
00476       call DenseMtx_print(S%mtx_bbs, noSize=1)
00477       write(stream,*) 'mtx_bbe = '
00478       call DenseMtx_print(S%mtx_bbe, noSize=1)
00479       write(stream,*)"-----------------------------"
00480 
00481     end subroutine SpMtx_printInfo
00482 
00483     !----------------------------------
00488         subroutine SpMtx_writeMatrix(A)
00489           implicit none
00490           
00491           type(SpMtx), intent(in) :: A !< matrix to be dumped to file
00492           
00493           logical   :: found
00494           integer   :: iounit, k, opened
00495           
00496           if (.not. ismaster()) &
00497             call DOUG_abort('[SpMtx_dumpMatrix] : Can only be called by master node.',-1)
00498           call FindFreeIOUnit(found, iounit)
00499           open(unit=iounit,iostat=opened,file=mctls%dump_matrix_file,status='replace', &
00500                                                         form='formatted',err=666)     !XXX TODO
00501           if (opened /= 0) &
00502                 call DOUG_abort('[SpMtx_dumpMatrix] : Could not open file.',-1)
00503           write(iounit, *) A%ncols, A%nnz
00504           do k = 1, A%nnz
00505             write(iounit, '(I10,I10,E24.16)') A%indi(k), A%indj(k), A%val(k)
00506           end do
00507           close(iounit)
00508           return
00509           
00510 666 call DOUG_abort('[SpMtx_dumpMatrix] : Error while writing to file.',-1)
00511         end subroutine SpMtx_writeMatrix
00512 
00513 
00515  subroutine SpMtx_writeLogicalValues(A, vals, fname)
00516    type(SpMtx) :: A !< Matrix
00517    logical :: vals(:) !< Values to write to file (does not have to be matrix values)
00518    character(*) :: fname !< Output file name.
00519 
00520    logical   :: found
00521    integer   :: iounit, k, opened
00522 
00523    call FindFreeIOUnit(found, iounit)
00524    open(unit=iounit,iostat=opened,file=fname,status='replace',err=666)
00525 
00526    write(iounit,*) max(A%nrows, A%ncols), size(vals)
00527    do k=1, size(vals)
00528       if(vals(k)) then
00529          write(iounit,*) A%indi(k), A%indj(k), 1
00530       else
00531          write(iounit,*) A%indi(k), A%indj(k), 0
00532       end if
00533    end do
00534 
00535    close(iounit)
00536    return
00537 
00538 666 call DOUG_abort('[SpMtx_dumpMatrix] : Error while writing to file.',-1)
00539  end subroutine SpMtx_writeLogicalValues
00540 
00541 End module SpMtx_util