1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8 9module class_Sparsity 10 !! Implements a reference-counted derived type ([bud](|page|/datastructures/1-buds.html)) 11 !! to hold the [sparsity pattern](|page|/datastructures/2-sparse.html) of the program's matrices. 12 !! Objects of this type (notably [[sparse_matrices:sparse_pattern]] can be passed around 13 !! and put in container types. 14 !! The legacy indexing arrays can be linked to this bud's data ([[class_Sparsity:attach]]) 15 !! when necessary. 16 use alloc, only: re_alloc, de_alloc 17 18 implicit none 19 20 public :: newSparsity 21 public :: print_type 22 public :: attach 23 public :: nrows, nrows_g 24 public :: ncols, ncols_g 25 public :: n_row, n_col 26 public :: nnzs, list_ptr, list_col 27 public :: equivalent 28 29 character(len=*), parameter :: mod_name= "class_Sparsity" 30 31 ! This is the "meat" of the type 32 ! 33 type Sparsity_ 34 integer :: refCount = 0 35 character(len=36) :: id = "null_id" 36 !---------------------- 37 character(len=256) :: name = "null_sparsity" 38 integer :: nrows = 0 ! Local number of rows 39 integer :: nrows_g = 0 ! Global number or rows 40 integer :: ncols = 0 ! Local number of columns 41 integer :: ncols_g = 0 ! Global number or columns 42 integer :: nnzs = 0 ! Local number of non-zeros 43 integer, pointer :: n_col(:) =>null() ! Nonzero cols of each row 44 integer, pointer :: list_col(:) =>null() ! Index of nonzero columns 45 integer, pointer :: list_ptr(:) =>null() ! First element of each row 46 end type Sparsity_ 47 48 ! This is a wrapper type to be passed around 49 type Sparsity 50 type(Sparsity_), pointer :: data => null() 51 end type Sparsity 52 53 interface nrows 54 module procedure nrowsSparsity 55 end interface 56 57 interface nrows_g 58 module procedure nrows_gSparsity 59 end interface 60 61 interface ncols 62 module procedure ncolsSparsity 63 end interface 64 65 interface ncols_g 66 module procedure ncols_gSparsity 67 end interface 68 69! ****************** 70! Specific routines to retrieve direct information 71! about the sparsity entries. 72 73 interface attach 74 module procedure attachSparsity 75 end interface 76 77 interface equivalent 78 module procedure equivalentSparsity 79 end interface 80 81 interface n_col 82 module procedure n_colSparsity 83 module procedure n_colSparsityI 84 end interface 85 86 interface n_row 87 ! This interface ensures that one can retrieve 88 ! all information from the sparse class. 89 module procedure n_rowSparsityI 90 end interface 91 92 interface nnzs 93 module procedure nnzsSparsity 94 end interface 95 96 interface list_ptr 97 module procedure list_ptrSparsity 98 module procedure list_ptrSparsityI 99 end interface 100 interface list_col 101 module procedure list_colSparsity 102 module procedure list_colSparsityI 103 end interface 104 105 interface print_type 106 module procedure printSparsity 107 end interface 108 109!=========================== 110#define TYPE_NAME Sparsity 111#include "basic_type.inc" 112!=========================== 113 114 subroutine delete_Data(spdata) 115 type(Sparsity_) :: spdata 116 call de_alloc( spdata%n_col, & 117 name="n_col " // trim(spdata%name),routine="Sparsity") 118 call de_alloc( spdata%list_ptr, & 119 name="list_ptr " // trim(spdata%name),routine="Sparsity") 120 call de_alloc( spdata%list_col, & 121 name="list_col " // trim(spdata%name),routine="Sparsity") 122 end subroutine delete_Data 123 124!-------------------------------------------------------------------- 125! For easy compatability we have added ncols and ncols_g as 126! optional arguments. 127! In case, not specified, values of nrows_g are used. (block-cyclic) 128 subroutine newSparsity(sp,nrows,nrows_g,nnzs,num,listptr,list,name, & 129 ncols,ncols_g) 130 131 type(Sparsity), intent(inout) :: sp 132 133 integer, intent(in) :: nrows, nrows_g, nnzs 134 integer, intent(in) :: num(:), listptr(:) 135 integer, intent(in), optional :: list(:) 136 character(len=*), intent(in), optional :: name 137 integer, intent(in), optional :: ncols, ncols_g 138 139 ! We release the previous incarnation 140 ! This means that we relinquish access to the previous 141 ! memory location. It will be deallocated when nobody 142 ! else is using it. 143 144 call init(sp) 145 146 sp%data%name = trim(name) 147 148 call re_alloc(sp%data%n_col,1,nrows, & 149 name="n_col " // trim(sp%data%name),routine="Sparsity") 150 call re_alloc(sp%data%list_ptr,1,nrows, & 151 name="list_ptr " // trim(sp%data%name),routine="Sparsity") 152 153 sp%data%nrows = nrows 154 sp%data%nrows_g = nrows_g 155 if ( present(ncols_g) ) then 156 sp%data%ncols_g = ncols_g 157 else 158 ! We default it to a Block-cyclic distribution, hence 159 ! the number of columns is the same as the number of 160 ! global rows (we cannot guess super-cells, that would indeed be amazing) 161 sp%data%ncols_g = nrows_g 162 end if 163 if ( present(ncols) ) then 164 sp%data%ncols = ncols 165 else 166 ! Again, block cyclic has the maximum number of columns 167 ! in each block 168 sp%data%ncols = sp%data%ncols_g 169 end if 170 sp%data%nnzs = nnzs 171 sp%data%n_col(1:nrows) = num(1:nrows) 172 sp%data%list_ptr(1:nrows) = listptr(1:nrows) 173 174 if (nnzs /= sum(num(1:nrows))) then 175 call die("nnzs mismatch in new_sparsity") 176 endif 177 178 call re_alloc(sp%data%list_col,1,nnzs, & 179 name="list_col " // trim(sp%data%name),routine="Sparsity") 180 if ( present(list) ) then 181 sp%data%list_col(1:nnzs) = list(1:nnzs) 182 else 183 sp%data%list_col(1:nnzs) = 0 184 end if 185 186 call tag_new_object(sp) 187 188 end subroutine newSparsity 189 190! This works beautifully with gfortran, but intel creeps out 191! when destroying the object. :( 192! subroutine pntSparsity(sp,nrows,nrows_g,num,listptr,list,name, & 193! ncols,ncols_g) 194! 195! use alloc, only : alloc_count 196! 197! type(Sparsity), intent(inout) :: sp 198! 199! integer, intent(in) :: nrows, nrows_g 200! integer, pointer :: num(:), listptr(:), list(:) 201! character(len=*), intent(in) :: name 202! integer, intent(in), optional :: ncols, ncols_g 203! 204! ! We release the previous incarnation 205! ! This means that we relinquish access to the previous 206! ! memory location. It will be deallocated when nobody 207! ! else is using it. 208! 209! call init(sp) 210! 211! sp%data%name = trim(name) 212! 213! if ( size(num) /= nrows ) then 214! call die('pntSparsity: n_col size does not match nrows') 215! end if 216! if ( size(listptr) /= nrows ) then 217! call die('pntSparsity: list_ptr size does not match nrows') 218! end if 219! 220! ! Make pointers 221! sp%data%n_col => num(:) 222! sp%data%list_ptr => listptr(:) 223! sp%data%list_col => list(:) 224! 225! sp%data%nrows = nrows 226! sp%data%nrows_g = nrows_g 227! if ( present(ncols_g) ) then 228! sp%data%ncols_g = ncols_g 229! else 230! ! We default it to a Block-cyclic distribution, hence 231! ! the number of columns is the same as the number of 232! ! global rows (we cannot guess super-cells, that would indeed be amazing) 233! sp%data%ncols_g = nrows_g 234! end if 235! if ( present(ncols) ) then 236! sp%data%ncols = ncols 237! else 238! ! Again, block cyclic has the maximum number of columns 239! ! in each block 240! sp%data%ncols = sp%data%ncols_g 241! end if 242! sp%data%nnzs = sum(num(1:nrows)) 243! if ( size(list) /= sp%data%nnzs ) then 244! call die('pntSparsity: list size does not match sum(n_col)') 245! end if 246! 247! ! Track the sparsity count, this *might* produce wrong 248! ! count of memory! 249! ! TODO CHECK MEMORY COUNT 250! call alloc_count( nrows, 'I', 'n_col '// trim(sp%data%name), "Sparsity" ) 251! call alloc_count( nrows, 'I', 'list_ptr '// trim(sp%data%name), "Sparsity" ) 252! call alloc_count( sp%data%nnzs, 'I', 'list_col '// trim(sp%data%name), "Sparsity" ) 253! 254! call tag_new_object(sp) 255! 256! end subroutine pntSparsity 257 258 pure function nrowsSparsity(this) result (n) 259 type(Sparsity), intent(in) :: this 260 integer :: n 261 n = this%data%nrows 262 end function nrowsSparsity 263 264 elemental function n_rowSparsityI(this,col) result (p) 265 type(Sparsity), intent(in) :: this 266 integer, intent(in) :: col 267 integer :: p 268 p = count(this%data%list_col == col) 269 end function n_rowSparsityI 270 271 pure function nrows_gSparsity(this) result (n) 272 type(Sparsity), intent(in) :: this 273 integer :: n 274 n = this%data%nrows_g 275 end function nrows_gSparsity 276 277 pure function ncolsSparsity(this) result (n) 278 type(Sparsity), intent(in) :: this 279 integer :: n 280 n = this%data%ncols 281 end function ncolsSparsity 282 283 pure function ncols_gSparsity(this) result (n) 284 type(Sparsity), intent(in) :: this 285 integer :: n 286 n = this%data%ncols_g 287 end function ncols_gSparsity 288 289 290 pure function nnzsSparsity(this) result (n) 291 type(Sparsity), intent(in) :: this 292 integer :: n 293 if ( initialized(this) ) then 294 n = this%data%nnzs 295 else 296 n = 0 297 end if 298 end function nnzsSparsity 299 300 function n_colSparsity(this) result (p) 301 type(Sparsity), intent(in) :: this 302 integer, pointer :: p(:) 303 p => this%data%n_col 304 end function n_colSparsity 305 elemental function n_colSparsityI(this,row) result (p) 306 type(Sparsity), intent(in) :: this 307 integer, intent(in) :: row 308 integer :: p 309 p = this%data%n_col(row) 310 end function n_colSparsityI 311 312 function list_ptrSparsity(this) result (p) 313 type(Sparsity), intent(in) :: this 314 integer, pointer :: p(:) 315 p => this%data%list_ptr 316 end function list_ptrSparsity 317 elemental function list_ptrSparsityI(this,i) result (p) 318 type(Sparsity), intent(in) :: this 319 integer, intent(in) :: i 320 integer :: p 321 p = this%data%list_ptr(i) 322 end function list_ptrSparsityI 323 324 function list_colSparsity(this) result (p) 325 type(Sparsity), intent(in) :: this 326 integer, pointer :: p(:) 327 p => this%data%list_col 328 end function list_colSparsity 329 elemental function list_colSparsityI(this,i) result (p) 330 type(Sparsity), intent(in) :: this 331 integer, intent(in) :: i 332 integer :: p 333 p = this%data%list_col(i) 334 end function list_colSparsityI 335 336 ! Generic routine for retrieval of information in one line... 337 subroutine attachSparsity(this,D,n_col,list_col,list_ptr, & 338 nrows,nrows_g,ncols,ncols_g, & 339 nnzs) 340 type(Sparsity), intent(inout) :: this 341 logical, intent(in) , optional :: D ! DUMMY, force named arguments 342 ! Optional arguments to retrieve all information with named arguments 343 integer, pointer , optional :: n_col(:), list_col(:), list_ptr(:) 344 integer, intent(out), optional :: nrows, nrows_g 345 integer, intent(out), optional :: ncols, ncols_g 346 integer, intent(out), optional :: nnzs 347 if ( present(D) ) call die('PROGRAMMING ERROR, named args please') 348 349 ! the arrays 350 if ( present(n_col) ) n_col => n_colSparsity(this) 351 if ( present(list_col) ) list_col => list_colSparsity(this) 352 if ( present(list_ptr) ) list_ptr => list_ptrSparsity(this) 353 ! the integers 354 if ( present(nrows) ) nrows = nrowsSparsity(this) 355 if ( present(nrows_g) ) nrows_g = nrows_gSparsity(this) 356 if ( present(ncols) ) ncols = ncolsSparsity(this) 357 if ( present(ncols_g) ) ncols_g = ncols_gSparsity(this) 358 if ( present(nnzs) ) nnzs = nnzsSparsity(this) 359 360 end subroutine attachSparsity 361 362 function equivalentSparsity(sp1,sp2) result(equivalent) 363 type(Sparsity), intent(inout) :: sp1, sp2 364 logical :: equivalent 365 integer, pointer :: ncol1(:), l_col1(:), ncol2(:), l_col2(:) 366 integer :: lno1, lno2, no1, no2 367 368 equivalent = same(sp1,sp2) 369 if ( equivalent ) return 370 if ( initialized(sp1) ) then 371 if (.not. initialized(sp2) ) return 372 else if ( initialized(sp2) ) then 373 if (.not. initialized(sp1) ) return 374 end if 375 ! If they are the same object, immediately return 376 ! with true 377 equivalent = sp1%data%id == sp2%data%id 378 if ( equivalent ) return 379 380 ! Both are initialized 381 call attach(sp1,nrows=lno1,nrows_g=no1, & 382 n_col = ncol1, list_col = l_col1 ) 383 call attach(sp2,nrows=lno2,nrows_g=no2, & 384 n_col = ncol2, list_col = l_col2 ) 385 equivalent = lno1 == lno2 386 if ( .not. equivalent ) return 387 equivalent = no1 == no2 388 if ( .not. equivalent ) return 389 equivalent = all(ncol1 == ncol2) 390 if ( .not. equivalent ) return 391 equivalent = all(l_col1 == l_col2) 392 393 end function equivalentSparsity 394 395 subroutine printSparsity(sp) 396 type(Sparsity), intent(in) :: sp 397 398 if (.not. initialized(sp) ) then 399 print "(a)", "Sparsity Not Associated" 400 RETURN 401 endif 402 403 print "(a,/,t4,2(a,i0),a,f0.4,2(a,i0),a)", & 404 " <sparsity:"//trim(sp%data%name), & 405 " nrows_g=", sp%data%nrows_g, & 406 " nrows=", sp%data%nrows, & 407 " sparsity=",real(sp%data%nnzs)/& 408 real(sp%data%nrows_g)/real(sp%data%ncols_g), & 409 " nnzs=",sp%data%nnzs,", refcount: ", & 410 refcount(sp), ">" 411 end subroutine printSparsity 412 413end module class_Sparsity 414