|
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 !!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
1.7.3-20110217