Branch data 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 : :
|
