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