DOUG 0.2

SpMtx_class.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 !!Sparse Matrix Class
00024 !!  Matrix constructor and destructor
00025 !!  Laplace Matrix Constructor
00026 !!----------------------------------------------------------
00027 module SpMtx_class
00028   use RealKind
00029   use globals
00030   use DOUG_utils
00031 
00032   implicit none
00033 
00034 #include<doug_config.h>
00035 
00036 #ifdef D_COMPLEX
00037 #define float complex
00038 #else
00039 #define float real
00040 #endif
00041 
00042   ! Sparse matrix arrangement modes:
00043   integer, parameter :: D_SpMtx_ARRNG_NO   = 0
00044   integer, parameter :: D_SpMtx_ARRNG_ROWS = 1
00045   integer, parameter :: D_SpMtx_ARRNG_COLS = 2
00046 
00047   ! Matrix shape:
00048   integer, parameter :: D_SpMtx_SHAPE_UNDEF    = -1
00049   integer, parameter :: D_SpMtx_SHAPE_SQUARE   =  1
00050   integer, parameter :: D_SpMtx_SHAPE_MOREROWS =  2
00051   integer, parameter :: D_SpMtx_SHAPE_MORECOLS =  3
00052 
00053   ! Matrix nonzero structure:
00054   logical :: D_SpMtx_STRUCTURE_SYMM = .true.
00055 
00056   ! Scaling of matrix values:
00057   integer, parameter :: D_SpMtx_SCALE_UNDEF         = -1
00058   integer, parameter :: D_SpMtx_SCALE_NO            =  0
00059   integer, parameter :: D_SpMtx_SCALE_DIAG          =  1
00060   integer, parameter :: D_SpMtx_SCALE_DIAG_FILTERED =  2
00061 
00062   ! Matrix file format
00063   integer, parameter :: D_SpMtx_FORMAT_TEXT   = 0
00064   integer, parameter :: D_SpMtx_FORMAT_BINARY = 1
00065   integer, parameter :: D_SpMtx_FORMAT_XDR    = 2
00066 
00067 
00069   type SpMtx
00071     integer                               :: nnz = -1
00072 
00073     integer                               :: nrows = -1, ncols = -1
00074 
00075     integer,        dimension(:), pointer :: indi, indj
00076 
00077     float(kind=rk), dimension(:), pointer :: val
00079     float(kind=rk), dimension(:), pointer :: val_intf_full
00081     float(kind=rk), dimension(:), pointer :: diag !< for scaled case
00082     
00083     logical,        dimension(:), pointer :: strong !< connections
00084     integer, dimension(:), pointer :: strong_rowstart,strong_colnrs !< !For strong connection reference:
00085     
00088 
00090     integer,        dimension(:), pointer :: M_bound
00091 
00096     integer                               :: arrange_type = -1
00097 
00098 
00100     integer :: shape = D_SpMtx_SHAPE_UNDEF
00101     logical :: symmstruct = .false.
00102 
00109     logical :: symmnumeric = .false.
00110 
00111     integer :: scaling = D_SpMtx_SCALE_UNDEF
00112 
00117 
00119     integer                          :: nblocks = -1
00120 
00121     integer                          :: mtx_inner_bound = -1
00122 
00125     integer, dimension(:,:), pointer :: mtx_bbs
00126 
00129     integer, dimension(:,:), pointer :: mtx_bbe
00130 
00131 
00138     integer                          :: ol0nnz = -1
00139 
00140 
00142     integer,   dimension(:), pointer :: perm_map
00143 
00145     integer :: subsolve_id
00146  end type SpMtx
00147 
00148 contains
00149 
00150 
00152   function SpMtx_New() result(M)
00153     implicit none
00154 
00155     type(SpMtx)    :: M !Sparse matrix (not allocated)
00156 
00157     M%indi => NULL()
00158     M%indj => NULL()
00159     M%val  => NULL()
00160     M%nnz = 0
00161 
00162     ! number of blocks
00163     M%nblocks = 0
00164 
00165     ! subblocks boundaries
00166     M%mtx_bbs => NULL()
00167     M%mtx_bbe => NULL()
00168     M%ol0nnz = -1
00169 
00170     ! rows & columns
00171     M%nrows = 0
00172     M%ncols = 0
00173 
00174     ! matrix shape
00175     M%shape = D_SpMtx_SHAPE_UNDEF
00176 
00177     ! symmetry of nonzero structure
00178     M%symmstruct = .false.
00179 
00180     ! numerical symmetry
00181     M%symmnumeric = .false.
00182 
00183     ! arrange type
00184     M%arrange_type = D_SpMtx_ARRNG_NO
00185     M%M_bound => NULL()
00186 
00187     ! permutation map
00188     M%perm_map => NULL()
00189 
00190     M%strong => NULL()
00191     M%strong_rowstart => NULL()
00192     M%strong_colnrs => NULL()
00193     M%diag => NULL()
00194 
00195     M%subsolve_id = -1
00196 
00197   end function SpMtx_New
00198 
00211   Function SpMtx_newInit(nnz,nblocks,nrows,ncols,symmstruct,symmnumeric,&
00212                       indi,indj,val,arrange_type,M_bound) result(M)
00213     Implicit None
00214     type(SpMtx)    :: M !Sparse matrix (not allocated)
00215     integer, intent(in):: nnz !Number of non-zero elements
00216     integer, intent(in), optional :: nblocks !Number of blocks
00217     integer, intent(in), optional :: nrows   !Number of rows
00218     integer, intent(in), optional :: ncols   !Number of columns
00219     logical, intent(in), optional :: symmstruct ! symmetry of nonzero structure
00220     logical, intent(in), optional :: symmnumeric! numerical symmetry
00221     integer, intent(in), dimension(:), optional :: indi
00222     integer, intent(in), dimension(:), optional :: indj
00223     float(kind=rk), intent(in), dimension(:), optional :: val
00224     integer, intent(in), optional               :: arrange_type
00225     integer, intent(in), dimension(:), optional :: M_bound
00226     !- - - - - - - - - - - - -
00227     M = SpMtx_New()
00228 
00229     allocate(M%indi(nnz), M%indj(nnz), M%val(nnz))
00230     M%nnz=nnz
00231 
00232     ! number of blocks
00233     if (present(nblocks)) then
00234        M%nblocks = nblocks
00235     else
00236        M%nblocks = 1
00237     end if
00238 
00239     ! subblocks boundaries
00240     allocate(M%mtx_bbs(2*M%nblocks,2*M%nblocks))
00241     allocate(M%mtx_bbe(2*M%nblocks,2*M%nblocks))
00242     ! default values in the case of unparallel case here:
00243     if (M%nblocks==1) then
00244       M%mtx_bbs = 1
00245       M%mtx_bbe = 0
00246       M%mtx_bbe(2,2) = nnz
00247     else
00248       M%mtx_bbs = -1
00249       M%mtx_bbe = -1
00250     endif
00251 
00252     ! rows & columns
00253     M%nrows=0
00254     M%ncols=0
00255     if (present(nrows)) then
00256        M%nrows = nrows
00257        if (.not.present(ncols)) M%ncols = nrows ! shape-symmetric matrix
00258     end if
00259     if (present(ncols)) M%ncols = ncols
00260 
00261     ! matrix shape
00262     if (M%nrows == M%ncols) then
00263        M%shape = D_SpMtx_SHAPE_SQUARE
00264     else if (M%nrows > M%ncols) then
00265        M%shape = D_SpMtx_SHAPE_MOREROWS
00266     else if (M%nrows < M%ncols) then
00267        M%shape = D_SpMtx_SHAPE_MORECOLS
00268     end if
00269 
00270     ! symmetry of nonzero structure
00271     M%symmstruct = .false.
00272     if (present(symmstruct).and.&
00273          M%shape == D_SpMtx_SHAPE_SQUARE) then
00274        M%symmstruct = symmstruct
00275     end if
00276 
00277     ! numerical symmetry
00278     M%symmnumeric = .false.
00279     if (present(symmnumeric).and.&
00280          (M%shape == D_SpMtx_SHAPE_SQUARE).and.&
00281          (M%symmstruct.eqv.(.true.))) then
00282        M%symmnumeric = symmnumeric
00283     end if
00284 
00285     if (present(indi)) then
00286       M%indi=indi(1:nnz)
00287     else
00288       M%indi=0
00289     endif
00290     if (present(indj)) then
00291       M%indj=indj(1:nnz)
00292     else
00293       M%indj=0
00294     endif
00295     if (present(val)) then
00296       M%val=val(1:nnz)
00297     else
00298       M%val=0.0_rk
00299     endif
00300 
00301     ! arrange type
00302     if (present(arrange_type)) then
00303       if (arrange_type==D_SpMtx_ARRNG_NO.or.&
00304           arrange_type==D_SpMtx_ARRNG_ROWS.or.&
00305           arrange_type==D_SpMtx_ARRNG_COLS) then
00306         M%arrange_type=arrange_type
00307       else
00308         print *,'WARNING: SpMtx_newInit arrange_type ',arrange_type,&
00309                 ' not defined!'
00310       endif
00311     else
00312       M%arrange_type=D_SpMtx_ARRNG_NO
00313     endif
00314 
00315     if (present(M_bound)) then
00316       allocate(M%M_bound(size(M_bound)))
00317       M%M_bound=M_bound
00318     endif
00319 
00320     ! scaling
00321     M%scaling=D_SpMtx_SCALE_NO
00322 
00323     ! permutation map
00324     M%perm_map => NULL()
00325   End Function SpMtx_newInit
00326 
00327 !!$  !--------------------------------------------
00328 !!$  ! Sets all blocks bounds
00329 !!$  !--------------------------------------------
00330 !!$  subroutine SpMtx_setBlocksBounds(A, bbs, bbe)
00331 !!$    implicit none
00332 !!$    type(SpMtx),             intent(in out) :: A   ! Sparse matrix
00333 !!$    integer, dimension(:,:), intent(in)     :: bbs ! Beginings of the blocks
00334 !!$    integer, dimension(:,:), intent(in)     :: bbe ! Ends of the blocks
00335 !!$
00336 !!$!    call SpMtx_setMtxInnerBound(A, nnz_intf+1)
00337 !!$
00338 !!$  end subroutine SpMtx_setBlocksBounds
00339 
00340 
00360   subroutine SpMtx_setMtxInnerBound(A, bound)
00361     implicit None
00362     type(SpMtx), intent(in out) :: A
00363     integer,     intent(in)     :: bound
00364 
00365     A%mtx_inner_bound = bound
00366   end subroutine SpMtx_setMtxInnerBound
00367 
00369   subroutine SpMtx_Resize(A, N)
00370     Implicit None
00371     type(SpMtx), intent(inout) :: A
00372     type(SpMtx)     :: M
00373     integer, intent(in):: N !Number of non-zero elements
00374     integer            :: minn
00375     !- - - - - - - - - - - - -
00376 
00377 !    reshape(A%indi(1:N))
00378 !    reshape(A%indj(1:N))
00379 !    reshape(A%val(1:N))
00380 
00381     minn = min(A%nnz, N)
00382 
00383     allocate(M%indi(N),M%indj(N))
00384     allocate(M%val(N))
00385     M%indi(1:minn) = A%indi(1:minn)
00386     M%indj(1:minn) = A%indj(1:minn)
00387     M%val(1:minn) = A%val(1:minn)
00388 
00389     if(associated(A%indi)) deallocate(A%indi)
00390     if(associated(A%indj)) deallocate(A%indj)
00391     if(associated(A%val)) deallocate(A%val)
00392 
00393     A%indi => M%indi
00394     A%indj => M%indj
00395     A%val => M%val
00396     A%nnz = N
00397 
00398     if (A%mtx_bbe(2,2)>N) A%mtx_bbe(2,2)=N
00399 
00400   End subroutine SpMtx_Resize
00401 
00410   subroutine ReadInSparseAssembled(A,filename)
00411     implicit none
00412     type(SpMtx),intent(in out) :: A   ! System matrix
00413     character*(*),intent(in) :: filename
00414     integer  :: fHandler  !< file Handler
00415     integer :: n,nnz,nsd,i,j
00416     float :: val
00417 
00418         call ReadInSparseAssembledHeader(filename, fHandler, n, nnz)
00419 
00420     nsd=2
00421     if (n>nnz) then
00422       i=nnz
00423       nnz=n
00424       n=i
00425     endif
00426     write(stream,*) 'n,nnz:',n,nnz
00427     A = SpMtx_newInit(nnz,nblocks=sctls%number_of_blocks, &
00428                          nrows=n,                      &
00429                          ncols=n,                      &
00430                     symmstruct=sctls%symmstruct,       &
00431                    symmnumeric=sctls%symmnumeric       &
00432                   )
00433     
00434     call ReadInSparseAssembledBulk(fHandler, A)
00435     
00436     do i=1,nnz
00437       if (i<3.or.i>nnz-2) then
00438         write(stream,'(I4,A,I2,I2,E24.16)') i,' mat:',A%indi(i),A%indj(i),A%val(i)
00439       endif
00440     enddo
00441     if (A%indi(1)==0) then
00442       A%indi=A%indi+1
00443       A%indj=A%indj+1
00444     endif
00445     call CloseSparseAssembledFile(fHandler)
00446     
00447     ! this is undistributed matrix:
00448     A%mtx_bbs(1,1)=1
00449     A%mtx_bbe(1,1)=0
00450     A%mtx_bbs(2,2)=1
00451     A%mtx_bbe(2,2)=A%nnz
00452     return
00453   end subroutine ReadInSparseAssembled
00454 
00455   !----------------------------------------------------------
00458   subroutine ReadInSparseAssembledHeader(filename, fHandler, n, nnz)
00459     implicit none
00460     
00461     character*(*),intent(in)  :: filename
00462     integer      ,intent(out) :: fHandler
00463     integer      ,intent(out) :: n
00464     integer      ,intent(out) :: nnz
00465     
00466     if (mctls%assembled_mtx_format == D_SpMtx_FORMAT_TEXT) then
00467       call ReadInSparseAssembledHeader_TextBinary(filename, fHandler, n, nnz, .false.)
00468     elseif (mctls%assembled_mtx_format == D_SpMtx_FORMAT_BINARY) then
00469       call ReadInSparseAssembledHeader_TextBinary(filename, fHandler, n, nnz, .true.)
00470 #ifdef HAVE_LIBFXDR
00471     elseif (mctls%assembled_mtx_format == D_SpMtx_FORMAT_XDR) then
00472       call ReadInSparseAssembledHeader_XDR(filename, fHandler, n, nnz)
00473 #endif
00474     else
00475       call DOUG_abort('[ReadInSparseAssembledHeader] Data format not recognized. 0=text, 2=XDR (if compiled in)', -1)
00476 
00477     endif
00478   end subroutine ReadInSparseAssembledHeader
00479   
00480   !----------------------------------------------------------
00483   subroutine ReadInSparseAssembledHeader_TextBinary(filename, fHandler, n, nnz, binary)
00484     implicit none
00485     
00486     character*(*),intent(in)  :: filename
00487     integer      ,intent(out) :: fHandler
00488     integer      ,intent(out) :: n
00489     integer      ,intent(out) :: nnz
00490     logical      ,intent(in)  :: binary
00491     
00492     logical :: found
00493 
00494     call FindFreeIOUnit(found, fHandler)
00495     if (.NOT. found) &
00496       call DOUG_abort('[ReadInSparseAssembledHeader_TextBinary] No free IO unit.', -1)
00497       
00498     if (binary) then
00499       open(fHandler,FILE=trim(filename),STATUS='OLD',FORM='UNFORMATTED', &
00500              ERR=444)
00501       read(fHandler,END=500) n,nnz
00502     else
00503       open(fHandler,FILE=trim(filename),STATUS='OLD',FORM='FORMATTED', &
00504              ERR=444)
00505       read(fHandler,FMT=*,END=500) n,nnz
00506     endif
00507     return
00508         
00509 444 call DOUG_abort('Unable to open assembled file: '//trim(filename)//' ', -1)
00510 500 continue ! End of file reached too soon
00511     call DOUG_abort('File '//trim(filename)//' too short! ', -1)
00512   end subroutine ReadInSparseAssembledHeader_TextBinary
00513 
00514 #ifdef HAVE_LIBFXDR  
00515   !----------------------------------------------------------
00518   subroutine ReadInSparseAssembledHeader_XDR(filename, fHandler, n, nnz)
00519     implicit none
00520     include 'fxdr.inc'
00521     
00522     character*(*),intent(in)  :: filename
00523     integer      ,intent(out) :: fHandler
00524     integer      ,intent(out) :: n
00525     integer      ,intent(out) :: nnz
00526     
00527     integer :: ierr
00528     
00529     fHandler = initxdr( trim(filename), 'r', .FALSE. )
00530     ierr  = ixdrint( fHandler, n )
00531     ierr  = ixdrint( fHandler, nnz )
00532 
00533   end subroutine ReadInSparseAssembledHeader_XDR
00534 
00535 #endif
00536 
00537   !----------------------------------------------------------
00540   subroutine ReadInSparseAssembledBulk(fHandler, A)
00541     implicit none
00542     
00543     integer        ,intent(in) :: fHandler
00544     type(SpMtx), intent(inout) :: A
00545     
00546     if (mctls%assembled_mtx_format == D_SpMtx_FORMAT_TEXT) then
00547       call ReadInSparseAssembledBulk_TextBinary(fHandler, A, .false.)
00548     elseif (mctls%assembled_mtx_format == D_SpMtx_FORMAT_BINARY) then
00549       call ReadInSparseAssembledBulk_TextBinary(fHandler, A, .true.)
00550 #ifdef HAVE_LIBFXDR
00551     elseif (mctls%assembled_mtx_format == D_SpMtx_FORMAT_XDR) then
00552       call ReadInSparseAssembledBulk_XDR(fHandler, A)
00553 #endif
00554     else
00555       call DOUG_abort('[ReadInSparseAssembledBulk] Data format not recognized. 0=text, 2=XDR (if compiled in)', -1)
00556     endif
00557     
00558   end subroutine ReadInSparseAssembledBulk
00559 
00560   !----------------------------------------------------------
00563   subroutine ReadInSparseAssembledBulk_TextBinary(fHandler, A, binary)
00564     implicit none
00565 
00566     type SpMtx_Triple
00567       integer :: i, j
00568       float(kind=rk) :: val
00569     end type SpMtx_Triple    
00570 
00571     integer        ,intent(in) :: fHandler
00572     type(SpMtx), intent(inout) :: A
00573     logical        ,intent(in) :: binary
00574 
00575     type(SpMtx_Triple), dimension(:), pointer :: triples
00576     integer :: i
00577     
00578     if (binary) then
00579       allocate(triples(A%nnz))
00580       read(fHandler, END=500) (triples(i), i=1,A%nnz)
00581       do i=1,A%nnz
00582         A%indi(i)=triples(i)%i
00583         A%indj(i)=triples(i)%j
00584         A%val(i)=triples(i)%val
00585       enddo
00586       deallocate(triples)
00587     else
00588       do i=1,A%nnz
00589         read(fHandler, FMT=*,END=500 ) A%indi(i),A%indj(i),A%val(i)
00590       enddo
00591     endif
00592     return
00593     
00594 500 continue ! End of file reached too soon
00595     call DOUG_abort('[ReadInSparseAssembledBulk_TextBinary] File too short! ', -1)
00596   
00597   end subroutine ReadInSparseAssembledBulk_TextBinary
00598 
00599 #ifdef HAVE_LIBFXDR
00600   !----------------------------------------------------------
00603   subroutine ReadInSparseAssembledBulk_XDR (fHandler, A)
00604     implicit none
00605     include 'fxdr.inc'
00606     
00607     integer        ,intent(in) :: fHandler
00608     type(SpMtx), intent(inout) :: A
00609     integer :: ierr, i
00610     
00611     do i=1,A%nnz
00612       ierr = ixdrint( fHandler, A%indi(i) )
00613           ierr = ixdrint( fHandler, A%indj(i) )
00614           ierr = ixdrdouble( fHandler, A%val(i) )
00615     enddo
00616     
00617   end subroutine ReadInSparseAssembledBulk_XDR
00618 
00619 #endif
00620   
00621   !----------------------------------------------------------
00624   subroutine CloseSparseAssembledFile(fHandler)
00625     implicit none
00626 #ifdef HAVE_LIBFXDR
00627     include 'fxdr.inc'
00628 #endif
00629     
00630     integer        ,intent(in) :: fHandler
00631     integer :: ierr
00632     
00633     if (mctls%assembled_mtx_format == D_SpMtx_FORMAT_TEXT) then
00634       close(fHandler)
00635     elseif (mctls%assembled_mtx_format == D_SpMtx_FORMAT_BINARY) then
00636       close(fHandler)
00637 #ifdef HAVE_LIBFXDR
00638     elseif (mctls%assembled_mtx_format == D_SpMtx_FORMAT_XDR) then
00639       ierr = ixdrclose(fHandler)
00640 #endif
00641     else
00642       call DOUG_abort('[CloseSparseAssembledFile] Data format not recognized. 0=text, 2=XDR (if compiled in)', -1)
00643     endif
00644     
00645   end subroutine CloseSparseAssembledFile
00646 
00654   Subroutine SpMtx_Destroy(M) ! SpMtx_Destroy
00655     Implicit None
00656     type(SpMtx), intent(in out):: M !Sparse matrix (in)
00657     !- - - - - - - - - - - - - - - - - - - - -
00658     M%nnz=0
00659     deallocate(M%indi,M%indj)
00660     deallocate(M%val)
00661     if (associated(M%perm_map)) deallocate(M%perm_map)
00662     if (associated(M%mtx_bbs))  deallocate(M%mtx_bbs)
00663     if (associated(M%mtx_bbe))  deallocate(M%mtx_bbe)
00664     if (associated(M%strong))   deallocate(M%strong)
00665     if (associated(M%strong_rowstart)) &
00666                                 deallocate(M%strong_rowstart)
00667     if (associated(M%strong_colnrs)) &
00668                                 deallocate(M%strong_colnrs)
00669     if (associated(M%diag))     deallocate(M%diag)
00670     !------
00671     !write(stream,*)'[SpMtx_Destroy]: PAY ATTENTION ON IT!'
00672     if (M%arrange_type/=D_SpMtx_ARRNG_NO) deallocate(M%M_bound)
00673     !------
00674     if (associated(M%M_bound)) deallocate(M%M_bound)
00675     M%arrange_type=0
00676 
00677   End Subroutine SpMtx_Destroy
00678 
00680   Function SpMtx_Copy(IM) result(FM)
00681     Implicit None
00682     Type(SpMtx), intent(in):: IM !Initial Sparse matrix(in)
00683     Type(SpMtx)            :: FM !copy of sparse matrix(out)
00684     !- - - - - - - - - - - - - - - - -
00685     FM=SpMtx_newInit(max(IM%nnz,IM%ol0nnz))
00686     FM%nnz=IM%nnz
00687     FM%ol0nnz=IM%ol0nnz
00688     FM%nrows=IM%nrows; FM%ncols=IM%ncols
00689     FM%Arrange_Type=IM%Arrange_Type
00690     if (FM%Arrange_Type /=0) then
00691       allocate(FM%M_bound(size(IM%M_bound)))
00692       FM%M_bound=IM%M_bound
00693     end if
00694     FM%indi=IM%indi; FM%indj=IM%indj
00695     FM%val=IM%val
00696     if (associated(IM%strong)) then
00697       allocate(FM%strong(size(IM%strong)))
00698     endif
00699   End Function SpMtx_Copy
00709   Function LaplSpMtx_New(N, special) result(M)
00710     Implicit None
00711     type(SpMtx)           :: M       !Sparse matrix (out:Laplace)
00712     integer, intent(in)   :: N       !Dimension of Laplace Matrix: N^2 x N^2
00713     !for elements:
00714     !.TRUE. : each elements are different from other
00715     !.FALSE.: classical Laplace Matrix
00716     logical, optional, intent(in):: special
00717     integer                      :: el,i,k  !counters
00718     integer                      :: Nsize   !number of non-zero elements
00719     real(kind=rk)                :: eps     !special: difference
00720     !- - - - - - - - - - - - - - - - - - - - -
00721     Nsize=N*(N+2*(N-1))+2*(N-1)*N              !non-zero elements
00722     M=SpMtx_newInit(Nsize)
00723     el=0                           !count non-zero elements
00724     if (present(special)) then     !initial eps value
00725       if (special) eps=0.001_rk
00726                           else
00727       eps=0._rk
00728     end if
00729     !!!Blocks in main diagonal
00730     !elements in main diagonal
00731     do i=1,N**2
00732       el=el+1
00733       M%indi(el)=i
00734       M%indj(el)=i
00735       M%val(el)=4.+el*eps
00736     end do
00737     !other non-zero elements in main diagonal blocks
00738     do k=1, N      !for blocks
00739       el=el+1    !first row
00740       M%indi(el)=(k-1)*N+1
00741       M%indj(el)=(k-1)*N+2
00742       M%val(el)=-1.-el*eps
00743       el=el+1    !last row
00744       M%indi(el)=k*N
00745       M%indj(el)=k*N-1
00746       M%val(el)=-1.-el*eps
00747       do i=(k-1)*N+2,k*N-1 !for rows; index: 2..N-1
00748         el=el+1
00749         M%indi(el)=i
00750         M%indj(el)=i-1
00751         M%val(el)=-1.-el*eps
00752         el=el+1
00753         M%indi(el)=i
00754         M%indj(el)=i+1
00755         M%val(el)=-1.-el*eps
00756       end do
00757     end do
00758     !!!Other blocks in Laplace matrix
00759     do i=N+1,N**2-N
00760       el=el+1
00761       M%indi(el)=i
00762       M%indj(el)=i-N
00763       M%val(el)=-1.-el*eps
00764       el=el+1
00765       M%indi(el)=i
00766       M%indj(el)=i+N
00767       M%val(el)=-1.-el*eps
00768     end do
00769     !first block
00770     do i=1, N
00771       el=el+1
00772       M%indi(el)=i
00773       M%indj(el)=i+N
00774       M%val(el)=-1.-el*eps
00775     end do
00776     !last block
00777     do i=N**2-N+1,N**2
00778       el=el+1
00779       M%indi(el)=i
00780       M%indj(el)=i-N
00781       M%val(el)=-1.-el*eps
00782     end do
00783     !Control
00784     if (el /= Nsize) print*, "ERROR: Laplace Matrix Constructor"
00785     M%nrows = N**2;
00786     M%ncols = N**2;
00787   End Function LaplSpMtx_New
00797   Function order_laplace_m(N) result(M)
00798     Implicit None
00799     type(SpMtx)           :: M     !Sparse matrix
00800     integer,intent(in)    :: N     !Laplace parameter
00801     integer               :: Nsize !number of non-zero elements
00802     integer               :: el, x, y, ind, k, nx, ny, nind
00803     integer,dimension(1:4):: neighx = (/0, 1, 0, -1/), neighy = (/1, 0, -1, 0/)
00804     Nsize=N*(N+2*(N-1))+2*(N-1)*N
00805     M=SpMtx_newInit(Nsize)
00806     el = 0
00807     do x = 1, N
00808       do y = 1, N
00809         ind = N*(x - 1) + y
00810         el = el + 1
00811         M%indi(el) = ind
00812         M%indj(el) = ind
00813         M%val(el) = 4.0_rk
00814         do k = 1, 4
00815           nx = x + neighx(k)
00816           ny = y + neighy(k)
00817           if ((nx >= 1) .and. (nx <= N) .and. (ny >= 1) .and. (ny <= n)) then
00818             nind = ind + neighx(k)*N + neighy(k)
00819             el = el + 1
00820             M%indi(el) = ind
00821             M%indj(el) = nind
00822             M%val(el) = -1
00823           end if
00824         end do
00825       end do
00826     end do
00827     M%nrows = N**2;
00828     M%ncols = N**2;
00829   end function order_laplace_m
00830 
00831 End Module SpMtx_class
00832 !----------------------------------------------------------------------
00833