|
DOUG 0.2
|
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
1.7.3-20110217