Line data Source code
1 : !!****m* ABINIT/m_spmat_ndcoo
2 : !! NAME
3 : !! m_spmat_ndcoo
4 : !!
5 : !! FUNCTION
6 : !! This module contains the a NDCOO (n-dimensional coordinate) format of sparse matrix.
7 : !! Datatypes:
8 : !! NDCOO_mat_t: ND COO matrix
9 : !!
10 : !! Subroutines:
11 : !! TODO: add this when F2003 doc style is determined.
12 : !!
13 : !!
14 : !! COPYRIGHT
15 : !! Copyright (C) 2001-2026 ABINIT group (hexu)
16 : !! This file is distributed under the terms of the
17 : !! GNU General Public License, see ~abinit/COPYING
18 : !! or http://www.gnu.org/copyleft/gpl.txt .
19 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
20 : !!
21 : !! SOURCE
22 :
23 : #if defined HAVE_CONFIG_H
24 : #include "config.h"
25 : #endif
26 : #include "abi_common.h"
27 :
28 : module m_spmat_NDCOO
29 : use defs_basis
30 : use m_abicore
31 : use m_errors
32 : use m_xmpi
33 : use m_spmat_base
34 : use m_dynamic_array, only: int2d_array_type, real_array_type, int_array_type
35 : implicit none
36 : !!***
37 : private
38 :
39 : !-----------------------------------------------------------------------
40 : !> @brief NDCOO_mat_t : N-dimensional COO matrix type
41 : ! The matrix is stored in two arrays, one for the indices, and
42 : ! the other for the values.
43 : ! e.g. a 3D matrix (M) with nnz non-zero entries will have an
44 : ! (3,nnz) array as indices, and a (nnz) array value.
45 : ! Both the index array and the value array are dynamic, so entries
46 : ! can be appended.
47 : ! let M be a 4*5*6 matrix,
48 : ! M(1,2,3)=1.9, M(1,3,4)=2.4, M(...) =0
49 : ! We'll have
50 : ! ndim=3, nnz=2, mshape=[4,5,6]
51 : ! ind = [ 1 2 3,
52 : ! 1 3 4 ]
53 : ! val = [1.9, 2.4]
54 : ! (ind and val are dynamic array.)
55 : !-----------------------------------------------------------------------
56 : type, public :: ndcoo_mat_t
57 : integer :: ndim=0 ! number of dimensions
58 : integer :: nnz=0 ! Number of (None-zero) entries.
59 : integer, allocatable :: mshape(:) ! the shape of the matrix. len(mshape)=ndim.
60 : ! Note that it is not checked if the index by the shaped. If the shape for some
61 : ! dimension is unkown, it can be set to -1.
62 : type(int2d_array_type) :: ind ! The index array
63 : type(real_array_type) :: val ! The value array
64 : logical :: is_sorted = .False. ! If the matrix is sorted by index
65 : logical :: is_unique = .False. ! If the matrix is made unique (no entry could have same index).
66 : logical :: is_pair_grouped = .False. ! If the matrix entries are grouped by first two indices
67 : type(int_array_type) :: pair_1list ! first index of the pair
68 : type(int_array_type) :: pair_2list ! second index of the pair
69 : type(int_array_type) :: pair_startend ! start and end matrix index of each group
70 :
71 : contains
72 : procedure :: initialize
73 : procedure :: finalize
74 : procedure :: add_entry ! add one entry
75 : procedure :: remove_zeros ! remove entries which are (or close to ) zero
76 : procedure :: sort_indices ! sort the matrix by indices
77 : procedure :: sum_duplicates ! remove duplicate indices by adding them up
78 : procedure :: get_val_inz ! get the z'th value.
79 : procedure :: get_ind_inz ! get the z'th indices.
80 : procedure :: get_ind ! get the indices for all in a dimension
81 : procedure :: group_by_1dim ! group the matrix by first dimension
82 : procedure :: group_by_pair ! group the matrix by first two dimensions
83 : procedure :: mv1vec ! multiply vector, return NDCOO entity with one dimension less
84 : procedure :: mv2vec ! multiply 2 vectors, return NDCOO entity with two dimensions less
85 : procedure :: vec_product2d ! multiply 2d matrix with one vector, return vector
86 : procedure :: vec_product ! multiply 3d matrix with two vectors, return vector
87 : procedure :: vec_product4d ! multiply 4d matrix with three vectors
88 : !procedure :: print
89 : end type ndcoo_mat_t
90 :
91 : public:: test_ndcoo
92 : contains
93 :
94 : !-------------------------------------------------------------------!
95 : ! ndcoo_mat_t initializer:
96 : ! Input:
97 : ! mshape: the shape of the N-dimension matrix. array(ndim)
98 : !-------------------------------------------------------------------!
99 27 : subroutine initialize(self, mshape)
100 : class(ndcoo_mat_t), intent(inout) :: self
101 : integer, intent(in) :: mshape(:)
102 27 : self%ndim=size(mshape)
103 81 : ABI_MALLOC(self%mshape, (self%ndim))
104 141 : self%mshape=mshape
105 27 : self%nnz=0
106 27 : self%is_sorted=.False.
107 27 : self%is_unique=.False.
108 27 : self%is_pair_grouped = .False.
109 :
110 27 : end subroutine initialize
111 :
112 : !-------------------------------------------------------------------!
113 : ! Finalizer of ndcoo_mat_t
114 : !-------------------------------------------------------------------!
115 38 : subroutine finalize(self)
116 : class(ndcoo_mat_t), intent(inout) :: self
117 38 : self%ndim=0
118 38 : self%nnz=0
119 38 : self%is_sorted=.False.
120 38 : self%is_unique=.False.
121 38 : if((self%is_pair_grouped)) then
122 2 : call self%pair_1list%finalize()
123 2 : call self%pair_2list%finalize()
124 2 : call self%pair_startend%finalize()
125 : endif
126 38 : self%is_pair_grouped = .False.
127 :
128 38 : if (allocated(self%mshape)) then
129 27 : ABI_FREE(self%mshape)
130 : endif
131 38 : call self%ind%finalize()
132 38 : call self%val%finalize()
133 38 : end subroutine finalize
134 :
135 : !-------------------------------------------------------------------!
136 : ! Add one entry to the ndcoo_mat_t
137 : ! Inputs:
138 : ! ind: indices of the matrix.
139 : ! val: value of matrix.
140 : ! Example:
141 : ! call m%add_entry([1,2,3], 0.5)
142 : !-------------------------------------------------------------------!
143 2553449 : subroutine add_entry(self, ind, val)
144 : class(ndcoo_mat_t), intent(inout) :: self
145 : integer, intent(in) :: ind(self%ndim)
146 : real(dp), intent(in) :: val
147 2553449 : self%nnz=self%nnz+1
148 2553449 : call self%ind%push(ind)
149 2553449 : call self%val%push(val)
150 2553449 : self%is_sorted=.False.
151 2553449 : self%is_unique=.False.
152 2553449 : self%is_pair_grouped=.False.
153 2553449 : end subroutine add_entry
154 :
155 :
156 : !-------------------------------------------------------------------!
157 : ! sort the entries by indices. (left to right)
158 : !-------------------------------------------------------------------!
159 12 : subroutine sort_indices(self)
160 : class(ndcoo_mat_t), intent(inout) :: self
161 24 : real(dp) :: tmp(self%nnz)
162 0 : integer :: reorder(self%nnz)
163 12 : if(self%is_sorted .or. self%nnz==0) return
164 12 : call self%ind%sort(order=reorder)
165 1278504 : tmp(:)=self%val%data(1:self%nnz)
166 1278504 : self%val%data(1:self%nnz)=tmp(reorder)
167 12 : self%is_sorted=.True.
168 12 : end subroutine sort_indices
169 :
170 :
171 : !-------------------------------------------------------------------!
172 : ! Remove zero entries in coo matrix.
173 : ! zero means abs(x)<eps
174 : !-------------------------------------------------------------------!
175 15 : subroutine remove_zeros(self, eps)
176 : class(ndcoo_mat_t), intent(inout) :: self
177 : real(dp), optional, intent(in) :: eps
178 : real(dp) :: eps1
179 : integer :: i, counter
180 : if (present(eps)) then
181 : eps1=eps
182 : else
183 : eps1=epsilon(1.0_dp)
184 : end if
185 15 : counter=0
186 1330623 : do i=1, self%nnz
187 1330623 : if (abs(self%val%data(i))> epsilon(1.0)) then
188 1283640 : counter=counter+1
189 5206800 : self%ind%data(:,counter) =self%ind%data(:, i)
190 1283640 : self%val%data(counter) = self%val%data(i)
191 : end if
192 : end do
193 15 : self%nnz=counter
194 15 : self%ind%size=counter
195 15 : self%val%size=counter
196 15 : end subroutine remove_zeros
197 :
198 : !-------------------------------------------------------------------!
199 : ! sum duplicate entries (also sort by indices)
200 : !-------------------------------------------------------------------!
201 15 : subroutine sum_duplicates(self)
202 : class(ndcoo_mat_t), intent(inout) :: self
203 30 : integer :: new_ind(self%ndim, self%nnz), i, counter
204 30 : real(dp) :: new_val(self%nnz)
205 15 : if (self%nnz==0) then
206 0 : self%is_unique=.True.
207 : return
208 : end if
209 15 : call self%remove_zeros()
210 15 : if (.not. self%is_sorted) then
211 12 : call self%sort_indices()
212 : end if
213 15 : counter=1
214 80 : new_ind(:, counter)= self%ind%data(:, 1)
215 15 : new_val(counter)=self%val%data(1)
216 1283640 : do i=2, self%nnz
217 3879670 : if (all(self%ind%data(:, i)==self%ind%data(:, i-1))) then
218 19440 : new_val(counter)=new_val(counter)+self%val%data(i)
219 : else
220 1264185 : counter=counter+1
221 5148400 : new_ind(:, counter)= self%ind%data(:, i)
222 1264185 : new_val(counter)=self%val%data(i)
223 : end if
224 : end do
225 15 : self%nnz=counter
226 5148495 : self%ind%data(:,1:counter)=new_ind(:,1:counter)
227 1264215 : self%val%data(1:counter)=new_val(1:counter)
228 15 : self%ind%size=self%nnz
229 15 : self%val%size=self%nnz
230 15 : self%is_unique=.True.
231 15 : end subroutine sum_duplicates
232 :
233 : !-------------------------------------------------------------------!
234 : ! Get the i'th value of the matrix.
235 : !-------------------------------------------------------------------!
236 0 : function get_val_inz(self, i) result(v)
237 : class(ndcoo_mat_t), intent(inout) :: self
238 : integer, intent(in) :: i
239 : real(dp) :: v
240 0 : v= self%val%data(i)
241 0 : end function get_val_inz
242 :
243 :
244 : !-------------------------------------------------------------------!
245 : ! get all the indices for the ith entry
246 : ! Input:
247 : ! i: ith entry
248 : ! Return:
249 : ! a integer array of indices.
250 : !-------------------------------------------------------------------!
251 5304 : function get_ind_inz(self, i) result(ind)
252 : class(ndcoo_mat_t), intent(inout) :: self
253 : integer, intent(in) :: i
254 : integer :: ind(self%ndim)
255 32160 : ind(:)=self%ind%data(:,i)
256 5304 : end function get_ind_inz
257 :
258 : !-------------------------------------------------------------------!
259 : ! Group the sparse matrix by the first dimension
260 : !> Output:
261 : !> ngroup: number of groups
262 : !> i1_list: list of 1st indices (array(ngroup))
263 : !> istartend: start and end of each group (array(ngroup+1))
264 : !> The starts will be istartend(1:ngroup)
265 : !> The ends will be istartend(2: ngroup+1)-1
266 : !-------------------------------------------------------------------!
267 3 : subroutine group_by_1dim(self, ngroup, i1_list, istartend)
268 : class(ndcoo_mat_t), intent(inout) :: self
269 : integer, intent(inout) :: ngroup
270 : integer, allocatable, intent(inout) :: i1_list(:), istartend(:)
271 :
272 : integer :: i, ii
273 3 : type(int_array_type) :: j1, jstartend
274 :
275 3 : if (.not. (self%is_unique)) then
276 3 : call self%sum_duplicates()
277 : end if
278 3 : if (self%nnz<1) then
279 : ngroup=0
280 3 : else if (self%nnz==1) then
281 0 : i=1
282 0 : ii=self%ind%data(1,i)
283 0 : call j1%push(ii)
284 0 : call jstartend%push(1)
285 0 : call jstartend%push(2)
286 : else
287 3 : i=1
288 3 : ii=self%ind%data(1,i)
289 3 : call j1%push(ii)
290 3 : call jstartend%push(i)
291 1109568 : do i=2, self%nnz
292 1109565 : ii=self%ind%data(1,i)
293 1109568 : if(ii == self%ind%data(1, i-1)) then
294 : cycle
295 : else
296 2317 : call j1%push(ii)
297 1109565 : call jstartend%push(i)
298 : end if
299 : end do
300 3 : call jstartend%push(self%nnz+1)
301 : end if
302 3 : ngroup=j1%size
303 3 : if(ngroup>0) then
304 9 : ABI_MALLOC(i1_list, (ngroup))
305 9 : ABI_MALLOC(istartend, (ngroup+1))
306 2323 : i1_list(:)=j1%data(1: j1%size)
307 2326 : istartend(:)=jstartend%data(1: jstartend%size)
308 : end if
309 3 : call j1%finalize()
310 3 : call jstartend%finalize()
311 3 : end subroutine group_by_1dim
312 :
313 : !-------------------------------------------------------------------!
314 : ! Group the sparse matrix by first two indices
315 : !> Output:
316 : !> ngroup: number of groups
317 : !> ilist: list of one index of the pair (array(ngroup))
318 : !> jlist: list of other index of the pair (array(ngroup))
319 : !> ijstartend: start and end of each group (array(ngroup+1))
320 : !> The starts will be ijstartend(1:ngroup)
321 : !> The ends will be ijstartend(2: ngroup+1)-1
322 : !-------------------------------------------------------------------!
323 :
324 2 : subroutine group_by_pair(self)
325 : class(ndcoo_mat_t), intent(inout) :: self
326 :
327 : integer :: i, ii, ij
328 : ! type(int_array_type) :: i1, j1, startend
329 :
330 2 : if((self%is_pair_grouped)) return
331 :
332 2 : if (.not. (self%is_unique)) then
333 2 : call self%sum_duplicates()
334 : end if
335 2 : if (self%nnz<1) then
336 :
337 2 : else if (self%nnz==1) then
338 : i=1
339 0 : ii=self%ind%data(1,i)
340 0 : ij=self%ind%data(2,i)
341 0 : call self%pair_1list%push(ii)
342 0 : call self%pair_2list%push(ij)
343 0 : call self%pair_startend%push(1)
344 0 : call self%pair_startend%push(2)
345 : else
346 2 : i=1
347 2 : ii=self%ind%data(1,i)
348 2 : ij=self%ind%data(2,i)
349 2 : call self%pair_1list%push(ii)
350 2 : call self%pair_2list%push(ij)
351 2 : call self%pair_startend%push(i)
352 139968 : do i=2, self%nnz
353 139966 : ii=self%ind%data(1,i)
354 139966 : ij=self%ind%data(2,i)
355 139968 : if(ii == self%ind%data(1, i-1) .and. ij == self%ind%data(2, i-1)) then
356 : cycle
357 : else
358 7774 : call self%pair_1list%push(ii)
359 7774 : call self%pair_2list%push(ij)
360 7774 : call self%pair_startend%push(i)
361 : end if
362 : end do
363 2 : call self%pair_startend%push(self%nnz+1)
364 : end if
365 :
366 2 : self%is_pair_grouped = .true.
367 :
368 : end subroutine group_by_pair
369 :
370 :
371 : !-------------------------------------------------------------------!
372 : ! Get the indices of the dim'th dimension
373 : ! Input:
374 : ! dim: dimension
375 : ! Returns:
376 : ! a integer array(nnz)
377 : !-------------------------------------------------------------------!
378 0 : function get_ind(self, dim) result(ilist)
379 : class(ndcoo_mat_t), intent(inout) :: self
380 : integer, intent(in) :: dim
381 : integer :: ilist(self%nnz)
382 0 : ilist(:)=self%ind%data(dim, 1:self%nnz)
383 0 : end function get_ind
384 :
385 :
386 : ! matrix vector product
387 0 : subroutine mv1vec(self, vec, iv, res)
388 : class(ndcoo_mat_t), intent(inout) :: self
389 : real(dp), intent(in) :: vec(:)
390 : integer, intent(in) :: iv ! which index is used for multiplication
391 : class(ndcoo_mat_t), intent(inout) :: res ! result
392 :
393 : integer :: iind, iiv, j, jv
394 0 : integer :: ind(1:res%ndim)
395 : real(dp) :: val
396 :
397 0 : if(self%ndim .ne. res%ndim+1) then
398 0 : ABI_ERROR('Dimension of resulting matrix is not equal to (dimension of initial matrix -1)')
399 : endif
400 :
401 0 : do iind =1 , self%nnz
402 0 : iiv=self%ind%data(iv, iind)
403 0 : jv=0
404 0 : do j=1, self%ndim
405 0 : if(j.ne.iv) then
406 0 : jv=jv+1
407 0 : ind(jv) = self%ind%data(j, iind)
408 : endif
409 : enddo
410 0 : val = self%val%data(iind)*vec(iiv)
411 0 : call res%add_entry(ind, val)
412 : end do
413 :
414 0 : call sum_duplicates(res)
415 :
416 0 : end subroutine mv1vec
417 :
418 :
419 1 : subroutine mv2vec(self, veci, vecj, iv, jv, res)
420 : class(ndcoo_mat_t), intent(inout) :: self
421 : real(dp), intent(in) :: veci(:), vecj(:) ! vectors to be multiplied with
422 : integer, intent(in) :: iv, jv ! which indices are used for multiplication
423 : class(ndcoo_mat_t), intent(inout) :: res ! result
424 :
425 : integer :: iind, iiv, ijv, jnew, j
426 2 : integer :: ind(1:res%ndim)
427 : real(dp) :: val
428 :
429 1 : if(self%ndim .ne. res%ndim+2) then
430 0 : ABI_ERROR('Dimension of resulting matrix is not equal to (dimension of initial matrix -2)')
431 : endif
432 69985 : do iind =1 , self%nnz
433 69984 : iiv=self%ind%data(iv, iind)
434 69984 : ijv=self%ind%data(jv, iind)
435 69984 : jnew=0
436 349920 : do j=1, self%ndim
437 349920 : if(j.ne.iv .and. j.ne.jv) then
438 139968 : jnew=jnew+1
439 139968 : ind(jnew) = self%ind%data(j, iind)
440 : endif
441 : enddo
442 69984 : val = self%val%data(iind)*veci(iiv)*vecj(ijv)
443 69985 : call res%add_entry(ind, val)
444 : end do
445 1 : call sum_duplicates(res)
446 :
447 1 : end subroutine mv2vec
448 :
449 :
450 4004 : subroutine vec_product2d(self, iv, veci, rv, res)
451 : class(ndcoo_mat_t), intent(inout) :: self
452 : real(dp), intent(in) :: veci(:)
453 : integer ,intent(in) :: iv, rv !
454 : real(dp), intent(inout) :: res(:)
455 : integer :: iind, iiv, irv
456 15571556 : do iind =1 , self%nnz
457 15567552 : iiv=self%ind%data(iv, iind)
458 15567552 : irv=self%ind%data(rv, iind)
459 15571556 : res(irv) = res(irv) + self%val%data(iind) * veci(iiv)
460 : end do
461 4004 : end subroutine vec_product2d
462 :
463 :
464 :
465 : ! matrix vector vector product. matrix should be dim3.
466 : ! n(vector)=ndim-1
467 : ! which returns a vecor
468 : ! res_r = \sum_ij M_{ijr} V_i V_j
469 : ! i, j, r can be in any order.
470 12014 : subroutine vec_product(self, iv, veci, jv, vecj, rv, res)
471 : class(ndcoo_mat_t), intent(inout) :: self
472 : real(dp), intent(in) :: veci(:), vecj(:)
473 : integer ,intent(in) :: iv, jv, rv !
474 : real(dp), intent(inout) :: res(:)
475 : integer :: iind, iiv, ijv, irv
476 6259209902 : do iind =1 , self%nnz
477 6259197888 : iiv=self%ind%data(iv, iind)
478 6259197888 : ijv=self%ind%data(jv, iind)
479 6259197888 : irv=self%ind%data(rv, iind)
480 6259209902 : res(irv) = res(irv) + self%val%data(iind) * veci(iiv)*vecj(ijv)
481 : end do
482 12014 : end subroutine vec_product
483 :
484 :
485 : ! matrix vector vector vector product. matrix should be dim4.
486 : ! n(vector)=ndim-1
487 : ! which returns a vecor
488 : ! res_r = \sum_ijk M_{ijkr} V_i V_j V_k
489 : ! i, j, k, r can be in any order.
490 6006 : subroutine vec_product4d(self, veci, vecj, kv, veck, rv, res)
491 : class(ndcoo_mat_t), intent(inout) :: self
492 : real(dp), intent(in) :: veci(:), vecj(:), veck(:)
493 : integer ,intent(in) :: kv, rv
494 : real(dp), intent(inout) :: res(:)
495 :
496 : integer :: iind, iiv, ijv, ikv, irv, igroup, istart, iend
497 : real(dp) :: scalprod
498 :
499 6006 : if(.not.(self%is_pair_grouped)) then
500 0 : call self%group_by_pair()
501 : endif
502 :
503 23357334 : do igroup = 1, self%pair_1list%size
504 : !precalculate scalar product of first and second columns for each group
505 23351328 : istart=self%pair_startend%data(igroup)
506 23351328 : iend=self%pair_startend%data(igroup+1)-1
507 23351328 : iiv=self%pair_1list%data(igroup)
508 23351328 : ijv=self%pair_2list%data(igroup)
509 23351328 : scalprod=veci(iiv)*vecj(ijv)
510 443681238 : do iind=istart, iend
511 420323904 : ikv=self%ind%data(kv, iind)
512 420323904 : irv=self%ind%data(rv, iind)
513 443675232 : res(irv) = res(irv) + self%val%data(iind) * veck(ikv)* scalprod
514 : end do
515 : enddo
516 6006 : end subroutine vec_product4d
517 :
518 :
519 :
520 0 : subroutine test_ndcoo()
521 0 : type(ndcoo_mat_t) :: m
522 : integer :: ngroup
523 0 : integer, allocatable :: i1list(:), ise(:)
524 0 : call m%initialize(mshape=[3,3,3])
525 0 : call m%add_entry(ind=[3, 2,1], val=0.3d0)
526 0 : call m%add_entry(ind=[1, 2,1], val=0.3d0)
527 0 : call m%add_entry(ind=[1, 2,1], val=0.4d0)
528 0 : call m%add_entry(ind=[3, 2,1], val=0.5d0)
529 0 : call m%add_entry(ind=[1, 1,2], val=0.5d0)
530 0 : call m%add_entry(ind=[2,5,1], val=0.0d0)
531 : !call m%print()
532 0 : call m%sort_indices()
533 0 : call m%sum_duplicates()
534 : !print *, "After sum"
535 : !call m%print()
536 : !print *, "Grouping"
537 0 : call m%group_by_1dim(ngroup, i1list, ise)
538 : !print *, "ngroup: ", ngroup
539 : !print *, "i1list: ", i1list
540 : !print *, "ise: ", ise
541 0 : ABI_SFREE(i1list)
542 0 : ABI_SFREE(ise)
543 0 : end subroutine test_ndcoo
544 :
545 0 : end module m_spmat_NDCOO
546 :
|
